基于MODIS數據和BFAST方法的植被變化監測

劉寶柱, 方秀琴, 何祺勝, 榮祁遠

(河海大學地球科學與工程學院，南京　210098)

?

基于MODIS數據和BFAST方法的植被變化監測

劉寶柱, 方秀琴, 何祺勝, 榮祁遠

(河海大學地球科學與工程學院，南京210098)

摘要：植被是聯結土壤、大氣和水分的自然“紐帶”，在全球氣候變化研究中具有“指示器”的作用。對歸一化植被指數(normalizeddifferencevegetationindex，NDVI)時間序列分析，可以為相關部門的工作和決策提供更好的支持。使用MODISNDVI數據結合BFAST(breaksforadditiveseasonalandtrend)方法實現對老哈河流域及周邊地區的植被變化監測，并確定其NDVI時間序列出現突變點的時間節點。結合氣象數據以及數據本身的質量作為影響因子，分析出現突變點的主要原因。研究結果表明，降水量、相對濕度、溫度、日照時數、流域蒸發量與NDVI變化趨勢呈正相關，風速與NDVI變化趨勢相關性很小。降水量對NDVI變化的影響具有滯后性，滯后時間與降水量大小有關。關鍵詞：NDVI； 時間序列；BFAST； 變化監測； 突變點

0引言

植被是地球系統生物圈的重要組成部分，包括森林、灌叢、草地、農田和果園等，既是生態環境的重要組成部分，又是維持生態環境、發揮有效生態效能的功能體，是衡量自然生態環境狀況和性質的主要指示物[1]。植被是聯結土壤、大氣和水分的自然“紐帶”，在全球氣候變化研究中具有“指示器”的作用[2]。植被覆蓋與氣候變化的關系是全球變化研究的重要內容。研究地表植被覆蓋變化和氣候因子之間的關系，對區域乃至全球的能量、生物化學循環和水循環以及氣候變化響應具有重要意義。

歸一化植被指數(normalizeddifferencevegetationindex，NDVI)能夠較好地反映植被的代謝強度及其季節性變化和年際變化特征，因而被廣泛應用于植被的監測與分類、物候分析、農作物估產、植被-氣候關系分析等方面。利用NDVI時間序列數據研究不同類型植被對氣候變化的響應，是當前植被生態系統研究的熱點之一。最近幾年提出的BFAST(breaksforadditiveseasonalandtrend)方法能夠很好地監測出NDVI時間序列中的季節和趨勢突變，特別是能夠指示突變發生的時空分布[3]。

老哈河流域是我國北方典型的半干旱區。本研究基于250m的MODISNDVI的時間序列數據，利用BFAST方法，監測該區域的NDVI變化狀況，檢測突變信息，結合氣象數據分析突變的原因，以及產生的結果，為老哈河流域的生態安全評價、區域環境規劃以及生態建設提供參考依據。

1研究區概況及數據源

本文選取老哈河流域作為研究區，其范圍如圖1所示。老哈河位于內蒙古自治區東部，是西拉沐淪河的一級支流，流域面積約為1.9×104km2，位于N41°～43°，E117°～120°，流域內地形較復雜，流域出口處高程約400m，源頭的七老圖山海拔2 000m左右。該流域屬中溫帶半干旱大陸性氣候，干燥少雨多風沙，年降雨量和年平均氣溫分別從西北部山地的300mm和2 ℃變化到東南部低地的600mm和7 ℃。該地區的農作物生長期為5—10月，農耕期為4—11月[4]。老哈河流域下墊面多為黃土丘陵區，地表植被覆蓋率較低，水土流失較為嚴重。

圖1　研究區范圍

研究選用的遙感數據是2000年2月—2004年12月間Terra衛星MOD13Q1的16d數據，空間分辨率為250m。這個數據可以頻繁地提供由人為因素造成的土地覆蓋變化信息[5]。MOD13Q1的16d數據是用受限視場角的最大NDVI值合成的[6]。影像在提取波段信息前，已經經過精幾何糾正和輻射校正。由于數據量較大，使用MRT(MODISreprojectiontool)軟件對MOD13Q1數據做拼接、裁剪、重投影以及NDVI波段提取的批處理。所用氣象數據來自區內6個氣象站點(圖1)。選取的氣象數據，包括日照時數、相對濕度、平均溫度、降水量、風速等主要用于與NDVI的突變結果做對比分析。其中NDVI值的提取是根據氣象站點的經緯度信息，以3像元×3像元范圍作為興趣區，取其平均NDVI值作為該氣象站點的NDVI值，以消除個別像元可能存在的偶然誤差。

2研究方法

2.1突變點檢測

流程如圖2所示。根據BFAST方法對老哈河流域6個氣象站點處于2000—2004年年底期間NDVI突變點的檢測，并結合氣象數據分析其突變點產生的原因及最大影響因子。BFAST方法具體使用的是R語言開發的BFAST程序包對得到的NDVI數據進行突變點檢測。BFAST程序包網址見http: //r-forge.r-project.org/R/?group_id=533。

圖2　突變點檢測流程圖

2.2分解模型

BFAST方法將集成的時間序列分解成趨勢、季節性和殘差3個部分[7]，根據一定的方法來檢測和表征趨勢和季節性部分中突然變化的部分。假設一個加法的分解模型可以迭代出與趨勢和季節性相匹配的分段線性模型。該模型的算法形式為[8]

Yt=Tt+St+et,t=1,…,n,

(1)

Tt=αj+βjt,

(2)

季節趨勢[9]可以表達為

(3)

式中： s為季節性的時長(如每年的觀察的次數)； γi,j為第i個時段的影響因子。當時間t為第i次觀測值時，dt,i=1，否則為0。因此，當t為季節性第0個時段時，dt,i-dt,0=-1。對于t為其他季節性時段時，dt,i-dt,0=1。dt,i通常被認為是季節性虛擬變量，它有2個可用的值為0和1來解釋回歸模型中的季節性時段。

2.3迭代算法

在擬合分段線性模型和估算突變點之前，先測試是否有突變點出現在時間序列中。基于移動總和MOSUM(movingsum)測試的普通最小二乘OLS(ordinaryleastsquares)殘差被用來測試是否存在一個或多個突變點。如果測試顯著性變化(P<0.05),則使用Bai和Perron的方法來估算突變點，并根據貝葉斯信息準則確定突變點的個數，以及突變點出現的時間和置信區間[10]。

2)趨勢系數αj,βj在區間j=1,…,m上，由基于M估算法的強健的回歸方程式(2)計算得到。趨勢估算為

(4)

4)季節性系數γi,j(i=1,…,s-1; j=1,…,m)可以從基于M估算方法的回歸方程中計算得出。季節性估算公式為

(5)

3結果與分析

3.1突變點的檢測結果

根據BFAST方法檢測突變點，6個氣象站點中發現赤峰氣象站和葉柏壽氣象站于2000—2004年年底期間出現季節性突變點，突變檢測結果分別如圖3所示。

(a) 赤峰氣象站 (b) 葉柏壽氣象站

圖3赤峰氣象站和葉柏壽氣象站NDVI時間序列分析

Fig.3ChifengmeteorologicalstationandYebaishoumeteorologicalstationNDVItimeseriesanalysis

輸出的突變點檢測結果見表1。表中括號內數值為第i次觀測時間節點1a內觀測的總次數。

表1　突變點檢測結果

從圖3(a)赤峰氣象站NDVI時間序列分析中可以看出，在季節性成分中出現了2個突變點，而在趨勢中沒有出現突變點，NDVI整體趨勢呈上升趨勢，但是上升趨勢并不明顯，其趨勢斜率為： 0.000 008 18。從整體數據可以看出，該點的植被生物量很低。從圖3(b)中可以看出，葉柏壽氣象站處的NDVI的季節性成分在2000年初至2004年底的時間序列中出現了2個突變點，趨勢成分中并沒有突變點的出現，且趨勢變化不大，其趨勢斜率為0.000 021 51，NDVI線性趨勢擬合始終維持在0.3左右。

3.2氣象數據與NDVI突變數據的對比分析

3.2.1赤峰氣象站數據對比分析

圖4為赤峰氣象站發生第一個季節性突變點處的NDVI趨勢圖和對應時間點的日照時數、相對濕度、氣溫、降水量以及風速圖的組合，便于進行赤峰氣象站季節性成分突變點位于時間節點2002(21)數據的對比分析。

(a) 赤峰氣象站第一個突變點 (b) 日照時數

(c) 相對濕度(d) 平均氣溫

(e) 降水量(f) 風速

圖4赤峰氣象站季節性成分突變點位于時間節點2002(21)的數據對比分析

Fig.4Chifengmeteorologicalstationseasonalcomponentbreakpointsdatacomparisonandanalysisat2002(21)

從圖4可以看出，風速、降水量、日照時數、相對濕度與NDVI之間沒有明確的相關關系，僅氣溫與NDVI變化呈正相關關系，所以氣溫是導致NDVI時間序列季節性成分上產生突變信息主要因素。

圖5為赤峰氣象站發生第二個季節性突變點處的NDVI趨勢圖和對應時間點的日照時數、相對濕度、平均氣溫、降水量以及風速圖的組合，便于進行赤峰氣象站季節性成分突變點位于時間節點2003(13)數據的對比分析。

(a) 赤峰氣象站第二個突變點 (b) 日照時數

(c) 相對濕度 (d) 平均氣溫

圖5-1赤峰氣象站季節性成分突變點位于時間節點2003(13)的數據對比分析

Fig.5-1Chifengmeteorologicalstationseasonalcomponentbreakpointsdatacomparisonandanalysisat2003(13)

(e) 降水量 (f) 風速

圖5-2赤峰氣象站季節性成分突變點位于時間節點2003(13)的數據對比分析

Fig.5-2Chifengmeteorologicalstationseasonalcomponentbreakpointsdatacomparisonandanalysisat2003(13)

從圖5中可以看出，突變點處的NDVI值與下一個時間點的NDVI相差約0.2，說明是由突變點處的氣象變化而產生的。同樣，風速與NDVI之間沒有明確的相關關系，從突變點的對應時間點看，相對濕度和降水量發生了較大的變化。因此，該處的季節性成分突變點主要是由相對濕度和降水量的變化造成的。因為降水量、相對濕度對植物的光合作用速率、生長發育速度、物候期出現的早晚都有很大的影響。而植被的生長快慢和物候期的提前到來都會對季節性成分中的突變做出貢獻。

3.2.2葉柏壽氣象站數據對比分析

圖6為葉柏壽氣象站發生第一個季節性突變點處的NDVI趨勢圖和對應時間點的日照時數、相對濕度、平均溫度、降水量以及風速圖的組合，便于進行葉柏壽氣象站季節性成分突變點位于時間節點2000(23)數據的對比分析。

(a) 葉柏壽氣象站第一個突變點 (b) 日照時數

(c) 相對濕度 (d) 平均氣溫

(e) 降水量 (f) 風速

圖6葉柏壽氣象站季節性成分突變點位于時間節點2000(23)的數據對比分析

Fig.6Yebaishoumeteorologicalstationseasonalcomponentbreakpointsdatacomparisonandanalysisat2000(23)

從圖6中可以得出，平均溫度、降水量與NDVI變化呈正相關關系，而日照時數、相對濕度、風速與NDVI相關性不明顯。從趨勢線中可以看出，平均溫度的趨勢線與NDVI趨勢線非常相似，則平均溫度為NDVI在20001001—20001230日期間變化的主要影響因子。植物蒸發蒸騰量通常主要包括植物蒸騰量和棵間土壤蒸發量。而平均溫度的降低導致植物的蒸騰量下降[11]。植物的蒸騰對植物的生長是有幫助的，蒸騰量越大，其生長速率越快。

圖7為葉柏壽氣象站發生第二個季節性突變點處的NDVI趨勢圖和對應時間點的日照時數、相對濕度、平均溫度、降水量以及風速圖的組合，便于進行葉柏壽氣象站季節性成分突變點位于時間節點2001(15)數據的對比分析。

(a) 葉柏壽氣象站第二個突變點 (b) 日照時數

(c) 相對濕度 (d) 平均氣溫

(e) 降水量 (f) 風速

圖7葉柏壽氣象站季節性成分突變點位于時間節點2001(15)的數據對比分析

Fig.7Yebaishoumeteorologicalstationseasonalcomponentbreakpointsdatacomparisonandanalysisat2001(15)

從圖7中可以看出，在2001年8月13日之前降水頻繁，降水量最高達到45mm，距突變時間點最近的一次降水發生2001年8月7日。降水量對NDVI的影響滯后了至少6d，降水量過大，導致物候期的提前到來，從而產生了季節性突變點。

3.3源數據精度分析

季節性突變點的產生很有可能是因為所采用的MODIS數據的本身精度就不夠高，或者源數據中有云的存在等等，導致突變點的產生不完全是因為氣候變化導致的。在ENVI中提取赤峰氣象站和葉柏壽氣象站處的數據精度信息，得到赤峰氣象站時間節點為2002(21)突變點的VI質量如表2所示。

表2　2002(21)突變點的VI質量

VI(2 189)： 很有可能有云像素，VI有用性為最低質量，氣溶膠含量低，不是云邊，沒有進行大氣二向性校正，不是混合云，檢測區域為陸地，無冰和雪覆蓋，無陰影。

VI(2 266)： 生產VI數據，但是由其他數據對其檢驗的，VI有用性不是很高，氣溶膠含量高，不是云邊，未進行大氣二向性校正，無混合云，檢測區域為陸地，無冰和雪覆蓋，無云陰影。

VI(18 508)： 產生VI，且質量好，VI有用性為最低質量，氣溶膠含量為中，不是云邊，未進行大氣二向性校正，不是混合云，檢測區域為暫時性有水區域，有冰和雪覆蓋，無云陰影。

從上面的數據質量結果可以看出，VI有用性普遍都很低，且不是云邊，未進行大氣二向性校正，不是混合云，在時間節點為2002(19)和2002(21)的檢測區域為陸地，而在2002(22)時間節點的檢測區域為暫時性水域，說明該處發生的突變情況主要由降雨造成； 此外，3個時間節點均未被冰或雪覆蓋。結合之前圖表對比分析，可以得出結論： 該季節性突變點出現的原因主要是由降雨造成的。

赤峰氣象站時間節點為2003(13)突變點的VI質量如表3所示。

表3　2003(13)突變點的VI質量

VI(2 257)： 產生的數據很有可能有云像素，VI有用性為質量較低，氣溶膠含量較低，不是云邊，未進行大氣二向性校正，無混合云，檢測區域為陸地，無冰和雪覆蓋，無陰影。

VI(34 842)： 產生VI，但是需要檢驗其他QA，VI有用性為質量較低，氣溶膠含量具有氣候性，不是云邊，未進行大氣二向性校正，無混合云，檢測區域為暫時性水域，無冰和雪覆蓋，有云陰影。

VI(2 189)： 很有可能有云像素，VI有用性為最低質量，氣溶膠含量低，不是云邊，沒有進行大氣二向性校正，不是混合云，檢測區域為陸地，無冰和雪覆蓋，無陰影。

從上面的數據質量結果可以看出，VI有用性普遍都很低，且不是云邊，未進行大氣二向性校正，不是混合云，在時間節點為2003(12)和2003(15)的檢測區域為陸地，而在2003(13)時間節點的檢測區域為暫時性水域，說明該處發生的突變情況主要是降雨造成，此外，3個時間節點均未被冰或雪覆蓋，且2003(15)時間節點處有云陰影，使得數據質量較低，但是突變點為2003(13)，影響不大。結合之前圖表對比分析，可以得出結論： 該季節性突變點出現的原因主要是由降雨造成的。

葉柏壽氣象站時間節點為2000(23)突變點的VI質量見表4。

表4　2000(23)突變點的VI質量

VI(2 120)： 產生VI，質量好，VI有用性為質量較低，氣溶膠含量為中，不是云邊，未進行大氣二向性校正，無混合云，檢測區域為陸地，無冰/雪覆蓋，無云陰影。

從上面質量檢查結果可以看出，該時間節點的數據質量較好，產生突變點的情況與數據本身無關，而從之前的圖表對比分析中得出，產生該突變的原因主要是與流域蒸發量有關。可以得出結論為： 2000(23)時間節點處的季節性突變主要是因為流域蒸發量的下降。

葉柏壽氣象站時間節點為2001(15)突變點的VI質量見表5。

表5　2001(15)突變點的VI質量

VI(2 116)： 產生VI，質量好，VI有用性為質量較低，氣溶膠含量為中，不是云邊，未進行大氣二向性校正，無混合云，檢測區域為陸地，無冰/雪覆蓋，無云陰影。

從上面質量檢查結果可以看出，該時間節點的數據質量較好，產生突變點的情況與數據本身無關，而從之前的圖表對比分析中得出，產生該突變的原因主要是與降水量有關。降水量最大值為45mm左右，且降水量對NDVI的影響具有滯后性，滯后時間為6d。所以可以得出結論為： 2001(15)時間節點處的季節性突變主要是因為降水量突然增多原因造成的。

4結論與展望

本文使用了BFAST方法將NDVI時間序列數據分成3個部分，分別為季節性成分、趨勢成分和殘差，并對季節性成分和趨勢成分做回歸分析，檢測并記錄突變點的信息； 結合氣象數據分析產生突變點的原因并得出結論。目前國內較少見做該方面研究的報道。本文結合了氣象數據對突變點產生的原因進行分析，更加詳細地解釋了產生突變的驅動因素和變化機制。研究得出以下結論：

1)在2000年初至2004年底，對6個氣象站點處的NDVI時間序列突變進行了監測，分別在赤峰氣象站監測到季節性突變點，時間節點為2002(21)和2003(13)以及葉柏壽氣象站監測季節性突變點，時間節點為2000(23)和2001(15)。

2)季節性突變點2002(21)產生的原因主要是由降雨造成的； 季節性突變點2003(13)產生的原因主要是由降雨造成的； 季節性突變點2000(23)產生的原因主要是由流域蒸發量下降造成； 季節性突變點2001(15)產生的原因主要是由降水量突然增多造成。

3)在2000年初至2004年底出現突變點的原因與數據本身質量無關，數據本身質量較好。

4)NDVI變化趨勢與降水量、流域蒸發量、氣溫、相對濕度、日照時數呈正相關趨勢，與風速相關性不大。降水量對NDVI變化的影響具有滯后性，滯后時間與降水量大小有關。

5)不同時期NDVI變化與降水量、流域蒸發量、氣溫、相對濕度、日照時數及風速等相關性都不一樣，影響程度也不同。

本文研究可以為自然和農業土地生產評估、農作物狀況以及物候變化研究提供參考依據。本研究還存在一些不足之處，只對氣象站點處的突變信息進行提取和分析，尚未實現整個面上的突變信息的提取和分析，在后續研究中，會結合全球同化數據對面上的突變點進行分析。

參考文獻(References)：

[1]張飛,塔西甫拉提·特依拜,丁建麗,等.新疆典型鹽漬區植被覆蓋度遙感動態監測——以渭干河-庫車河三角洲綠洲為例[J].林業科學,2011,47(7):27-35.

ZhangF,TashpolatT,DingJL,etal.DynamicallymonitoringvegetationcoverbyremotesensinginthetypicalsalinizationregionofXinjiang:AcasestudyindeltaoasisofWeiganandKuqaRivers[J].ScientiaSilvaeSinicae,2011,47(7):27-35.

[2]劉蕾.2001—2005年新疆植被覆蓋動態變化及原因分析[J].干旱環境監測,2007,21(3):146-148.

LiuL.AnalysisonXinjiangvegetationcoveragedynamicchangesanditscausesfrom2001to2005[J].AridEnvironmentalMonitoring,2007,21(3):146-148.

[3]VerbesseltJ,HyndmanR,NewnhamG,etal.Detectingtrendandseasonalchangesinsatelliteimagetimeseries[J].RemoteSensingofEnvironment,2010,114(1):106-115.

[4]方秀琴,任立良.西遼河的老哈河流域土地利用遙感動態監測[J].地球信息科學學報,2009,11(1):125-131.

FangXQ,RenLL.DetectionoflandusechangeinLaohaheRiverBasinofWestLiaoheRiverbasedonremote-sensingimages[J].JournalofGeo-InformationScience,2009,11(1):125-131.

[5]TownshendJRG,JusticeCO.Selectingthespatialresolutionofsatellitesensorsrequiredforglobalmonitoringoflandtransformations[J].InternationalJournalofRemoteSensing,1988,9(2):187-236.

[6]HueteA,DidanK,MiuraT,etal.OverviewoftheradiometricandbiophysicalperformanceoftheMODISvegetationindices[J].RemoteSensingofEnvironment,2002,83(1/2):195-213.

[7]HueteA,JusticeC,vanLeeuwenW.ModisVegetationIndex(MOD13):AlgorithmTheoreticalBasisDocument[EB/OL].(1999-04-30).http://eospso.gsfc.nasa.gov/ftp_ATBD/REVIEW/MODIS/ATBD-MOD-13/atbd-mod-13.pdf.

[8]VerbesseltJ,HyndmanR,ZeileisA,etal.Phenologicalchangedetectionwhileaccountingforabruptandgradualtrendsinsatelliteimagetimeseries[J].RemoteSensingofEnvironment,2010,114(12):2970-2980.

[9]LinderholmHW.Growingseasonchangesinthelastcentury[J].AgriculturalandForestMeteorology,2006,137(1/2):1-14.

[10]BaiJS,PerronP.Computationandanalysisofmultiplestructuralchangemodels[J].JournalofAppliedEconometrics,2003,18(1):1-22.

[11]屈艷萍,康紹忠,張曉濤,等.植物蒸發蒸騰量測定方法述評[J].水利水電科技進展,2006,26(3):72-77.

QuYP,KangSZ,ZhangXT,etal.Areviewofmethodsformeasurementofevaportranspirationfromplants[J].AdvancesinScienceandTechnologyofWaterResources,2006,26(3):72-77.

(責任編輯： 李瑜)

Monitoring the changes of vegetation based on MODIS data and BFAST methods

LIU Baozhu, FANG Xiuqin, HE Qisheng, RONG Qiyuan

( School of Earth Science and Engineering, Hohai University, Nanjing 210098, China)

Abstract:Vegetationisanatural“link”whichlinkssoil,airandwaterandan“indicator”inglobalclimatechangeresearch.Usingnormalizeddifferencevegetationindex(NDVI)time-seriesanalyses,wecanprovidebettersupportfortherelevantresearchesanddecision-making.UsingMODISNDVIdatabindingwithBFAST(breaksforadditiveseasonalandtrend)method,theauthorsimplementedmonitoringvegetationdynamicsintheLaohaheRiverBasinandthesurroundingareas,andidentifieditsNDVItime-seriesabruptchangepointsoccurringintime.Themeteorologicaldataandthequalityofthedataitselfwerealsousedasaninfluencefactoranalysisofthemainreasonforthebreakpoints.Itisfoundthatprecipitation,relativehumidity,temperature,sunshineandwaterevaporationarepositivelycorrelatedwithNDVItrends,whilewindspeedislesscorrelatedwithNDVItrends.What’smore,theprecipitationandsunshinehourimpactonNDVIchangehasacertainlag.

Keywords:NDVI;time-series;BFAST;changemonitoring;breakpoints

doi:10.6046/gtzyyg.2016.03.23

收稿日期：2015-04-14；

修訂日期：2015-08-28

基金項目：水利部公益性行業科研專項經費項目“東北灌區節水灌溉生態與增產效應評估研究”(編號： 201401001)、國家自然科學基金項目“復雜地形條件下多源遙感數據森林生物量協同反演研究”(編號： 41101308)及“北方半干旱區典型土地利用變化的水文效應”(編號： 41201027)共同資助。

中圖法分類號：TP79

文獻標志碼：A

文章編號：1001-070X(2016)03-0146-08

第一作者簡介：劉寶柱(1993-)，男，碩士研究生，主要從事定量遙感研究。Email:bzliu@hhu.edu.cn。

通信作者：何祺勝(1981-)，男，博士，副教授，主要從事定量遙感、GIS建模與開發等方面的研究。Email:heqis@hhu.edu.cn。

引用格式： 劉寶柱,方秀琴,何祺勝,等.基于MODIS數據和BFAST方法的植被變化監測[J].國土資源遙感,2016,28(3):146-153.(LiuBZ,FangXQ,HeQS,etal.MonitoringthechangesofvegetationbasedonMODISdataandBFASTmethods[J].RemoteSensingforLandandResources,2016,28(3):146-153.)