閆俊杰， 張 靜， 雷 雨， 沙吾麗·達吾提拜， 呂光輝

(1. 新疆大學資源與環(huán)境科學學院， 新疆 烏魯木齊 830046； 2. 中國科學院新疆生態(tài)與地理研究所， 新疆 烏魯木齊 835000；3. 中國科學院大學， 北京 100049；4. 伊犁師范學院生物與地理科學學院， 新疆 伊寧 835000；5. 新疆水利廳， 新疆 烏魯木齊， 830000)

我國天然草原面積有4億hm2[1]，廣泛分布于我國西南及西北地區(qū)。草地不僅是我國陸地生態(tài)的重要組成，具有重要的生態(tài)地位，同時草地還是我國畜牧業(yè)發(fā)展的基本生產(chǎn)資料，對我國西南及西北經(jīng)濟發(fā)展和社會穩(wěn)定意義重大。然而據(jù)資料顯示，我國90%的天然草原出現(xiàn)不同程度的退化，其中30%出現(xiàn)嚴重退化[2]。伊犁河谷位于我國西北邊陲，河谷內草地發(fā)育良好，是我國優(yōu)質的畜牧基地[3]。目前，隨著伊犁河谷經(jīng)濟發(fā)展和人口增多，人類對草地生態(tài)的干擾逐步加劇，以生產(chǎn)力降低為典型特征的草地退化問題日趨嚴重[4]，影響當?shù)厣鷳B(tài)穩(wěn)定，制約畜牧業(yè)健康發(fā)展。開展草地生態(tài)動態(tài)監(jiān)測，加強草地退化規(guī)律研究，對該區(qū)草地生態(tài)系統(tǒng)恢復與重建具有重要意義。

近30年來，遙感數(shù)據(jù)以其區(qū)域性及時效性等優(yōu)勢日漸成為區(qū)域植被動態(tài)監(jiān)測的重要數(shù)據(jù)源[5]。由遙感數(shù)據(jù)反演的植被指數(shù)(vegetation index，VI))可以定量表達植被的生長狀況，已被廣泛用于植被動態(tài)的參數(shù)化表達[6]。在眾多VI之中，歸一化植被指數(shù)(normalized difference vegetation index，NDVI)對植被類型、生長狀態(tài)及生產(chǎn)力反應敏感，是目前應用最廣的VI[5-6]。MODIS NDVI數(shù)據(jù)起始于2000年，時間序列完整，且相對AVHRR NDVI及SPOT-VGT NDVI時間序列數(shù)據(jù)，其在空間分辨率上也具有一定優(yōu)勢[7-8]，已被廣泛應用于不同區(qū)域植被動態(tài)監(jiān)測[7-11]。

伊犁河谷不僅是新疆西部重要生態(tài)屏障，還是重要畜牧基地，該區(qū)草地變化一直被廣泛關注，如：周鑫及張旭琛等人[12-13]利用MODIS數(shù)據(jù)對其生物量反演進行了研究；尚二萍等人[14]對其碳存儲動態(tài)進行了評估；閆俊杰等人[15-16]對其草地退化及草地退化對生態(tài)服務價值的影響進行了分析；劉芳等人[17]研究了伊犁河谷草地NDVI變化及其對降水的敏感性。前人的研究為今后伊犁河谷草地保護和恢復提供了重要參考，但目前，針對伊犁河谷草地植被覆蓋變化趨勢的研究仍相對缺乏，尤其是對變化趨勢持續(xù)性的研究，基于此，本文利用MODIS NDVI遙感數(shù)據(jù)，借助于GIS空間分析技術及Mann-Kendall趨勢分析和檢驗方法，對伊犁河谷草地植被覆蓋的變化趨勢進行分析，并利用變標度極差分析法(rescaled range analysis，R/S分析)，計算其Hurst指數(shù)的空間分布，探討伊犁河谷草地植被覆蓋變化趨勢的持續(xù)性特征，以期為該區(qū)草地生態(tài)恢復與重建提供參考依據(jù)和科學指導。

1 材料與方法

1.1 研究區(qū)概況

伊犁河谷位于新疆西天山區(qū)，行政上包括有伊寧市、霍爾果斯市、可克達拉布、伊寧縣、察布查爾縣、昭蘇縣、特克斯縣、鞏留縣、新源縣、尼勒克縣及霍城縣等3市8縣。伊犁河谷地形特殊，總體像一個向西敞開的喇叭，利于西風帶濕潤水汽抬升凝結而形成降雨，成就了伊犁河谷“西域濕島”及“塞外江南”的美譽。氣候特性方面，伊犁河谷處于溫帶大陸性氣候區(qū)內，但河谷內地形復雜，高山氣候特征明顯，從平原到山區(qū)多年平均降水量的變化幅度達200～1 000 mm，多年平均氣溫變化幅度達9.2～2.8℃，受此影響，河谷內植被垂直分異規(guī)律明顯，其中草地植被分布廣泛，面積最大，是河谷內生態(tài)的主體，分布的草地類型主要有高寒草甸、山地草甸、溫性草甸草原、溫性草原、溫性荒漠草原、溫性荒漠及低平地草甸[12]。

1.2 數(shù)據(jù)來源與預處理

本文用到的NDVI數(shù)據(jù)為MODIS MOD13Q1產(chǎn)品，其空間和時間分辨率為分別為250 m和16 d，時間序列為2000年1月-2016年12月，每年23期數(shù)據(jù)，共計391期。草地類型分布數(shù)據(jù)來自中國科學院中國植被圖編輯委員會繪制的中國1:100萬植被類型分布圖，但本文通過對2015年Landsat 8 OLI影像解譯，剔除了研究區(qū)的非草地植被類型，并對草地類型進行了組合歸類(圖1)，之后對解譯數(shù)據(jù)進行了柵格化處理。對獲得的遙感數(shù)據(jù)除進行了數(shù)據(jù)格式轉換、鑲嵌、投影轉換及研究區(qū)提取等預處理外，為降低噪音信息對影像數(shù)據(jù)質量的影響，還對每年23期NDVI數(shù)據(jù)進行了Savitzky-Golay濾波處理，之后為獲得年NDVI數(shù)據(jù)，對濾波處理后的數(shù)據(jù)進行了最大值合成(max value composite，MVC)。為保證NDVI及草地類型柵格數(shù)據(jù)的空間匹配，兩種數(shù)據(jù)的像元大小均重采樣為50 m×50 m。氣象數(shù)據(jù)來自中國氣象局氣象數(shù)據(jù)中心，包括伊寧市、伊寧縣、霍爾果斯、霍城、察布查爾、尼勒克、鞏留、新源、昭蘇和特克斯10個氣象站年平均氣溫和年累積降水量數(shù)據(jù)。

圖1 伊犁河谷草地類型分布圖Fig.1 Distribution of grassland types in Ili Valley

1.3 趨勢檢驗

用Mann-Kendall檢驗與Theil-Sen median趨勢分析判斷NDVI時間序列的變化趨勢，并量化其變化率[18]。Mann-Kendall檢驗統(tǒng)計量計算過程為：

(1)

(2)

(3)

(4)

式中，sign為符號函數(shù)，NDVIk及NDVIi為時序數(shù)據(jù)集合，n為集合長度。在α顯著性水平下，當|ZC|>U1-α/2時，則時間序列在α水平上變化趨勢顯著。若ZC>0，則變化趨勢上升，若ZC<0，則變化趨勢為下降。

Theil-Sen median趨勢分析用于量化NDVI的變化趨勢，表示單位時間內NDVI的變化量，計算公式為：

(5)

式中10，則NDVI時間序列呈上升趨勢；反之，則呈下降趨勢。

1.4 R/S分析

R/S分析為英國水文學家Hurst在研究尼羅河水文時間序列資料時所提出[19]。其計算過程如下：

對于時間序列NDVI(t)，t=1，2，3…，n，以及任意正整數(shù)p≥1，定義均值序列NDVI(p)為：

(6)

定義累積離差序列NDVI(t,p)為：

(7)

定義極差序列R(p)為：

(8)

定義標準差序列S(p)為：

(9)

Hurst通過長時間的實踐總結，建立了如下關系：

(10)

式中：H為Hurst指數(shù)，通過利用雙對數(shù)坐標{ln(P)，ln(R(p)/S(p))}擬合求的。根據(jù)R/S的原理，當H=0.5時，NDVI時間序列為隨機序列，具有隨機游走的特性；當0.5

2 結果與分析

2.1 NDVI變化趨勢

借助于IDL8.5軟件平臺，計算2000-2016年伊犁河谷全區(qū)及不同類型草地NDVI年平均值，以此繪制全區(qū)和各草地類型年平均NDVI的變化曲線(圖2，圖3)，判斷和計算其變化趨勢和變化率(表1)，分析伊犁河谷草地NDVI總體變化趨勢。

圖2 伊犁河谷全區(qū)草地年均NDVI變化曲線Fig.2 Average annual NDVI variation curve of grassland in Ili Valley

經(jīng)計算，伊犁河谷全區(qū)草地NDVI多年平均值為0.68，表明伊犁河谷草地植被覆蓋總體為中高水平，但是從圖2可以看出，2000-2016年伊犁河谷草地NDVI總體呈下降趨勢，根據(jù)趨勢檢驗結果，其ZC值為-2.35，|ZC|>1.98，下降趨勢達到顯著水平，其變化率為-0.0038單位·a-1。對于NDVI的年際波動，伊犁河谷全區(qū)平均NDVI在2000-2002年持續(xù)升高，于2002年和2003年達到最高值0.74，之后波動下降，于2014年達到最低值0.59，2014-2016又持續(xù)升高，其中2003-2014年間NDVI降低明顯，11年內降低20.27%。

圖3 伊犁河谷各草地類型年均NDVI變化曲線Fig.3 Average annual NDVI variation curve of each grassland type in Ili Valley

對于不同草地類型，由圖3和表1可知，2000-2016年伊犁河谷各草地類型NDVI均呈現(xiàn)下降趨勢，且年間波動相似。但各種植被類型NDVI下降的速率和顯著性均存在一定的差異性。山地草甸和溫性草原NDVI下降速率最大，分別為0.0047單位·a-1和0.0041單位·a-1(表1)，且山地草甸NDVI下降趨勢達到極顯著水平，而溫性草原則未達到顯著水平；其次為高寒草甸、溫性草甸草原和溫性荒漠草原，其NDVI的減小速率分別為0.0030單位·a-1、0.0035單位·a-1和0.0034單位·a-1，其中，高寒草甸NDVI呈極顯著下降，而溫性草甸草原和溫性荒漠草原下降不顯著；溫性荒漠和低平地草甸NDVI的變化速率最小，分別為0.0018單位·a-1和0.0028單位·a-1，下降趨勢均未達到顯著水平。

表1 伊犁河谷全區(qū)及各草地類型NDVI變化趨勢及變化率Table 1 The NDVI changing trend and rate of the whole area and each grassland type in Ili Valley

注：**表示達到極顯著水平，p<0.01;*表示達到顯著水平，p<0.05

Note：**indicates significance at the 0.01 level;*indicates significance at the 0.05 level

2.2 NDVI變化趨勢的空間分異

根據(jù)Mann-Kendall檢驗計算方法，對伊犁河谷草地NDVI進行逐像元計算，獲得其ZC值的空間分布數(shù)據(jù)，以此制作NDVI變化趨勢空間分布圖(圖4)，并統(tǒng)計其面積比例(表2)。

圖4 伊犁河谷草地NDVI變化趨勢空間分布圖Fig.4 Spatial distribution of the NDVI changing trend of the grassland in Ili Valley

由圖4可知，2000-2016年伊犁河谷草地NDVI變化趨勢為下降的區(qū)域占據(jù)河谷草地面積的絕大比例。由表2可知，伊犁河谷88.07%的草地NDVI出現(xiàn)不同程度降低，且49.62%的下降趨勢達到顯著水平，空間上總體呈帶狀分布于河谷北、南及西部烏孫山等區(qū)域的中高山區(qū)，而高海拔的高山區(qū)以及海拔較低的低山、山前洪-沖積扇和丘陵區(qū)草地NDVI雖也以下降為主，但其下降趨勢卻未達到顯著水平。NDVI變化趨勢為升高的草地分布零散，主要分布在河谷南部和東部高海拔的高山區(qū)及伊犁河下游北部山前洪-沖積扇和丘陵區(qū)、尼勒克縣東部丘陵區(qū)、昭蘇盆地東部邊緣和鞏乃斯河上游南岸河谷平原區(qū)等區(qū)域，其面積比例為11.92%，僅有1.05%達到了顯著水平，集中分布在伊犁河出國境口北部部分區(qū)域。

表2 伊犁河谷全區(qū)及各草地類型NDVI不同變化趨勢面積比例Table 2 The area ratio of different NDVI changing trends of the whole area and each grassland type in Ili Valley /%

不同草地類型中，山地草甸是伊犁河谷主要的優(yōu)質草原，分布最為廣泛(圖1)，但，其NDVI下降最為嚴重，其面積的95.76%的NDVI出現(xiàn)不同程度下降，下降趨勢達到顯著水平的比例也高達73.10%(表2)。溫性荒漠NDVI出現(xiàn)下降的面積比例最小，為70.88%，且其中僅有16.00%達到了顯著水平，同時其NDVI呈顯著升高的比例高達4.81%(表2)，遠高于其它類型，這主要是由于整個河谷內NDVI有所升高的草地主要分布于溫性荒漠的分布區(qū)內(圖1和圖4)。剩余草地類型NDVI下降的面積比例相當，比例均介于79%～87%之間。

2.3 持續(xù)性特征

由圖5可知，伊犁河谷草地NDVI Hurst指數(shù)空間差異明顯，但全區(qū)絕大部分區(qū)域H>0.50，其面積比例高達97.77%(表3)，而全區(qū)平均H值也達到0.71(表4)，呈現(xiàn)強的正持續(xù)性，表明2016年之后一定時期內伊犁河谷絕大部分區(qū)域草地NDVI仍將持續(xù)下降。而全區(qū)H<0.50的面積比例僅為2.23%，空間上主要分布在NDVI升高的分布區(qū)內(圖4和圖5)，意味著河谷內面積有限的NDVI升高的區(qū)域在2016年后的一定時期內其NDVI的變化趨勢將由升高轉變?yōu)榻档汀?/p>

具體來看(表3)，全區(qū)56.37%草地NDVI的H>0.70，表現(xiàn)出相對較強的正持續(xù)性，其NDVI為下降趨勢且達到顯著水平的比例為39.73%，而NDVI為升高趨勢且達到顯著水平的比例僅為0.47%；全區(qū)草地NDVI正持續(xù)性相對較弱(0.5

表3 伊犁河谷草地NDVI變化趨勢與其Hurst指數(shù)矩陣表Table 3 Matrix table of the changing trend and Hurst index of the grassland NDVI in Ili Valley/%

圖5 伊犁河谷2000-2016年草地NDVI Hurst指數(shù)空間分布圖Fig.5 Spatial distribution of the NDVI Hurst index of the grassland in Ili Valley from 2000-2016

對于不同草地類型(表4)，全區(qū)各草地類型NDVI變化的H值均>0.5，具有正持續(xù)性特征，表明未來一定時期各草地類型的NDVI均將持續(xù)下降。各草地類型中，山地草甸H值最大，為0.74，同時山地草地還是全區(qū)面積最大、平均NDVI值最大，NDVI下降速率最快的草地類型，可見山地草甸NDVI的下降不僅是全區(qū)草地NDVI下降的主要原因，還是全區(qū)NDVI變化呈現(xiàn)正持續(xù)性的主要原因；溫性草甸草原和溫性草原H值均為0.70，其NDVI的變化也具有較強的正持續(xù)性；低平地草甸的H值最小，為0.64，其NDVI變化的正持續(xù)最弱；其他草地類型H值也均介于0.6～0.7之間。

表4伊犁河谷不同草地類型NDVI Hurst指數(shù)Table 4 NDVI Hurst index of different grassland types in Ili Valley

3 討論

3.1R/S分析結果檢驗

R/S分析提出時是用于水文時間序列變化趨勢的預測[19]，近年來該方法也逐步被應用于MODIS NDVI等遙感時間序列數(shù)據(jù)變化趨勢的預測[6-7]，水文數(shù)據(jù)時間序列長度均較長，而MODIS NDVI時間序列長度有限，僅有17年，對于短時間序列數(shù)據(jù)R/S分析結果是否可靠性，目前尚缺乏統(tǒng)一的檢驗標準或方法。Peng等[8]指出，根據(jù)R/S分析原理，當時間序列為隨機序列時，其H= 0.50[19]，基于此，對NDVI時間序列數(shù)據(jù)進行隨機打亂，然后計算其H值，可以根據(jù)隨機處理后NDVI時間序列的H值是否接近0.50來對R/S分析進行檢驗[8]。因此，利用伊犁河谷草地正常NDVI序列H值與隨機處理后NDVI序列H值空間分布數(shù)據(jù)，提取兩者所有像元點的H值，建立其散點圖(圖6)，對所有像元點NDVI序列隨機處理前后H值進行對比。

圖6 正常NDVI序列Hurst指數(shù)與隨機處理后NDVI序列Hurst指數(shù)散點圖Fig.6 Scatter plot of the Hurst index of the Normal NDVI sequence and randomly processed NDVI sequence

由圖6可知和表3可知，對于正常NDVI序列，全區(qū)97.77%像元的H>0.50，且絕大部分像元H值介于0.45～0.90之間。而對于隨機處理后的NDVI序列，全區(qū)很大一部分像元H<0.50，經(jīng)統(tǒng)計其比例達到了38.18%，遠高于正常NDVI序列的該比例，且絕大部分像元的H值介于0.35～0.75，H值域范圍整體有所減少。此外，由圖6還可以看出，全區(qū)絕大部分像元位于圖6中斜率為1的線的右側，表明絕大部分像元隨機處理之后其NDVI序列的H均有所減小，經(jīng)過統(tǒng)計，其比例達到了89.86%。可見，雖然隨機處理后的NDVI序列雖沒有轉變?yōu)闃藴孰S機序列(H=0.50)，但相對于正常NDVI序列，其隨機特性明顯增高，驗證了本文中R/S分析結果的可靠性。

3.2 NDVI變化原因

根據(jù)上文分析，伊犁河谷草地NDVI呈現(xiàn)下降的趨勢，結果與劉芳[17]、劉洋[18]及李凈[20]等人的研究結果一致。對于伊犁河谷草地植被NDVI的變化原因，利用伊犁河谷10個氣象站點的年平均氣溫數(shù)據(jù)和年累積降水數(shù)據(jù)(圖7)分別與NDVI數(shù)據(jù)進行相關分析，結果顯示NDVI與氣溫的相關系數(shù)為0.14(p=0.61)，與降水的相關系數(shù)為0.50(p=0.04)，表明NDVI年際變化與年累積降水的相關性要高于年平均氣溫，結果與李凈[20]和祝穩(wěn)[21]等人對于我國西北植被NDVI變化與氣溫和降水關系的分析結果一致。

圖7 2000-2016年伊犁河谷降水和氣溫年際變化Fig.7 Inter-annual variation of precipitation and temperature in Ili Valley from 2000-2016

對于植被變化的原因，已有大量學者利用反映植被生長狀況的植被指數(shù)數(shù)據(jù)，以及反映氣候條件的降水數(shù)據(jù)或氣溫數(shù)據(jù)等數(shù)據(jù)，對氣候變化與植被變化相關性進行了廣泛的研究[22]。如Kawabata等[23]的研究結果表明北半球氣候變暖使得植被活動增強，生長季延長[24]。而對于人類活動對植被變化的影響，人類活動是影響植被變化的關鍵因素的這一觀點已被廣泛接受[25]。例如，針對我國新疆、內蒙古以及青藏高原等區(qū)域的草地植被退化，大量研究證明過度放牧是其發(fā)生和不斷惡化的主要因素[26-28]。雖然人類活動是植被變化關鍵影響因素[29]，但由于受到數(shù)據(jù)源以及數(shù)據(jù)處理技術和方法等多方面因素的限制，在區(qū)域尺度上，對人類活動以及人類活動與植被變化關系的量化仍是有待深入研究的難題[30]。Evans和Geerken[31]通過預先建立降水量與植被NDVI之間的關系模型，之后在降水數(shù)據(jù)的支持下預測出植被的NDVI，該預測值被稱為氣候影響條件下的NDVI，而實際NDVI與預測NDVI之間的殘差被認為是由人類活動所引起，以此實現(xiàn)植被變化自然因素和人為因素的分離，該方法被稱為殘差分析法?；跉埐罘治龅幕舅枷?，Wang[32]、Jiang[29]及Sun[30]等對中亞、我國北方及西北榆林地區(qū)植被動態(tài)的氣候與人為影響進行了區(qū)分和量化，獲得了較好的結果。然而，殘差分析在建立氣候和植被關系模型的過程中所用植被數(shù)據(jù)往往是已包含了人類活動對植被的影響，使得植被變化的人為因素和氣候并不能被完全分離，成為殘差分析法所存在的一個重要缺陷。盡管如此，該法方的科學性仍具有較高的信服力，并已經(jīng)逐步成為分離植被變化自然與人為影響因素的重要方法，而被廣泛所采用[32]，這為伊犁河谷草地植被變化原因的深入研究指出了方向，提供了思路。

4 結論

2000-2016年伊犁河谷草地NDVI總體呈顯著下降趨勢，其平均變化率為-0.0038單位·a-1。各草地類型的NDVI均呈下降趨勢，但山地草甸及高寒草甸下降趨勢達到顯著水平，呈極顯著下降，下降速率分別為0.0047單位·a-1和0.0030單位·a-1。

2000-2016年伊犁河谷全區(qū)88.07%草地NDVI變化趨勢為下降，而49.62%達到顯著水平，空間上呈帶狀分布于河谷北、南及西部烏孫山等區(qū)域的中高山區(qū)；伊犁河出國境口北部部分區(qū)域草地的NDVI變化呈顯著增加趨勢，但面積非常有限。各草地類型中山地草甸分布最為廣泛，其面積的95.76%的NDVI出現(xiàn)下降，且73.10%達到顯著水平。

R/S分析表明，伊犁河谷草地97.77%草地NDVI的H> 0.50，且全區(qū)56.37%草地NDVI的H> 0.70，呈現(xiàn)較強的正持續(xù)性，未來一定時期伊犁河谷絕大部分草地的NDVI將繼續(xù)呈下降趨勢。各草地類型平均NDVI的H值均> 0.50，山地草甸H值最大，為0.74，是全區(qū)NDVI變化呈現(xiàn)正持續(xù)性的主要原因。

2000-2016年伊犁河谷草地年NDVI與年平均氣溫的相關系數(shù)為0.14(p=0.61)，與年累積降水的相關系數(shù)為0.50(p=0.04)，與氣溫相比，伊犁河谷降水的年際變化與草地NDVI年變化的相關性更高。