基于LandTrendr模型的亚热带森林干扰与恢复动态变化分析

刘姗姗，黄鑫毅，赵帅，陈耀亮*

(福建师范大学 a.地理学国家级实验教学示范中心,b.地理科学学院，福州 350007)

自20世纪80年代以来，中国亚热带地区森林覆盖率显著增加，较小的林龄结构和充沛的雨热条件使该地区森林有可能成为全球较大的碳汇区[1-2]；然而，持续的森林干扰及其引起的年龄效应可能对森林碳收支产生重要影响。森林资源状况是生态系统的重要评价指标[3]，森林生态系统通常会受到诸如森林干旱、火灾和砍伐等来自自然或人类的干扰，对森林生态系统的结构和功能产生重要影响[4]。揭示和掌握森林扰动及恢复的时空变化特征，可为准确评估区域森林发展态势、实施有针对性的森林经营规划和管理提供相应的理论依据和技术支撑，对提高亚热带森林的抗干扰能力、增强亚热带森林的自然恢复能力具有重要的借鉴意义。传统的森林资源及其变化调查监测以地面测量为主，存在工作量大、劳动强度大、成本高、周期长、效率低和时效性差等问题[5]，且其调查精度低，难以满足当今森林资源变化监测的需要。遥感以其大区域、时效性好等特点成为森林干扰监测的重要技术手段。主要方法有：分类后比较法[6]、影像差异法[7]、分类及统计分析法等[8]。近年来，有学者基于Landsat长时间序列数据提出了时间序列分析法进行森林干扰监测，包括基于检测时序光谱指数相对于稳定状态的突变方法和基于拟合算法区分森林干扰信号与背景噪声方法[9-11]。与传统干扰监测方法相比，时间序列分析法能够确定森林干扰发生的年份、持续时间、干扰强度等信息，有效地监测森林的长期变化状况[12]。时间轨迹拟合干扰算法主要包括植被变化追踪VCT[13]、季节与趋势断点算法BFAST[14]和基于光谱轨迹的LandTrendr[15]。VCT能够较好地检测出森林转变为其他地类的变化，但不能有效地检测间伐、森林退化等干扰变化；BFAST算法是基于植被的物候模型，需要对每年四季多期影像进行曲线拟合，对影像的要求较高，因此在云覆盖度高的亚热带区域难以实现；LandTrendr方法能够很好地对森林轻微干扰进行检测，这些森林干扰监测方法主要基于Landsat和MODIS时间序列数据。

本研究基于长时间序列的Landsat TM/ETM+/OLI影像数据集，采取LandTrendr时间分割算法[15-17]对研究区内的森林干扰及恢复情况进行监测，并对研究区内从沿海至内陆地区的森林干扰及恢复动态变化进行分析。

1 研究区概况

研究区域地处浙江省中东部沿海，囊括磐安、天台、三门、仙居和临海等5个县的县域范围(图1)。区域地属中亚热带季风气候，冬夏季风交替显著，四季冷暖干湿分明，光照适中，雨量充足，主要气象灾害有春季低温阴雨、五月寒潮，夏季的台风、冰雹、高温干旱，雨季洪涝等。区域内森林植被类型丰富多样，主要分布着阔叶林、针叶林、针阔混交林、毛竹林等。森林生态系统易受多种因素影响，浙江省中东部是中国典型的亚热带森林区域，地形复杂多变，濒临东海，受到当地气候要素及海陆位置的深刻影响，夏季高温、台风等自然灾害频发，同时人类经济活动对森林的干扰也比较频繁，因而其森林干扰特性较为复杂。

图1 研究区：浙江省中东部沿海区域 Figure 1 Study area: Eastern coastal areas of Zhejiang province

2 研究材料和方法

2.1 研究数据

选取1984—2015年间共15期轨道号(path/row)为118/040的Landsat TM/ETM+/OLI影像数据为研究数据。影像选取原则是：(1)尽量获取在植被生长茂盛期(6—9月)的影像，以减少物候对植被光谱识别的干扰；(2)尽量选取云量少的影像，以保证时间序列内有相对较高的影像质量。

收集1994、1999和2004年共3期森林资源清查数据，并对样地数据进行筛选及分类：(1)筛选3期均有的固定样地；(2)根据任意2期样地的植被类型、优势树种和林龄等信息将样地分为干扰、恢复及稳定3类。共筛选研究区内94个样地，并作为验证数据。

2.2 影像预处理

影像预处理主要包括辐射定标、大气校正、云掩膜、森林信息提取等过程。

2.2.1辐射定标及大气校正

在Linux环境下，运行LEDAPS模型算法完成影像的辐射定标及大气校正。LEDAPS[18]是由美国航空航天局(NASA)针对长时间序列Landsat影像数据集的预处理而开发的模型。该模型对Landsat数据的预处理包括辐射定标和6S辐射传输模型大气校正[18]。模型的输出是地表反射率数据。

2.2.2云及阴影的掩膜

采用Fmask算法[20-21]对时间序列影像进行云及云影的掩膜，以消除云及阴影对时间分割算法的影响。

2.2.3森林信息提取

通过剔除水体、不透水地表及农田，确定森林区域。综合考虑31年间森林扰动或恢复等对森林面积的影响，分别对1984、2015年的影像作森林信息提取，最后取二者并集作为研究的森林区域。

表1 主要研究数据Table 1 Data sets used in research

2.3 研究方法

采用LandTrendr时间分割算法[15]对影像进行处理，得到研究区31年间的森林植被干扰与恢复产品。选用归一化燃烧比(NBR)指数为算法监测指数，其计算公式为：

(1)

式(1)中：NIR为近红外波段(850～880 nm，TM/ETM+影像中的Band4，OLI中为Band5)，SWIR为短波中红外波段(1 570～1 650 nm，Landsat影像中为Band7)。

在对Landsat时间序列影像掩膜及标准化处理后，被掩膜的像元将不参加时间序列分析，轨迹拟合时将根据分割前后的指数值进行插值填补。LandTrendr算法得到的结果见表2。

表2 LandTrendr主要输出结果说明Table 2 The main output of LandTrendr algorithms

像元的整体变化趋势包括恢复、干扰、稳定以及平衡4种状态。恢复指变化轨迹趋势是整体向上的；干扰是指变化轨迹是整体下降的；稳定则指整个时间长度NBR未发生变化或变化微小；平衡指变化轨迹波动大，上升分割段与下降分割段呈现抵消的态势。

2.4 精度验证

采用森林资源连续清查数据进行验证。选取1994、1999、2004年3期清查中均有调查的共94个固定样地，根据样地的地类、优势树种及林龄等信息，判断每2期清查样点的变化情况(干扰、恢复或稳定)，再与LandTrendr时间序列分割结果做验证分析。

2.5 分类统计

2.5.1森林干扰与恢复等级划分

根据对研究区林业与气候资料的收集与分析，该区域内的森林扰动因素主要有人为砍伐、火灾、极端高温干旱、极端低温冰冻、台风等。为进一步分析该区域内的干扰程度，按茂密森林变为非森林(人为砍伐，林地变为裸地、建筑物等)、森林植被严重损害(高温干旱等)、森林退化(长期缓慢变化)等变化分为严重干扰、中度干扰、轻微干扰3个等级。同理，对研究区内的恢复变化划分为重度恢复(非森林-茂密森林)、中度恢复(疏林地-茂密森林)和轻微恢复(林分自然生长)等3个等级。

以NBR指数作为等级划分的指标，基于2015年夏季影像在研究区内选取大量森林(茂密森林、疏林地等)与非森林(裸地、建筑物、水体等)样本，计算其NBR值域区间，并结合不同干扰事件(如人为砍伐、极端高温干旱等)所造成的NBR值变化的范围，最终确定干扰等级划分的阈值区间。同理可得到恢复等级划分的阈值区间(表3)。

表3 森林干扰与恢复等级定义Table 3 Class definition of forest disturbance and recovery level

2.5.2森林干扰与恢复总体分布

在LandTrendr算法中，根据1984年和2015年像元的NBR值的变化情况，将其分为干扰、恢复、无变化和平衡(干扰后又恢复到原始状态)4种状态，并计算出像元在整个变化轨迹中处于该状态的总年数。由此可得到1984—2015年间研究区内的森林干扰与恢复总体时空分布格局。

2.5.3森林最大干扰量与恢复量的计算

在LandTrendr算法中，对像元31年间的变化轨迹进行分割，每个分割均代表像元的某一种变化，即像元可能有多个或重复变化，如干扰-恢复-干扰等。因此，为了更直观地展现干扰与恢复的时空分布格局，本研究中筛选出每个像元变化轨迹中NBR变化量(mags)最大的分割，并以此对像元的干扰或恢复属性进行划分。最后按照表3中的阈值对其森林扰动和恢复程度进行划分。

此外，找出像元每个影像年份所在的分割，以此分割的属性判断其处于干扰或恢复状态，并对其干扰或恢复等级进行划分。由此得到研究区每个影像年份的森林干扰与恢复面积并分析其变化特征。

3 结果与分析

3.1 森林扰动及恢复结果精度评价

对94个森林清查固定样地进行验证，得到精度评价表(表4)。清查数据中共有20个发生过干扰的样地，69个恢复样地及5个稳定样地。从表4中可看出动态检测的总体精度为74.47%，干扰的生产精度为50.00%，恢复的为79.71%，稳定的为100%。其中，干扰的生产精度与其他二者相比偏低，结合清查数据与LandTrendr处理结果分析可以判断，可能是样地遭到轻微干扰(如间伐等)但林分整体还是处于恢复状态只引起其恢复速率变缓，因而造成错分。恢复与稳定混淆严重，监测的稳定结果中有近60%的稳定样点在清查数据中表现为恢复，这是由于在实际判读的森林清查数据中，林分的林龄持续增加的这一状态是被认为处于森林恢复状态，但林分在光谱信息中可能未有明显表现，因此在LandTrendr中会被监测为稳定。

表4 森林清查数据的验证精度 Table 4 Validation of the forest inventory data accuracy

3.2 森林总体变化空间格局

图2为研究区的森林总体变化空间格局。经统计，研究区内74.62%的森林均处于恢复状态，1.42%的森林总体呈现出干扰变化，23.81%的森林一直未有变化，其余森林处于干扰与恢复达到平衡的状态。恢复的森林中，以20～31年的恢复为主，占86.6%，说明大部分的森林均未受到干扰或受到了短期轻微干扰；其次是1～10年的恢复，占9%，表明森林为2003年后发生的干扰后恢复的过程。由于浙江省内城镇建设扩张，用地需求大增，仙居县等多县均出现了“毁林造地”等现象，根据监测结果可知森林干扰多分布在不透水地表周边，呈破碎化分布，且以短期(1～5年)干扰为主，占38%；其次是中长期(10～20年)干扰，占29.31%，说明经济发展、城镇扩张等人为因素是导致森林干扰的重要因素。

图2 研究区的森林总体变化空间格局Figure 2 Overall changes of forest spatial pattern in the study area

3.3 森林干扰量级时空格局

研究区的森林干扰量级时空分布格局如图3。经统计，研究区内的最大森林干扰主要发生在1993—2004年间，占50.7%;其次是发生在2004—2015年间，占30.55%;1984—1993年间的干扰最少。该区域内的森林干扰以中度干扰为主，占47.68%，严重干扰占19.38%，且以短期干扰(1～5年)为主，占比达92.80%。

图3 研究区的森林干扰时空格局Figure 3 The spatial and temporal pattern of forest disturbance in the study area

各县区的最大干扰事件以中度干扰为主，仙居县和临海市的严重干扰面积最大，临海市的中度干扰面积最大。最大干扰事件中，超过90%的干扰事件持续时间均集中在1～5年，5～10年的干扰事件中以仙居县最多，占比4%左右，10～20年的干扰事件中以临海市为最多，占该市森林干扰事件中的4.5%。

图4为研究区各年份的干扰面积。1992、2003、2007、2013等年份各县区的干扰面积最大，且在2000年后，干扰面积逐年增大。其中，仙居、临海、天台、磐安等均在1992年达到第一个峰值，在2000年后均大幅度增长，在2013年及2015年达到最高峰。结合台州市统计年鉴信息可知，在1991年、2002—2004年间、2007年、2013年前后等均有夏季极端高温干旱气象灾害发生，同时研究区内干扰面积大幅度升高，可认为这些年份的森林干扰主要是由高温干旱引起的。2000年以后受严重高温干旱影响的森林干扰面积远大于20世纪90年代；同时2000年以后无高温干旱灾害的年份的森林干扰面积大于20世纪90年代的森林干扰面积，这表明人为干扰因素已成为森林干扰的重要因素。

图4 研究区各年份森林干扰面积Figure 4 The forest disturbance area of each year in the study area

从图4中可以看出，研究区内的森林扰动面积在2013—2015年呈现跳跃式增长，这是由于2013年影像为当年7月至9月所有无云区域拼接而成的影像，影像的拼接处由于光谱信息差异较大，出现大面积的干扰噪音信息，导致结果中出现2013年的干扰面积大幅度增加。

3.4 森林恢复量级时空格局

图5为研究区的森林恢复量级时空分布格局。该区域的最大森林恢复主要发生在1984—1993年间，占96.99%;其次发生在1993—2004年间，占2.54%;2004—2015年间的恢复最少。该区域内的森林恢复以中度恢复为主，占49.32%，轻微恢复占46.8%，且以长期恢复(20～31年)为主，占比达91.06%，其次为中长期(10～20年)恢复，占6.66%，表明该区域内的森林大部分未受到干扰或干扰后恢复较好，处于良好的恢复状态。

图5 研究区的森林恢复时空格局Figure 5 The spatial and temporal pattern of forest recovery in the study area

各县区的恢复状态以中度恢复为主，仙居县和临海市的面积最大。此外，磐安和三门县超过92%的森林恢复持续时间集中在20～31年，说明磐安和三门县的森林生长恢复状态良好；仙居县和临海市超过2.5%的森林恢复持续时间集中在1～10年，说明仙居和临海的森林生长恢复状态不佳。在森林恢复等级上，磐安县近60.2%的森林恢复处于轻微恢复状态，即自然生长状态，三门县的中度恢复和最大恢复比例均为五县区中最大，分别为58.77%和5%，表明三门县的森林恢复力度大，人工造林成效显著。

研究区各个年份的恢复面积如图6。在1984—1998年间，各区县森林恢复面积基本保持稳定水平，其中天台县、磐安县、仙居县森林恢复总面积最低值出现在1992年，而三门县和临海市最低值则出现在1987年，但其在1987年与1992年的森林恢复总面积相差不大；而在2000年后，森林恢复面积相较于20世纪90年代有一定程度下降，森林恢复面积变化状况基本与森林干扰面积的变化状况相对应。

图6 研究区各年份森林恢复面积Figure 6 The forest recovery area of each year in the study area

4 讨论与结论

4.1 讨论

LandTrendr算法主要思想是从时间序列数据中提取地表光谱变化轨迹，对时间序列轨迹做适当分割及线性拟合，去除影像光谱噪音信息，将复杂的变化特征简化为相连的几段直线以突出要点，从而捕获时间序列数据的干扰特征，得到干扰年份、干扰量及恢复信息。LandTrendr算法在监测干扰年份、干扰量的基础之上, 能够监测特定年份内发生的干扰以及多次的扰动变化，从而能够识别主要与次要干扰的分布[22-23]。LandTrendr算法相较于其他时间轨迹拟合干扰算法可操作性更强，同时也更易识别出轻微干扰。

本次研究利用LandTrendr算法对浙江省中东部地区进行森林干扰和恢复监测，采用的是NBR指数，该指数常用于北美等区域的森林干扰监测研究中，此次监测结果的总体精度为74.47%，恢复与稳定监测的生产精度较高，分别为79.71%和100%，然而干扰的生产精度仅有50%，监测精度低。LandTrendr算法中提供多个控制参数的阈值可调试以达到更好的森林干扰监测效果，本次研究只对其中一部分参数进行了调整，这可能会影响森林干扰与恢复的监测精度。此外，基于遥感影像数据构建的时序地表光谱变化轨迹受到影像质量的制约，研究区遥感影像云量多，一是无法满足一年一景或隔年一景的高质量影像的需求，研究所采用的影像多数隔2～3年一景，从而导致森林干扰信息被过滤，监测精度降低；二是2013年影像为当年7月至9月所有无云区域拼接而成的影像，物候差异和拼接边界云及阴影遗漏严重，给监测带来许多噪音，导致监测结果中2013年森林干扰面积异常增高。在精度验证方面，国家森林资源连续清查数据只有1994—2004年间的数据，每隔5年一期，这也导致许多森林干扰信息被过滤，难以与时许遥感数据的森林干扰监测结果相互验证，从而影响监测结果精度。

在今后的研究中，可以采用波段、其他指数或多种指数进行监测，对比得到最适用于监测中亚热带地区森林干扰与恢复的方法，此外可将各年份的夏季影像无云部分进行组合，保证每年一景或多景影像，进行森林干扰监测，并对LandTrendr算法监测指数及控制参数阈值等进行调试，以提高中亚热带森林干扰监测精度。同时，结合原始时间序列影像和时序高分影像对多年的森林扰动与恢复监测结果进行人为判读，建立混淆矩阵，弥补国家森林资源连续清查数据的时间滞后性[24]。

4.2 结论

研究表明，虽然临海市与三门县均为沿海区域，自然灾害频发，但近31年来，由于三门县沿海防护林建设等工程的积极开展[25]，有效减少了人为森林干扰事件的发生，自20世纪90年代起，三门县相较于其他县市，森林干扰面积较低；临海市由于人口密度大、经济发展相对发达，其森林干扰频繁且干扰程度重；仙居县与磐安县离海域较远，但仙居县人口密度比磐安县大，且经济发展较为发达，近31年来森林干扰面积大且干扰程度重，磐安县的森林干扰相对较少；5个县区的森林恢复面积基本处于稳定状态，波动幅度小，绝大部分的森林均处于自然恢复状态。

由以上结果可看出，沿海区域的森林虽比内陆更易遭受台风等自然灾害侵袭，但随着森林保护力度的加强，森林抗干扰的能力更强，如三门县整体呈现出的森林扰动变化较少，在内陆森林区域，由于森林保护力度低，其森林干扰程度反而比沿海区域的大。因此，除自然因素外，经济发展导致的森林扰动或是政策驱动下的森林恢复等人为因素是森林干扰与恢复的重要驱动力。在后续的研究中，可对研究区每年的森林干扰量及森林干扰驱动力等进行分析。