基于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.)