基于时序突变检测的植被空间变化特征识别方法——以海河北部山区为例

海 月,杨广斌,李若男,郑 华

1 贵州师范大学地理与环境科学学院, 贵阳 550025 2 贵州省山地资源与环境遥感应用重点实验室, 贵阳 550001 3 中国科学院生态环境研究中心城市与区域生态国家重点实验室, 北京 100085

植被是覆盖地表植物群落的总称,是生态系统的重要组成部分。植被既受大气、土壤、水分等环境要素的综合影响,又能反作用于环境[1-2]。研究植被长期变化不仅能够揭示生态系的演变过程并指示生态系统健康状况；更有利于人类对植被进行合理的保护、利用及管理[3-4]。随着遥感技术的发展、数据获取的便捷及分析手段的提升,基于植被长时序趋势特征的研究越来越多,趋势分析成为常用的分析方法。归一化植被指数(Normalized Vegetation Index,NDVI)能够反映植被的生长状态、覆盖情况、生物量及环境要素影响[5-8],并作为植被质量评估指标应用十分广泛。长时序NDVI分析的优势不仅在于能有效的反应植被变化的情况,还在于不同数据源能够混合使用,使得大空间尺度的植被质量评估成为可能[9-11]。

在全球气候变化的大背景下,利用长时序NDVI的空间特征评估气候变化带来的影响是关注的重点。许多研究已从不同空间尺度上(全球、区域、流域等)研究了气候变化对植被及土地利用变化的影响[12-17]。从研究结果看,NDVI空间特征与气候变化密切相关,气候的时空异质性对NDVI的变化趋势有十分显著的影响[18-23]。此外,长时序NDVI分析还被广泛用于分析土地利用演变,包括生态恢复工程带来的植被变化[24-27],以及城市扩张带来的植被退化等[28-33]。从时间尺度上看,除年际间长时序分析外[34],部分研究还分析了NDVI在季节和年内的变化趋势[35],并在小空间尺度上,结合不同物种的物候信息揭示气候和人类活动对植被生长变化的影响[36-37]。

现有的大多数研究中,对NDVI的趋势分析多针对时间和空间总体变化展开[34,38],无法描述精细变化特征,不能反映长时序数据中的波动特征,从而无法满足植被保护及恢复差异化决策和管理需求,使得长时序NDVI分析结果受到质疑[34, 39-40]。基于此,本研究建立了一种新的空间特征识别方法,利用2000—2018年MODIS NDVI产品数据集,将时序突变检测与趋势分析相结合,基于栅格尺度识别NDVI空间突变拐点,分析了突变拐点前后NDVI变化趋势与常规全时序变化趋势的异同,分析了两种趋势分析的空间特征差异,并探讨了趋势变化的影响因素。

1 研究区域

海河流域北部(东经115°30′至119°45′,北纬39°10′至42°40′)由滦河及潮白河流域组成,北接内蒙古坝上高原,南连京津冀城市群。区域内年均温1℃至11℃,年降水量400—700 mm,属暖温带半干旱-半湿润季风气候。滦河及潮白河流域面积89870 km2,其中上游山区面积占73.9%(66376 km2),跨河北、内蒙古、北京、天津、辽宁五省(市)。区域内西北部以草地覆盖为主,中部及东部以林灌覆盖为主,其中林灌草的覆盖面积占78.1%(图1)。研究区域内分布着众多针对森林、草地、地质、生物等国家、省、县级保护区,总面积4093 km2,占区域面积的6.2%。从行政区域来看,海河流域北部山区主要涵盖承德市大部分区域,以及北京、唐山北部,张家口西部,秦皇岛东部,是京津冀都市圈重要的生态涵养区；此外,潮白河及滦河还是首都水源地的源头和首都应急水源地所在地。区域内的植被具有十分重要的生态屏障作用,其空间变化特征不仅关系到本地生态系统健康,同时还会对下游都市圈的生态安全产生影响。

图1 研究区域土地利用及气象站点分布示意图Fig.1 Distribution of land use and meteorological stations

2 研究方法

为了揭示植被在不同时段内的空变化特征,从而精确评估变化趋势,为明确区域人类活动对植被的影响以及差异化的植被恢复策略提供支持。本研究构建了基于栅格尺度的NDVI时序突变检测方法,采用Pettit突变分析方法识别NDVI空间突变的拐点,采用最小二乘拟合计算突变拐点前后NDVI的变化趋势,分析各时段趋势变化空间特征,并与常规全时序NDVI趋势分析及时空变化特征对比,提出NDVI空间变化特征差异,并分析差异的影响因素(图2)。

图2 基于时序突变检测的NDVI空间特征识别框架Fig.2 Framework of NDVI spatial feature recognition based on temporal mutation detection

2.1 数据来源与处理方法

NDVI数据来自Google earth engine中MODIS13Q1数据集,空间分辨率为250 m,时间分辨率为16 d。本研究中使用数据的时间段为2000年至2018年。在提取数据集中的各期NDVI后,使用Savitzky-Golay滤波方法平滑数据并排除异常值,识别研究区域多年NDVI季节变化特征,采用每年4—9月生长季节均值作为NDVI的表征。

气象数据来自资源环境数据云平台,采用2000—2018年研究区域及周边59个国家级气象站的年均值,利用ANUSPLIN插值软件对年均气温和降水进行插值,空间分辨率为250 m,进而提取海河北部山区气象因子的空间分布。

2.2 突变分析及突变拐点识别

突变分析用于识别基于栅格尺度的NDVI在2000—2018年间的突变时间,采用Pettit突变分析方法。该方法不仅能够识别长时间序列中突变点的位置及数量,还能够判断突变点的统计意义是否显著,方法如下：

(1)

(2)

K=max|Ut,n|,t=1,2,…,n

(3)

(4)

式中,Ut,n为定义的统计量；n为时间序列长度；xj,xi为检验样本中的随机变量；K为突变年份；P为显著性水平检验值；在给定显著性水平α情况下,若P≤α,则认为该突变点满足显著性水平。

确定各栅格的突变年份后,统计各年突变栅格的数量,得到累积突变栅格比例,将累计曲线拐点对应年份作为时序分析的分界点。

2.3 趋势分析及分级

NDVI趋势分析采用最小二乘拟合计算植被变化趋势的斜率,计算方法如下：

(5)

式中,SNDVI为变化趋势,SNDVI>0表明植被变化呈现增加趋势,SNDVI<0表明植被变化呈现降低趋势,SNDVI=0表明植被变化无明显变化趋势；n为时间序列长度；i为研究年份序号；NDVIi为第i年的生长季平均NDVI值。从显著性水平来看,选取0.05显著性水平将NDVI变化分为3个等级(表1)。

表1 NDVI变化趋势等级

3 结果

3.1 NDVI突变分析及空间变化

图3 年度突变栅格及累计突变比例Fig.3 Annual mutation grids and cumulative percentage

基于栅格尺度的突变空间分析结果表明,显著突变(P<0.05)年份的空间分布呈现出较为明显的空间异质性：初期(2004—2008年)突变区域较小,主要集中在区域东南部；中期(2009—2011年)突变范围急剧扩大,由东南向西北覆盖绝大部分区域；后期(2012—2014年)突变区域减少且向西北部集中。从面积比例上看,NDVI存在显著突变出现的时间为2004—2014年,其中突变集中的时间为2009—2011年(图3)。从突变累积曲线看,2011年显著突变的栅格出现拐点。至2011年,研究区域内累积共有69.1%的面积发生了突变,占突变面积的88.7%,基于以上结果,本研究以2011年作为分界点,分析2000—2011年和2011—2018年NDVI的变化趋势,并将其与2000—2018年的全时序趋势分析结果进行对比。

3.2 NDVI变化趋势及其与气象因子相关性分析

图4为海河北部山区不同时间段NDVI的变化趋势。对比全时序和突变拐点前后的结果可以看出,全时序变化趋势中,区域整体呈现改善趋势(87.2%),退化区域集中中部以及东南部(1.7%)。突变拐点前的趋势变化中,呈现改善的区域分布在中部和东部(占54.3%),退化区(占1.6%,)主要分布在西北部草原区(正蓝旗、克什克腾旗、多伦县),中部承德市周边(双滦区、双桥区),东南部唐山北部(遵化市、迁西县、迁安市、滦平县)。突变拐点后,改善区域主要分布在北部和西部(占52.7%),退化区(占1.2%)主要分布在东部宽城县、青龙县、以及迁西县和迁安市北部。

图4 海河北部山区NDVI变化趋势Fig.4 Classification of NDVI variation trend in the mountain area of Northern Haihe Basin

对比不同时间段NDVI变化趋势的重叠率可以看出(表2),突变拐点前后显著改善和基本不变的区域与全时序分析的重叠率较高(均大于50%),显著退化区的重叠率较低,特别是2011—2018年间,显著退化区域与全时序分析的重叠区仅为15.8%。这表明基于时序突变识别的变化趋势与全时序分析的识别结果存在空间差异,特别是显著退化区的识别结果差异最大,而退化区正是植被恢复并需要人类正向干预的重点区域。

表2 突变前后NDVI变化趋势占全时序分析比例

图5 海河北部山区NDVI与气象因子相关性Fig.5 Correlation of NDVI to meteorological factors in the mountain area of Northern Haihe Basin

对比全时序和时序突变拐点前后时间段,NDVI与气象因子相关性也存在明显的空间异质性(图5)。气温与NDVI总体上呈现负相关,其中突变拐点前呈现负相关的面积最大,占77.9%,突变后负相关面积降至50.6%,集中在西部和东部；相比之下,全时序分析的负相关面积比例为72.3%,显著负相关占比(11.7%)大于突变拐点前后(1.7%和3.4%)。三种时段中,降水与NDVI的正相关面积比例均在75%以上；突变后,呈现负相关的面积有所增加,占总面积的24.46%,主要分布在西部赤城、丰宁县和东南部遵化、青龙、宽城等地。结合NDVI变化趋势来看,不同时段内,气象因子对NDVI显著改善的作用并不大,但却是突变后中东部局部呈现退化的原因之一。

3.3 典型退化区时序突变与传统检测比较

从识别结果来看,退化主要集中在三个区域(图6),分别为西北部草原区(包括多伦县、围场县、以及正蓝旗和克什克腾旗部分区域),中部林灌区(包括隆化县、滦平县、承德县、双滦区和双桥区),和东南林农区(包括宽城县、兴隆县、遵化市、迁西县、迁安市、以及滦县部分区域)。

图6 典型退化区空间分布Fig.6 Spatial distribution of typical degradation area

从面积上看(表3),三个典型退化区的面积占全时序分析结果的67.0%,其中中部林灌区的占比最大,为43.1%,其次为东南林农区(15.1%)和西北草原区(8.8%)。对比两种趋势分析结果,退化区的分布存在明显的空间异质性。突变拐点前后显著退化区的比例与全时序突变趋势的重叠比例均低于50%,特别是2011—2018年,平均重叠比例仅为16.3%。结合区域气候变化及经济发展政策来看,气候变化带来的退化影响在2000—2011年主要影响西北部草原区,在2011—2018年主要影响东南林农区；城镇扩张导致的退化主要是在2011—2018年间影响中部林灌区和东南部林农区。

表3 典型退化区面积比例

4 讨论

从面积上看,基于时序突变分析与常规全时序分析在海河北部山区NDVI趋势空间特征上存在明显差异。常规全时序趋势分析方法中,显著改善的区域占80%以上,退化区域占1.7%；基于时序突变分析方法中,突变拐点前后时段内,改善区面积均下降了30%左右,基本不变的区域增加了45%左右,退化区占1.2%—1.5%左右。从空间分布上看,2000—2018年区域NDVI整体呈现显著上升的趋势,植被处于改善状态,研究结果与海河及我国北方2000年以后NDVI的变化趋势一致。何航等[38]对中国北方NDVI变化特征的研究中表明,在气候变化背景下,1982—2015年间华北地区植被程良性改善状态,显著增加面积大于40%。王永财等[41]针对海河流域1998—2011年不同土地覆被类型的NDVI变化趋势进行了分析,显示各类土地覆被的NDVI均呈现上升趋势,其中森林、农田、湿地、草地的上升趋势明显。基于时序突变趋势分析结果表明,虽然改善区仍占主导(>50%),但面积显著减少,突变拐点前后空间格局差异明显。突变前,改善区主要分布在中部和南部,突变拐点后主要分布在北部和西部；退化区则由突变前的西北和东南部变为突变后的东部地区。对比前人的研究成果,基于时序突变的分析方法准确的识别了突变拐点前锡林郭勒东部草原草甸区的退化[42],以及突变拐点后唐山北部林农区的退化[43]。

NDVI的变化受到多种因素的共同影响,主要包括气候变化和人类活动两方面。在全球暖湿化的总体趋势下,气象因子对NDVI的影响一直是研究的热点。结合区域气象因子变化,海河北部山区NDVI变化趋势与降水的时空变化密切相关。突变拐点前,降水下降区分布在西北部草原区,突变拐点后则转移至东部林农区。从NDVI与气温的相关性变化趋势看,呈现负相关的区域主要分布在承德周边的林灌覆盖区,表明植被覆盖对局地气温起到了一定调节作用。相比突变前,突变后NDVI与气温呈负相关的区域减少(由77.9%变为50.6%),但呈现显著负相关的区域增加(由1.7%变为3.4%),表明虽然植被对气候调节作用的面积有所减少,但强度依然明显。从相关性分析来看,气象因子对NDVI的变化的空间异质性十分明显,呈现促进作用的区域集中在承德西部、北京北部和张家口赤城县,主要源于该区域在2011年后降水及气温呈现整体增加趋势,有利于植被的生长。呈现负面作用的区域集中在东部(迁西、遵化、迁安、青龙和宽城),其原因可能主要由于东部在2011年后降水呈现下降趋势,其中迁西、迁安和青龙等地呈现出显著下降趋势。从不同植被类型的面积来看,森林仍然是受影响最大的区域,其次为灌丛和耕地；从植被类型的比例看,55%以上的森林和灌丛、47.0%的耕地及27.7%的草地受到了气象因子的正向影响,仅有4.5%以上的森林灌丛、9.9%的耕地和10.7%的草地受到负向影响。相比气象因子,人类活动是海河北部山区NDVI整体呈上升趋势的主要因素,这主要是源于近几十年来在我国北方开展的生态造林及森林保育工程。从研究区域整体来看,基于京津风沙源和退耕还林还草等工程,2000年至2015年间,区域内林灌由41868 km2增加至42086 km2,草地由9598 km2增加至9,752km2。统计表明,承德市自2002年京津风沙源治理工程建设以来,共造林61.3万hm2,截至2015年,承德市的森林面积占京津冀三省市林地面积的36%[44]。张家口市辖赤城县在2001—2013年间通过飞播或人工种植共造林4.0万hm2[45]。随着工程后期植被恢复增速的降低以及区域发展等因素,人类活动对NDVI带来的负面影响也逐步体现,特别是2011年“首都经济圈”概念提出后,唐山、承德等城市周边的城市化进程也直接导致了2011—2018年中部承德周边和东南部唐山-秦皇岛周边NDVI的退化。综合来看,海河北部山区的植被整体改善主要是由于人类活动的影响,而突变点前后局部植被退化的空间异质性则是气象因素和人类活动共同作用的结果。

从研究结果看,常规全时序分析能够识别总体趋势变化,但却无法识别由于数据区间波动带来局部时空特征,无法准确评估生态系统的退化状况和识别退化区域,从而不能准确判别退化的影响因素,无法为生态工程效果评估和后续政策的调整提供科学依据。基于时序突变分析方法在识别NDVI分时段趋势的基础上,能够揭示植被空间变化的动态趋势,为明确气象因素和人类活动对植被的影响(特别是负面影响),并为制定差异化修复策略提供了有效的决策信息。

5 结论

(1)结合归一化植被指数NDVI,建立了基于时序突变检测的植被空间变化特征识别方法,并与常规长时序趋势检测方法及结果对比；结果表明,基于时序突变检测的方法能够识别NDVI精细的空间变化特征,从而为后续修复提供完整有效的信息。

(2)以海河北部山区为例,通过时序突变分析表明海河北部山区NDVI在2011年存在突变拐点；突变拐点前后NDVI变化趋势均以改善为主,但退化区域呈现出空异质性,突变拐点前的退化区主要位于西北草原区,突变拐点后的退化区主要位于中东部林农区。

(3)分析表明,区域内植被整体改善主要源于人类活动的影响(生态工程等),气候因素的作用有限；但区域内突变拐点前后植被局部退化的空间变化特征则受到气候因素和人类活动共同影响。

(4)基于时序突变分析方法与常规长时序分析方法的评估结果存在空间差异,特别是退化区域的识别结果差异显著；对比前人研究成果,基于时序突变分析结果能够更为准确的识别区域生态系统的变化及影响因素。