干扰是森林生态系统动态变化的主要驱动力, 干扰的历史影响林分的生长状态, 不同干扰的类别、强度和大小能改变林分物种组成和林分结构^[1-3]。典型的自然干扰 (雨雪灾害、病虫害等) 和人为干扰 (采伐、土地利用变化等) 以及干扰后更新事件将影响森林碳汇^[4]。不同景观变化下 (如城市化、气候或灾害变化及森林管理活动等), 干扰的状态受不同程度的影响^[5]。森林采伐作为人为干扰及景观动态的重要部分, 一方面影响着森林碳储量, 另一方面它与森林覆盖度、物种分布、地形及社会经济因素等密切关联^[6]。自从1980年开始, 大尺度的造林和再造林计划的实施使得中国具有了最大的人工林分布面积, 占全球人工林面积的1/4^[7]。人工林的存在使得森林碳储量得以增加^[8]。在南方区域, 造林和再造林活动是区域碳汇变化的主要原因之一^[9]。人工速生丰产林建设承担着缓解木材需求, 发挥改善经济、社会、生态效益的重要作用, 被列为国家产业政策, 具有举足轻重的作用^[10]。但是人工林普遍存在树种单一、结构和功能较差、病虫害频繁发生的特点。另外, 人工林营造时树种选择主要考虑速生、丰产、短轮伐期的特性^[11-12], 这使得中国南方区域的森林存在快速、广泛的时空变化。发展科学可靠的方法来监测森林干扰或者恢复时空模式对于人工林健康经营管理及其碳核算意义深远^{[9, 13]}。

然而, 仅凭借两期或者三期长时间间隔的影像^[14-15]往往不能充分捕捉南方人工林频繁发生的采伐及伐后更新活动。近年来基于30m分辨率的Landsat卫星遥感技术的发展为森林干扰的监测提供了数据支撑, 特别是时间序列的遥感数据被成功用于森林干扰或恢复变化监测^[16]。但是如何能快速有效地识别出不同持续时间和地域分布的森林干扰或恢复现象仍然是项技术难题。尽管一些森林变化自动化分析算法^[17-21]为森林生态系统的近实时监测提供了强有力的工具, 但是很少有方法既能监测短期的森林变化又能探测持续的森林变化现象。基于Landsat时间序列数据^[22-23]并充分利用不同干扰类型的时间和强度特性识别干扰或恢复的技术, 如LandTrendr算法^[20], 已被成功用于西北太平洋区域森林干扰的监测^{[19-20, 24]}。而且利用1986年起至当前, 近乎连续的时间记录南方人工林干扰或者恢复分布的研究也未见报道。所以, 基于监测所得综合性和连续性的人工林干扰和恢复数据记录对于南方区域森林管理和政策制定、生态系统保护和重建、生物多样性保护和碳评价有重大价值。

本研究主要目的基于不同来源的Landsat长时间序列遥感影像数据堆栈p122r043(条带号), 利用LandTrendr算法分析广东省佛冈县的森林干扰与恢复的时空模式, 包括森林干扰年份、森林干扰量、森林干扰持续时间等, 并借助现有影像以及历史数据与实地调查记录, 验证森林干扰与恢复的制图精度, 分析森林变化的驱动力及识别典型干扰类型时空特性及其关联要素。

1 研究区域与数据

佛冈县 (图 1) 位于广东省中部, 珠江三角洲边缘, 清远市东南角, 境内低山、丘陵、谷地、平原交错, 大多在海拔300m以下。全县山地丘陵区大部分是人工林, 主要树种为桉树、马尾松、湿地松、杉木、荷木、藜蒴、毛竹、荔枝和龙眼等果树, 目前全县森林以速生丰产的商品林和经济林占主导 (图 1)。本研究选取了1986—2011年间24期Landsat TM/ETM+遥感影像 (红色方形框为p122r043一景数据边界) 为研究数据, 该一景数据仅覆盖佛冈县大部分区域 (以下均简称佛冈县) 作为后续研究区 (图 1)。

本研究选取的24期Landsat TM/ETM+遥感影像来源于美国地质调查局的地球资源观测与科学 (USGS/EROS) 中心 (以下简称EROS), 其余购买于中国科学院对地观测与数字地球中心 (中科院地面站, 以下简称BJGS) 或者从地理空间数据云处免费获得。影像获取日期的选择必须满足植物生长旺盛的6—9月份, 在一些低纬度地区可以顺延至5—10月份, 以减少物候对于影像光谱识别的影响。在云量选择上, 由于获取完全无云影像有一定难度, 尽可能保证较少云量, 在一个时间序列内要保证一定比例的质量较好的数据。表 1描述了研究所用的Landsat TM/ETM+遥感数据的相关信息。广东省林业生态规划设计院提供了佛冈县2005年—2011年森林资源数据库和2005年—2006年的小班数据用于识别森林类型和不同地表类型 (图 1)。借助Google Earth提供的高清晰图像了解现地情形, 并根据现地调查和查阅当地历史资料及向有关人士咨询发现该区域存在大量森林采伐, 曾出现一些灾害现象。

2 研究方法 2.1 预处理与光谱指数设定 2.1.1 大气校正

将以上两种来源数据基于LEDAPS进行标准化处理构建LTSS, 该方法于沈文娟等^[25]中被具体描述。

2.1.2 云、阴影掩膜

为了进一步保证影像数据光谱识别的有效性, 基于需要修复的目标影像和无云污染的参照影像, 通过参照影像整合创建一个平均的参照影像 (参照影像的数量至少3—5幅, 确保影像清晰或者云掩膜效果好), 然后将目标影像与参照影像进行比较, 通过人为设定的阈值进行云和阴影的掩膜。云掩膜主要借助波段1进行阈值设定, 当像元值小于阈值在影像上显示为黄色, 即为需要掩膜的云覆盖区域；阴影掩膜主要利用波段4和波段5进行阈值设定, 黄色显示区域表明需要掩膜的区域 (图 2)。

2.1.3 光谱指数选取

在Kennedy等^[20]的研究中, 归一化燃烧指数 (NBR) 对于捕捉干扰事件具备最大敏感性, 比起归一化植被指数 (NDVI) 和缨帽变化湿度指数 (Wetness) 而言有较高干扰匹配的可能性。另在Cohen等^[19](2010) 研究中, 表明其对于干扰具备较好的解释能力。本研究中使用NBR^[26]作为算法监测的指数, 该指数通过近红外 (NIR) 和短波红外 (SWIR) 的差值与和值的比值用于监测干扰。

(1)

该比值中, NIR反映健康绿色植被, SWIR反映岩石和裸土。健康的森林有高的NIR值和低的SWIR值, 并具备高的NBR值。森林干扰后, SWIR值增加, NBR值表现为更低。然后, 时间序列堆栈中的影像数据被转化成NBR值, 与云掩膜数据进行像元与像元之间的匹配。并将这些时间序列源数据用于算法分割。

2.2 时间序列分割算法与滤波

LandTrendr算法的核心为时间序列分割。分割算法使用时间序列中的每个像元来识别变化的NBR值, 并保留时间序列中NBR值的最大变化, 对应的小的变化进行滤波操作, 目的是增加信噪比。具体参数设置如表 2所示。分割的第一步是由顶点的年份决定, 用于定义分割的结束点, 第二步是基于这些顶点, 使用点与点或者回归线的灵活混合得到最好的直线轨迹, 从分割算法得到的值有：每年的源数据 (每年对应像元最好的无云NBR值)、顶点对应的年份、每年的拟合NBR数据 (分割中每个点的NBR值被用于描述轨迹), 这些数据被用于后续的制图算法。整个过程识别出NBR值的3种特性, 即增加 (恢复)、下降 (干扰) 和保持不变。其中基于地表数据使用回归模型转换NBR值来评估植被覆盖的比例生成干扰量值^[20]。考虑到该回归模型是基于西北太平洋区域不同植被覆盖与密度构建, 本研究区植被组合和密度与其的不一致性, 在算法参数设置时进行了相应的调整以适应目前的研究区的植被特征 (表 2)。干扰持续时间值为1—25a (取决于影像堆栈的年份数), 由下降的NBR值对应的年份数决定。

2.3 森林干扰制图

基于年份干扰量和干扰持续时间影像创建干扰图 (图 3)。最小制图单位用于确保干扰的斑块至少有11个像元, 约为1 hm²的大小, 此阈值可以捕捉到大多数的森林采伐活动, 也避免了干扰事件制图面积太小而使得大范围内验证困难的问题 (表 2)。生成的干扰图有3个图层, 即干扰发生年份、相对干扰量、干扰持续时间。当出现超过一种干扰的情形, 根据斑块面积与相对干扰量的乘积之和进行干扰斑块的分类, 来平衡不同干扰量事件。获得最高值的斑块即为主要干扰, 第2高值的为次要干扰, 第3则为第3干扰, 但第3种情况较少而不被考虑。最后使用人工解译获取的森林训练区样本对影像进行随机森林分类^[27]。掩膜中识别的森林像元在整个时间序列中始终保持森林的状态, 最终的干扰图滤波去除了非森林像元中的任何干扰。本研究获取了佛冈县研究区域变化最大的干扰模式、持续4a之内变化最大的干扰模式、持续4a之内变化次要大的干扰模式、持续超过4a的干扰模式。由于该算法能监测急剧或持续干扰 (恢复) 事件, 将发生于4a内的干扰定义为急剧干扰或恢复, 超过4a则定义为发生了持续干扰或恢复。另外整个算法处理的是整幅影像, 最后用研究区域边界切割干扰结果影像。

2.4 森林恢复制图

为了获取干扰后的森林恢复, 本研究计算了植被恢复的两种指标 (绝对和相对指标)。绝对指标与干扰后的恢复有关, 表示干扰发生后5a之间NBR值的变化。计算公式如下：

(2)

式中, NBR_{fitted, t5}和NBR_{fitted, t0}分别表示干扰发生后拟合的5a时间内NBR值和顶点定义恢复分割开始的NBR值 (末端为干扰的分割)。ΔNBR_disturbance定义为时间序列开始之前干扰引起的森林损失。干扰后恢复值RI, 称为恢复指标, 计算公式为：

(3)

式中, 分母在干扰分割时随着NBR值变化, 并先于恢复发生。

2.5 森林干扰产品验证

基于LandTrendr算法得到佛冈县森林变化最大的干扰年份图, 由于变化最大的干扰年份图记录的是像元发生多次干扰, 只记录其中变化最大的一次干扰, 即记录在该区域发生的最主要的干扰, 也为典型森林变化区域。典型区域一般分布着主要干扰和次要干扰。主要干扰, 如皆伐、林分替换事件 (火), 干扰后的像元在影像上显示不再为森林像元；次要干扰其干扰后像元显示的亮度大, 比起干扰前像元绿色程度下降, 但是始终保持为森林像元。通常情况下主要干扰能够反映区域干扰事件的发生, 影响次要干扰的主要为一些低干扰量事件, 一般发生在主要干扰事件之后, 与主要干扰保持紧密的关系, 在一定程度上, 次要干扰事件可以忽略。再结合主要干扰和次要干扰图的分布, 通过光谱可视化方法对于主要干扰和次要干扰的发生进行目视识别。由于森林恢复一般需要经过较长的过程, 获取当地验证资料不容易, 本研究暂不考虑森林恢复的验证, 主要验证佛冈县变化最大的干扰年份图。

本研究选取佛冈县2个典型的区域, 使用受干扰年份对应的Landsat影像验证干扰制图 (图 4), 使用2km×2km的样方, 将其叠加到时间序列堆栈中的每个影像上, 并在每个影像对应的区域提取干扰区, 将目视解译得到的干扰区叠加到主要干扰年份图上, 对比真实解译的干扰面积与算法分割所得的面积, 确认正确识别的精度。参与验证的资料主要来自于对现有24期Landsat影像的目视解译及含小班区域的野外实地调查资料, 同时结合Google Earth进行辅助判别。

3 结果与分析 3.1 验证森林干扰产品

基于影像光谱目视识别及对照进行干扰制图精度的验证, 很容易地获取到干扰发生的典型区域。通过解译典型区域1每年的Landsat影像, 可知图 4典型区域1的主要干扰发生在2006年, 选取2002年和2005年的影像进行比对, 发现其他年份也有一定量的干扰, 但都不是主要干扰。对比目视解译的面积和干扰年份图中的面积, 目视解译的干扰面积为0.8768km², 主要干扰年份图中自动干扰制图面积为0.7659 km²。通过解译典型区域2每年的Landsat影像, 可知图 4典型区域1的主要干扰发生在2008年, 选取2002、2005、2006、2007年的影像进行比对, 发现各年份都有一定量的干扰发生, 2008年的干扰为主要干扰, 与算法所得的主要干扰的定义相符合。对比目视解译面积和干扰年份图中的面积, 目视解译的面积为0.3165 km², 主要干扰年份图中的自动干扰制图面积为0.2511 km²。经过以上的研究, 目视解译的结果与LandTrendr算法所得滤波干扰图像的结果有一定的误差, 约为0.1km²。此外, 通过影像光谱的识别, 也能辨别干扰分布的区域。

3.2 森林最大干扰和恢复变化的模式分析 3.2.1 森林最大干扰和恢复变化时间分布特征

图 5为森林变化最大的干扰年份图, 最大的干扰即像元发生多次干扰, 只显示从分割开始到结束最大的一次干扰, 图中水体、非森林、持续森林已被掩膜, 其余不同颜色表示不同的干扰年份分布。本研究基于像元尺度, 刻画了不同年份的森林干扰面积以及森林恢复面积 (图 6)。发现佛冈县的森林干扰较为剧烈, 一般都在1000 hm²左右。而1987、2002、2004、2005、2006、2007和2009年的干扰面积均超过2000 hm², 其中1987、2007年两年的干扰面积达到6000 hm²以上。相比森林干扰的变化, 佛冈县的森林恢复面积随时间的变化相对平稳, 以1987年为最大, 达到了19000 hm²。通过对佛冈县森林干扰和恢复面积的趋势分析, 发现20世纪80年代末到90年代森林干扰和恢复的面积基本少于2000年以后的变化面积, 变化趋势比2000年以后的平缓；从2000年开始, 森林干扰面积逐渐上升, 变化趋势高于森林的恢复, 但森林的恢复面积仍有所提升。

将干扰的像元进行邻近斑块组合, 并确保干扰斑块的最小面积为1hm², 用于计算不同干扰持续时间的面积百分比。图 7显示佛冈县的森林干扰面积随着持续时间总体呈下降趋势, 即干扰事件的发生时间保持较短的年份, 主要发生在1—5a范围内, 其中1—2a的干扰持续时间最多, 发生在1a的干扰面积所占比例最高, 约占26%, 表明佛冈县的森林干扰主要为急剧的干扰事件。

3.2.2 森林最大干扰变化空间分布特征

佛冈县的森林干扰存在较强的空间分布特征, 图 8为该县森林干扰面积随海拔高度变化的分布图。从图中可以看出随着海拔的升高, 研究区的干扰面积百分比呈现明显的下降趋势, 其中, 海拔分布于125—179m之间的干扰分布面积最大。总的来说, 超过50%的森林干扰面积发生在海拔小于179m的低海拔区域, 而当海拔大于787m之后, 森林干扰面积占研究区森林干扰总面积的比例仅为0.96%。

3.3 森林急剧干扰模式分析

佛冈县主要的急剧干扰事件 (主要干扰) 和次要的急剧干扰事件 (次要干扰) 均得到识别 (表 3)。研究表明20世纪80年代末到90年代末森林的主要干扰面积明显小于2000年以后, 说明急剧干扰在2000年以后发生较为频繁, 2005、2007、2009三年主要干扰面积在2000 hm²以上, 2004、2006、2008年主要干扰面积也在1000 hm²左右, 到2010年以后逐渐下降, 其中2005年主要干扰面积达到2700hm²。而且2000年以后主要干扰的发生有明显的波动趋势, 与佛冈县速生丰产林、人工林短轮伐期的特性相符, 表现出每隔两年基本出现的一次波动, 这种情况在2010年之前特别明显。图 9为急剧干扰持续时间分布, 主要为发生在1—4a内的干扰事件, 每年的干扰面积比例各不相同, 其中, 约38%的干扰持续1a的时间, 约28%的干扰持续2a时间, 约25%的干扰持续3a时间, 约7%的干扰持续4a时间, 说明短期的干扰事件发生较多。

3.4 森林持续干扰和恢复模式分析

由于1987年的相关数据异常, 移除该年的数据。另外持续干扰和恢复的定义而言, 2007年之后没有统计数据。研究发现持续的森林干扰和恢复在佛冈县有一定量的分布 (图 10), 但总的面积不大, 不超过100 hm², 2000年之前的森林持续恢复面积较持续森林干扰大, 到2000年之后持续干扰的面积呈现波动上涨的趋势, 2007年的时候达到持续干扰和恢复的最大值, 分别为380 hm²、220 hm²。此外, 从2005年开始, 森林的持续恢复面积增加, 表明在这段时间森林覆盖保持生长与更新的状态, 而在1989年也有森林增加的状态。另外, 通过比较佛冈县森林持续干扰和急剧干扰的面积, 发现2000年之前持续干扰和急剧干扰面积相当, 变化比较平缓, 到2000年之后, 急剧干扰的面积远大于持续干扰, 最高约达2800 hm², 但两者的变化趋势都呈现波动上升 (图 11)。

4 讨论与结论 4.1 驱动力分析

佛冈县属于人工林分布典型区域, 分布着大量的人工林、速生丰产林, 其森林干扰剧烈的原因主要由于一些急剧的干扰事件引起, 如人工林的采伐、炼山还有一些沿道路和河流的城市建设、恶劣气候影响等, 与已有森林资源调查和小班资料的记载较为符合。20世纪80—90年代广东省的林业政策使得大片造林的出现, 干扰相对较少^[28]。然而, 1995年佛冈县实施速生丰产桉树营造计划, 轮伐期为6a^[29]。而2000年以后随着人工林、速生丰产林达到成熟期, 大量的林木被砍伐。同时, 为了进一步的造林, 造林前炼山整地行为也催生大量的森林干扰事件^[30-31]。并且随着2004年经济开始增速, 城市化建设不断发展, 靠近道路和河流的森林不断受影响, 使得森林面积急剧减少, 干扰的面积呈现明显波动式增长^[32]。2005年到2007年分布着较为广泛的火烧迹地, 说明人为炼山在这些年份密集出现。2005年至2011年间, 用材林、速生丰产林、薪炭林采伐加剧, 形成了众多采伐迹地。其中2008年的冰冻雨雪灾害造成了严重的森林损毁, 随后的林业快速恢复政策有助于森林的恢复更新。同时, 人工林、速生丰产林比较适宜种植于中低海拔区域, 一方面易于种植管理, 另一方面易于运输管理, 海拔越高, 森林受干扰可能性越小。1987年和1989年干扰相关数据出现异常值, 由于图像的全云量覆盖, 数据质量较差, 不能正确地识别地表地物, 对于干扰和恢复的错分现象明显。

然而, 在驱动力分析中, 并没有考虑到影响人工林森林干扰和恢复所有可能的自然、人为驱动力因子, 而且没有建立它们之间定量的关系, 仅从定性角度进行了分析。

4.2 算法与技术分析 4.2.1 预处理与算法分析

LandTrendr算法在监测干扰年份、干扰量的基础之上, 能够监测特定年份内发生的干扰, 并能够识别主要与次要干扰的分布, 特别是能够监测急剧干扰和持续干扰的情况。LandTrendr获取到的是根据干扰量和面积确定的最大干扰发生的情况, 已经略去了一些细微干扰的发生。此外, LandTrendr算法在像元尺度基础上进行不同类型的滤波, 保证了干扰斑块面积的可识别性, 有助于监测一些典型干扰事件。该算法需要建立在Landsat时间序列数据堆栈基础之上, 影像数据配准、质量的好坏以及数据的格式都将会影响后续算法的精度。本研究中原始数据有不同的数据来源, 不同数据来源数据的处理与存档方式不一致, 在进行预处理之前必须进行标准化处理操作, 而非标准化数据的标准化过程必将存在误差；云掩膜、算法参数设置时人为阈值设定使得结果误差存在；算法要求的影像日期分布在7、8月份, 由于不可能获取到如此确切日期的数据并且要保证好的数据质量, 月份上的偏差产生物候的影响；影像1到2年的间隔使得不能很好捕捉缺失年份的信息, 影响干扰年际监测结果。本研究中含有若干数据质量不好的年份, 致使结果的误差存在, 如出现1987年异常高的森林恢复面积值。基于已有研究中所证明的最优指数基础上选取了单一指数用于算法监测, 不同指数对于算法运行好坏的贡献未知, 或者不同指数对于不同地域的监测效果可能存在差异性, 可以尝试多种指数的比较, 获取最优指数的监测。

4.2.2 验证方法分析

国外的验证技术较为成熟^{[17, 19, 33]}。如基于Timesync^{[19, 34]}工具的LandTrendr算法验证能够自动化获取解译结果, 并与算法分割结果做比对, 考虑该方法在国内使用缺乏历史地面数据和森林资源数据以及高清分辨率的影像, 未将其在本研究中推广。但基于本研究中的验证方法所得的结果, 发现干扰监测的结果与目视解译的结果相比呈偏小的趋势, 主要考虑到该算法对高密度针叶林覆盖区域的森林干扰监测更为敏感, 对于本研究中非完全针叶覆盖区域存在一定的限制, 并导致结果的错判和误判^[35]。综合现有分析结果, 尽管存在识别上的不足, 监测结果仍具备可参照性。在以后国内推广使用中需将本研究的方法形成具备相应标准的体系, 归档所需要的森林干扰历史资料, 收集丰富的地面数据辅助验证, 使得方法应用所得结果更为可靠。

4.2.3 算法与技术应用分析

森林干扰的4个特征, 即持续时间、强度、频率和面积^[36]。本研究未能实现森林干扰各项特征的量化与分析。由于连续年份的干扰类型与生物量历史统计资料的缺乏, 仅能监测出急剧和持续干扰发生的持续时间和面积, 对于具体的干扰类型、强度及其发生原因定位不足^{[1, 5-6]}。

4.3 结论

本研究尝试探索南方人工林区域森林干扰与恢复的影像自动化监测技术。尽管LandTrendr算法在国内的初次使用存有一些技术与数据选择困难, 但是通过整合与改进预处理技术能够成功的进行干扰制图, 并基于量化的森林变化进行驱动力定性分析, 不仅能够提供森林干扰与恢复的历史, 也组成了未来评价森林碳储量和气候变化下生物多样性保护的关键部分, 并为国内关于森林干扰与恢复长时间变化的自动化分析提供了参考。

致谢: 美国俄勒冈州立大学杨志强博士在算法使用中提供技术支持，美国马里兰大学黄成全博士帮助写作，特此致谢。