基于ANUSPLIN插值的星载LiDAR（GEDI与ICESat-2/ATLAS）提升巨竹地上碳储量估算精度研究

本研究创新性地将ANUSPLIN插值技术应用于GEDI与ICESat-2/ATLAS星载LiDAR数据，显著提升了巨竹（Dendrocalamus giganteus）地上碳储量（AGC）的区域尺度估算精度。通过对比协同克里金（CK）方法，证明ANUSPLIN在空间扩展与可视化表现上均更优；结合XGBoost模型（R2=0.93，RMSE=5.89 Mg/ha），实现了高精度碳储量制图，为竹林碳监测与全球气候变化研究提供了可靠方法支撑。

　　

1 引言

地上碳储量（AGC）是森林生态系统中最基本的数量特征之一，反映了物质循环、能量流动及植物与环境相互作用的复杂关系。竹林作为“世界第二大森林”，是森林生态系统的重要组成部分。中国作为“竹子王国”，竹林资源丰富，2021年竹林碳储量达2.10×10⁸ Mg，占全国森林碳储量的1.96%。准确估算竹林AGC对理解全球气候变化机制、制定碳排放政策及减缓全球变暖具有重要意义。

传统森林AGC估算方法主要依赖野外调查，耗时耗力，难以满足快速、大尺度获取AGC分布数据的需求。遥感技术的发展为森林碳储量监测提供了快速高效的手段。当前竹林生物量与碳储量遥感研究多采用被动光学遥感数据，但由于竹林生长密度高、枝干交错，光学遥感估算精度较低。例如，基于Landsat TM数据估算毛竹地上生物量（AGB），最大相关系数仅为0.48；利用Sentinel-2数据和随机森林回归（RFR）估算浙江省竹林AGB，R2为0.46。机载激光雷达（LiDAR）基于主动遥感技术，可提高估算精度（R2=0.64），但成本高昂，不适于大尺度应用。星载LiDAR数据（如GEDI和ICESat-2/ATLAS）为区域尺度竹林AGC估算提供了低成本、高精度的解决方案，但其离散采样特性无法提供全覆盖数据，需借助外推方法。空间插值作为获取连续空间分布信息的有效手段，已广泛应用于星载LiDAR数据处理。

目前，地统计方法如协同克里金（CK）和序贯高斯模拟已广泛用于星载LiDAR数据的空间插值，但ANUSPLIN方法的应用相对较少。插值方法的有效性受研究区域影响，且选择标准不统一。ANUSPLIN基于薄板样条理论，最初用于气象数据插值，可引入高程数据作为协变量。与地统计插值方法相比，ANUSPLIN无需事先校正半变异函数参数，操作更简便，插值精度更高。在气温与降水研究中，ANUSPLIN表现优于普通克里金和反距离加权法，尤其在复杂地形区域。

2 材料与方法

2.1 研究区

研究区位于云南省新平县，地处哀牢山东坡，地势西北高、东南低，地形复杂，土壤肥沃。海拔从低谷到中山变化，生境多样。气候属亚热带至温带，年均温约18°C，年均降水量869 mm，年日照2839 h，无霜期约316天。林地面积23520.1 ha，占全县面积的55.8%，森林覆盖率70.99%。植被资源高度多样，包括1402种高等植物。竹林面积2417.75 ha，涉及12个乡镇，以巨竹为主，广泛分布于450–1800 m海拔范围。

2.2 巨竹AGC测量

野外调查于2024年1月进行，依据代表性和可及性原则，预设80个巨竹样地，最终建立51个圆形样地（半径12.5 m，面积490.63 m2）。记录每样地的胸径（DBH）、地理位置和植株数。AGC估算分三步：确定代表平均木的AGB；通过碳含量系数将个体AGB转化为AGC；根据样地内总植株数计算样地水平AGC。巨竹AGC计算公式为C_t=W×f_c，其中碳系数竹杆、竹枝和竹叶分别为0.45、0.45和0.43。51个样地的巨竹AGC最小值、最大值、平均值和标准差分别为4.08 Mg/ha、101.78 Mg/ha、41.63 Mg/ha和20.55 Mg/ha。

2.3 全波形数据

GEDI是首个星载全波形LiDAR系统，搭载于国际空间站，采样范围51.6°N至51.6°S。本研究使用GEDI L2B数据产品，空间分辨率25 m，覆盖研究区2022年1月至2023年2月。共有57217个激光脚印样本，系统采样（间隔24）后得到2384个预设样本。

2.4 光子计数

ICESat-2数据产品分为4级（ATL00–ATL22），本研究主要使用ATL03和ATL08产品。ATL08由ATL03经DDBSCAN和KNNB算法处理，再通过改进的PTD分类方法生成，提供植被冠层和地形参数。研究区共21080个激光脚印样本，系统采样（间隔9）后得到2342个预设样本。

2.5 辅助数据

数字高程模型（DEM）采用ALOS PALSAR RTC 12.5 m DEM，通过NASA Earthdata Search获取，重采样至25 m空间分辨率以匹配脚印和样地面积。

3 研究方法

3.1 LiDAR数据面生成

3.1.1 ANUSPLIN插值方法

ANUSPLIN基于普通和局部薄板样条插值理论，利用平滑样条函数对多元数据进行插值，确保插值曲面的平滑性和准确性。局部薄板平滑样条插值允许引入自变量和协变量（如高程）。计算涉及最小二乘估计，通过广义交叉验证（GCV）确定平滑参数。具体步骤包括：将插值数据导入SPSS并导出批处理文件；将数据文件放入ANUSPLIN插件文件夹，根据模型修改sp.txt文件中的样条度数、自变量和协变量参数；运行run.cmd进行插值；转换生成的grd文件格式，处理异常值，得到最终结果。

ANUSPLIN允许多种自变量、协变量和样条度数组合，共18种模型。本研究以GEDI和ICESat-2/ATLAS参数为自变量，高程为协变量，样条度数设为2、3和4，进行最优空间插值。基于原始数据和ANUSPL用户指南，编写了运行SPLINA和LAPGRD程序模块的批处理命令文件。SPLINA生成曲面拟合系数和误差统计文件，LAPGRD利用该系数文件获取插值曲面。ANUSPLIN在日志文件中提供一系列统计参数，用于诊断误差来源和评估插值质量，包括有效拟合曲面参数数量、信号自由度（Signal）、残差自由度（Error）和信噪比（SNR）。

最佳模型选择标准：基于GCV或最大似然法，SNR应最小化；Signal应小于观测站点数的一半；SPLINA生成的误差文件中不应出现“”符号。Signal代表拟合曲面复杂度，若超过观测站数的一半，表明数据误差大或存在不适于曲面模型的数据，此时拟合过程无法找到最优平滑参数，这些问题在误差统计文件中以“”标记。

本研究采用ANUSPLIN提供的GCV方法验证插值结果。评估标准为：GCV、均方残差（MSR）、方差（VAR）、样条均方根误差（RTGCV）、样条残差均方根误差（RTMSR）和样条方差（RTVAR）越接近0，且RTGCV值越接近GCV值、RTMSR接近MSR、RTVAR接近VAR，拟合精度越高。

3.1.2 CK插值方法

CK是普通克里金的扩展，不仅利用预测点的空间自相关性，还利用主变量与多个协变量之间的关系。

3.2 插值结果精度验证

采用交叉验证评估插值方法准确性，使用四个指标：平均误差（ME）、均方误差（MSE）、均方根误差（RMSE）和平均标准误差（ASE）。ME和MSE越接近0，预测值偏差越小；RMSE最小化且ASE接近RMSE时，精度提高。

3.3 AGC估算模型及精度评估

选择三种机器学习回归模型：LightGBM、RFR和XGBoost。LightGBM在R 4.3.1环境中通过RStudio调用“LightGBM”包，网格搜索优化超参数（num_leaves=20, max_depth=5, learning_rate=0.1）。RFR使用“randomForest”包，设置importance=TRUE获取变量重要性，采用%IncMSE进行特征选择，网格搜索优化超参数（mtry=9, ntree=200）。XGBoost使用“XGBoost”包，网格搜索优化超参数（Nrounds=45, max_depth=3, eta=0.1）。采用十倍交叉验证，使用四个精度指标评估模型性能：决定系数（R2）、RMSE、总体估算精度（P）和相对RMSE（rRMSE）。

4 结果

4.1 特征变量选择

基于RFR模型的特征重要性排序，按%IncMSE降序排列变量。逐步添加变量后，当变量数增至9时，RMSE最小（8.23 Mg/ha）。最终选择的9个变量包括4个GEDI参数：cover（13.35%）、digital_elevation_model（9.88%）、pai（11.54%）和sensitivity（13.48%），以及5个ICESat-2/ATLAS参数：toc_roughness（10.80%）、h_median_canopy_abs（11.22%）、asr（14.53%）、h_canopy_abs（8.86%）和h_max_canopy_abs（8.80%）。

4.2 LiDAR数据面生成与比较

4.2.1 ANUSPLIN

GEDI参数cover、pai、sensitivity和ICESat-2/ATLAS参数toc_roughness、h_median_canopy_abs、h_canopy_abs、h_max_canopy_abs在样条函数数为4时模型精度最高，显著高于样条数2和3。而GEDI参数digital_elevation_model和ICESat-2/ATLAS参数asr在样条数为2时精度更高。除digital_elevation_model外，其余8个变量的评估指标（GCV、MSR、VAR、RTGCV、RTMSR、RTVAR）均趋近于0，表明ANUSPLIN插值适用于星载LiDAR空间扩展，实验结果可靠。

4.2.2 CK

基于半变异函数最优原则选择插值模型。GEDI四个参数（cover、digital_elevation_model、pai、sensitivity）最优模型为高斯模型；ICESat-2/ATLAS五个参数（toc_roughness、h_median_canopy_abs、asr、h_canopy_abs、h_max_canopy_abs）最优模型为球状模型。九个参数中，仅GEDI的digital_elevation_model和ICESat-2/ATLAS的h_median_canopy_abs在25%–75%范围内存在块金效应，表现为中等空间自相关性，其余七个参数表现出强空间自相关性。

4.2.3 ANUSPLIN与CK插值精度比较

九个参数中，除digital_elevation_model的MSE和ASE在ANUSPLIN高于CK外，其余参数ANUSPLIN精度均优于CK。插值结果显示，ANUSPLIN生成的曲面更平滑，更好地保留了局部细节，一定程度上减轻了条带效应。CK由于经验协方差函数仅近似空间协方差变化，拟合过程导致局部细节丢失，曲面呈现机械梯度或不均匀模式。ANUSPLIN插值表现出强烈的自然平滑性，有效捕捉地形特征，在复杂地形区域清晰勾勒出山地轮廓。

4.3 区域尺度模型估算与制图结果

将选择的9个GEDI和ICESat-2/ATLAS自变量与51个样地数据输入LightGBM、RFR和XGBoost模型。XGBoost模型（R2=0.93, RMSE=5.89 Mg/ha, P=85.84%, rRMSE=14.16%）优于LightGBM（R2=0.52, RMSE=14.61 Mg/ha, P=64.84%, rRMSE=35.16%）和RFR（R2=0.90, RMSE=8.23 Mg/ha, P=79.79%, rRMSE=20.21%），相对LightGBM提升78.85%、59.67%、32.36%和59.71%，相对RFR提升3.33%、28.41%、7.58%和29.94%。

基于XGBoost模型，绘制新平县巨竹AGC分布图。51个样地预测平均AGC为41.59 Mg/ha。全县平均AGC为40.62 Mg/ha，范围15.70–81.96 Mg/ha，总储量1.14×10⁷ Mg。空间上，巨竹主要集中分布在县西北部，其他区域分布较少。巨竹主要分布于373–1755 m海拔，仅少量出现在1755 m以上，与野外观测一致。

5 讨论

5.1 不同传感器对竹林AGC估算的影响

基于光谱数据的AGC预测精度通常低于LiDAR数据，因多光谱数据仅提供冠层结构水平信息，无法捕获垂直结构信息。而垂直结构参数（如树高和冠层高度）与森林AGC密切相关。LiDAR能快速获取三维空间信息，穿透能力强，可精确监测垂直冠层结构参数，在森林参数反演、生物量和碳储量估算中优势显著，大大提高了估算精度。与CK插值相比，ANUSPLIN插值进一步提升了XGBoost和RFR模型的精度，表明ANUSPLIN插值能有效增强星载LiDAR数据的空间外推能力，提高模型预测精度，为竹林碳储量监测和区域尺度森林碳汇评估提供了更稳健的技术支持。

5.2 插值方法对地图可视化效果的影响

本研究首次将ANUSPLIN插值方法用于星载LiDAR空间扩展，插值精度验证与其他插值精度结论一致，结果可靠。DEM分辨率影响插值精度，本研究使用12.5 m DEM，重采样至25 m，有利于提高ANUSPLIN插值精度。样条数量显著影响结果精度，ANUSPLIN插值高度依赖模型选择。通过ANUSPLIN插值星载LiDAR数据并基于既定原则选择最优模型，交叉验证结果表明该方法优于地统计插值方法。从可视化角度看，ANUSPLIN结果优于CK，因CK结合已知点和半方差分析构建统计模型，易受统计规律影响，产生条带效应，尤其在斑点众多、空间异质性高的区域。ANUSPLIN则利用控制点拟合最小曲率曲面，结果更平滑。相比CK，ANUSPLIN插值在插值曲面准确性、平滑性和细节保持上更优，且效率更高，需较少人工参数调整，降低了插值过程的不确定性。

5.3 星载LiDAR不同点密度对插值的影响

样本量影响模型精度，通常样本量越大，模型可靠性越好，但过大样本量浪费资源。星载LiDAR点密度亦然：点越多，插值时间越长，条带效应越明显。因星载LiDAR脚印点沿轨道均匀分布，插值前需根据选择标准进行质量筛选，过滤同一或相邻条带的低质量点，增强脚印点的空间随机性。即使预处理后，GEDI和ICESat-2/ATLAS的剩余点仍不满足空间随机性要求。一些研究者在插值前进行点稀释以减轻条带效应，且不影响精度。未来研究可探索插值的最优点数。当样本量小时，CK插值表现相对较好，图谱更平滑，因样本与网格点距离增加，对应于经验半方差函数的平滑部分。因此，未来研究可探讨点密度对插值的影响，寻求最优点数以减少工作量。

5.4 变量选择对模型精度的影响

变量选择是模型开发的前提，建模变量的选择直接影响反演模型的预测性能。通过选择相关变量，可剔除无贡献变量，降低数据维度，简化模型，提高模型精度和泛化能力。估算森林结构参数时，不同变量选择方法导致模型预测精度变化，取决于数据组合，进而影响估算结果。由于巨竹AGC与星载LiDAR参数关系复杂，变量选择至关重要。本研究采用随机森林特征重要性方法进行变量选择，确保优先选择解释性强、影响力大的变量。通过该特征重要性排序，逐步将变量添加至模型，当变量数达9时，RMSE最小（RMSE=8.23 Mg/ha）。未来研究可探索更多方法，如SVM-RFE、Boruta和KNN-FIFS，进行特征选择，以确定更优参数，提高星载LiDAR估算模型的预测精度。

5.5 未来研究展望

本研究的方法学亮点是首次将ANUSPLIN插值方法用于星载LiDAR数据的空间外推，证明了其相对于传统地统计CK方法的优越性。此外，开发了基于XGBoost的模型，以更高精度估算竹林AGC。结果对竹林碳储量监测、森林资源调查和区域碳汇评估具有重要意义。但本研究存在一些局限性，如协变量考虑不全、仅关注单一物种（巨竹）、样地数量有限。未来研究可引入更多地形和环境因子，增加样地数和LiDAR脚印数，并在不同森林类型和更大区域验证该方法，为森林碳储量估算和碳中和目标提供更可靠的技术支持。

6 结论

本研究系统探讨了利用GEDI和ICESat-2/ATLAS数据的空间插值方法，并评估了其在巨竹AGC估算中的性能。结果表明，ANUSPLIN插值优于CK插值。当基于随机森林特征重要性选择9个关键变量（GEDI：cover、digital_elevation_model、pai、sensitivity；ICESat-2/ATLAS：toc_roughness、h_median_canopy_abs、asr、h_canopy_abs、h_max_canopy_abs）时，RMSE降至8.23 Mg/ha。机器学习模型中，XGBoost表现最佳（R2=0.93, RMSE=5.89 Mg/ha, P=85.84%, rRMSE=14.16%），预测新平县巨竹平均AGC为40.62 Mg/ha，总AGC为1.14×10⁷ Mg。本研究表明，结合ANUSPLIN与星载LiDAR的框架在大尺度森林碳储量监测中潜力巨大，为气候变化研究提供了宝贵参考。