基于GEE的諾木洪洪積扇植被時空變化特征、成因及趨勢分析

姚金璽， 張志， 張焜

(1.中國地質大學(武漢)地球物理與空間信息學院，武漢 430074； 2.青海省青藏高原北部地質過程與礦產資源重點實驗室，西寧 810300)

0 引言

植被是陸地生態系統中的重要組分。理解植被對氣候變化的響應對認識與預測未來生態系統的演化具有重要意義[1]，是干旱、半干旱區生態環境監測的內容[2]。近幾十年來，受到全球氣候變暖的影響，青海省諾木洪地區氣候向暖濕化發展，徑流有所增加。盡管對諾木洪河的管理得到加強，但由于農業用水的增加，對下游生態用水、植被發育產生一定的影響。

利用時間序列植被指數監測陸地植被生長變化規律并探討氣候要素的驅動作用已成為植被-氣候相互作用研究中的重要方向之一[3-4]。氣候變化和人類活動可能導致植被出現新的變化，需要對其重新認識[5]。基于傳統數據平臺和低空間分辨率影像分析大區域植被變化和驅動因素已開展較多研究[6-9]。基于遙感和地理信息系統(geographic information system,GIS)技術,徐浩杰等[6]和王林林等[7]利用生長季MODIS數據分析植被季節時空變化及其驅動因子，發現柴達木盆地植被覆蓋呈逐步改善趨勢; 張斯琦等[8]和李艷麗等[9]利用MODIS數據分析植被覆蓋度及人類活動與植被演化、徑流改變之間的關系。但前人研究多利用低空間分辨率遙感數據監測大區域植被宏觀變化，并采用傳統方法下載、存儲與處理遙感影像，研究成果精度與效率低。

Google Earth Engine(GEE)平臺中數據包括近40 a的全球衛星影像，且提供超過800 種功能函數[10]。基于云計算的 GEE 平臺可在線高效處理大規模地理數據集，有效地解決了遙感大數據中的處理難題[11-12]。

基于GEE平臺，利用2000—2017年MODIS增強植被指數(enhanced vegetation index，EVI)數據和Landsat歸一化差異植被指數(normal difference vegetation index，NDVI)數據，分析長時間序列植被變化特征，重點研究不同年代諾木洪洪積扇上枸杞種植區和鹽堿化區植被變化特征以及氣候對植被變化的影響。在此基礎上分析人為因素在植被時序演化中所起到的作用，并利用重標極差分析方法研究諾木洪洪積扇植被變化趨勢的可延續性，以期為諾木洪河水資源綜合管理和可持續發展決策提供參考。

1 研究區概況及數據源

研究區主體位于青海省海西蒙古族藏族自治州都蘭縣諾木洪地區，地理坐標為N36°10′～36°35′，E96°15′～96°40′，海拔為2 767～3 191 m(圖1)，屬高原干旱大陸性氣候，干燥少雨，年均氣溫為1.2～4.3 ℃，降水量為17.8～177.5 mm，集中在6—9月，表現出雨熱同季[13]。區內主要植被類型為梭梭、檉柳和蘆葦等植物，覆蓋度不高，在晚期洪積扇扇緣及河流沿岸等含水量高的地區，植被長勢良好。

圖1 諾木洪洪積扇位置示意圖

本文MODIS，Landsat5與Landsat8衛星數據均經過了輻射定標和大氣校正，使用的EVI數據是250 m空間分辨率的MODIS衛星16 d合成的EVI產品MOD13Q1[14]，NDVI數據是基于Landsat5與Landsat8影像計算所得。通過GEE算法去云得到高質量的影像數據后還需進行空間和時間雙維度的數據篩選，空間篩選利用clip或clipToCollection函數，時間跨度為2000—2017年。預處理過程和數值計算等操作在GEE平臺上由JavaScript語言編程實現。

氣象數據為諾木洪氣象觀測站的年均溫度及總降水量，時間跨度為2000—2017年，其來源于中國氣象科學數據共享服務網(http: //cdc.cma.gov.cn)。

2 研究方法

研究主要內容包括： 遙感影像的收集與預處理，獲取降水量和溫度數據； 計算植被覆蓋特征因子和自然環境因子； 變化特征的分析及結果驗證； 植被變化影響因素分析和未來變化趨勢分析4個部分。技術路線如圖2所示。

圖2 技術路線

2.1 最大化合成與最大值計算法

對2000—2017年研究區Landsat衛星數據進行計算NDVI值，并篩選每個像素的NDVI最大值，予以合成，針對MODIS衛星EVI產品進行篩選最大的EVI值。利用2組數據以及不同的計算方法來表達整個區域當年植被生長最繁茂的時期，分析該時間段內月際、年際植被變化特征。

2.2 Pearson相關分析與一階偏相關分析法

對溫度、降水量和各個時間點的最大EVI值數據做相關分析與偏相關分析，由此解釋溫度、降水量與最大化EVI數據的相關程度[15]。2個變量之間Pearson相關系數ρXY定義為這2個變量的協方差與二者標準差積的商，即

(1)

式中：X和Y為2個變量；cov(X,Y)為X與Y的協方差；σXσY為X與Y的標準差之積；E(X)為X的數學期望。

當相關系數處于[0.8,1]時，二者高度相關； 處于[0.5,0.8)時，二者中度相關； 處于[0.3,0.5)時，二者低度相關； 小于0.3時則為弱度相關。

在3個變量中，任意2個變量的偏相關系數是在排除其余一個變量影響后計算得到的，稱為一階偏相關系數，公式為：

(2)

式中：rij·h為排除變量Xh影響后的Xi和Xj的偏相關系數；rij為變量Xi與Xj的簡單相關系數；rih為變量Xi與Xh的簡單相關系數；rjh為變量Xj與Xh的簡單相關系數。

偏相關系數檢驗的零假設為： 總體中2個變量間的偏相關系數為0。使用t檢驗方法，公式為：

(3)

式中：r為相應的偏相關系數；n為樣本觀測數；k為可控制變量的數目；n-k-2為自由度。當t>t0.05(n-k-2)或p<0.05時，拒絕原假設。

2.3 Hurst指數

Hurst指數常用于定量描述時間序列變化趨勢的可延續性[16]，可以更好地分析植被的年際變化特征。Hurst指數由英國水文學家 Hurst[17]提出，現已應用于地質、遙感和水文等領域中。

本文中利用重標極差分析方法計算Hurst指數，分析研究區域內未來短期的植被變化趨勢。Hurst指數值H總體處于[0,1]區間之內，將此區間進行劃分代表著不同的含義。當0≤H<0.5時，表明時間序列具有長期相關性，但將來的總體趨勢和過去相反，即反持續[18]； 當H≥0.5時，表明時間序列具有長期相關的特征，即此過程具有持續性，且H越接近1，持續性越強。

3 結果與分析

3.1 植被時序變化特征分析

3.1.1 植被覆蓋度

利用像元二分模型結合NDVI數據獲取研究區2000年和2017年的植被覆蓋度(圖3)。參考相關植被覆蓋度劃分的文獻[19],結合該研究區的氣象數據、植被生長狀態以及實地勘察情況，綜合分析可知，在戈壁以及鹽堿化嚴重的區域，往往只有少量梭梭樹以及檉柳，覆蓋度比較低，而在枸杞種植園地由于人為的作用，植被發育較好。因此，將植被覆蓋度以0.5為節點，處于[0,0.3)，[0.3,0.5)，[0.5,0.8)和[0.8,1]分別劃分為極低植被覆蓋度、低等植被覆蓋度、中等植被覆蓋度與高等植被覆蓋度。在整個區域18 a的時間內，極低植被覆蓋度由66.04%上升到67.71%，低等植被覆蓋度由12.84%上升到13.63%，中等植被覆蓋度由12.33%下降到11.28%，高等植被覆蓋度由8.79%下降到7.38%，總體上來看植被變化并不大，整體植被發育比較弱。

(a) 2000年(b) 2017年

3.1.2 植被指數

MODIS衛星和Landsat衛星的空間分辨率分別為250 m和30 m，時間分辨率分別為1 d和16 d。采用不同空間、時間分辨率的衛星產品，在不同空間、時間尺度上分析植被的變化特征，以對比分析使結果更具代表性。

利用最大化合成法得到的NDVI年均值和最大EVI年均值來表達每年植被最茂盛的時期(圖4)。最大化合成NDVI年均值由0.029上升到0.054，增幅為0.025，平均年增速為0.208%(k=0.002 19,R2=0.806 73)； 最大EVI年均值從0.633上升到0.771，增幅為0.138，平均年增速為1.15%(k=0.007 73,R2=0.816 56)。結果顯示在2000—2008年內2組數據均處于相對平穩狀態，在此時間段內植被變化不大，在2008—2017年內2組數據顯示的植被變化趨勢處于大幅度的上升狀態。

(a) 最大化合成NDVI年均值(b) 最大EVI年均值

利用時間分辨率為1 d的MODIS數據分析月際植被變化特征，最大EVI月均值的結果高值集中在每年5—10月(圖5(a))，以此判斷植被比較茂盛的時間段是在5—10月份，植被增長最快的時間約為每年6—7月。

(a) 最大EVI月均值(b) 月平均降水量和溫度

植被變化時間特征與降水量、溫度在每年的月際分布基本一致，平均降水量和溫度在1月份出現最小值(0.3 mm，-8.4 ℃)，在6—7月份出現最大值(17.1 mm，19.0 ℃)。從5月份開始，溫度逐漸上升，降水量也有所增加(圖5(b))，且最大EVI月均值對降水、溫度的變化亦有所響應，為了進一步反映它們之間的關聯，則引入相關與偏相關分析方法。植被生長季開始前溫度多在0 ℃以下，溫度升高有利于植物萌發。在植物生長季期間，平均溫度約為10 ℃，熱量充足且降水量的增加有利于提高土壤的含水量，進而促進植被發育[20]。

3.2 植被變化影響因素分析

基于18 a間諾木洪洪積扇的MODIS衛星16 d合成的EVI產品MOD13Q1，1 a內會有比較多的合成時間點，對每2個合成時間點間的溫度數據取平均，降水量取總和，初步形成對應時間點的變量集合(時間點也稱為樣本點，根據時間點篩選自然因素數據更有利于研究其與植被變化之間的關系)。為了使每個月的數據點具有代表性，對每個月的最大化EVI和溫度數據取平均值，降水數據取累計便得到需要進行分析的最大化EVI、溫度與降水變量集合(3個數據集合的變量個數均為411個)。

對最大化EVI值的平均數據集分別與平均溫度數據集、降水數據集進行相關分析，結果顯示EVI數據集與溫度數據集之間的Pearson相關系數是0.839，顯著性0.000<0.01,說明二者之間高度相關，即溫度的變化在植被變化過程中起到很大的作用。EVI數據集與降水數據集之間的Pearson相關系數是0.457，顯著性0.000<0.01,說明二者之間具有低度相關性，則該研究區域的植被對溫度的敏感度要大于對降水的敏感度。

對3個變量集合進行一階偏相關分析，EVI數據集與溫度、降水數據集的一階偏相關系數分別為0.816和0.327，顯著性0.000<0.01。2個一階偏相關系數均大于0，即EVI數據與溫度、降水均呈現正相關的相關關系，且EVI數據和溫度、降水之間的正相關具有顯著性。

綜合分析結果可知： 在2000—2017年間，EVI數據與溫度之間具有顯著的強正相關性，與降水量則是呈現出弱正相關關系，說明溫度與降水是影響植被發育的重要因素，且影響該區域植被變化的主要因素為溫度，而降水對植被變化的影響效果比較弱。原因是諾木洪洪積扇降水總體稀少，年降水量僅17.8～177.5 mm，植被多為農作物與人工防護林，生產用水多依靠地下水與遠程高山融雪[21]。

3.3 枸杞種植區和鹽堿化區植被變化與制約關系

3.3.1 重點區植被覆蓋度變化趨勢分析

諾木洪洪積扇植被在18 a間內整體變化并不明顯，但枸杞園地和鹽堿化區變化則是相對突出，故在研究區內圈定2個重點研究區——枸杞種植區和鹽堿化區(圖1斜線區域)，劃分依據是以研究區2000年和2017年植被覆蓋度等級圖為基礎，結合水系支流數量和方向，將諾木洪河左側支流(哈西瓦河)、青藏公路和新扇扇緣部分劃為良好綠洲——枸杞種植區。水是影響鹽堿化區植被發育好壞的最重要的因素，一般河流尾部是植被變化相對劇烈的位置[22]，故將早期洪積扇的扇緣位置劃為鹽堿化區。

根據像元二分模型對枸杞種植區和鹽堿化區近18 a內的植被覆蓋度進行計算和擬合(圖6)。

(a) 枸杞種植區(b) 鹽堿化區

枸杞種植區的極低植被覆蓋度等級在時間段內有很大程度降低(k=-14.41)，中等植被覆蓋度有所降低，但變化程度不大(k=-2.02)； 相反高等植被覆蓋度上升較多(k=10.90)，低等植被覆蓋度也有所上升(k=5.58)。綜合分析枸杞種植區極低覆蓋度降低而高等植被覆蓋度上升現象，即在2000—2017年內枸杞種植區植被發育較好。而鹽堿化區的各個植被覆蓋度等級變化均不大，高等植被覆蓋度有降低的趨勢(k=-0.45)，低等和中等植被覆蓋度雖有略微上升趨勢，結合覆蓋度等級之間轉變關系及該區域主要受高等植被覆蓋度影響分析，鹽堿化區域在近18 a內植被是有所退化的。

3.3.2 重點區植被覆蓋度變化趨勢驗證

利用GEE平臺在枸杞種植區和鹽堿化區均勻生成驗證點，計算近18 a的NDVI進行一元線性回歸以獲取NDVI值和時間回歸方程的斜率值k，如圖7所示。枸杞種植區的k值較大，即植被改善速率較快，而樣本點10，15，23和55斜率值為負值，說明有些地塊植被出現衰減，考慮其是不適合農作物生長的地塊，通過人為因素的干預導致植被衰減。鹽堿化區k值比較小，多數值處于[0,5]，說明鹽堿化區域植被在此時間段內并沒有很大的改善。相反，樣本點1，8，16，17，37，45，51和54斜率值出現負值，說明鹽堿化區域部分植被退化與消亡。

圖7 2000—2017年間諾木洪洪積扇枸杞種植區和鹽堿化區NDVI變化k值

綜合分析18 a整個區域內植被逐漸發育，但枸杞種植區和鹽堿化區植被改善速率差異較大，枸杞種植區植被變化速率k均值為13.23，下游植被變化速率均值為3.47。故考慮是由于枸杞種植區域植被受到人為作用的干擾[23]，導致用水量逐漸增大，導致鹽堿化區域水資源減少，進而使得鹽堿化區植被改善速率遠低于枸杞種植區，反而有部分區域植被逐漸退化。

3.4 不同空間尺度的植被未來變化趨勢分析

3.4.1 研究區植被未來變化趨勢

利用Hurst指數對研究區未來變化趨勢分析(圖8)，研究區像元的H主要集中0.5～0.8區間之內，H≥0.5的像元占比為95.89%，H<0.5的像元占比為4.11%。故在未來一段時間內該區域植被變化與18 a來植被變化具有極大的一致性，洪積扇扇緣位置的持續性要低于扇中和扇根位置，原因是扇緣有發育較好的植被，植被發育受多方面因素影響，而扇中和扇根主要是戈壁和巖體，植被發育很低。

(a) Hurst指數分布(b) Hurst指數分布折線圖

3.4.2 重點區植被未來變化趨勢

枸杞種植區植被變化具有強持續性(圖9(a))，H均值為0.759，其中：H≥0.5，即持續性比重為99.169%，H<0.5，即反持續性比重為0.764%，即表示未來一段時間內植被變化與18 a間植被變化趨勢有很大程度的一致性。

(a) 枸杞種植區Hurst指數分布 (b) 鹽堿化區Hurst指數分布

鹽堿化區植被變化持續性相對于枸杞種植區要弱，但仍具有強持續性(圖9(b))，H均值為0.624，H≥0.5，即持續性比重為95.788%，H<0.5，即反持續性比重為3.605%，即表示未來一段時間內鹽堿化區植被變化與18 a間的植被變化趨勢有較大程度一致性。

分析可知鹽堿化區的持續性相對于枸杞種植區要低，由于鹽堿化區除了受到溫度和降水等氣象因素的影響之外，更重要的還受上游枸杞地種植區植被變化影響，上游人為因素干擾使得上游植被進一步發育時，則會導致下游區域植被持續退化。

4 結論

通過對研究區與重點區進行綜合分析其植被變化特征、植被變化成因以及未來趨勢，可以得到以下結論：

1)2000—2017年間，諾木洪洪積扇植被向好的方向轉變，但改善速率不大，主要原因是受到氣象因子的影響，且溫度的影響要遠大于降水的影響。

2)枸杞種植區受到人為作用使得植被持續發育，枸杞種植區的植被增長對鹽堿化區植被有著一定的制約效應，且在未來一段時間內會一直存在。

3)未來枸杞種植區與鹽堿化區均具有明顯的強持續性。

本文是依托于自然因素對植被變化的影響上來反映人類活動的作用，雖然采用18 a間的遙感和氣候數據進行分析，時間上滿足長時間序列分析的要求，但研究區內地表植被的變化十分復雜，物候信息在植被變化中也十分重要，若增加每個時間段內植被物候信息，分析結果精度將會進一步提高，也可以更好地說明人類活動的作用。