Distribution of forest carbon stocks in “Three Parallel Rivers” region based on GTWR-GBRT model

“三江并流”保护区位于滇西北高原横断山脉纵谷地区(25°30′~29°00′N，98°00′~100°30′E)，是指怒江、澜沧江、金沙江3条大江在滇西地区平行并流的最近区域，包括云南省怒江市、迪庆州、丽江市的3个市州及8个市县(图1)；拥有丰富多样的地貌类型，是多种地质演化的代表区域；拥有众多孑遗物种和珍稀濒危物种，是中国重要的水能资源、矿产资源、旅游资源和生物资源的集聚地^[13]。因多样的自然地貌和气候条件(涵盖南亚热带至寒带多种气候类型)而形成丰富植被类型，主要包括常绿阔叶林、硬叶常绿阔叶林、针阔混交林以及暖温性、温性针叶林等，体现了显著的生态多样性^[14]。

Figure 1.  Schematic map of the study area

实测数据为2012、2017、2021年的云南省森林一类清查数据。本研究样地为方形样地，边长为25.82 m。数据包括样地号、优势树种、树龄、经纬度坐标、郁闭度等信息。林木蓄积量通过样木胸径和树高的实测数据，利用二元立木材积公式计算获得。

在美国地质调查局(USGS，https://earthexplorer.usgs.gov)网站获取研究区Landsat卫星遥感影像，分辨率为30 m，其中2012年数据来自Landsat 7，2017与2021年数据来自Landsat 8。考虑到云覆盖产生的影像质量问题，选取每年6—12月云量低于5%所有可用的卫星影像数据，每年取8景(序列号均为130/42、131/41、131/42、132/40、132/41、132/42、133/40、133/41)，共为24景影像，且经查验无坏带，进行大气校正、辐射定标^[15]。在地理空间数据云平台(https://www.gscloud.cn/)下载研究区30 m的DEM高程数据。

基于蓄积量-生物量转换模型^[16−17]，估算样点的生物量并计算碳储量，公式如下：

式(1)和(2)中：B为单位面积生物量(t)，V为单位面积蓄积量(m³)， a、b为参数；C为碳储量(t·hm⁻²)；F为含碳系数。各优势树种的参数a、b及含碳系数F均参考文献[18]，研究区森林碳储量见表1。

利用光谱、地形及纹理特征等信息对森林碳储量进行建模，可以提高估测精度^[19]。基于遥感影像提取研究区的地形和遥感因子，分别提取6个单波段因子及对应的48个纹理特征(3 × 3窗口)、4个植被指数因子，以及坡度、坡向、海拔共61个特征因子(表2)，这些因子涵盖了遥感影像中的主要信息类型。

由于变量间的线性关系会掩盖其统计显著性，增大参数估计误差，导致模型预测性能降低，因此需要检验变量间是否存在共线性问题，而方差膨胀因子(VIF)是诊断多重共线性的关键指标^[20]。为探究研究区地形因子和遥感因子与碳储量之间的关系，利用SPSS软件绘制碳储量与各变量的散点图，剔除异常值后与变量进行逐步回归分析和共线性检验，一般认为VIF＜10代表不存在多重共线性^[21]，本研究中VIF＜2，最终确定6个因子作为GTWR模型的输入变量(表3)。

在第2阶段对碳储量以变量相关性较大的因子进行RF的特征参数重要性评价，GTWR依赖共线性筛选的变量保证参数可解释性，GBRT引入重要性变量增强预测能力。而机器学习中，变量过多会对模型性能、训练效率及泛化能力产生多方面的负面影响。因此本研究选取重要性贡献度大于80%的变量，即第6波段差异性纹理(b6diss)、第4波段方差纹理(b4var)、第6波段(b6)、第6波段均值纹理(b4mean)，以及第1阶段逐步回归变量和GTWR模型估算的碳储量(GTWR_tcl)作为机器学习的输入变量。

空间自相关性检验用于评估研究区内变量的空间关联性，其系数能够量化变量在空间分布上的相互依赖程度^[22]。通常采用全局莫兰指数(I)来判别^[23]，公式如下：

式中：w_ij为空间权重，x_i、x_j分别为第i个样点和第j个样点的数据观察值，$ \bar{x} $为数据平均值。一般，I取值为[−1,1]，越接近1表示空间自相关越强，越接近−1表示差异性越强，0表示在空间随机分布。

由表4可得：2012—2021年莫兰指数均大于0，P均小于0.01，表明森林碳储量在空间上存在显著聚集。

以云南森林调查的样地数据为主要数据源，结合地形、遥感因子，并与碳储量建立单一GTWR、RF、GBRT及组合模型GTWR-RF、GTWR-GBRT。GTWR模型在考虑影响因子时空异质性的情况下，具有反演碳储量分布的优势。GBRT模型能够整合驱动因子间复杂的非线性关系，而不考虑变量之间的多重共线性问题，具有较高的预测性能。GTWR模型难以描述复杂的非线性关系，GBRT模型缺乏对异质性的考虑，组合模型可以有效避免单一模型的缺陷。用80%的样地数据进行模型构建，20%的样地数据进行精度验证。组合模型主要分为2个阶段，第1阶段通过多元线性回归选出适合GTWR模型的变量，得到具有时空特征的GTWR_tcl；第2阶段将第1阶段剔除的变量进行特征重要性分析，选择排前4的变量及第1阶段的变量纳入非线性模型中，构建研究区的最佳拟合模型并进行反演，以提高对森林碳储量的估测精度。

为评价模型的可靠性与估计性能，采用10折交叉验证法，通过比较模型预测值与实际测量值来评估其预测精度。评价指标包括决定系数(R²)、相对均方根误差(r_RMSE)及预测精度(A)，R²越高表明模型拟合越好，而r_RMSE越低，A越高，则表明预测精度越高。公式如下：

式中：$ y_{i} $为实测值，$ {\hat{y}}_{i} $为模型预测值，i为第i个样点，n为总样点数。

由表5可知：GTWR模型的性能(R²=0.73)相对有限，这是由于GTWR为统计回归模型，虽能体现碳储量的时空相关性，但难以充分捕捉数据的复杂非线性特征。RF模型性能(R²=0.92)优于GTWR，能更稳定地捕捉碳储量与其影响因子间的复杂非线性关系，从而获得更稳健的预测结果。GBRT模型性能(R²=0.95)高于RF，关键在于它能通过加法模型与前向分步算法，进一步优化拟合。GTWR-RF模型(R²=0.95)的拟合效果优于单一模型，而GTWR-GBRT模型表现则更为突出，其R²达0.98，预测精度较GTWR-RF提升0.06，且r_RMSE最低。这表明GTWR-GBRT模型在拟合与预测中均具有优势，尤其适用于兼具时空依赖性和非线性特征的复杂数据集。

根据GTWR-GBRT模型得到的估测结果，绘制2012—2021年“三江并流”区域的碳储量分布图(图2)。2012—2021年，研究区森林碳储量的高值区范围及整体空间格局呈现动态演变特征。从空间分布来看，碳储量整体呈现“西高东低，局部偏高”的格局，高值区主要分布在迪庆州，这与其地处横断山区、地形复杂、人类活动干扰少、原始森林保存完好密切相关；北部高碳储量区的连片性虽随年份有变化，但核心范围始终稳定，生态系统的基础稳定性形成了高效碳汇结构。与之形成对比的是，丽江碳储量总体偏低，其可能原因在于该地区旅游业发达，人类活动干扰对森林生态系统产生了潜在的影响；南部区域碳储量整体偏低且空间异质性较小，同样受到地形限制和人类活动(如开发、产业活动等)的双重影响。从时间演变看，空间格局变化明显。2012年高碳储量区范围相对广，说明当时区域生态系统碳汇能力较强；2017年，高碳储量区面积有所缩减，低碳储量区范围扩张，说明生态系统碳储存受干扰，可能是干旱、森林火灾或人类活动加剧，破坏植被、降低碳固定；2021年，高碳储量区有一定的恢复，主要得益于生态修复工程和天然林保护。

Figure 2.  Spatial distribution map of predicted carbon stocks based on the GTWR-GBRT model

由表6可知：研究区碳储量呈现明显的海拔梯度分布特征，不同海拔梯度的碳储量占比在2012—2021年保持相对稳定，没有出现剧烈波动，表明其空间分布格局具有较高的稳定性。海拔梯度的变化会通过温度、光照和水分等环境因子的协同作用，直接或间接地影响森林碳储量。总体碳储量主要分布在海拔2 500 m以上，海拔高于3 500 m的碳储量占比最大，为33％。一般而言，高海拔地区因植被减少和生长缓慢，碳储量通常较低，而“三江并流”区域位于中国南部热带地区，雨量充沛，气温适宜，有利于林木生长，生态系统稳定，碳汇能力强。此外，这些地区的人为干扰相对较少，一定程度上降低了人类活动对林木生长的影响，碳储存潜力大，因此碳储量相对较高。

根据模型预测碳储量均值变化(图3)，2012年碳储量均值为12.91 t·hm⁻²，2017年为13.14 t·hm⁻²，2021年为13.63 t·hm⁻²，说明该区域碳储量呈现增长趋势。这得益于生态保护相关政策的持续推进，区域内可能有序开展植树造林、退化林修复、退耕还林等生态工程，过度砍伐、开垦等破坏植被的行为得到遏制，恢复随时间推移，树木持续生长发育，生物量不断累积，碳储存能力逐步增强，使得“三江并流”区域碳储量随年份呈现增长态势。

Figure 3.  Trend of changes in mean carbon stocks estimation from 2012 to 2021

传统GWR仅能刻画参数在空间维度的异质性。为同时捕捉时空动态，GTWR引入时间维度，使空间权重不仅依赖地理距离，还考虑时间序列上的关联性，进一步允许参数随时间变化，反映了时空异质性。然而，GTWR本质上仍是线性模型，对于复杂的非线性关系刻画能力有限^[24]。为此，本研究引入GBRT，该算法通过集成多棵回归树，凭借高精度、灵活性和鲁棒性，成为机器学习中的核心算法之一^[25]，能更好拟合碳储量与变量之间的非线性关系。基于此，构建了GTWR-GBRT组合模型，其核心优势在于兼具时空动态建模的地理可解释性，与非线性关系的高精度拟合能力。本研究对比了GTWR、RF、GBRT及组合模型GTWR-RF、GTWR-GBRT，结果显示GTWR-GBRT组合模型的效果最佳，说明时空异质性模型结合非线性模型估测森林碳储量具有优势。

研究区80%的碳储量普遍分布在海拔2 500 m以上，其中海拔3 500 m以上的碳储量占比最高，高达33%。西部的碳储量最高，而研究区属于云南西北部，由于该地区海拔高，人口稀疏且人为干扰小，植被保存较好，森林碳储量普遍高于东部地区，这表明地质环境条件和人类活动干扰是造成云南省碳储量差异的重要因素^[26]。本研究结果与CHENG等^[27]和梁瑄等^[28]的估算结果一致：碳储量主要集中在滇西部的普洱市、迪庆州、西双版纳州和怒江州等地。此外，丽江市碳储量较低，这与范华鹏等^[29]采用RF模型的预测结果相符。整体来看，在2012—2021年期间，2017年森林碳储量的高值分布较多，说明天然林资源保护，退耕还林、还草措施和重点防护林建设对森林碳储量的增加有一定的效果。总体而言，模型预测的碳储量均值说明该区域碳储量呈现增长趋势。结合空间分布来看，西部地区的碳储量最高，应优先在碳储量丰富的地区进行保护，以实现最佳的固碳效益，东部地区碳储量相对较少，可通过实施可持续森林管理做法来增加东部的碳储存。后续可加强低海拔地区植被恢复，同时重点保护高海拔地区脆弱生态系统(应对气候变化、减少人为干扰)，巩固区域碳汇能力。

在森林碳储量估算的过程中，传统光学遥感影像受饱和效应的影响，常存在低值被高估而高值被低估的情况，从而影响估算精度和反演结果。本研究利用Landsat卫星数据获取了遥感因子，虽然有长时间的遥感影像，但其分辨率为30 m，对森林碳储量的精准估测有一定的影响。穆喜云等^[30]利用机载LiDAR遥感数据，使用RF构建预测模型及经典的统计方法-多元线性回归，对内蒙古大兴安岭生态站的森林碳储量进行估测，获得了比较理想的结果。此外，在模型构建方面，GTWR需要选择合适的带宽来保证时空加权估计的准确性，GBRT需调整树的数量、深度等参数。若联合调参不当可能导致过拟合或陷入局部最优。后续研究可以在优化模型参数的基础上，结合高分辨率数据，以获得最佳预测性能。

为克服单一模型在捕捉时空异质性与非线性关系上的局限，本研究构建了两阶段GTWR-GBRT混合模型，融合GTWR的时空动态能力与GBRT的非线性拟合优势。该模型估算精度优于GTWR、GBRT、RF及GTWR-RF，说明同时考虑时空异质性与非线性特征是提升森林碳储量估算可靠性的关键。基于该模型估算结果，2012—2021年“三江并流”区域森林碳储量的空间格局呈“西高东低、局部集聚”，高海拔区域碳储量普遍较高，显示出明显的时空异质性。