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

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

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

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

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

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

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

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

1 材料與方法

1.1 研究區(qū)概況

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

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

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

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

1.3 數(shù)據(jù)處理與分析

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

(1)

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

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

圖2 NDVI最大值合成流程圖

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

(2)

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

參考中華人民共和國水利部2008年發(fā)布的《土壤侵蝕分類分級標準》，結(jié)合本研究區(qū)實際情況和相關(guān)研究[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圖像差值法 為反映研究區(qū)植被覆蓋空間變化特征，利用圖像差值法計算不同影像之間的植被覆蓋變化ΔFg，差值大于0表示植被覆蓋增加，小于0表示植被覆蓋減少[25]。其公式為：

ΔFg=Fyears2-Fyear1

(3)

式中，F(xiàn)year1為前一時期植被覆蓋度等級，F(xiàn)year2為后一時期植被覆蓋度等級。

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

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

(4)

(5)

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

2 結(jié)果與分析

2.1 NDVI年內(nèi)合成最大值驗證和植被覆蓋精度檢驗

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

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

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

表1 阿拉爾墾區(qū)解譯標志

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

圖4 每類樣本點個數(shù)和分布圖

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

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

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

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

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

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

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

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

2.3 植被覆蓋空間變化特征

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

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

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

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

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

表4 1990—2019年阿拉爾墾區(qū)植被覆蓋度等級轉(zhuǎn)移矩陣

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

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

2.4 植被覆蓋重心轉(zhuǎn)移

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

圖9 1990-2019年阿拉爾墾區(qū)植被覆蓋重心轉(zhuǎn)移變化

3 討論

3.1 氣候變化影響

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

氣候變化對植被覆蓋的影響，還間接表現(xiàn)在對流域徑流的影響[33]。阿拉爾墾區(qū)年均氣溫在1990—2019年呈現(xiàn)出增加趨勢，特別是在2004年后，平均升高了1.09℃[34]。塔里木河水源主要來自天山山脈的冰雪融化，溫度的升高將導(dǎo)致流域徑流量的增加[35]，從新疆塔里木河流域管理局提供的徑流量數(shù)據(jù)顯示，2005年塔里木河年均徑流量比2004年增加了27.6×108m3，為墾區(qū)植被擴張?zhí)峁┝艘欢ǖ乃幢Ｕ?，有利于植被覆蓋面積的增加，同時也是墾區(qū)近30年平均植被覆蓋度增加的主要原因。

3.2 人類活動影響

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

表5 阿拉爾墾區(qū)不同時期土地利用結(jié)構(gòu)

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

4 結(jié)論

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