貴州烏江流域NDVI降尺度適應性評價

楊江州，周旭，熊軍，周志凱，韋小茶

(1.貴州省地質礦產勘查開發局 117地質大隊，貴陽 550018；2.貴州師范大學 地理與環境科學學院，貴陽 550025)

0 引言

植被是陸地生態系統的主體，具有廣泛性、復雜性、豐富性的物種資源，是地表水循環、生物循環的重要途徑[1-2]。隨著遙感技術的發展，其觀測面積廣、信息量大、分辨率高、時間序列長等優勢，在植被監測中得到廣泛應用[3-5]。歸一化植被指數(normalized difference vegetation index，NDVI)是目前公認的表征植被變化的最有效參數之一，已廣泛用于植被生長及其動態變化研究中[5]。王茜等[6]應用NDVI數據研究我國30年NDVI時空變化，并對氣候因素影響進行分析；蒙吉軍等[7]基于NDVI、NPP數據研究西南喀斯特植被變化對氣候的響應；朱林富等[8]對山地城市進行植被變化時空分異特征分析；張繼等[9]基于NDVI對生態工程建設背景下貴州高原植被變化及因素分析；王冰等[10]基于NDVI數據分析貴州喀斯特地區植被覆蓋變化趨勢。在降尺度研究中，統計降尺度方法得到廣泛運用。如Huth[11]運用神經網絡與線性回歸方法對歐洲日氣溫進行降尺度研究，發現較為簡單的線性方法也可以取得較好的降尺度效果；Immerzeel等[12]通過建立TRMM降水數據與NDVI 關系，對伊比利亞半島的TRMM年降水數據進行降尺度研究，將RMM 數據的空間分辨率從 0.25°提高至 1 km。之后Duan等[13]、李凈等[14]、范雪薇等[15]在此基礎上引入地形因子進行多元回歸，實現了TRMM降水數據降尺度，且能夠體現更細致的降水分布特征。然而在植被變化的研究中，高精度、小尺度的時間序列均較短，多為2000年至今的研究，對植被降尺度研究較為缺乏，且1982—2000年間NDVI產品多是8 km低分辨率數據。表明利用該尺度數據揭示植被變化僅能反映大范圍的總體特征，小尺度探究植被變化存在不足。

貴州烏江流域地理環境空間異質性大，高分辨率遙感數據揭示中小流域NDVI變化及其驅動因素尤其重要[16-17]。烏江流域地處中國西南喀斯特腹地，巖溶發育典型，土壤瘠薄，水資源賦存條件極弱[18-19]。流域特殊的地質背景，再加上人類活動影響，生態環境不斷退化，使其演化為典型生態脆弱區[20-21]。因此，研究喀斯特流域NDVI長時間變化趨勢，有助于揭示生態環境演化特征，減弱不合理人類活動對其的影響，對山地生態文明建設具有積極意義。本文基于空間統計降尺度方法，探究喀斯特流域NDVI的降尺度與時空變化。研究成果對低分辨率NDVI數據降尺度獲取,喀斯特地區石漠化治理與生態建設具有重要指導意義。

1 研究區概況、數據與方法

1.1 研究區概況

貴州烏江流域位于中國西南腹地，自西向東橫貫整個貴州省北部，向西與云南省接壤，向北與四川省和重慶市接壤。流域干流全長847 km，是貴州省最大流域，支流較多，呈羽狀水系分布。流域屬亞熱帶季風性濕潤氣候，氣候溫和，降水豐沛，雨熱同季，年平均氣溫15 ℃，年平均降水量約1 200 mm。流域地質構造上主要屬于揚子準地臺中的黔北臺隆，在地域獨特的水熱條件的驅動下，碳酸鹽巖被流水侵蝕，使流域內巖溶地貌發育強烈，75.6%的地區是石灰巖、白云巖等碳酸鹽巖[22]。由于碳酸鹽巖抗侵蝕能力較強，母質造土能力差，成土過程緩慢，不利于土壤的形成與植被生長，使其基巖裸露率高，石漠化敏感性強。地貌類型主要分為上游的巖溶峽谷地貌，中游的巖溶高原地貌，下游的巖溶槽谷地貌，平均海拔約1 100 m。

1.2 數據來源

GIMMS NDVI 數據來自NOAA網站(https：//www.noaa.gov)網站下載，空間分辨率為8 km，時間分辨率為15 d，時間序列為1982—2015年；MODIS NDVI數據來自美國NASA網站(https：//modis.gsfc.nasa.gov)網站下載MODIS13 Q1，空間分辨率為250 m，時間分辨率為16 d，時間序列為2000—2016年。

1.3 研究方法

本文采用的NDVI數據是基于GIMMS和MODIS13 Q1 2種遙感產品獲取，都源于美國國家航空航天局(National Aeronautics and Space Administration，NASA)的一系列遙感衛星。GIMMS NDVI數據集是使用時間較長的植被數據，該數據空間分辨率為8 km，時間分辨率是15 d，下載1982—2015年數據，總共816幅影像，利用Matlab編程語言轉換GIMMS NDVI3g.v1原始產品為Geottif格式，采用WGS_1984地理坐標和Albers Equal-Area Conic投影坐標。MODIS NDVI數據是植被數據中具有較高分辨率且應用廣泛的數據，空間分辨率為250 m，時間分辨率為16 d，下載2000—2016年數據，總共408幅影像，利用MRT(MODIS reprojection tools)軟件對數據進行投影轉換、數據拼接等預處理，采用與GIMMS NDVI數據一致的地理坐標和投影坐標。再運用ENVI/ArcGIS對數據集進行裁剪、拼接、幾何校正等預處理工作，同時剔除NDVI產品的無效填充值。采用最大值合成法(maximum value composites，MVC)將每月中2期NDVI影像提取最大值，作為每月的NDVI值，此方法能進一步消除大氣、衛星軌道飄移、太陽高度角等的干擾，在地形復雜的山地具有較好的效果。其計算如式(1)所示。

NDVIi=Max(NDVIij)

(1)

式中：NDVI為歸一化植被指數；NDVIi指第i個月或者第i年的NDVI；NDVIij指第i年的第j個半月的NDVI或第i年第j月的NDVI。

回歸模式法是應用最早和最廣泛的統計降尺度技術，是通過建立預報變量與預測因子之間線性或非線性關系的一種概念性方式[11]。應用IDL(interactive data langguage)編寫線性回歸方程代碼，實現NDVI數據降尺度，獲取更早時間序列及較高空間分辨率的數據。其原理利用2001—2015年的GIMMS NDVI與MODIS NDVI數據進行15年逐像元回歸分析，且該GIMMS NDVI數據的空間分辨率都重采樣為250 m分辨率，與MODIS NDVI分辨率一致。得出2001—2015年的GIMMS NDVI與MODIS NDVI數據像元回歸方程的參數后，基于1982—2000年GIMMS NDVI數據，反推出各像元上的NDVI值，從而得到1982—2000年250 m的NDVI數據，具體流程如圖1所示。式(2)為線性回歸方程。

GNi=a+b×MNi+ci

(2)

式中：GNi代表第i月、年的GIMMS NDVI；MNi代表第i月、年的MODIS NDVI；a和b為常數；c為回歸方程的殘差。

圖1 降尺度流程圖

采用平均絕對誤差(M1)、平均相對誤差(M2)、均方根誤差(M3)、Theil不等系數(M4)探究GIMMS NDVI數據空間降尺度得到的NDVI和MODIS NDVI數據的精度，對其做時間序列與空間的可靠性分析(式(3)～式(6))[23-24]。

(3)

(4)

(5)

(6)

式中：ai實際值；bi為模擬值；n為檢驗點的個數；M4介于0和1之間，M4越趨近0，則代表實際值與模擬值誤差越小，模型模擬的精度越高。

利用線性趨勢法進行NDVI時間趨勢的分析，模擬出每個柵格的變化趨勢，以單個像元時間變化特征綜合反映出整個區域的時空演變格局和空間變化規律。隨著時間序列的改變，NDVI表現為序列整體的上升或下降變化趨勢，呈現出空間分布格局的變化和在某時刻出現突變或轉折，NDVI變量可以看作是時間的一元線性回歸，利用最小二乘法逐像元計算NDVI變化斜率如式(7)所示[25-26]。

(7)

式中：B為某像元的趨勢線斜率；t為年份；n為研究年段35(1982—2016年)。當B值為正值時，表現為隨時間t變化NDVI呈上升趨勢；反之，當B為負值時，隨時間t變化NDVI呈下降趨勢。其絕對值越大表示研究區域NDVI上升或下降的趨勢越顯著。

2 統計降尺度模型結果驗證

2.1 基于降尺度數據的可靠性分析

由于GIMMS與MODIS數據在2000年開始重疊，以2000年GIMMS NDVI降尺度為例，更能說明降尺度的效果特征，如圖2所示。二者對比，整體來看，分辨率為8 km的GIMMS NDVI降尺度后得到分辨率250 m的NDVI影像更加清晰，更有利于揭示流域NDVI變化；雖然在GIMMS NDVI中，烏江流域西部和中南部的低值在統計降尺度后有一定的變化，但整體的趨勢都呈現西北低、東南高的變化趨勢。表明降尺度的NDVI整體分布狀況與GIMMS NDVI在空間上有著較好的一致性。

圖2 2000年貴州烏江流域GIMMS NDVI降尺度效果圖

2.2 降尺度數據的時間驗證

利用1982—2000年的GIMMS NDVI數據和降尺度NDVI數據的年際平均值統計尺度和多年各月平均值統計尺度進行時間相關性檢驗。年際平均值統計尺度相關性驗證為1982—2000年全省平均的年時間序列，有19個時間點；多年各月平均值統計尺度相關性驗證為1982—2000年多年各月平均值組成的站點序列，有12個時間點。

以年尺度時間序列來看，運用1982—2000年的GIMMS NDVI數據和降尺度NDVI數據做誤差檢驗分析(表1)，結果表明：從實際誤差(M)來看，誤差值的正數值與負數值的個數各在一半，且多年總和為0.007 11，表明GIMMS數據降尺度得到NDVI數據與實際波動較為合理，誤差較小。在絕對誤差與相對誤差中，1993年絕對誤差(M1)與相對誤差(M2)達到最高，分別為0.013 02、0.022 95；1992年絕對誤差(M1)與相對誤差(M2)值最低，分別為0.000 22、0.000 37。從數值趨勢來看，90年代與80年代相比，具有上升的趨勢。19年間，二者均方根誤差(M3)為0.006 52，Theil不等系數(M4)為0.005 45，雖各年份的誤差均有波動，但整體較小。

表1 年際平均值誤差檢驗

從月尺度時間序列來看，運用1982—2000年各月均值的GIMMS NDVI數據和降尺度NDVI數據做誤差檢驗分析(表2)，結果表明：實際誤差(M)中，有7個正值，5個負值，大致各在一半，同樣表明降尺度數據的合理性。從絕對誤差(M1)與相對誤差(M2)來看，6月的多年各月平均值絕對誤差(M1)與相對誤差(M2)達到最高，分別為0.084 32、0.126 32；原因可能在于6月是植被茂盛時期，不同植被每年長勢不一樣，而降尺度數據的NDVI被多年數值回歸后導致誤差偏大；9月絕對誤差(M1)與相對誤差(M2)值最低，分別為0.000 70、0.001 00。多年各月平均值均方根誤差(M3)為0.029 56，Theil不等系數(M4)為0.029 45，相對于年尺度變化可發現，誤差有所增大，但整體降尺度的月份趨勢與GIMMS的月份趨勢大致一樣，可以反映月份序列的變化。

表2 多年(1982—2000年)各月平均值誤差檢驗

圖3 年際平均值統計尺度和多年各月平均值統計尺度擬合分析

根據GIMMS NDVI數據和降尺度NDVI數據做年際平均值統計尺度和多年各月平均值統計尺度相關性分析(圖3)。年際平均值統計尺度相關分析的線性方程為y=0.635 6x+0.215，R2=0.998 9；多年各月平均值統計尺度相關分析的線性方程為y=1.238 6x-0.139 8，R2=0.949 6，都且在0.05顯著性水平上顯著相關。表明降尺度得到的NDVI數據與GIMMS NDVI數據擬合度較好。

綜上所述，年際平均值統計尺度和多年各月平均值統計尺度都具有擬合精度較高，良好地反映時間變化趨勢，說明降尺度得到的NDVI數據在時間序列上具有一定的代表性。但結合表1和表2比較來看，年際平均值統計尺度模擬更優。原因可能在于年際平均值統計尺度模型中，2000—2015年進行回歸數據較少，一年1個數據，總共15個數據；而多年月均值統計尺度，一年模擬12個數據，總共180數據。表明多數據線性回歸后，得到的參數，模擬NDVI具有均一性，而各年月氣候環境的不確定性，使之降尺度得到的NDVI值誤差精度增大，但整體誤差值較小，誤差精度滿足NDVI時間變化的研究需要。

2.3 降尺度數據的空間驗證

以2000年MODIS NDVI數據與2000年降尺度NDVI數據進行空間上的相關性檢驗，進行空間誤差檢驗分析(圖4)。從圖4可知，2000年降尺度的NDVI值相對于MODIS數據來說，最大值從0.82變為0.84，相差0.02，整體值有所增加。但整體均值相差不大，差值為0.003 96，變化趨勢保持較好的一致性。

圖4 2000年貴州烏江流域MODIS NDVI與降尺度NDVI分布圖

根據以上原理方法，利用ENVI對柵格數據進行波段運算，求出降尺度NDVI值與MODIS NDVI值的平均絕對誤差(M1)、平均相對誤差(M2)、均方根誤差(M3)、Theil不等系數(M4)，進行統計分析得出表3，再運用ArcGIS進行數據處理得出空間圖(圖4)。由于NDVI絕對值處于0～1之間，采用理論值為小于等于5%的誤差，視為精度效果好；理論值為5%～10%之間的誤差，視為精度效果較好；大于等于10%的誤差，說明誤差較大。因此，本文按照實際情況誤差值的大小分為3級：小于等于5%，精度良好；5%～10%，精度較好；大于等于10%，精度較差。

表3 2000年MODIS NDVI值與降尺度NDVI值誤差占比 %

結合表3與圖5來看，可知平均絕對誤差(M1)、均方根誤差(M3)、Theil不等系數(M4)空間占比小于等于5%誤差值的大部分集中于中部與東北部，誤差占比分別為80.82%、80.82%、86.20%；5%～10%誤差值集中于西部，誤差占比分別為18.11%、18.11%、13.02%；較少出現大于等于10%誤差值，占比分別為1.07%、1.07%、0.80%。雖然平均相對誤差(M2)大于等于10%的誤差值有所增加，占比為15.50%，集中分布西部，在中部和東部都有零散的分布，但整體小于等于10%的占比為85.85%，整體趨勢精度較好。表明降尺度的NDVI值可以良好地反映在空間的動態趨勢。

圖5 2000年貴州烏江流域MODIS NDVI值與降尺度NDVI值誤差空間檢驗

3 植被時空變化特征及其趨勢變化分析

3.1 植被年際時空變化特征

以1982—2000年降尺度NDVI數據與2001—2016年MODIS NDVI數據相結合，統計出烏江流域1982—2016年NDVI的年際變化，如圖6所示。從整體趨勢來看，多年的線性方程為y=0.001 4x-2.185 6，R2=0.321 6，貴州烏江河流域NDVI變化呈明顯增加波動變化趨勢。根據烏江流域NDVI的年際變化，可劃分為2個階段：第1個階段為1982—2000年，流域NDVI值從0.590 1下降到0.564 2，線性方程為y=-9E-0.5x+0.764 3，R2=0.001 6，有向下的波動變化趨勢，但并不明顯；第2個階段為2000—2016年，流域NDVI值從0.564 2上升到0.673 1，線性方程為y=0.006x-11.508，R2=0.827 7，NDVI值有明顯的增加趨勢。其原因在于，20世紀80年代處于改革開發初期，烏江流域居民對植被保護的意識不高，再加上人類活動對植被的亂砍濫伐，造成第一個階段不明顯的下降波動趨勢。90年代末期，國家加大西南喀斯特地區的石漠化治理，使烏江流域植被逐步恢復改善，因此形成第2階段明顯增加的波動趨勢。

圖6 1982—2016年貴州烏江流域植被年際變化

應用ArcGIS對每年的NDVI影像進行柵格運算，最后得出多年平均NDVI空間分布，如圖7所示。貴州烏江河流域逐年NDVI值具有明顯的空間分布差異，由于受水平地帶性、氣候和地形地貌等因素的綜合影響，呈現出東高西低的分布格局，NDVI值為0.06～0.80。NDVI低值區域主要分布于流域中西部，從西部區域看，形成低值區域原因為該區域處于貴州地勢較高山區，熱量較低，年降雨量偏少，不利于植被生長；從中部區域來看，低值區域主要集中于貴陽市與遵義市周邊，其中貴陽市周邊的低值區域最大，主要原因是城市的開發增強，土地利用格局發生深刻改變，影響了植被生長。這表明人類活動頻繁區域對植被生長起到抑制作用，因此在開發建成區中，應加強城市建成區的綠化，改善人居生態環境。

圖7 1982—2016年貴州烏江流域NDVI多年平均空間分布

3.2 植被趨勢變化分析

基于利用線性趨勢法，在像元尺度上分析了貴州烏江河流域NDVI年際變化趨勢(圖8)。按照實際情況趨勢線斜率值的大小分為5級：B<-0.001為顯著下降；-0.001≤B<-0.000 5為輕微下降；-0.000 5≤B≤0.000 5為基本不變；0.000 50.001為顯著上升。從1982—1999年趨勢來看，基本不變的趨勢占主體，占總面積的43%，而上升趨勢與下降趨勢呈零散分布狀，其中上升趨勢占總面積的42%，下降趨勢占總面積的15%，整體呈增加趨勢。原因在于20世紀90年代石漠化治理在21世紀初期取得了初步成效，生態壞境得到了改善，植被覆蓋度增加。

圖8 貴州烏江流域NDVI年際變化趨勢

4 討論與結論

4.1 討論

烏江流域地處貴州喀斯特高原，多云雨霧天氣，覆蓋研究區的Landsat等高分辨率遙感數據云量較多，數據質量較差，難以得到時間序列連續數據，而SPOT、AVHRR等數據滿足長時間序列，但是空間分辨率極低，再加上喀斯特山區地形特征差異顯著，可進入性極差，地表實測數據等基礎資料相對缺乏。本文基于統計降尺度模型，將1982—2000年GIMMS NDVI數據的8 km分辨率提升到250 m，且與2000—2016年MODIS NDVI數據(分辨率為250 m)相結合，從而得到1982—2016年250 m的長時間序列高分辨率NDVI數據。通過該方法獲得研究區長時間序列高分辨率的NDVI數據，高分辨率的NDVI數據更能提高喀斯特小空間尺度的植被變化參數特征，有利于揭示喀斯特山區流域植被變化及其驅動因素[14]。因此，基于NDVI數據建立2種各像元間的統計關系，實現數據統計降尺度，從而延長插補時間序列獲取高分辨率長時間序列數據的方法。該方法也可對今后在難以獲取長時間序列數據或小空間尺度需求更高分辨率的相關研究提供參考。

王秀春等[27]基于1998—2006年SPOT生長季(4—9月)逐旬最大合成NDVI數據，分析了烏江流域近10年的植被覆蓋動態變化，結果數據顯示 NDVI年際變化呈增加趨勢，年增長率為0.009/a，流域內地級市及省會城市呈減少的趨勢。在文獻[9-10]基于NDVI數據對貴州喀斯特地區植被覆蓋變化研究中，均發現貴州喀斯特地區總體處于上升趨勢，降水對植被變化影響減弱。相對于前人的研究，本文以烏江流域為研究區，得出1982—2016年研究區的NDVI年際變化同樣呈增加趨勢，年增長率為0.006/a；流域植被年際趨勢分布變化中，人類活動較大的地方植被變化呈減少的趨勢；從整體來看，本研究的結果與前人研究基本一致。

在1982—2016年貴州烏江流域植被年際變化中(圖6)，NDVI在2000年出現明顯下降，原因可能為該年降水量偏小，年均溫度較低，不利植被生長。該年現象與文獻[28-29]研究NDVI變化趨勢基本一致。但本文方法只是數據的交叉檢驗，沒有輔助數據的佐證，缺乏一定的物理機制，在未來模型可以通過加入第3個輔助數據來進一步的探討。同時也存在一定的不確定性，如 NDVI 在高植被覆蓋區易出現飽和的問題；不同傳遞函數與降尺度結果精度的關系等問題。

4.2 結論

通過IDL編程，實現植被遙感數據的空間統計降尺度，對降尺度數據進行可靠性分析后將得到降尺度NDVI數據與MODIS數據相結合，得出貴州烏江流域NDVI時空變化，結論如下。

1)通過數據的處理與空間統計降尺度模型應用，將1982—2000年GIMMS NDVI數據的8 km分辨率降尺度得到250 m NDVI數據，且降尺度的NDVI數據通過了GIMMS NDVI數據和MODIS數據在空間尺度與時間尺度的誤差檢驗，能夠滿足研究區NDVI空間變化的研究需要。

2)基于IDL編程，實現了植被遙感數據的空間統計降尺度，降尺度的NDVI數據良好反映時空變化的趨勢，表明該方法具有一定的適用性及科學性。

3)1982—2016年，貴州烏江河流域NDVI變化呈明顯增加波動變化趨勢，在逐年NDVI空間分布中，呈現東高西低的分布格局。

4)近35年來，貴州烏江河流域NDVI年際變化趨勢呈上升趨勢。表明20世紀90年代石漠化治理在21世紀初期取得初步成效，生態壞境得到了改善，植被覆蓋度也有所增加。