基于Landsat时序数据的森林干扰监测*

2020-06-15陈芸芝汪小钦

林业科学 2020年5期

钟莉陈芸芝汪小钦

(福州大学空间数据挖掘与信息共享教育部重点实验室卫星空间信息技术综合应用国家地方联合工程研究中心数字中国研究院(福建) 福州 350108)

森林是陆地生态系统的重要组成部分，是地球上最大的碳储存库，其动态变化直接影响全球陆地森林生态系统的碳循环和碳蓄积(葛全胜等， 2008；赵金龙等， 2013)。森林动态变化包括森林突发性变化和植被时序趋势变化。森林干扰是森林突发性变化，根据干扰的起因，通常可分为人为干扰和自然干扰，其形式主要表现为森林火灾、森林采伐以及因受雨雪和病虫灾害等所引起的森林改变(Gowardetal., 2013)。近年来，随着经济发展，全球各地森林植被频繁遭到严重破坏，人为干扰(包括森林采伐、城市化、改林为耕)也逐渐增加(李维长， 2000)。

遥感技术的不断进步，卫星遥感数据作为时间尺度变化监测的重要手段越来越受到研究者重视(Coppinetal.， 1996)。森林干扰研究的数据源通常分为基于区域和全球尺度遥感数据2类，其中，AVHRR、MODIS和SPOT VEG等数据空间分辨率较低，多用于大尺度的森林变化监测(宋富强等， 2011；王强等， 2012)；相比而言，Landsat系列影像数据具有较高精度，多用于区域尺度的森林干扰监测及其造成的森林碳通量变化研究(杨辰等， 2015)。美国地质调查局于2008年免费开放了Landsat堆栈数据(Woodcocketal.， 2008)，该数据具有良好的稳定性和持续性，基于Landsat时序数据的研究取得了长足发展(Townshendetal.， 2012； Wulderetal.， 2015)。以往研究中，森林变化监测通常是对影像分类后进行两两比较(Cohenetal.， 2010)，不仅费时费力导致效率低下，而且针对长时序分析精度不高，不能满足应用需求(杨辰等， 2013)；影像差异法利用差值、比值、典型相关性变化等方法构建差异影像，虽然消除了分类后带来的误差积累，但是受辐射差异影响较大，除进行严格的辐射归一化外，只能通过选择同一传感器和同一季相数据以减少噪声(Deroseetal.， 2011； Maseketal.， 2008)。干扰监测中，运用几期影像很难区分突发干扰和背景噪声，而时间序列分析中，可利用背景噪声和突发变化的时间序列轨迹不同将其分开(Hostertetal.， 2003)。国外学者利用时间序列影像联合分析，通过提取时序特征，进而识别森林干扰，大大提高了工作效率和监测精度(Huangetal.， 2009； Jamalietal.， 2015； Kennedyetal.， 2010； Sulla-Menasheetal.， 2014； Verbesseltetal.， 2010)。

福建是我国南方重点林区之一，森林覆盖率约66%，居全国首位。森林在全球碳循环和碳管理中占有重要位置(Faheyetal.， 2010)。本研究以福建省长汀县为研究区，基于Landsat时序数据，采用时间序列轨迹分析(LandTrendr)方法进行森林干扰监测，以期为陆地森林生态系统碳循环和碳蓄积以及气候变化研究提供参考。

1 材料与方法

1.1 研究区概况

长汀县地处福建省西部，25°18′40″—26°02′05″N，116°00′45″—116°39′20″E，南与广东省近邻，西与江西省接壤，呈四周凸起、中间下凹的盆地特征，地形破碎，岭谷相间。该县处中亚热带季风气候区，平均气温18.3 ℃，无霜期年均260天，四季分明。花岗岩发育的红壤抗蚀能力弱，原生植被多遭破坏，现有植被主要包括马尾松(Pinusmassoniana)、灌丛和荒草坡在内的次生植被和人工植被。研究区地理位置见图1。

1.2 数据源与数据预处理

选用长汀县2000—2016年15期Landsat影像数据(表1)，条带号121/42，数据来源于美国地质调查局(USGS)。所选影像尽量保持同一季相，以减少植被物候差异对光谱识别的影响。长汀县夏季雨水丰富，为获取高质量影像，将时间集中在9、10月，个别年份由于云覆盖等因素不满足要求，将时间间隔延长至2年。数字高程资料来源于地理空间数据云平台提供的全球30 m分辨率GDEM数字高程产品。

对每年影像进行辐射定标和集成于ENVI的基于MODTRAN4+辐射传输模型的FLAASH大气校正，得到长汀县2000—2016年时序地表真实反射率影像。为消除云和云阴影对结果的干扰，利用Fmask算法(Zhuetal.， 2012)逐一去除云像元。对掩膜后留下的空值采用最邻近年份的无云像元时间内插方法进行填充，对于时间序列影像的第1年和最后1年，使用最邻近年份的无云像元替代：

xi=xp+(i-p)(xn-xp)/(n-p)。

(1)

式中：xi为第i年的填充值；xp和xn分别为第p年和第n年的指数值。

图1 研究区地理位置

表1 遥感影像数据介绍

1.3 光谱指数选取

选用归一化燃烧率(normalized burn ratio，NBR)构建时序数据进行干扰监测。国外学者(Kennedyetal.， 2010； Cohenetal.， 2010)基于大量数据研究表明，NBR对多种类型干扰事件具有最大的敏感性。NBR由对植被叶绿素含量较敏感的近红外(NIR)和对水汽含量较敏感的短波红外(SWIR)组合构成，健康植被因高NIR值和低SWIR值，使其具有较高的NBR值，而受干扰的像元由于呈现土壤特征，SWIR值增大，使NBR值降低。据此，可根据时序轨迹中NBR值的突变识别森林变化信息：

(2)

1.4 时间序列轨迹拟合算法

LandTrendr算法(Kennedyetal.， 2010)的核心是将复杂的时序轨迹简化为一组相连的线段，并对分割后的分段进行线性拟合，旨在消除噪声，突出重要的信息。算法主要参数设置如表2所示，关键流程如下：

1) 时序轨迹提取及噪声消除逐像元提取时序轨迹(图2a)，并通过设置噪声值进一步去除预处理过程中残留的细小噪声(图2b)。

2) 时序轨迹分割对整个时间序列，首先将时序中的第1年和最后1年作为第1次分割的起始点，得到观测年份与指数值的一阶回归方程【图2c(1)】，并计算每个观测年指数值与预测值的绝对差，选取绝对差最大的一年作为下次分段的分段点，据此将时间序列分为2部分【图2c(2)】。然后对2部分分别计算回归方程和每段的均方误差(MSE)，选取MSE较大的一段参与下次分割，以此类推，直到分段数达到设定的初次分段数【图2c(3)～(5)】。为防止分段数过多造成过拟合现象，选用角度阈值进一步判断，使分段数达到最大分段数。通过计算对比原始时序轨迹中段与段的夹角，去除最“浅”角对应的分段【图2c(6)】。

3) 时序轨迹拟合联合基于点与点连线及基于回归连线2种方法进行时序轨迹拟合。对于第1个分段，分别计算2种方法的MSE，选取MSE较小的方法作为分段的连接线。从下一分段开始，每一个连接线必须与前一个相连，同样计算2种方法的MSE，选取MSE较小的方法连接，最后得到一条完整的拟合轨迹。

4) 模型简化采用恢复率阈值(1/n)剔除不合理的分段点，对上一过程拟合的时序轨迹通过迭代方式简化模型并重新计算拟合。在研究区内，绝大多数森林干扰为突发性事件，但植被恢复是个缓慢的过程,通过设置恢复率阈值(1/n)，将植被恢复时间小于n年的森林干扰视为噪声。根据长汀县植被恢复状况，将恢复率阈值设为0.5，恢复时间小于2年，即干扰发生1年后，植被恢复到干扰发生前状态的区域视为噪声。

5) 光谱滤波通过回归模型将NBR值转化为植被覆盖率进而生成干扰量，以植被变化强度阈值(表2滤波过程)为滤波条件，剔除原始分割算法中小的光谱分段点，以减少光谱异常信号或植被物候差异造成的时序轨迹过度拟合。

基于以上过程，可识别时间序列中NBR下降、上升和保持不变3种特性。经过噪声去除及对由光谱异常信号或植被物候差异造成的伪变化进行滤波后，NBR值下降即由森林干扰导致。

图2 LandTrendr分割流程

表2 LandTrend参数设置

1.5 森林干扰制图

制图前需先排除水体、裸地等非森林地物。基于NBR指数的时序分割算法获取顶点的湿度、亮度和绿度影像，计算每年每个像元的缨帽变换角(TCA)。TCA值越大，对应的森林覆盖率越高，选择平均值最高的TCA值，其对应的像元最像森林像元(Huangetal.， 2002)。通过TCA阈值提取森林像元，提取的像元在整个时序中被识别为持续森林。由于研究所用时间序列较短，可监测到的二次干扰区域很少，故只提取最大干扰作为监测结果。TCA计算公式为：

式中： TCG和TCB分为缨帽变换绿度和湿度分量。

2 结果与分析

2.1 森林干扰监测结果

图3为长汀县2000—2016年森林干扰分布，图中非森林和持续森林已被腌膜，以白色背景表示，其余不同颜色代表不同年份的干扰分布。

2.2 森林干扰时空分布特征

图4为长汀县森林干扰面积随时间的变化特征。经统计分析， 2000—2016年长汀县森林干扰总面积192.49 km2，平均每年受扰动森林12.83 km2； 2001年干扰量最小，不足1 km2； 2004、2008和2009年森林受干扰较为严重，均在30 km2以上，约占当年森林面积的1.3%，共占干扰总面积的50%，其中2004年干扰面积高达32.85 km2； 2003、2006、2007和2010—2011年干扰面积略大于10 km2，其余年份干扰面积小于或远小于10 km2。

图4 长汀县森林干扰面积随时间的变化特征

从长汀县气象局获取的数据得知， 2004年是长汀县近20年来森林火灾次数最多的一年， 2008和2009年也爆发了大面积森林火灾(江帆等， 2015)。2008年受“拉尼娜”特大雨雪冰冻灾害侵袭，南方森林资源损失惨重(曹坤芳等， 2010)。2002—2004年、2006—2008年干扰面积呈持续上升趋势， 2001和2005年干扰降低，特别是从2009年起干扰面积持续减少直到2016年近乎最低。2000—2016年森林干扰面积在个别年份波动较大，但总体上随时间呈下降趋势。

对每年影像进行森林提取，获得长汀县历年森林面积及所占土地总面积比例。由图5可知， 2004、2008和2009年森林面积下跌较大， 2009年起呈缓慢增长趋势，与扰动情况相呼应。对比图4和图5，2004、2008和2009年的森林干扰面积占当年森林面积的1.3%左右，其余各年均小于当年森林面积的0.6%。

图5 长汀县历年森林面积及所占总面积比例

由图6可知，长汀县森林干扰面积随持续时间呈极速下降趋势，后保持低值平缓，说明干扰持续时间保持较短的年数，主要为1～2年，少数为3年，其中发生在1年的干扰面积比例最大，达82%；发生在2和3年的干扰面积分别占13%和4%。结合Google影像目视解译，长汀县森林干扰主要是由森林火灾和人工砍伐造成的急剧干扰事件。

图6 长汀县森林干扰面积随持续时间变化分布

由图3各年份干扰分布可知，长汀县森林干扰面积具有较强的空间分布特征。2004、2008和2009年森林干扰斑块面积较大，主要集中在长汀县中东部非森林区域附近，受人为干扰影响。2008年除分布在非森林区域附近的干扰外，有大量因受雨雪灾害影响的细小斑块分布在长汀县西部。其他年份亦如此，主要分布在长汀县非森林区域附近，四周分布零散细小的斑块。图7清晰展现了长汀县森林干扰面积在各海拔段的分布情况，随海拔升高，干扰呈明显下降趋势，其中海拔380～460 m之间的干扰面积最大，超过60%的干扰发生在中低海拔地区，海拔大于940 m的干扰仅占0.45%。

图7 长汀县森林干扰面积随海拔变化分布

2.3 基于地面调查的精度评估

图8 长汀县森林干扰实地调查

结合森林干扰监测结果和地面调查资料对4个典型样区(样区分布见图3)进行验证，获取干扰发生的时间和位置。图8a-d从左到右分别为样区干扰年TM影像、监测结果、实地调查照片和NBR指数时序轨迹，实地调查照片均拍摄于2013年7月。当干扰发生时，NBR指数大幅度降低，之后逐渐增加，植被慢慢恢复。经目视解译并对照提取结果可知，干扰斑块可被完整提出，且边界准确清晰，细小干扰也能逐一识别。图8a、b干扰分别发生在2010和2008年，到2013年植被有所恢复但不是特别理想，从实地调查照片可知植被主要以灌木为主；从图8c、d可以看出，植被恢复效果良好。实地调查结果与研究得到的干扰监测结果一致，表明该方法适用于长汀县森林干扰监测，且具有较好的监测能力。

2.4 基于像元尺度的精度验证

采用分层随机抽样方法，以像元为单位提取森林干扰样本进行精度验证。选用每年作为一层(包括15期森林干扰区域，森林和非森林各为1层)，共17层。以最少像素层30为基准，对森林干扰层进行等比例随机抽样，15期干扰样本数量共9 075个，森林和非森林像元数偏多，适当减少，分别为3 200和1 000个，即共13 275个像元样本参与验证。验证结果统计为混淆矩阵，见表3。

由表3可知，时间序列轨迹算法可以很好提取长汀县2000—2016年森林干扰，总体精度达96.26%，Kappa系数为0.92，具有较高的监测精度。各年份用户精度均在80%以上， 2014和2016年相对其他年份出现一定漏提。2016年处于时序轨迹的末端且前1年影像缺失，加之目视解译发现该年干扰范围和斑块面积较小，分布零星，导致2016年出现45%左右的漏提。

3 讨论

LandTrendr算法根据数据本身决定变化轨迹，利用干扰与背景噪声不同，通过参数控制进行分割拟合，消除时序中的噪声和小的光谱变化点，以突出重要的信息，最终得以识别干扰。该方法不仅可提取多种类型的干扰，还可监测同一区域多次扰动变化，同时获取植被干扰和恢复相关信息。

由于云和云阴影的影响，不能保证每年都有满足要求的高质量影像。某些年份影像缺失，可能导致干扰监测时间延误1年，或某些区域出现植被恢复过快现象导致轻微漏提。后续工作可重构高时间分辨率影像，降低因影像缺失造成的漏提、误提现象。经目视解译可知，长汀县森林干扰类型主要为森林火烧和人工采伐，本研究所用方法虽可提取多种类型的干扰，但未能划分出具体的干扰类型，后续工作可结合其他技术方法分别提取出火烧迹地和采伐迹地。