基于MODIS數據的新疆草地物候提取方法及變化趨勢分析

張仁平，郭靖，馬曉芳，郭偉勇

（1.新疆大學資源與環境科學學院，新疆 烏魯木齊 830046；2.新疆大學綠洲生態教育部重點實驗室，新疆 烏魯木齊 830046；3.新疆林業科學院，新疆 烏魯木齊 830000；4.中國科學院西北生態環境資源研究院，甘肅 蘭州 730000；5.察布查爾縣農業農村局，新疆 伊犁 835300）

植被物候是指受氣候和環境因素的影響而出現的以年為單位的周期現象，它是指示氣候與自然環境變化的重要指標，對于研究全球生態環境變化至關重要［1—3］。以往的植被物候觀測是通過定點定時的野外目視觀察法，雖然能夠準確地獲取觀測點周圍的植被物候信息，但無法應用推廣到大區域［4］。遙感監測因具有長時間序列、覆蓋范圍廣等特點，目前成為物候觀測的主要手段［5—6］。歸一化植被指數（normalized difference vegetation index，NDVI）能夠反映植被光合有效輻射的變化情況，被廣泛應用于多尺度下植被動態變化監測、土地覆被變化、植物生物物理參數反演及區域環境變化等方面［7—9］，但由于傳感器本身及其受云、大氣氣溶膠、冰雪等因素干擾，從衛星數據中獲得的NDVI時間序列數據不可避免地受到不同程度的噪聲污染，無法真實的反映植被生長狀況，因而需要對其進行重建［10—12］。NDVI時序數據重構方法有很多，其中雙邏輯斯蒂函數擬合法（double logistic，D-L）僅依賴NDVI曲線的變化特征，即可客觀地反映植被生長過程的細節特征，且該方法對噪聲的處理及物候特征提取能力較為突出［13］。S-G濾波法（Savitzky-Golay）通過NDVI數據的上包絡線來擬合NDVI時序數列的變化趨勢，保留了NDVI曲線的特征，通過一個迭代的過程可消除偏離正常生長趨勢線的噪聲，而使平滑達到最好的效果。同時該方法不僅不受數據時間空間尺度和傳感器的限制，還能清晰地描述NDVI時間序列的長期變化趨勢以及局部的突變信息［14］。非對稱性高斯函數擬合法（asymmetric gaussian，A-G）能夠較好地描述NDVI時序數據的總體變化趨勢和全局特征，在噪聲點優化及保留原始數據質量方面，具有更好的擬合效果［15］。基于此，本研究選用A-G、D-L和S-G算法重構NDVI時序數據。美國NASA根據Zhang等［16］基于Logistic函數對美國東北地區的MODIS EVI（enhanced vegetation index）進行擬合，形成物候產品數據（MCD12Q2）。因為此種算法不需要界定植物物候的閾值，普遍適用于大范圍的物候監測［17—18］。因此，選擇最佳的遙感物候提取方式，對于特定區域植被物候監測及準確進行物候監測具有一定的意義。

目前，已有大量關于草地物候時空變化的研究。在國外，Delbart等［19］利用Spot-vegetation（VGT）和NOAA advance very high resolution radiometer（AVHRR）數據分析了1982—2004年歐亞大陸植物物候變化趨勢，結果表明該地區植物返青日期有所提前。Zeng等［20］利用MODIS NDVI數據對2000—2010年北半球的植物物候研究顯示，植被返青仍然保持提前的趨勢。在國內，王宏等［21］基于NOAA NDVI和MSAVI數據對中國北方植物物候變化趨勢的分析表明，荒漠草原和典型草原的返青期（start of growing season，SOS）提前，枯黃期推遲。倪璐等［22］基于1986—2015年的GIMMS NDVI數據，定量分析了中國天然草地物候的時空變化特征，結果表明1998年以前全國草地SOS平均提前速率為0.37 d·a—1，1998年以后草地物候變化趨勢相反，草地SOS平均推遲速率為0.28 d·a—1。何寶忠等［23］通過MODIS NDVI數據提取了新疆地區2001—2016年典型植被物候期，發現新疆植被物候具有明顯的緯向分布和垂直地帶性分布特征，且與全球大背景下典型植被物候特征變化趨勢相反，新疆植被SOS呈推遲趨勢，推遲幅度為1.9 d·10 a—1。然而，張仁平等［24］基于MCD12Q2數據，對新疆地區2001—2014年草地物候的時空變化進行了研究，前者的研究表明新疆草地返青期整體呈現提前趨勢，提前速率為0.11 d·a—1。通過以上研究可以看出，不同區域和不同時間段物候的變化區域存在著一定的差異，進一步研究新疆地區草地物候變化趨勢，有助于更深刻地理解氣候變化帶來的影響。

1 材料與方法

1.1 研究區概況

新疆維吾爾自治區地處我國西北隅（34°22′—49°33′N，73°22′—96°21′E），總面積為1.66×106km2，約占全國總面積的1/6。新疆具有“三山夾兩盆”地理環境特征，自南向北依次為昆侖山脈地—塔里木盆地—天山山脈—準格爾盆地—阿爾泰山脈，平均海拔2000 m以上。由于遠離海洋，且四周高山阻隔，該區形成了典型的溫帶大陸性氣候，具有氣溫溫差較大、降水量少，光熱資源充足等特點。新疆草地資源豐富，主要包括低地草甸、高寒草原、高寒草甸、溫性荒漠、溫性草原、高寒荒漠等6種草地類型，草地總面積約57.26萬km2，居全國第3位。該區草地資源不僅是當地農牧民重要的生產資料，而且草地生態系統關系著當地和周邊地區的生態與環境狀況以及國家生態安全（圖1）。

圖1 研究區草地類型及地面觀測點分布Fig.1 Distributions of the grasslands and ground observation points in Xinjiang

1.2 數據源

NDVI數據來源于美國地質勘探局（https：//earthdata.nasa.gov/eosdis/daacs/lpdaac）的MOD13Q1 V 006，該數據是由16 d合成的250 m×250 m空間分辨率的植被指數產品，本研究下載了2001—2019年共計437景影像為數據源，對新疆地區的草地返青期進行提取和分析。物候產品源于美國國家航空航天局（https：//lpdaac.usgs.gov/）的MODIS土地覆蓋動態產品（MCD12Q2），時間為2001—2019年，空間分辨率為500 m×500 m。為使遙感反演物候數據與后期地面物候觀測資料的對比分析結果較為可靠，本研究將MOD13Q1 NDVI重采樣至500 m。

地面草地返青數據采用2012—2014年2—5月的物候觀測數據。該數據共計123個樣地，樣地大小為500 m×500 m，每個樣地設5～8個樣方，每個樣方大小為1 m×1 m，每隔5～7 d對采樣點完成一次觀測。當樣地中植物開始展葉并正常發育的株（叢）或返青蓋度占總株（叢）的百分率為50%時，即認為該樣地返青［6］。原始地面物候觀測資料為日期型格式，需將其轉換為數值型格式以便后期研究分析。

1.3 研究方法

1.3.1NDVI時間序列重構 1）Savitzky-Golay濾波法（S-G）是一種基于最小二乘法的卷積算法，由Savitzky等［25］于1964年最早提出。其基本原理：用一定長度的過濾器和待處理數據作卷積，對待處理的數據作加權多項式擬合，擬合的目標是求得最小均方根誤差，而一些遠離大多數點的邊沿點不會參與擬合。這樣擬合時過于偏離正常生長趨勢線的噪聲部分會被丟棄。其公式可表示為：

式中：Y*j為平滑后的NDVI數據；Y j+i代表原始NDVI數據；Ci為擬合系數，表示從濾波器首部開始第i個NDVI值對應的權重；N為濾波器的長度，其值越大濾波結果越平滑，但也會遺漏能反映真實NDVI的細節數據，反之易產生大量的重復數據。因此本研究經過反復的實驗，將Savitzky-Golay濾波窗口的值設為4。

顱內靜脈竇血栓屬于特殊性的腦靜脈系統缺血性腦血管病，患者的主要臨床癥狀包括意識障礙、視乳頭水腫以及局灶性神經功能障礙等[1]。顱內靜脈竇血栓具有發病形勢多樣、發病原因復雜以及臨床癥狀、體征缺乏特異性等特點，且患者病情進展迅速，增加了誤診以及漏診率[2-3]。因此，本研究通過分析診斷時間與顱內靜脈竇血栓患者臨床特征及預后的關系，旨在為顱內靜脈竇血栓的早期診斷提供參考依據。現報道如下。

2）非對稱性高斯函數擬合法（A-G）是基于不對稱高斯函數的非線性最小二乘擬合算法，由Jonsson等［26］于2002年提出。其原理是對植被生長過程用分段高斯函數模擬，然后通過對各高斯擬合曲線進行平滑連接實現時間序列重建，公式可表示如下：

式中：f（t）為局部擬合函數；t為時間；c1和c2控制曲線的基準和幅度；a1決定峰值和谷值的位置；a2，a3，a4，a5控制曲線的寬度和陡峭度；g（a1，…，a5）為時間t的高斯函數；局部擬合后通過整體函數F（t）連接。

式中：F（t）為整體函數；tl，tr為時間序列數據中未擬合部分的區間；（t）為［0，1］之間的剪切數；fl（t）、fc（t）、fr（t）分別為區間內左邊最小值、中間最高值和右邊最小值的局部擬合函數。

3）雙邏輯斯蒂函數擬合法（D-L）是一種半局部擬合的經典方法，由Beck等［27］在2006年提出。其原理是首先將整個時間序列中時間點對應值按極大、極小分成多個區間，然后分別對該區間進行D-L函數局部擬合，其局部擬合方式與非對稱高斯擬合方法（A-G）類似。該方法具有無需確定閾值，僅依賴NDVI曲線的變化特征，即可客觀地反映植被生長過程的特征。其公式表示為：

式中：a1，a2，a3，a4控制左、右半部分的寬度和陡峭度。其整體連接函數F（t）與非對稱性高斯函數擬合法（A-G）相同，見公式4。

1.3.2物候期的提取 植被物候信息提取的方法較多，常見的有擬合法、最大斜率法、主成分分析法、滑動平均法、閾值法等［28］，雖然以上方法都是基于植物的生理特征和遙感機理，但對于不同的應用區域、數據質量等方面有所差異。本研究基于前人研究經驗及新疆地區實際情況［23，29］，采用Jonsson等［26］提出的動態閾值法，將NDVI增長達到當年NDVI振幅20%的時刻定義為返青期的閾值，基于TIMESAT軟件逐年提取新疆地區的草地返青期。

MCD12Q2產品是以MODIS增強型植被指數（enhanced vegetation index，EVI）為數據源，先對EVI時序數據進行預處理，消除氣溶膠、云和觀測角度等影響，然后用滑動窗口（連續5個時相EVI值）來判斷持續升高和持續降低區間，當所在區間極大值和振幅滿足給定閾值條件后，判斷該升高—降低區間為一個生長周期［30］，最后對于每一個生長周期采用分段Logistic函數擬合，并利用其曲線變化率確定植被生長始期、成熟期、末期等。其中，將EVI增長的節點日期定義為返青期。

1.3.3對比驗證 地面實測數據是植被物候遙感監測精度評價的最有效手段。基于新疆地區2012—2014年共計123個樣地觀測數據作為真值，采用皮爾森相關系數（R），偏差（Bias）和均方根誤差（root mean square error，RMSE）對4種方法提取的返青期進行對比分析。

皮爾森相關系數［31］是用來描述兩個隨機變量線性相關的統計量，其值介于—1與1之間。

式中：n為用于驗證的樣本數據；x i為遙感提取的返青期；xˉ為由遙感提取的2012—2014年返青期平均值；yi為地面觀測的返青期；為由地面觀測的2012—2014年返青期平均值。當R＞0時，表明遙感提取的返青期與地面觀測的返青期呈現正相關；當R＜0時，結果相反；當R＝0時，表示遙感提取的返青期與地面觀測的返青期不存在線性相關關系。

Bias是指遙感返青期與地面實測值之間的差別，其描述的是模型本身的準確度［32］。其計算公式如下：

RMSE同樣用來表示遙感返青期與地面實測值的偏差，從而反映不同方法提取返青期的模擬精度［33］。

式中：N代表地面觀測點的個數；Ei為第i個遙感提取的返青期；Oi為第i個地面觀測點的返青期。Bias與RMSE越接近0，表明遙感提取的返青期與地面觀測的返青期越接近。

1.3.4趨勢分析 本研究采用線性回歸模型，從像元尺度上來定量描述新疆地區草地返青期（start of growing season，SOS）的變化趨勢，以此來分析研究區草地物候的時空變化特征［34］。其公式表示如下：

式中：b代表斜率，其值大于0，說明變化趨勢呈現推遲的趨勢，相反則是提前的趨勢；n為年數；i為年序號，Pi表示第i年對應的植被返青期。

同時，采用F檢驗對變化趨勢的顯著性進行評價。

式中：k＝1，n＝19。通過F檢驗表查不同水平對應的R值，得到F0.01（1，18）＝8.53，F0.05（1，18）＝4.94，F0.1（1，18）＝3.05。當|R|≥0.59時，通過置信度0.01水平的檢驗（變化趨勢極顯著）；當0.49＜|R|＜0.59時，通過置信度0.05水平的檢驗（變化趨勢顯著）；當0.40＜|R|＜0.49時，通過置信度0.1水平的檢驗（年際變化趨勢略顯著）。

2 結果與分析

2.1 對比驗證

研究中用于4種方法對比分析一致的地面實測點，均為123個。受MCD12Q2產品數據缺失較多的影響，研究區部分地面實測點對應的MCD12Q2產品返青期為空值，最終可用于分析的僅有26個地面實測點。

從圖2可以看出，2012—2014年新疆地面實測點返青期主要集中在每年第68～122天，A-G、D-L、S-G反演的返青期主要集中在第60～120天，而MCD12Q2產品的返青期在第100～130天。可見，A-G、D-L和S-G反演的返青期比地面實測返青期提前，而MCD12Q2產品返青期推遲。從返青期提取結果與地面觀測數據的相關系數R來看，A-G、D-L返青結果和地面實測數據的相關性較高，分別達到0.879、0.880，而S-G、MCD12Q2產品與地面實測數據的相關性較低，分別為0.870、0.626。從返青期提取結果與地面觀測數據的誤差指標RMSE來看，MCD12Q2最小，RMSE為12.692 d，其次是A-G和D-L監測結果，RMSE分別為16.395 d、18.754 d，而S-G方法監測結果與實測返青期誤差量較大，RMSE大于20 d。從返青期提取結果與地面觀測數據的擬合度Bias來看，AG和D-L提取結果與地面觀測數據的擬合度較優，分別為—0.058、—0.059，MCD12Q2提取結果與地面觀測數據的擬合度較差，為—0.163。上述分析結果表明，4種方法中，基于A-G方法返青期提取結果精度較高。

圖2 返青期觀測數據與遙感反演數據的對比分析Fig.2 The comparison between phenology observed in the ground stations and phenology extracted by asymmetric gaussian，double logistic，Savitzky-Golay and MCD12Q2 retrieval

由于遙感提取的物候參數與地面實測存在尺度不統一等問題［35］，因而進一步從空間層面上探究4種方法提取新疆草地返青的最優模型。標準差衡量的是樣本數據的離散程度。一般來說，標準差越小，表明數據越聚集；標準差越大，則表明數據越分散。

從圖3可以看出，4種方法反演的返青期在北疆地區的標準差總體上較低，而南疆地區返青期的標準差較高。結合表1，4種方法中返青期標準差（＜15 d）和標準差（15～30 d）的面積比例較大。A-G、D-L、S-G及MCD12Q2產品的返青期標準差（＜30 d）的面積比例分別為82.19%、76.53%、67.98%、96.81%，以A-G方法提取的返青期標準差（＜30 d）的面積比例最多。MCD12Q2產品的返青期標準差（＜30 d）的面積比例最高，但是其反演的返青期在空間上空值較多。

表1 不同方法反演的返青期標準差面積比例Table1 Standard deviation area ratio of SOSinverted by different methods（%）

圖3 不同方法反演的返青期（SOS）標準差分布Fig.3 Distr ibution of the standar d deviation of SOSin Xinjiang inverted by differ ent methods

綜上所述，4種方法中A-G方法反演的返青期的結果最佳。因此，后面關于新疆草地返青期的空間分布及時間變化趨勢等的研究分析都是基于A-G提取的返青期來展開的。

2.2 草地物候空間分布特征

圖4為2001—2019年新疆草地平均返青期的空間分布格局。新疆地區草地返青期分布具有明顯的地域差異，呈現出自北向南逐漸推遲的分布特征。2001—2019年該地區返青期主要集中在第60～140天，多年均值約為第100天。其中，北部準噶爾盆地和伊犁河谷區域的草地返青時間最早，早于第80天，而阿爾泰山、天山中部及昆侖山等區域的草地返青時間最晚，晚于第140天。結合圖5得到，新疆不同草地類型返青期存在一定的差異。其中，高寒草甸與高寒草原的返青期最晚，約為第130天，其次是高寒荒漠與低地草甸，分別為第127、124天，再次為溫性草原，為第103天，溫性荒漠返青時間最早，約為第86天。

圖4 基于A-G方法新疆草地返青期（SOS）的空間分布Fig.4 Spatial distribution of SOSin Xinjiang based on A-G method

圖5 基于A-G方法不同草地類型返青期（SOS）的多年均值Fig.5 Multi-year mean values of SOS for different grassland types based on the A-G method

2.3 草地物候年際變化特征

基于像元尺度對2001—2019年新疆草地的返青期與年份之間進行了線性擬合，得到過去19年間該區返青期的變化趨勢。從圖6可以看出，北部準噶爾盆地邊緣地帶、天山南部及昆侖山東段少數區域的草地返青期呈現提前的趨勢，提前的面積比例約為46.93%，提前幅度多為1 d，而草地返青期推遲的區域主要集中于天山東部和昆侖山北部邊緣地區，推遲的面積比例約為53.07%，推遲幅度約為1 d。總體上來看，草地返青期通過P＜0.1的有效像元比例約為10.00%，其中2.0%像元的SOS在0.05的顯著性水平上呈顯著的提前趨勢，平均每年提前1～2 d，主要分布在阿爾泰山以南、準噶爾盆地西部邊緣區，而2.7%像元的SOS在0.05的顯著性水平上呈顯著的推遲趨勢，平均每年約推遲2 d，主要在塔里木盆地西北部邊緣及準噶爾盆地東部少數區域。

圖6 基于A-G方法新疆草地返青期（SOS）的空間變化及顯著性檢驗Fig.6 Var iation tr end and significance test of grassland SOSin Xinjiang based on A-G method

由圖7可知，低地草甸、溫性荒漠和高寒荒漠類型下，A-G方法提取的返青期均呈現推遲的變化趨勢，溫性草原、高寒草原及高寒草甸的返青期表現為提前。從返青期推遲變化來看，溫性荒漠的SOS推遲最為明顯，推遲速率為0.456 d·a—1，其次是低地草甸，推遲速率為0.379 d·a—1，而高寒荒漠的返青期稍有推遲，推遲速率為0.313 d·a—1。從返青期提前變化來看，溫性草原的返青期明顯提前，提前速率為0.384 d·a—1，高寒草原次之，返青期的提前速率為0.196 d·a—1，而高寒草甸的返青期略有提前，提前速率為0.127 d·a—1。分時間段來看，2001—2010年，低地草甸、溫性荒漠、溫性草原、高寒草原及高寒草甸的返青期均呈現推遲趨勢，分別推遲2、13、13、1、11 d，而2010—2019年溫性荒漠、溫性草原、高寒草原及高寒草甸的返青期變化趨勢相反，分別提前8、14、3、5 d。

圖7 基于A-G方法不同草地類型返青期（SOS）的年際變化趨勢Fig.7 The interannual variation trend of SOSin different grassland types based on A-G method

3 討論

不同的遙感物候提取方法各有優缺點，本研究采用4種不同構建原理的植被物候提取算法，提取了新疆地區草地植被返青期。通過123個地面實測數據對不同方法的提取結果進行比較和分析，發現A-G方法反演的返青期結果和地面實測數據的相關性較高，擬合度較優，誤差相對較小。雖然4種數據中MCD12Q2產品的返青期標準差（＜30 d）的面積比較最大，但由于該產品數據缺失較多，地面有效實測點僅有26個，因而該方法驗證結果誤差相對較大。

新疆屬于典型的干旱半干旱地區，生態系統脆弱，研究該區域的植被物候特征變化趨勢，對深入認識氣候變化對陸地生態系統的影響具有典型性。安佑志［36］的研究發現，受海拔因素的影響，新疆北部山區的返青時間要比同緯度地區植被稍晚些，其中阿爾泰山的返青時間最晚，準噶爾盆地和伊犁河谷地區的返青時間最早。同時得到20世紀90年代之前新疆北部地區的返青期以提前為主，南部部分地區的返青期顯著延遲，本研究估算的返青期具有明顯的垂向分布差異結論與此基本一致。已有研究表明，氣候變暖導致了21世紀前物候返青期的提前，但隨著升溫程度的減弱，這種趨勢在近幾十年可能會發生變化［3，22］。Piao等［37］認為，自1982年至今青藏高原植被返青期并沒有明顯的變化趨勢，而是在不同階段返青期提前和推遲有所不同。張雯等［38］的研究表明，1982—2013年內蒙古地區植被返青期總體呈現推遲趨勢，尤以2010—2013年推遲趨勢較為明顯。從新疆草地返青期整體變化趨勢來看，本研究與何寶忠等［23］的研究結果較為一致，均認為新疆草地的返青期呈現推遲趨勢。從以上分析可以看出，植被物候的變化與植被類型、研究區域、研究時段以及遙感數據等有關系。

4 結論

本研究通過4種不同物候提取方法，進行新疆草地返青期的提取，并基于地面實測數據，對遙感提取返青期進行精度對比驗證，同時基于最優模型分析了新疆草地返青期時空變化。主要結論如下：

1）4種方法中A-G方法反演的新疆草地返青期的結果最佳。基于A-G方法提取的新疆草地返青期和實測點的草地返青期相關性較高（0.879），RMSE相對較小（16.395 d）。同時A-G方法提取的新疆草地返青期標準差（＜30 d）的面積比例最多，達到82.19%。

2）新疆地區草地返青期主要集中在第60～140天，具有自北向南逐漸推遲的明顯地域差異。其中，北部準噶爾盆地和伊犁河谷區域的草地返青時間最早，而阿爾泰山、天山中部及昆侖山等區域的草地返青時間最晚。不同草地類型返青期存在一定的差異，以高寒草甸與高寒草原的返青時間（第130天）最晚，而溫性荒漠返青時間（第86天）最早。

3）新疆地區的草地返青期總體呈現微弱的推遲趨勢，其中阿爾泰山以南、準噶爾盆地西部邊緣區（2.0%像元）的SOS在0.05的顯著性水平上呈顯著的提前趨勢，平均每年約提前1～2 d，而塔里木盆地西北部邊緣及準噶爾盆地東部少數區域（2.7%像元）的SOS在0.05的顯著性水平上呈顯著的推遲趨勢，平均每年約推遲2 d。不同草地類型中，低地草甸、溫性荒漠和高寒荒漠的返青期均呈現推遲的變化趨勢，而溫性草原、高寒草原及高寒草甸的返青期表現為提前趨勢。此外，溫性荒漠、溫性草原、高寒草原及高寒草甸的返青期以2010年為界，呈現“先推遲后提前”的變化趨勢。