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

姚金璽， 張志， 張焜

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

0 引言

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

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

Google Earth Engine(GEE)平臺中數(shù)據(jù)包括近40 a的全球衛(wèi)星影像，且提供超過800 種功能函數(shù)[10]?；谠朴嬎愕?GEE 平臺可在線高效處理大規(guī)模地理數(shù)據(jù)集，有效地解決了遙感大數(shù)據(jù)中的處理難題[11-12]。

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

1 研究區(qū)概況及數(shù)據(jù)源

研究區(qū)主體位于青海省海西蒙古族藏族自治州都蘭縣諾木洪地區(qū)，地理坐標為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月，表現(xiàn)出雨熱同季[13]。區(qū)內(nèi)主要植被類型為梭梭、檉柳和蘆葦?shù)戎参?，覆蓋度不高，在晚期洪積扇扇緣及河流沿岸等含水量高的地區(qū)，植被長勢良好。

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

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

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

2 研究方法

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

圖2 技術(shù)路線

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

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

2.2 Pearson相關(guān)分析與一階偏相關(guān)分析法

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

(1)

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

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

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

(2)

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

偏相關(guān)系數(shù)檢驗的零假設(shè)為： 總體中2個變量間的偏相關(guān)系數(shù)為0。使用t檢驗方法，公式為：

(3)

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

2.3 Hurst指數(shù)

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

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

3 結(jié)果與分析

3.1 植被時序變化特征分析

3.1.1 植被覆蓋度

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

(a) 2000年(b) 2017年

3.1.2 植被指數(shù)

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

利用最大化合成法得到的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)。結(jié)果顯示在2000—2008年內(nèi)2組數(shù)據(jù)均處于相對平穩(wěn)狀態(tài)，在此時間段內(nèi)植被變化不大，在2008—2017年內(nèi)2組數(shù)據(jù)顯示的植被變化趨勢處于大幅度的上升狀態(tài)。

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

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

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

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

3.2 植被變化影響因素分析

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

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

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

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

3.3 枸杞種植區(qū)和鹽堿化區(qū)植被變化與制約關(guān)系

3.3.1 重點區(qū)植被覆蓋度變化趨勢分析

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

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

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

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

3.3.2 重點區(qū)植被覆蓋度變化趨勢驗證

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

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

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

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

3.4.1 研究區(qū)植被未來變化趨勢

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

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

3.4.2 重點區(qū)植被未來變化趨勢

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

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

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

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

4 結(jié)論

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

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

2)枸杞種植區(qū)受到人為作用使得植被持續(xù)發(fā)育，枸杞種植區(qū)的植被增長對鹽堿化區(qū)植被有著一定的制約效應(yīng)，且在未來一段時間內(nèi)會一直存在。

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

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