基于TINDVI的伊犁河谷草地覆蓋變化

閆俊杰, 陳 晨, 趙 陽, 沙吾麗·達吾提拜, 劉海軍

(伊犁師范大學 資源與生態研究所, 新疆 伊寧 835000)

草地占據了全球將近20%的的陸表面積[1]，它不僅是陸地生態系統重要組成部分[2]，而且在生物多樣性保護、水源涵養、大氣調節及原材料供給等生態系統維護與調節方面占有重要地位[3]。全球變化背景下，草地植被動態變化及其對生態環境的影響被廣泛關注[4-7]，而對草地植被生長科學而精確的量化表達是研究草地動態變化的基礎。

遙感技術在數據獲取的時效性、成本及數據豐富度等方面具有難以比擬的優勢[8-9]，是目前區域植被動態監測的主要方法及數據來源。基于遙感技術，不同學者將植被反射光譜中的不同波段進行線性或非線性組合來突出地表植被信息，提出了多種植被指數來對植被的生長狀況進行量化表達[10]，為植被動態變化研究提供了有效且計算簡便的指標。在眾多植被指數中，歸一化植被指數(Normalized difference vegetation index，NDVI)應用最為廣泛，被廣泛應用于全球、區域及流域等不同尺度的植被動態研究[4-5,9,11-14]。MODIS NDVI，NOAA-AVHRR NDVI及SPOT NDVI時間序列數據是主要數據源[15]。在對NDVI時間序列數據預處理過程中，最大值合成法(maximum Value Composites，MVC)常被用于計算年內NDVI最大值(NDVImax)，來表征單個生長季內植被生長達到最好時的狀況[10,16]，從而為植被的年際變化分析提供基礎。然而植被的年際變化并不僅僅體現在不同年份植被生長可達到的最高水平的差異上，植被在不同年份的年內生長過程的差異也是植被動態變化的重要方面。時間累積NDVI(Time integered NDVI，TINDVI)被定義為年內NDVI變化曲線與時間軸圍城的面積或一定時期內NDVI的累積值[17-19]。相對于NDVImax，TINDVI包含了植被生長過程的信息，且已有研究證明TINDVI與植被干物質量具有很好線性關系[20-21]，然而以TINDVI為衡量指標來分析植被變化的研究卻鮮有報道。對比分析TNDVI與NDVImax在反映植被時空變化方面的比較優勢，探討TINDVI的應用價值，將為全面和深入認識植被的動態變化提供有利幫助。

西北地區是我國草地資源的重要分布區，而位于天山西段的伊犁河谷更是我國久負盛譽的優質牧場。已有學者以NDVImax或覆蓋度為參考，指明伊犁河谷草地退化嚴重，且退化仍在持續惡化[22-25]。本研究利用包含生長過程信息的TINDVI作為表征草地植被的參數指標，分析伊犁河谷草地植被TINDVI時空變化，并探討TINDVI與NDVImax的差異，以期進一步探明伊犁河谷草地植被的變化特征，并為相關研究提供借鑒。

1 研究區概況

伊犁河谷位于80°09′42″—84°56′50″E，42°14′16″—44°53′30″N，地處天山西段，整個河谷地勢西低東高，地形呈向西敞開的“V”字形(圖1)，西風濕潤氣流受此特殊地勢地形影響，容易抬升凝結，給河谷帶來豐沛降水，成就了其“塞外江南”的美譽。伊犁河谷氣候雖屬于溫帶大陸性氣候，但卻表現出明顯的高山氣候特征[26]，平原區與山區氣候明顯不同，平原區干旱少雨、蒸發能力強，而山區降水豐沛、天氣冷涼[27]，平原區降水200 mm左右，但山區降水可達1 000 mm，氣溫變幅約為9.2～2.8 ℃[26]。河谷內植被類型以草地主，且類型豐富多樣，是新疆乃至我國著名的優質牧場，分布有高寒草甸、山地草甸、溫性草甸草原、溫性草原、溫性荒漠草原、溫性荒漠及低平地草甸等多種類型草地[24,28]。伊犁河谷行政區劃上屬于伊犁哈薩克自治州，畜牧業在國民經濟中占有重要比例，且畜牧方式仍然以傳統放牧為主。相關研究表明，近年來受氣候變化及人為影響，草地植被覆蓋度持續降低[29]，草地生產力降低及毒害草蔓延等草地退化問題嚴重[22,30]。

圖1 研究區示意圖

2 材料與方法

2.1 數據來源與預處理

文章用到的NDVI數據為美國NASA EOS數據中心發布的MODIS MOD13Q1產品。MOD13Q1數據為16 d合成的時間序列數據，空間分辨率為250 m，每年23期，本研究用到2000年2月—2018年12月共19 a數據。在數據預處理過程中，對NDVI時間序列數據進行了鑲嵌、轉投影、重采樣、研究區裁剪以及Savitzky-Golay濾波處理。用于年際變化分析的NDVI年數據通過對每年23期數據的最大值合成(MVC)處理獲得[10,16]，即用年內植被生長達到最好狀況時的NDVI(NDVImax)作為年尺度上的NDVI。

草地分布邊界的矢量數據是通過對2018年6—7月份sentinel-2假彩色合成影像目視解疑獲得。對解疑獲得的矢量數據進行了柵格化處理，作為NDVI數據研究區裁剪的掩膜。為保障數據的空間匹配，草地分布數據的柵格化處理及NDVI數據重采樣處理過程中，像元大小均設置為250 m×250 m。

氣溫和降水數據來自中國氣象科學數據共享服務網(http:∥data.cma.gov.cn/)，包括伊犁河谷內伊寧市、霍爾果斯市、霍城等10個氣象站和其周邊巴音布魯克在內的11個氣象站的月數據。在11個氣象站點數據的基礎上，借助于Anusplin4.2軟件及DEM數據，進行空間插值，獲得伊犁河谷月氣溫和月降水空間分布數據，之后合成生長季(3—10月份)內的平均氣溫和累積降水空間分布數據，插值過程中柵格大小設置為250 m，與NDVI數據柵格大小一致。

2.2 研究方法

2.2.1 TINDVI計算 TINDVI通常被定義為生長季內NDVI變化曲線與時間軸圍城的面積或NDVI的累積值[17-19]。累積值計算方便且簡單，因此TINDVI利用生長季內NDVI的累積值表示，即：

(1)

式中:i即代表NDVI數據的日期。MOD13Q1 NDVI產品每年23期數據,從每年的第1天到第353天每16 d一期數據。根據研究區草地植被返青期和枯萎日期實際情況及其空間差異[31];TINDVI計算取3月初(第65天)到11月初(第305天)時段內共16期NDVI數據,并在計算過程中對3月初和11月初NDVI<0數據進行了剔除。

2.2.2 單調趨勢Mann-Kendall非參數檢驗 Mann-Kendall非參數檢驗被廣泛應用于時間序列數據趨勢檢驗，其統計量Zc計算過程如下[32-33]：

(2)

(3)

(4)

(5)

式中:xi為xk為時間序列數據的集合。在α置信水平上,若|Zc|≥Z1-α/2,即時間序列在α置信水平上存在顯著上升或下降趨勢。Zc>0,則變化趨勢為上升趨勢;Zc<0,則變化趨勢為下降趨勢。

時間序列變化速率可用Kendall傾斜度來表示，其計算公式如下：

(6)

式中:1

3 結果與分析

3.1 TINDVI空間分布特征

利用2000—2018年伊犁河谷草地TINDVI的年空間分布數據，計算其多年平均值。基于TINDVI的多年平均數據分析伊犁河谷草地TINDVI的空間分布特征(圖2)。

圖2 草地TINDVI空間分布(A)及其直方圖(B)

由圖2可知，伊犁河谷草地多年平均TINDVI的值域約為0～11.0。空間上伊犁河谷草地TINDVI呈現明顯的垂直分異特征。具體表現為TINDVI先隨海拔的升高而逐步增加，河谷西部海拔較低的平原及山前洪積沖積扇區TINDVI較低(2.08.0)，但空間上主要集中在河谷中南部及東部的中山區域，之后隨海拔升高TINDVI變化轉為逐步減小，河谷周邊積雪時間較長的高山區TINDVI最小值(TINDVI<2.0)。

經統計，整個河谷6.08.0的比例達到了21.67%；TINDVI<2.0的比例最少，為9.34%；2.0

3.2 2000-2018年伊犁河谷草地TINDVI時間動態

對2000—2018年伊犁河谷草地TINDVI的年空間分布數據求逐年平均，繪制TINDVI的年際變化曲線(圖3)，分析伊犁河谷草地TINDVI總體的年際變化。

圖3 2000-2018年TINDVI年際變化曲線

由圖3可知，2000—2018年伊犁河谷草地TINDVI總體呈波動增加的變化趨勢，但根據Mann—Kendall趨勢檢驗，其Zc值為1.05，未達到顯著水平，其年變化率β值為0.015。在年際波動上，TINDVI于2008年和2014年出現較大谷值，值分別為5.34,5.36；于2007年、2013年和2016年出現較大峰值，值分別為6.26，6.40,6.39；2008年谷值前后的2000—2007年及20008—2013年分別出現兩個明顯的波動增加時段。

3.3 2000-2018年伊犁河谷草地TINDVI空間動態

基于2000—2018年伊犁河谷草地TINDVI的逐年空間分布數據，根據Mann-Kendall檢驗的計算步驟，進行逐像元計算，獲得TINDVI變化趨勢Zc值和變化率β值的空間分布數據(圖4)，基于此分析TINDVI變化的空間特征。

注：1.96為Mann—Kendall檢驗p=0.05顯著水平上Zc臨界值。

由圖4可知，伊犁河谷絕大部分草地TINDVI呈增加趨勢，經統計，其面積比例達到了73.95%，但其中僅有11.56%達到了顯著水平，在空間上離散分布在河谷南部、東部及中西部的部分區域；整個區域有26.05%的區域TINDVI呈降低變化趨勢，空間上主要分布在昭蘇盆地、特克斯河下游周圍及河谷北部的中山區域，而TINDVI降低達到顯著水平的比例僅為0.67%。

對于TINDVI的變化速率，全區大部分區域的變化率為-0.025～0.050(圖4B)，其比例達到了85.91%，其中TINDVI變化率為0～0.025的分布比較廣泛，比例達到了41.29%，TINDVI變化率為-0.025～0的比例較小，為19.95%，空間上主要離散分布于河谷中部及北部的中山和低山區域，TINDVI變化率為0.025～0.050的比例為24.67%，空間上離散分布于河谷的不同區域；TINDVI變化率<-0.025和>0.050的比例均小，比例分別為5.99%和8.10%，分布也均較為離散。

3.4 TINDVI與NDVImax年際變化的時空差異

相關研究表明，近年來伊犁河谷草地的NDVImax及覆蓋度變化總體均呈現減小趨勢[24-25,29]，但根據本文的分析，伊犁河谷草地TINDVI卻總體有所增大，空間上73.95%草地的TINDVI呈現不同程度增加。鑒于本文研究與前人研究得出相反的結論，有必要對TINDVI及NDVImax變化的差異性進行探討。

對比TINDVI與NDVImax年際變化曲線(圖5A)，2000—2018年NDVImax變化總體呈減小趨勢，據Mann-Kendall趨勢檢驗，其Zc值為-0.28，未達到顯著水平，變化速率為-0.000 1；而TINDVI總體呈增加趨勢，Zc值為1.05，也未達到顯著水平。在年際波動上，TINDVI以及NDVImax均在2008年及2014年出現較大谷值，在2011年出現較大峰值，但在2005—2006年、2009—2010年及2016—2017年，NDVImax呈現增大變化，而TINDVI則呈現減小變化。

變化的空間分異方面，利用TINDVI變化趨勢Zc值圖與NDVImax變化趨勢Zc圖進行疊加，對比TINDVI及NDVImax變化趨勢的空間差異(圖5B)。由圖5B可知，伊犁河谷57.28%的草地NDVImax有所減小，42.58%草地NDVImax和TINDVI的變化趨勢相反，其中NDVImax減小而TINDVI增加的比例達到了36.91%，空間上主要位于NDVImax及TINDVI值較高的中山區域，而NDVImax和TINDVI呈相同變化趨勢的面積比例雖然也達到了57.42%，但在空間上主要分布在NDVImax及TINDVI值相對較低的山前洪積沖積扇、低山丘陵及積雪時間較長的高山區域。

圖5 伊犁河谷草地TINDVI與NDVImax年際變化曲線(A)及變化趨勢的疊加圖(B)

4 討 論

4.1 TINDVI與NDVImax年際變化時空差異的原因

由上文分析可知，2005—2006年、2009—2010年以及2016—2017年，全區平均TINDVI和NDVImax出現相反的波動變化(圖5A)，同時不同區域內，TINDVI和NDVImax的年際變化也不盡相同(圖5B)。NDVImax代表了生長季內植被覆蓋到達最大時的NDVI值，而TINDVI代表生長季內NDVI的累積值。生長季內NDVI變化曲線即能直接展示NDVImax值的大小，也能間接展示TINDVI值的大小，因此，對比不同年份及不同區生長季內NDVI變化曲線，能為解答NDVImax及TINDVI時空變化產生差異的原因提供答案。

以16 d為間隔，計算生長季內(每年第65—305天)TINDVI與NDVImax變化趨勢相反年份(2005—2006年、2009—2010年以及2016—2017年)，以及研究期前期(2000—2005年)和末期(2014—2018年)TINDVI與NDVImax變化趨勢不同疊加類型區的區域NDVI平均值，制作生長季內TINDVI與NDVImax變化趨勢相反年份(圖6A)及研究期前期和末期兩者變化趨勢不同疊加類型(圖6B)的NDVI變化曲線，探討NDVImax及TINDVI時空變化產生差異的原因。

對于在2005—2006年、2009—2010年及2016—2017年NDVImax與TINDVI呈現截然相反的變化，由圖6A可知，2005年草地NDVI在達到最大值之前春夏季的第65—161天和之后夏秋季的第209—289天內其值均小于2006年同時間的值，致使NDVImax雖有所增加，但TINDVI卻發生減小；2009—2010年主要是在春夏季的第65—161天內2009年NDVI值小于2010年；2016—2017年則主要是在夏秋季的第193—257天內2016年NDVI值小于2017年。

圖6 TINDVI與NDVImax變化趨勢相反年份及兩者變化趨勢不同疊加類型在研究期前期和末期NDVI生長季內的變化

而對于NDVImax及TINDVI變化趨勢在空間上的差異，圖6B分別展示了NDVImax及TINDVI變化趨勢不同交互類型(圖5B)中2000—2018年的前5 a(2000—2004年)及末端5 a(2014—2018年)的NDVI平均年內變化曲線。由圖可知，對于NDVImax及TINDVI均呈減小趨勢的區域，主要是由于NDVImax減小的同時，在夏秋季的絕大部分時間內(第129—289天)NDVI均有所減小(圖6B1)；對于NDVImax減小而TINDVI增大的區域，主要是由于在春夏季的第65—145天內及秋季的第257—305天的NDVI均有所增大(圖6B2)；對于NDVImax增大TINDVI減小的區域以及NDVImax增大TINDVI增大的區域，分別是由于夏秋季的第209—305天的NDVI均有所減小(圖6B3)或整個生長季NDVI幾乎都有所增大(圖6B4)。

通過上述討論和分析可以看出，生長季內NDVI隨植被從返青到枯萎的變化過程也隨區域的不同和年份的不同而存在明顯差異，單利用植被生長達到最大時的NDVImax作為基礎指標來評估草地的年際變化明顯存在一定局限，用TINDVI作為參考指標來衡量草地植被變化具有一定優勢。然而，雖然相關研究證實了TINDVI與玉米和大麥等農作的干物質量具有線性關系[20-21]，且TINDVI也常被作為植被生物量的相似或替代指標[34-35]，但產草量及凈初等生產力是衡量草地光合作用所積累有機質或干物質的關鍵指標，也是評價草地退化與否的重要參考[4-5,28]，而TINDVI與其關系如何仍需深入探究。

4.2 TINDVI及NDVImax與降水和氣溫關系及其對比

氣候變化是驅動植被變化的關鍵因素。分別計算TINDVI及NDVImax與降水和氣溫的偏相關系數，探討TINDVI及NDVImax與降水和氣溫相關性及其差異(圖7—8)。

由圖7可知，無論是對于TINDVI還是NDVImax，降水與草地植被變化的相關性明顯高于氣溫，且與降水以正相關為主，這與前人關于天山及中亞區域的研究結果一致[13-15]。而由圖8可知，2000—2018年生長季內伊犁河谷草地降水以減少為主，而全區氣溫均有所升高，氣溫升高促使水分蒸發增強，加劇了水分散失，而降水的減少又減少了水分來源，加劇水分對植被生長的脅迫，這可能是大面積草地NDVImax減小的原因。但在水分脅迫加劇的條件下，TINDVI卻發生增加，而由上文分析可知，TINDVI增加主要得益于春季(3—4月份)植被NDVI較大幅度的增加(圖6B2)，經分析，2000—2018年3—4月平均氣溫和累積降水均有所增加(氣溫和降水Mann—Kendall檢驗Zc值分別為0.28及0.98)，為春季(3—4月份)NDVI的增加提供了有利條件。

注：UU為與降水及氣溫相關性均非顯著；PU為與降水顯著正相關；UP為與氣溫顯著正相關；NU為與降水顯著負相關；UN為與氣溫顯著負相關；PP為與降水及氣溫均顯著正相關；NN為與降水和氣溫均顯著負相關；PN為與降水顯著正相關，與氣溫顯著負相關；NP為與降水顯著負相關，與氣溫顯著正相關。

氣候是驅動植被變化的關鍵因素，而海拔又使氣候對植被的影響產生明顯的海拔分異，因而植被與氣候變化的相關性也常存在明顯的海拔分異[4-5]。根據圖7A，在伊犁河谷中部低海拔區和中海拔區，植被與降水呈正相關而與氣溫相關性較弱，而在河谷周邊的高海拔區，植被氣溫呈正相關而與降水呈負相關。伊犁河谷地處亞歐大陸干旱與半干旱區的腹地，水分是植被生長的關鍵限值因子，尤其在海拔較低的區域，因而低海拔區植被與降水變化呈正相關[14]；而高海拔區，降水雖然相對豐沛，但冬季漫長，全年氣溫明顯低于平原區，氣溫則是限制植被生長的關鍵因子，因而植被變化氣溫呈正相關，與降水呈負相關[14,36]。植被與氣候相關性的海拔分異不僅在伊犁河谷有所表現，在整個天山、祁連山和青藏高原等區域均有所表現[14,4-5,37]。

對比TINDVI及NDVImax與降水和氣溫的相關性(圖7)，經計算，TINDVI與降水和氣溫呈顯著相關的比例分別達到了38.51%和12.22%，而NDVImax的該比例則分別為22.33%和5.27%，TINDVI與降水和氣溫變化相關性明顯高于NDVImax，尤其是與氣溫的相關性，可見TINDVI對氣候變化響應更為敏感；同時，在空間上，相對與NDVImax，TINDVI與降水和氣溫相關性的海拔分異規律也更為明顯，這也間接表明TINDVI對氣候變化響應更為敏感。相對于植被生長達到覆蓋最高時的NDVImax，TINDVI包含了整個生長季內植被生長過程中NDVI變化的信息，而累積降水和平均氣溫也包含有整個生長季內各自信息，這應該是TINDVI對氣溫和降水變化響應更敏感的原因。

5 結 論

(1) 根據多年平均數據，NDVI累加計算得到的伊犁河谷草地TINDVI的取值范圍約在0～11；空間上TINDVI海拔分異明顯，TINDVI的高值(>8.0)區域主要分布在中山區域，河谷平原區及冰雪覆蓋時間較長的高山區為TINDVI低值(<4.0)區；全區75.69%的區域TINDVI>4.0。

(2) 2000—2018年伊犁河谷草地TINDVI總體呈增大趨勢，但未達到顯著水平，其年變化速率為0.015。空間上，全區73.95%的區域TINDVI有所增大，其中11.56%達到了顯著水平，而TINDVI變化呈減小趨勢的面積比例僅為0.67%。

(3) 2000—2018年伊犁河谷草地NDVImax與總體呈非顯著減小趨勢，其年變化速率為0.000 1，空間上57.28%的NDVImax發生減小，42.58%的區域NDVImax與TINDVI變化趨勢相反，而生長季內草地生長過程所存在的年際變化是兩者變化趨勢相反的原因。

(5) 相對與NDVImax，伊犁河谷草地TINDVI對降水和氣溫變化響應更敏感。河谷低海拔和中海拔區草地覆蓋與降水呈顯著正相關，而高海拔區與氣溫呈顯著負相關，降水與草地覆蓋變化的相關性明顯高于氣溫。

(4) 伊犁河谷草地植被TINDVI變化的時間和空間差異是因為生長季內草地NDVI變化過程在不同年份及不同區域存在差異，利用TINDVI作為衡量草地植被變化的參考指標具有一定優勢。