1990—2019年阿拉爾墾區植被覆蓋時空變化特征分析

宋 奇，高 琪，馬自強，王 楠，王明玥，彭 杰*

(1.塔里木大學植物科學學院，新疆 阿拉爾 843300；2.北京大學地球與空間科學學院，遙感與地理信息系統研究所，北京 100871；3.浙江大學環境與資源學院，浙江 杭州 310058)

植被作為全球陸地生態系統的重要組成部分，其通過光合作用、呼吸作用與生物圈其他自然要素形成緊密聯系，在生態系統的能量流動和物質循環中扮演著重要角色[1-2]。作為連接土壤、水、生物和大氣等要素的重要紐帶，植被不僅具有維持區域氣候穩定和生態平衡的功能，而且能夠反映區域綠洲和荒漠2大生態系統的相互作用過程，為區域生態研究提供重要信息[3-4]。此外，植被還具有防風固沙、涵養水源、調節區域小氣候、治理大氣污染和美化環境等眾多生態功能[5-6]。明確植被覆蓋時空變化規律，探明氣候變化和人類活動對植被覆蓋變化的影響，對區域環境質量評價和生態系統維護具有重要意義[7]。

目前，研究植被覆蓋變化的主要方法有地面調查法和遙感監測法[8]。傳統地面調查法需投入大量人力、物力和財力且不適合于大面積研究。隨著遙感技術的不斷發展，遙感監測法因其具有信息更新快、覆蓋范圍廣、獲取便捷、成本低、利用價值高等優勢，成為監測植被覆蓋變化的重要手段[9-10]。利用遙感手段分析植被覆蓋空間分布特征已得到廣泛運用[11]。在植被覆蓋變化分析中，歸一化植被指數(Normalized difference vegetation index,NDVI)是植被覆蓋度遙感監測方法中最常用的植被指數，綜合反映了區域植被覆蓋密度和植被生長狀況，是描述植被空間分布密度的最佳指標。此外，植被覆蓋度(Vegetation Fractional Coverage,VFC)也是對區域植被長勢進行量化的重要指標[12]，它是指單位面積內植被樹冠在地面形成的垂直投影面積占單位面積的百分比。通過像元二分模型[13]計算歸一化植被指數并對植被覆蓋度進行提取的方法被廣泛運用于植被覆蓋動態監測研究[14]。

多數植被覆蓋變化研究的處理方法是對單一的遙感影像進行NDVI提取展開相關研究[15-17]。然而，這樣的研究未能捕捉連續性信息，容易造成關鍵拐點的遺漏，而且由于云覆蓋、植被物候差異、耕地輪作休耕和病蟲害的影響，單一遙感影像無法反映當年植被覆蓋最全面的情景，難以準確分析研究區的植被覆蓋情況。因此，基于長時序多時相合成數據進行的植被覆蓋變化研究仍需進一步加強。

阿拉爾墾區地處我國天山山脈和塔克拉瑪干沙漠的過渡帶，屬于西北干旱區典型的人工綠洲。阿拉爾墾區不僅是新疆主要的棉花(Gossypiumspp.)及特色林果生產基地，更是防止塔克拉瑪干沙漠北移的重要生態屏障，在維持整個新疆地區生態系統平衡中發揮著重要作用。近些年，墾區人口不斷增加，機械化作業水平逐步提高，人們為滿足生活所需，盲目開墾土地資源，人工植被面積(農作物和經濟林等)不斷增加，自然植被面積(天然林地和草地等)不斷減少，生態危機日益嚴重。

本研究收集了阿拉爾墾區1990—2019年共405景Landsat影像，并對其進行年內NDVI最大值合成，利用像元二分模型計算逐年的植被覆蓋度，運用多年平均、圖像差值法和面積加權重心轉移模型分析墾區植被覆蓋的空間分布特征、空間變化特征和重心轉移情況，并探討了氣候變化和人類活動對墾區植被覆蓋變化的影響。旨在探明人工綠洲區域的植被覆蓋變化過程，為實現墾區水土資源合理利用、植被保護、生態安全和可持續發展提供理論支撐。

1 材料與方法

1.1 研究區概況

阿拉爾墾區(80°30′0″～81°58′0″ E，40°22′0″～40°57′0″ N)，地處天山南麓，塔克拉瑪干沙漠北緣，內含勝利、上游和多浪3大水庫。墾區面積為4 105.91 km2，主要土地利用類型為草地和林地，分別占墾區面積的25.67%和14.83%。墾區主體部分以塔里木河為中軸線，次級河流和灌溉渠道為紐帶，構成密集網狀水體廊道，形成了較為成熟的人工農田綠洲景觀。墾區氣候具有典型暖溫帶極端大陸性干旱荒漠氣候特征，溫暖干燥、光熱資源豐富，日照率58%～69%、年均日照時數2 556.3～2 991.8 h、年均蒸發量1 876.6～2 558.9 mm、年均降水量40.1～82.5 mm。墾區多受揚塵和沙暴等惡劣天氣影響。

1.2 數據來源與預處理

從美國地質勘探局USGS網站(http://giovis.usgs.gov)中獲取阿拉爾墾區1990—2019年間所有Landsat遙感影像(影像軌道行/列號：146/32，分辨率：30 m)，從中篩選出云覆蓋度低于100%的可用影像共計405景，圖1為各月份的遙感影像云覆蓋及每年的影像數量情況，其中1990—1998年選用Landsat 5 TM影像，1999—2012年選用Landsat 7 ETM+影像，2013—2019年選用Landsat 8 OLI影像。利用ENVI 5.3對影像逐一進行輻射定標、大氣校正、幾何校正、裁剪和圖像增強等預處理操作。1990—2019年阿拉爾墾區(區站號51730)氣象數據來自《地面氣象記錄月報表》，土地利用現狀遙感監測數據是我國精度最高的土地利用遙感監測產品，此類遙感監測數據來自中國科學院資源環境科學數據中心的共享平臺。

圖1 遙感影像的云覆蓋度及影像數量情況

1.3 數據處理與分析

1.3.1NDVI最大值合成 對篩選出的405景可用影像進行預處理后，逐景提取NDVI[18]，計算公式為：

(1)

其中，NIR近紅外波段反照率(Near infrared spectrum,NIR)，R為紅波段反照率(Spectrum,R)。

進而將每年各月份影像進行NDVI周年最大值合成，得到1990—2019年間30景NDVI合成后影像(圖2)。由圖2所示，合成后影像的植被覆蓋度更高，更能全面反映當年植被覆蓋最全面的情景，有利于植被覆蓋度的計算，從而得到更精細的解譯結果。

圖2 NDVI最大值合成流程圖

1.3.2植被覆蓋度計算 基于NDVI和植被覆蓋度之間存在的極顯著線性相關關系，運用像元二分法模型對植被覆蓋度進行計算[19]。其公式為：

(2)

式中，NDVIsoil為無植被覆蓋的裸土像元NDVI值，NDVIveg為全植被覆蓋像元NDVI值。取累計直方圖的 5%和 95%為置信度區間，確定研究區的NDVIsoil和NDVIveg值[20-21]。

參考中華人民共和國水利部2008年發布的《土壤侵蝕分類分級標準》，結合本研究區實際情況和相關研究[22-24]，將植被覆蓋度閾值按照等間距法劃分為5個等級：裸地或極低植被覆蓋(0～0.15)，等級為1；低植被覆蓋(0.15～0.3)，等級為2；中植被覆蓋(0.3～0.45)，等級為3；高植被覆蓋(0.45～0.6)，等級為4；極高植被覆蓋(0.6～1)，等級為5。

1.3.3圖像差值法 為反映研究區植被覆蓋空間變化特征，利用圖像差值法計算不同影像之間的植被覆蓋變化ΔFg，差值大于0表示植被覆蓋增加，小于0表示植被覆蓋減少[25]。其公式為：

ΔFg=Fyears2-Fyear1

(3)

式中，Fyear1為前一時期植被覆蓋度等級，Fyear2為后一時期植被覆蓋度等級。

為突出1990—2019年植被覆蓋空間變化情況，對1990年和2019年植被覆蓋度提取結果進行差值量化分析：其中，ΔFg=3為極度改善；ΔFg=2為中度改善；ΔFg=1為輕微改善；ΔFg=-3為極度退化；ΔFg=-2為中度退化；ΔFg=-1為輕微退化；ΔFg=0為穩定。

1.3.4植被覆蓋面積加權重心轉移分析 重心是地理學中描述地理要素或空間對象轉變的重要指標。重心的動態轉移變化反映了地理要素空間分布的轉移軌跡，加權重心則是通過對重心坐標賦予權重來表示地理現象分布的不均勻性，選用面積加權重心模型計算不同等級植被覆蓋度的重心，并反映研究時段內墾區植被空間演變過程[26]。其公式為：

(4)

(5)

式中，X，Y分別為墾區植被分布重心的經緯度坐標；Ci表示第i個墾區植被分布斑塊的面積；Xi，Yi分別表示第i個墾區植被分布斑塊分布重心的經緯度坐標。

2 結果與分析

2.1 NDVI年內合成最大值驗證和植被覆蓋精度檢驗

考慮到每個年份使用的遙感影像數量有限，有可能造成年內NDVI被低估，因此對本文的NDVI年內合成最大值進行驗證。將本文的NDVI年內最大值與中分辨率成像光譜儀(Moderate resolution imaging spectroradiometer,MODIS)的NDVI產品計算的年內最大值進行對比分析，確保本文計算的年內NDVI最大值未被低估。獲取了阿拉爾墾區2010年間所有MODIS的NDVI產品進行最大值合成，隨后與Landsat影像的NDVI最大值進行對比，如圖3所示，分別將兩產品的NDVI最大值合成后結果進行兩處細節放大，從中可以明顯看出Landsat產品的NDVI最大值合成后結果更清晰，NDVI值更高，效果最佳。并且Landsat和MODIS產品都是16 d一期，兩產品的空間分辨率不同，Landsat是30 m，MODIS是250 m，明顯Landsat影像的成像效果更佳，因此本文選擇Landsat影像進行NDVI最大值合成。

圖3 中分辨率成像光譜儀和遙感影像產品的歸一化植被指數最大值合成后比較

為驗證所得植被覆蓋情況的可靠性，根據阿拉爾墾區的土地利用/覆被現狀及相關參考文獻[27-28]，經綜合對比分析，建立了一套適合于阿拉爾墾區的解譯標志。以2010年為例，各植被覆蓋類型的地表特征、影像特征和實地調研情況如表1所示。

表1 阿拉爾墾區解譯標志

根據該解譯標志建立各植被覆蓋類型所對應的驗證樣本，樣本點個數及分布情況如圖4所示，各樣本點的分布遵循均勻、全方位覆蓋原則，最終各樣本所占的像素點個數分別為：極高覆蓋100 603個、高覆蓋60 387個、中覆蓋10 394個、低覆蓋53 570個、極低覆蓋240 571個，同理將其余年份逐一建立驗證樣本進行精度驗證。

圖4 每類樣本點個數和分布圖

進而對其建立混淆矩陣進行各植被覆蓋類型的精度評價。由表2可知，阿拉爾墾區的總體精度(Overall accuracy,OA)為88.50 %，Kappa系數為0.83。相關文獻表明[29-30]，總體精度OA值和Kappa系數均大于0.8時，所得的結果可用于相關研究，因此本研究的結果可信。

表2 阿拉爾墾區2010年植被覆蓋精度評價

3.2 植被覆蓋空間分布特征及變化趨勢

基于NDVI像元二分模型，計算出墾區1990—2019年不同等級的植被覆蓋度，以此為依據進行等級劃分，得到阿拉爾墾區1990—2019年的植被覆蓋等級分布情況(圖5)。由圖5可知，阿拉爾墾區的植被空間分布總體以塔里木河為軸線，從內向外擴展。極高和高植被覆蓋區主要集中于墾區中部，呈現出大片塊狀集中分布的趨勢；中植被覆蓋區主要集中于高植被覆蓋區外圍，植被類型多以人工防護林為主；低和極低植被覆蓋區集中于墾區南、北兩端，主要為稀疏的天然植被。

圖5 阿拉爾墾區不同時期植被覆蓋等級分布

統計得到的阿拉爾墾區平均植被覆蓋度變化和不同等級植被覆蓋面積變化情況(圖6)。1990—2019年，阿拉爾墾區平均植被覆蓋度呈波動式增加的趨勢，擬合優度R2為0.9042，平均植被覆蓋度在0.1和0.5之間變化，但存在年際差異：1991最低，約為0.1286；2018年達到峰值，約為0.4656。

極低植被覆蓋區面積呈線性減少趨勢明顯(R2=0.9639)，面積由1990年的3 050 km2減至2019年的1 226 km2，降幅59.8 %。低植被覆蓋區面積呈增加趨勢顯著(R2=0.9729)，增幅450%。30年間中植被覆蓋區面積的變化波動較大，整體呈減少趨勢(R2=0.1475)。高植被覆蓋區面積的變化呈明顯增加趨勢(R2=0.9858)，增幅162%。極高植被覆蓋區面積的年際間變化波動性減少(R2=0.9514)，面積變化呈增加趨勢，增幅1 880 %(圖6)。對比可知，在不同覆蓋等級的面積變化中，極低植被覆蓋區的面積減少最顯著，高植被覆蓋區的面積增加最顯著。

圖6 阿拉爾墾區平均植被覆蓋度和不同等級植被覆蓋面積變化

2.3 植被覆蓋空間變化特征

由于每年的植被覆蓋波動較大，因此用多年平均的方式來表征不同階段阿拉爾墾區的指標覆蓋變化和不同等級的植被覆蓋面積情況。在ArcGIS軟件中對阿拉爾墾區植被覆蓋進行多年平均，結果如圖7所示。

圖7 阿拉爾墾區植被覆蓋的多年平均變化情況

統計多年平均中不同等級的面積變化(表3)可知，1990—1995年，以極低覆蓋為主(2 431.12 km2)，主要分布在墾區南北兩端的偏遠地區；1995—2000年，極低覆蓋面積減少，面積為1 918.70 km2；2000—2005年，極低覆蓋面積在所有覆蓋等級中依然為最大值(1 694.10 km2)；2005—2010年，墾區西北部開始出現植被，以中覆蓋類型為主；2010—2015年，極高覆蓋面積高達772.73 km2；2015—2019年，極高覆蓋面積成為所有覆蓋等級中的最大值(1 113.94 km2)；1990—2019年，整體來講，墾區還是以極低覆蓋為主，而其他密度較高的覆蓋類型主要分布在塔里木河區域。

表3 阿拉爾墾區植被覆蓋多年平均中不同等級的面積變化情況 單位：km2

綜上可知，阿拉爾墾區植被變化不僅存在時間階段性差異，同時存在區域性差異。為進一步分析阿拉爾墾區植被覆蓋面積的變化情況，運用轉移矩陣對1990年和2019年不同等級的植被覆蓋面積進行計算。1990—2019年，阿拉爾墾區植被覆蓋主要以極低覆蓋向高覆蓋轉移和高覆蓋向極高覆蓋轉移為主，具體轉移量為：極低植被覆蓋轉向高植被覆蓋的面積為1 119.48 km2，高植被覆蓋轉向極高植被覆蓋的面積為792.49 km2。從最終變化量來看，極低植被覆蓋的面積減少量最大，減少了1 861.54 km2，極高和高植被覆蓋面積的增加量較大，分別增加了1 541.54 km2和1 376.37 km2(表4)，墾區的植被覆蓋情況整體上表現出極低植被覆蓋區的面積逐漸減少，高和極高植被覆蓋區的面積逐漸增加趨勢。

表4 1990—2019年阿拉爾墾區植被覆蓋度等級轉移矩陣

通過差值量化提取的阿拉爾墾區植被覆蓋面積變化看出(圖8)，阿拉爾墾區中心區域周圍成為墾區植被覆蓋面積退化的熱點區域。同時墾區中心區域因經濟相對發達，周邊鄉鎮和公路干線集中，所以植被覆蓋面積的退化程度更為顯著。相對而言，離塔里木河流域較遠的區域成為墾區植被覆蓋面積改善的熱點區域。

圖8 1990-2019年阿拉爾墾區植被覆蓋變化等級量化分布

2.4 植被覆蓋重心轉移

重心分布變化能夠從空間上反映墾區植被分布的時空演變特征[31]。阿拉爾墾區植被覆蓋重心變化結果顯示(圖9)，近30年，阿拉爾墾區植被覆蓋重心不斷朝東北方向轉移，但在2005年后開始往相反的西南方向來回移動。具體來說，1990—2005年向東北方向轉移；2005—2010年向西南方向轉移；2010—2015再次向東北方向轉移；2015—2019再次向西南方向轉移。在各時段中，1990—1995年重心轉移最明顯，直線轉移距離為4.12 km。這主要是因為1994年在阿拉爾墾區的東北區域大面積開墾耕地，致使墾區植被覆蓋重心向東北方向移動，而在2006年國家在墾區西南區域大面積增設團場，使得墾區西南方向的覆被面積迅速增加，重心向西南方向轉移，在2005年后墾區重心在東北和西南方向之間來回轉移。

圖9 1990-2019年阿拉爾墾區植被覆蓋重心轉移變化

3 討論

3.1 氣候變化影響

阿拉爾墾區人為種植作物多以高和極高植被覆蓋的形式呈現，天然疏林灌叢和草地多以低植被覆蓋為主，而自然氣候的變化主要影響天然植被的生長[32]。因此，阿拉爾墾區年均降水量的增加對墾區天然疏林灌叢和草地產生了積極作用，這也正是墾區低植被覆蓋向中植被覆蓋轉移的原因。

氣候變化對植被覆蓋的影響，還間接表現在對流域徑流的影響[33]。阿拉爾墾區年均氣溫在1990—2019年呈現出增加趨勢，特別是在2004年后，平均升高了1.09℃[34]。塔里木河水源主要來自天山山脈的冰雪融化，溫度的升高將導致流域徑流量的增加[35]，從新疆塔里木河流域管理局提供的徑流量數據顯示，2005年塔里木河年均徑流量比2004年增加了27.6×108m3，為墾區植被擴張提供了一定的水源保障，有利于植被覆蓋面積的增加，同時也是墾區近30年平均植被覆蓋度增加的主要原因。

3.2 人類活動影響

人類活動對生態環境造成的影響與土地利用有關[36]。從中國科學院資源環境科學數據中心獲取了阿拉爾墾區1990，2000，2010和2018年的土地利用數據(表5)。1990—2018年阿拉爾墾區耕地和城鄉工礦居民用地的面積增加幅度最大，其中耕地面積從1990年的61.59 km2(1.5%)增至2018年的1 180.04 km2(28.74 %)，凈增加1 151.3 km2，說明阿拉爾墾區植被覆蓋面積的增加主要得益于耕地的擴張，擴張最明顯的區域位于墾區西北部(圖8)。從表5中還可看出，墾區未利用土地面積減少。據1990—2018年阿拉爾墾區土地利用變化轉移矩陣計算，阿拉爾墾區1 750 km2的自然草地和2 667 km2的未利用土地被轉移成為耕地。

表5 阿拉爾墾區不同時期土地利用結構

人類活動除通過開墾農田、種植作物影響墾區植被覆蓋面積變化外，城鄉基建對墾區植被覆蓋面積變化的影響也十分顯著[37-38]。耕地和城鄉基建的擴張促進了墾區經濟發展，也帶來了巨大的用水壓力，造成墾區生態用水不足，在墾區人工植被面積不斷增加的同時天然植被面積不斷減少。這種以犧牲天然植被(草地)換取人工植被(耕地)的方式，破壞了墾區的正常生態平衡。干旱區天然植被作為維持墾區生態穩定的重要屏障，其面積的減少勢必加快了墾區荒漠化進程，威脅到墾區的長期發展。同時，塔里木河的過度引流造成下游輸水減少，嚴重威脅下游的生態安全和可持續發展。

4 結論

本研究對阿拉爾墾區植被覆蓋時空變化特征進行分析，主要結論為：1990—2019年，阿拉爾墾區平均植被覆蓋度呈增加趨勢。極低植被覆蓋區的面積呈明顯的減少趨勢。植被覆蓋空間變化存在時序性和區域性差異。整體上，墾區植被覆蓋重心整體向東北轉移。氣候變化對阿拉爾墾區植被覆蓋有一定的促進作用，但人類活動的影響最為直接。大面積耕地開墾是墾區植被覆蓋面積增加的主要原因，這種以犧牲自然草地換取耕地的方式破壞了墾區的生態穩定和可持續發展。