1. 摘要
人类活动改变了世界许多地区的干扰模式，但关于中国森林干扰格局和趋势的定量信息尚无记载。我们采用了光谱-时间分割方法 ，利用所有可用的Landsat数据绘制了单个干扰斑块，并表征了中国森林中干扰率、规模、频率和严重程度的模式和趋势。从1986年到2020年，大约39.7%的中国森林受到了干扰，年均干扰率为1.16 ± 0.41%。干扰率以每年−390 ± 142 km²的速度下降，主要由于自2000年代以来有效实施的森林保护政策。干扰率、频率和规模在东南地区普遍增强，但在东北地区减弱。我们的高质量、空间明确的干扰图为理解重要但不太为人所知的泛温带森林地区的景观尺度驱动因素和森林动态与功能提供了基本的数据层。
2. 引言
3. 结果
3.1 干扰机制
卫星观测显示，2000年中国的总森林面积为2.098亿公顷，这与国家森林清查数据、全球森林资源评估以及其他卫星观测估计结果一致。森林面积主要分布在东北（NE，20%）、西南（SW，30%）和东南（SE，37%）地区。我们的结果显示，1986-2020年间约有8330万公顷森林受到干扰，占中国森林面积的39.7%，年干扰率为1.16 ± 0.41%（平均值 ± 标准差，以下简称相同）（图1）。全国范围内，森林火灾仅占干扰的一小部分（约5%），区域热点（>30%）主要位于东北地区的大兴安岭落叶针叶林和西南地区的藏东南山区针叶林（补充图3a）。然而，气候变暖增加了东北地区大兴安岭的森林火灾比例（补充图3b）。
我们使用用户准确率和生产者准确率作为主要准确性评估指标。视觉导出参考数据与随机森林（RF）分类干扰图的比较显示，用户准确率和生产者准确率分别为0.76和0.96，总体准确率为0.96（CI95%=（0.955, 0.963），Kappa值=0.93）（补充表1）。时间准确性（即干扰年份）在3年内，与其他研究一致11。我们的空间明确干扰图还显示，与2001-2019年间基于Landsat的全球森林覆盖变化（r2=0.44，补充图3和补充图4）和2003-2009年间国家森林火灾记录（r2=0.5，补充图5）显著正相关。这些结果共同表明，本研究使用的光谱-时间分割方法大致捕捉了中国森林干扰的时空特征。
中国森林干扰模式存在显著的空间差异（图2，补充表2）。在0.5°网格单元中，年均干扰森林面积为10.17 ± 11.78 km²，呈高度左偏分布（补充图6），其中75%的干扰斑块年均干扰面积小于17.95 km²（补充表2）。少数地区（如西南山区针叶林）因较高的火灾活动水平，年均干扰森林面积相对较大（补充图3）。按每个网格单元的森林面积归一化的年均干扰森林比例大致显示出与年均干扰森林面积相似的模式（图2a和b）。年均干扰频率为每平方公里森林面积1.70 ± 1.67个斑块，东南地区的频率较高，那里自2000年代以来的森林管理实践（如人工林）提供了中国大部分木材。干扰斑块的平均大小为1.35 ± 0.45公顷，东南地区由于集约的人工林种植和东北地区频繁的大型森林火灾导致斑块较大。斑块大小分布也呈左偏，其中75%的网格单元大小为1.56公顷或更小。干扰严重程度，计算为干扰前一年和干扰后一年光谱植被绿度指数的变化（即ΔNDVI=NDVIt-1−NDVIt+1），为0.36 ± 0.51，长江以北的干扰严重程度普遍较高，可能由于火灾发生比例较高。
3.2 扰乱制度的趋势
1986年至2020年间，干扰模式发生了深刻变化，但其变化趋势在不同地区有所不同（图3，补充表3）。在全国范围内，年干扰面积以每年390 ± 142 km²的速度减少（p < 0.01）（图1b），相对于年干扰森林面积（24489 ± 8707 km²/年）减少了-1.6 ± 0.58%。全国干扰率的趋势主要由东北地区的减少趋势驱动（图1c和图3a），而西南和东南地区的整体干扰率没有显著变化（图1）。干扰模式的趋势还表现出南北梯度（图3）。总体而言，东南地区的干扰模式（率、频率和规模）加剧，而东北地区则减弱。干扰严重程度在中国的森林中总体上呈下降趋势，只有东南地区的少数区域显示出上升趋势。北方和南方中国的干扰模式不同趋势，广泛反映了区域性恢复策略，以促进生长较慢地区（如东北地区）的森林生态功能，并维持生长较快地区（如东南、西南地区）的木材产品。
3.3 森林恢复政策的作用
为了理解干扰模式变化的潜在驱动因素，我们使用广义线性模型（GLM）研究了社会经济因素、气候变化、植被类型、土壤特性以及Landsat时间序列变异性对每0.5°网格单元年干扰森林面积百分比的影响。GLM分析解释了年干扰面积变化的约42%，其中森林恢复比其他因素具有更大的影响。尽管Landsat图像变异性的非平稳性可能影响干扰趋势的检测，但我们发现其影响相对较小（图4）。蒸汽压亏缺、归一化植被指数（NDVI）和土壤有机质含量对年干扰森林面积有负面影响，这表明生长条件较好的森林经历了较低的干扰率。相反，海拔较高、人口密度较大、距离道路较近的地区更可能经历较高的干扰率（图4）。
4. 讨论
5. 方法
5.1 干扰检测
我们使用了一种成熟的光谱-时间分割方法，即基于Landsat的干扰和恢复趋势检测（LandTrendr），在Google Earth Engine（GEE）云计算平台上检测干扰。LandTrendr的核心是通过Landsat观测将像素级年度时间序列光谱轨迹分解为线性特征，从而提取一组与干扰相关的指标。LandTrendr已广泛应用于各种森林环境中的变化检测，关于算法的详细信息可见之前的出版物。这里我们简要描述了在GEE平台上使用LandTrendr生成中国森林干扰年份和类型的关键步骤。总体分析流程见补充图10。
首先，我们通过使用所有可用的、质量良好的（云量≤20%）第1类Landsat 5（主题制图仪）、Landsat 7（增强型主题制图仪）和Landsat 8（操作陆地成像仪）图像生成年度光谱一致的时间序列数据，这些图像是在1986年至2020年生长旺季（6月1日至9月30日）期间获取的。选择生长旺季以排除冰雪和土壤的复合影响，并最大限度地利用森林干扰后的光谱变化。为了应对Landsat传感器之间的光谱不一致性，我们根据带相应的系数通过线性变换调整光谱值。使用Fmask算法去除云、云影、雪和水。采用Medoid方法计算1986年至2020年间30米空间分辨率的年度带复合图。
其次，我们使用五个光谱指数运行LandTrendr，包括两个光谱带（短波红外I和II，即B5和B7）、塔塞尔帽湿度（TCW）、归一化燃烧比（NBR）和归一化植被指数（NDVI）。这五个指数是表征植被绿度和结构的有效指标，常用于检测森林干扰和恢复的变化62。对于每个光谱指数，LandTrendr生成一组参数来描述像素级的可能干扰事件，包括干扰前水平的光谱值（preval）、变化幅度（mag）、持续时间（dur）和变化率（rate）以及信噪比（dsnr）（n=5）。利用这五个光谱指数，我们生成了一组干扰相关参数层（n=25，5个光谱指数×5个参数），随后使用从参考数据衍生的机器学习模型检测和分类干扰。
5.2 干扰分类
5.2.1 参考数据
高质量的一致性参考数据是训练和分类干扰类型的关键。为此，我们采用分层方法生成了总共31225个参考点。我们首先使用2001年至2020年的森林损失数据生成了大量潜在干扰点。然后，通过将MODIS烧毁区（BA）与潜在干扰点叠加，按照Liu等人（2019）使用的程序将火灾干扰与非火灾干扰分开。具体而言，如果MODIS BA数据与Landsat衍生的火灾年份及火灾后2年内（即t+0, t+1, t+2）的森林损失数据一致，则确定为火灾干扰。接着，我们提取作为潜在干扰的点，包括火灾和非火灾干扰（包括森林转化为其他土地利用类型和不同强度的林业实践）。我们还生成了经历了无干扰（如持续森林）的相同数量的点，这些点是通过选择光谱指数变化很少的像素确定的。这些参考点，包括火灾、非火灾干扰和持续森林，随后用于采样1986年至2020年的时间序列光谱数据。最后，逐个检查每个参考点的时间序列光谱数据，确保它们准确代表干扰事件。此过程最终生成了总共31225个参考数据点，包括2356个火灾干扰点、13242个非火灾干扰点和15627个无干扰点（持续森林）（补充图2）。
5.2.2 随机森林分类
我们使用机器学习模型将每个像素分类为火灾干扰、非火灾干扰或无干扰。参考数据点用于采样上述LandTrendr衍生的干扰相关参数层，从而生成包含干扰类型的数据集。我们将数据集分为70%的训练数据和30%的验证数据。利用训练数据，训练了一个随机森林（RF）模型，将每个参考点分类为火灾、非火灾干扰或无干扰。我们的RF方法显示，基于短波红外（SWIR）的湿度指数（如B7，TCW）是检测森林干扰的强预测因子（补充图11），可能是因为它们对植被水分含量和冠层结构的敏感性。最后，我们将训练好的RF模型应用于整个分类堆栈，以一致地绘制1986年至2020年中国森林的干扰类型，假设参考数据时期2001-2020年的光谱轨迹可以外推到整个1986-2020年的映射期。然而，请注意，我们的方法旨在检测引起冠层开口的相对急剧和离散的干扰，而不是低强度林业实践和慢性干扰导致的森林结构或组成的细微变化。
5.2.3 干扰年份
我们使用LandTrendr通过光谱变化幅度的开始时间确定干扰年份。由于我们在五个光谱指数上运行LandTrendr，每个像素有五个可能的干扰年份。因此，我们使用至少三个不同指数的中值来确定干扰年份（即NDVI、NBR、TCW、B5、B7）。通过这种方式，我们仅保留了使用至少三个指数检测到的干扰像素，从而减少了假阳性错误。LandTrendr生成的光谱变化最大的年份通常在3年内准确11。根据显示可能干扰事件的指数数量，还为每个干扰像素分配了置信水平。具体来说，如果干扰由三个、四个或五个光谱指数检测到，则分别分配低、中、高置信度。
5.2.4 精度验证
我们在像素和国家层面对干扰图进行了验证。在像素级别，我们使用前一部分描述的验证子样本对最终图进行验证。我们通过混淆矩阵报告用户准确率和生产者准确率（补充表1）作为主要的准确性评估指标。在国家层面，我们将本研究检测到的森林干扰与现有数据集进行了比较。具体来说，我们比较了本研究与2003年至2009年国家火灾记录之间的森林火灾干扰面积（补充图5）。我们还比较了本研究与2001年至2019年基于Landsat的全球森林覆盖变化之间的干扰率（补充图4）。
5.2.5 数据后处理
我们应用了一系列空间滤波器，以最小化来自时间序列光谱轨迹或图像之间误配引起的两个潜在不确定性来源的不现实离群值。这可能导致（a）个别像素或包含少量像素的小斑块被检测为干扰，从而增加了假阳性错误，或（b）未检测到干扰，而其周围像素大部分受到干扰，从而增加了假阴性错误。为了解决问题（a），我们通过设置最小映射单位为两个30×30 m²像素（0.18 ha）来删除所有单像素干扰斑块。为了解决问题（b），我们应用了3×3移动窗口，通过根据周围像素的年份分配干扰年份来填补空白。最后，我们通过在3×3窗口内使用多数规则分配中心像素的干扰年份，从而解决了与正确识别干扰年份的不确定性相关的伪影。
5.3 表征干扰模式及其趋势
基于前一步生成的年度森林干扰图，我们在中国森林的每个0.5°网格单元（n=1946）内使用五个指标表征干扰模式。在每个网格单元内，我们计算了：（1）年均干扰森林面积（km²/年），（2）年均干扰森林比例（%/年），即年干扰森林面积除以总森林面积，（3）干扰斑块大小（公顷），即使用八邻规则计算每个斑块的干扰像素数，（4）干扰频率（每年每1000 km²森林面积的斑块数），即每年干扰斑块数除以总森林面积，（5）干扰严重程度（ΔNDVI = NDVIt−1 − NDVIt+1），即LandTrendr分析中干扰前一年和干扰后一年NDVI变化的幅度。我们使用（1）和（2）来表征干扰率，使用（3）-（5）描述斑块特征。将（2）和（4）标准化为每个网格单元的森林面积，使其在网格单元之间具有可比性。对于（3）-（5），我们仅计算大于0.45 ha（五个30×30 m²像素）的斑块，因为小于0.45 ha的斑块仅贡献了不到10%的干扰面积，同时大大增加了计算需求。较大的斑块也更有可能被更准确地检测到，对森林动态和功能有更大的影响。在我们的分析中，面积通过像素数乘以像素大小（0.09 ha）计算。
我们报告了每个网格单元的五个干扰模式指标，并在每个网格单元内汇总斑块基于斑块的指标（干扰大小、频率和严重程度）。我们使用算术平均值可视化了森林干扰模式的空间模式，并使用算术平均值和中位数报告了统计数据。使用Mann Kendall趋势检验量化了每个网格单元中五个干扰模式指标的趋势，这是一种对时间序列中的单调趋势不敏感的非参数检验。
我们还评估了卫星图像可用性变化对干扰趋势的影响。随着Landsat任务的发展，每年可用的Landsat图像数量在增加（补充图1）。尽管这为我们提供了更多的机会来检测地面变化，但这可能会导致假阳性趋势。1986年至1999年间可用的Landsat图像数量较少，但随时间增加。因此，年干扰面积存在很大的年际变化，且呈不显著的下降趋势。然而，自2000年以来，可用的Landsat图像数量翻了一番，除了2000年和2012年由于TM和ETM + 传感器退役而导致的可用Landsat图像数量特别低（补充图1）。这导致了年干扰面积的相对较小的年际变化和下降趋势。尽管更详细的分析（见下文）表明Landsat观测的增加趋势有助于1986-2020年间干扰率的增加，但我们的结果显示，与其他因素相比，其影响相对较小（图4）。此外，如果1986-1999年间Landsat观测的密度保持一致，干扰率的下降幅度将更高，因为后期观察期的更密集的Landsat观测可能会产生假阳性趋势。
5.4 森林干扰模式变化的潜在驱动因素
为了理解森林干扰模式变化的潜在驱动因素，我们使用广义线性模型（GLM）探讨了气候因素、生物物理因素（地形、植被、土壤）、社会经济因素（如管理、人口密度、道路密度和距离最近的道路的距离）以及Landsat观测的变化对每个0.5度网格单元年均干扰森林面积的影响。气候因素用于研究背景气候是否对干扰模式有影响。我们总结了1985年至2020年期间使用4公里月度网格化TerraClimate数据的年均气温（MAT）、年均降水量（MAP）和蒸汽压亏缺（VPD）。由于VPD与MAT和MAP存在共线性，故只将其纳入分析。森林管理变量包括森林管理面积的百分比和生态项目的状态。森林管理面积的百分比计算为具有森林管理活动的像素数除以0.5°网格单元内的所有森林像素数，使用2015年的新衍生森林管理数据集。生态项目的边界用于确定干扰斑块是否位于生态项目区域内（即1-生态项目区域内，0-生态项目区域外，生态项目的详细描述见下文）。人口密度（人/km²）、道路密度（km/km²）和距离最近道路的距离（km²）在0.5°网格单元内总结。Landsat观测的趋势通过所有Landsat观测的年可用性计算（补充图1）。生物物理因素包括地形（如海拔、坡度）、植被（如NDVI作为生长率的代理）和土壤（如土壤有机质含量、土壤容重）。
GLM使用年均干扰森林面积作为因变量，气候因素、生物物理因素、管理和Landsat观测的可用性作为自变量进行拟合。GLM使用R的stats包中的glm函数进行拟合。基于Akaike信息准则（AIC）使用逐步算法选择最简模型。GLM的标准化回归系数可以解释为自变量的相对影响。部分图用于指示自变量的方向（图4）。
5.5 森林恢复政策及其对干扰的影响
自1949年中华人民共和国成立以来，森林管理政策经历了重大变化。1970年代之前，过度的木材采伐减少了森林覆盖，导致了森林条件的退化，并引发了严重的生态和环境后果，如水土流失、沙尘暴和洪水。为了恢复退化的森林生态系统，自1987年以来，中国政府实施了几项国家重点生态恢复项目，包括三北防护林工程（TNSFP，也称为绿色长城）、长江和珠江防护林工程（RSFP）、天然林保护工程（NFPP）、退耕还林工程（GGP，也称为中国坡地转换工程）。TNSFP于1970年代末启动，是一项造林计划，覆盖了中国北部的一半地区。其目的是改善环境条件，减少水土流失，并促进三北地区的森林/草地相关产品。RSFP于1989年启动，是一项森林恢复计划，覆盖了长江和珠江流域。其目的是减少南方中国的水土流失、岩石荒漠化和洪水发生。NFPP于1999年启动，是一项森林保护/保护计划，旨在保护中国的天然林，目标是保护生物多样性、保护水质、防止水土流失和荒漠化，并减少与砍伐森林相关的洪水和其他自然灾害的可能性。在NFPP地区，主要的恢复策略包括减少和调整采伐强度，并加速在生态敏感地区，特别是长江和黄河流域的造林。GGP于2000年启动，是一项全国性的生态补偿计划，鼓励农民将坡地（>20度）的耕地转为森林。除上述森林相关计划外，还有其他几个生态项目，包括京津风沙源治理工程和退牧还草工程。然而，这两个项目侧重于干旱生态系统，因此未纳入我们的分析。
为了进一步评估森林管理政策对干扰模式的影响，我们比较了生态恢复项目边界内外的年均干扰森林面积在0.05显著性水平（即p < 0.05）下正趋势和负趋势的比例（补充图8）。生态恢复项目边界内负趋势比例较高表明项目的有效实施导致了干扰率的下降。这些项目的边界是通过从行政机构获得的地图在省级进行数字化的。
5.6 绘制森林扩展图
除了绘制森林内部的干扰图外，我们还通过森林恢复/造林绘制了森林外的扩展模式和速率，以全面捕捉调节中国森林变化的干扰模式。由于LT-GEE无法绘制森林范围动态，我们通过合并1990-2019年的30米年度土地覆盖图30和1985-2016年的AVHRR、2000-2020年的MODIS67和2000-2020年的Landsat3的平均树冠覆盖图，衍生了森林范围（定义为树冠覆盖>20%）。年度森林扩展定义为先前非森林土地利用类型上的新获得森林，并在研究期间保持森林状态。森林扩展率（km²/年）在中国森林的每个0.5°网格单元（n = 1946）内绘制（补充图7）。
5.7 森林掩膜
1986年的森林掩膜通过以下程序衍生。首先，通过合并高级甚高分辨率辐射计（AVHRR）的1986年树冠覆盖图（TC1986）和Landsat传感器3的2000年树冠覆盖图（TC2000）获得1986年复合树冠覆盖图（TCcompostite1986）。TCcompostite1986的树冠覆盖百分比确定为TC1986或TC2000中的最大树冠覆盖（即TCcompostite1986 = max(TC1986, TC2000)）。通过这种方式，我们可以捕捉到2000年之前可能发生的森林干扰（如TC1986 > TC2000），以及1986年至2000年期间森林面积的扩展（如TC1986 < TC2000）。然后，按照Liu等人（2019）的方法，将TCcompostite1986 > 20%定义为1986年的森林掩膜。我们应该注意到，我们的研究区域不包括2000年后的新造林区。所有分析均在森林掩膜内进行，从而排除了其他土地覆盖类型的潜在混杂因素。有关TC1986和TC2000的描述可见于Hansen（2013）和Song（2018）。
6. 数据可用性
7. 代码可用性
8. 补充资料
引用格式：Liu, Z., Wang, W.J., Ballantyne, A. et al. Forest disturbance decreased in China from 1986 to 2020 despite regional variations. Commun Earth Environ 4, 15 (2023). https://doi.org/10.1038/s43247-023-00676-x
