通过捕捉生物量密度的时间变化来估算森林退化程度：来自热带干旱森林的见解

《Remote Sensing Applications: Society and Environment》：Estimating forest degradation by capturing temporal change in biomass density: insights from a tropical dry forest

本研究通过结合2019年和2022年的多时相地面实测数据与Sentinel-1/2遥感影像，采用随机森林回归模型和纹理特征分析，量化了墨西哥Ayuquila河流域热带干旱森林的AGB变化。结果表明，尽管平均AGB未显著变化，但冠层覆盖和高度下降，茎密度增加，表明存在渐进式退化。经严格验证，退化面积占4.3%，恢复1.3%，94.4%区域稳定。验证了多源数据融合在森林退化评估中的有效性。

　　### 森林退化与恢复的估算：结合多时相地面数据与遥感技术

森林退化通常表现为森林冠层覆盖度和生物量密度随时间的减少，这种变化会削弱森林生态系统功能和所提供的服务。在许多热带干林地区，森林退化常常出现在森林的局部区域，因此仅依靠地面数据进行检测变得困难，因为地面调查通常是在广泛分布的采样点进行的。本研究利用多时相的地面生物量数据和遥感技术，对一个长期受到轮耕影响的热带干林进行了森林退化的估算。通过随机森林回归模型和Sentinel-1和Sentinel-2影像，我们对2019年和2022年的生物量进行了建模，并计算了这两个年份之间的生物量变化（ΔAGB = AGB₂₀₂₂ - AGB₂₀₁₉），并将ΔAGB图重新分类为三个类别：生物量增加、无变化和生物量减少。通过随机分配1257个采样点到这三个类别，我们进行了严格的精度评估，以获得无偏的面积估算。

尽管2019年和2022年的地面数据平均生物量密度在统计上没有显著差异，但观察到冠层覆盖度和冠层高度略有下降，而树干密度有所上升，这与轮耕模式中常见的萌发现象相符。2019年的生物量模型RMSE为14.5 ± 4.8 Mg ha^-1，而2022年的模型RMSE为14.0 ± 4.52 Mg ha^-1，表明两个年份的生物量估算具有可比性。我们的生物量损失估算精度高于生物量增加，这说明使用两个时相的地面数据能够实现比单纯依赖遥感更高的估算精度。生物量损失占总森林面积变化的4.3%，而生物量增加则仅占1.3%，大部分森林区域（94.4%）保持稳定，因此整个研究区域的生物量呈现出轻微的负趋势。

### 研究区域与数据方法

本研究在墨西哥哈利斯科州的Ayuquila河流域（ARB）进行，该地区因在生态系统服务方面的重要性而被纳入REDD+的早期行动计划。ARB的面积约为4000平方公里，位于19°25'至20°10'N和104°3'至103°3'W之间。该地区具有亚热带气候，年平均温度为18-22°C，年平均降水量为800-1200毫米，主要集中在6月至10月。热带干林（TDF）是该地区的主要森林类型，覆盖了ARB的约24%；其他主要森林类型包括温带森林，如松林和橡树林，合计约占12%。ARB地形复杂，海拔范围从260米到2500米；平地主要用于永久农业，而TDF则分布在坡地上，受轮耕影响较大。

在2019年和2022年，我们分别采集了41个和22个圆形500平方米的样地数据。2019年的样地设计遵循随机抽样，并按森林类型（保护森林和退化森林）进行分层。2022年的样地设计也采用随机抽样，但考虑到可及性因素，如距离公路的距离（>500米且<5000米）。样地的生物量密度因人为退化和不同的生物物理条件而有所变化。除了胸径（DBH）外，我们还收集了冠层高度、冠层覆盖度、树木数量和树干数（>2.5厘米）等数据。由于在轮耕制度下，大多数树木会迅速萌发，因此在2022年的调查中，我们记录了树种的分类信息，并利用2022年所有树木的平均木材密度（WSG）来计算2019年树木的WSG。

### 遥感数据与预处理

我们使用Google Earth Engine（GEE）平台上的JavaScript代码下载Sentinel-1和Sentinel-2影像。这些影像的获取时间与样地数据采集时间相隔约三个月。我们选择了Sentinel-2的处理级别2A，该影像经过正射校正和大气校正，具有接近零云量的观测条件。我们选择了具有10米空间分辨率的光谱波段，包括蓝、绿、红和近红外波段，并计算了归一化植被指数（NDVI），作为植被生长的代理指标。

Sentinel-1影像通过Sentinel-1工具箱进行校准和正射校正，包括热噪声去除、辐射校正和地形校正。地形校正后的数据通过对数缩放转换为分贝（dB）（10*log10(x)）。Sentinel-1影像注册在干涉式宽波束（Interferometric Wide Swath）模式下，包含垂直-垂直（VV）和垂直-水平（VH）的发射接收极化波段。影像来自下降轨道，分辨率为10米。我们应用了一个30米半径的圆形平滑滤波器以去除斑点噪声。

我们还从NASADEM_HGT v001下载了数字高程模型（DEM）数据，并使用QGIS 3.28中的栅格分析估计了坡度。坡度以弧度测量，是景观在最陡方向上的高程变化率。这些数据用于结合其他遥感变量进行生物量预测。

### 纹理特征的提取与应用

纹理特征可以衡量图像中灰度级强度的空间排列。这些特征为生物量估算提供了额外的信息，并被发现可以提高预测准确性。灰度级共生矩阵（GLCM）包含图像中固定相对位置的像素数字对的相对出现频率，是遥感中常用的一种纹理特征提取方法。GLCM纹理特征由纹理统计、窗口大小、量化级别和像素间距离与方向决定。我们使用R中的GLCM包提取了Sentinel-2A影像的绿色波段和Sentinel-1影像的VV和VH波段的纹理统计，在四个窗口大小（3x3、5x5、7x7和9x9像素）和四个方向上进行提取。我们进行了初步实验，并确定了从绿色波段和VV、VH波段计算出的平均值、方差和熵作为生物量预测中最有用的纹理特征。我们考虑9x9像素的窗口大小为合理的对应尺度，使光谱信号与地面数据之间有合理的匹配。

### 预测变量的准备与地面数据对比

预处理后的数据使用Sentinel-2栅格作为模板，通过默认的双线性算法进行重采样，并重新投影到UTM WGS84 14N坐标系，空间分辨率为10米。在每个采样年份，我们根据纹理数据的窗口大小构建了四个数据集。每个数据集中的所有栅格层通过R中的“raster”包命令stack和brick堆叠成多层栅格砖。我们还将2019年和2022年的地面测量点重新投影到相同的投影。使用R中的extract命令，从栅格堆叠中提取出落在样地范围内的均值。我们使用“rstatix”库中的wilcox_test命令进行了Wilcoxon符号秩检验。

### 随机森林回归模型用于生物量预测

我们应用了随机森林（RF）回归模型来估算2019年和2022年的生物量。在RF中，树的构建方式类似于构建一个自助样本（bootstrap）。对于RF参数，我们采用了默认的ntree值500，以及解释变量数量的平方根作为mtry。节点大小设置为1，随机种子设置为4。对于每个采样年份，我们使用不同的预测变量（即不同窗口大小的纹理指标）训练和评估随机森林模型。由于我们拥有相对较小的数据集（2019年41个样地，2022年22个样地），我们使用了一个包含100次随机化的自助数据集来评估不同模型的性能。在每次随机化中，不同的数据被包含在训练集和测试集中。我们使用RMSE（均方根误差）和相对RMSE（RMSEr）来选择最佳的生物量预测模型。RMSEr是通过将RMSE除以测量数据的平均值计算得到的。对于每组模型，我们使用每个随机化测试集的均值和标准差计算RMSE。最终，我们选择了在测试集中获得最低平均RMSE的模型作为最佳模型。一旦确定了预测变量，我们使用随机选择的90%数据训练最终模型，并在剩余的10%数据上进行评估，以预测整个研究区域2019年和2022年的生物量。

### 生物量模型验证与误差图计算

一旦选择了最佳模型，我们通过使用自助法（bootstrapping）生成10个随机训练和测试集，来估算2019年和2022年生物量预测的不确定性。在每个年份，我们预测了10个生物量图，并计算了这些图的均值和标准差。基于这些图，我们通过将标准差图除以均值图，计算了相对误差的变异系数（CV）。CV图可以被解释为生物量预测不确定性的指标，突出了均值预测中的差异。

### 森林退化与恢复的估算

我们对2019年和2022年的森林生物量进行了地图绘制，并计算了ΔAGB（AGB₂₀₂₂ - AGB₂₀₁₉）。我们将ΔAGB图在整个研究区域重新分类为三个类别：生物量增加区域、无变化区域和生物量减少区域。由于我们发现ΔAGB的标准差为7.6，我们定义了-7.6 Mg ha^-1至7.6 Mg ha^-1之间的区域为无变化区域。我们采用分层随机抽样程序评估了ΔAGB图的准确性，并使用方程（2）计算了样本数量。我们分配了300个点用于生物量损失区域，300个点用于生物量增加区域，其余661个点用于无变化区域。我们假设森林覆盖与生物量具有等价关系，因此可以利用森林覆盖变化来验证ΔAGB图。为此，我们在每个点上放置了10米网格，并通过Google Earth的视觉解释对网格内的土地覆盖进行了分类，参考历史高分辨率影像作为基准。在需要时，我们还参考了Sentinel-2影像。在参考数据标记过程中，我们删除了4个最初分配给生物量损失的点，由于缺乏解释信息，最终保留了1257个点。我们使用FAO开发的Open Foris精度评估工具（通过R实现）进行了无偏精度评估和面积估计，并使用方程（3）和（4）计算了95%置信区间（CI）。

### 研究结果

森林结构变量的统计结果见表1。基面积、生物量和树木数量在2022年的数值高于2019年，但这些增加在Wilcoxon符号秩检验中并未显著。另一方面，冠层覆盖度和冠层高度在2022年显著降低，而树干密度显著增加。这些发现表明，由于轮耕和树干提取，森林正在经历持续的退化。在AGB变化图中，ΔAGB损失的范围为-32.91至-7.6 Mg ha^-1，而ΔAGB增加的范围为7.6至22.8 Mg ha^-1。ΔAGB值在-7.6至7.6之间的区域被分类为无变化区域。AGB损失和增加的箱线图见图6。

### 森林退化与碳损失的面积估算

通过使用地面数据进行校正，我们估计了AGB损失的面积为4696.99 ± 1165.08公顷，占总面积的4.3%。值得注意的是，校正后该比例从19.9%降低到4.3%，这表明仅依赖遥感数据绘制的ΔAGB图大大高估了AGB损失。同样，AGB增加的面积校正后从14.9%增加到61.3%。根据2019年的平均AGB密度26.83 Mg ha^-1，我们可以计算出AGB损失总量为126,020.24 ± 31,259.10 Mg。然而，由于我们缺乏再生植被的AGB密度数据，AGB增加的估算较为困难。再生植被的AGB密度可能远低于2022年热带干林的平均AGB密度。

### 讨论

在本研究中，我们使用了两个时间点（2019和2022）的地面测量数据和同一时期的卫星影像，以估算热带干林的退化情况。通过全区域的AGB地图绘制，我们获得了AGB增加、无变化和AGB损失的区域。通过使用1257个采样点进行严格的精度评估，我们发现AGB损失的估算精度高于AGB增加。校正后的AGB损失区域（4,697 ± 1,165公顷）占总面积的4.3%，而AGB增加区域（1,435 ± 622公顷）仅占1.3%，大部分森林区域（94.4%）保持稳定，表明尽管存在一定的恢复，整个研究区域的森林覆盖呈现出轻微的负趋势。

森林退化通常发生在相对稳定的森林区域中，尤其是在热带干林这种通常具有较低森林覆盖率的地区。如果没有在整个研究区域进行评估，很难在退化变得严重之前检测到这些退化区域。此外，如果没有全区域的绘制，高度不可达地区的退化可能会被遗漏。因此，像REDD+这样的项目特别关注开发全区域的森林退化和砍伐估算方法。我们的案例表明，多时相的AGB地图绘制能够有效估算森林退化和恢复。同时，我们强调了严格的精度评估和面积校正的重要性，因为单独依赖遥感数据往往高估了相对罕见的森林退化和再生现象。这一方法有潜力减少碳核算中的不确定性，并提高REDD+实施的效率。

### 结论

森林退化降低了森林对生物多样性的保护能力，并影响了生态系统功能和服务。本研究利用多时相的地面测量数据和Sentinel-1和Sentinel-2影像，通过全区域的AGB地图绘制，估算了一个长期受到轮耕影响的热带干林的退化情况。我们获得了AGB增加、无变化和AGB损失的区域，并通过严格的精度评估方法，发现AGB损失的估算精度高于AGB增加。校正后的AGB损失区域（4,697 ± 1,165公顷）占总面积的4.3%，而AGB增加区域（1,435 ± 622公顷）仅占1.3%，大部分森林区域（94.4%）保持稳定，表明尽管存在一定的恢复，整个研究区域的森林覆盖呈现出轻微的负趋势。本研究为利用多时相的生物量地图来检测森林退化提供了一种改进的方法，并表明严格的精度评估和面积校正对于准确估算森林碳损失和恢复至关重要。这有助于减少碳核算中的不确定性，并提高REDD+项目的实施效率。