基于无人机多光谱的高山松地上生物量估测

杨正道，舒清态，黄金君，周文武，胥 丽，罗绍龙，王书伟

（1.西南林业大学 林学院，云南昆明 650224；2.广西壮族自治区中国科学院广西植物研究所，广西桂林 541006）

森林是陆地重要的生态系统之一，发挥重要作用，不仅具有丰富的生物多样性和提供众多生态服务，还是地球上最大的碳库[1]。森林生物量是森林生态系统的重要组成部分，包括地上生物量和地下生物量；其中，地上生物量由树干生物量、树枝生物量和树叶生物量组成，一般用单位时间或单位面积内累积的干物质或能量表示[2]，在全球森林生态系统效益评估中发挥主导作用[3-5]。对森林地上生物量进行估测对于了解碳循环、评估碳汇潜力及制定森林资源管理和林业可持续发展策略有重要价值和意义[6-7]。

森林地上生物量估测方法主要包括野外实测和遥感估测[8]。基于样地调查的野外实测生物量估测方法准确性高，但存在一些缺点，如消耗大量时间、人力和资金，破坏性较大，在偏远和崎岖地区实地调查困难等，不适用于大范围森林地上生物量估测和动态变化监测[9-10]。随着地理信息技术不断发展，遥感（Remote Sensing，RS）数据在森林地上生物量估测方面得到广泛应用。遥感连续动态监测和地理信息系统（Geographic Information System，GIS）空间数据分析能力不断增强，学者们结合地面数据和卫星遥感数据进行森林地上生物量估测，主要利用的卫星遥感数据为Landsat 8光学遥感影像[11]。尽管卫星遥感技术在森林地上生物量估测方面取得了一定进展，但仍存在一些局限。首先，获取生物量像元的真实值并不容易；森林地上生物量是通过间接推断得出的，通常需依赖实地调查数据进行校正和验证，这一过程需耗费大量时间和人力资源。其次，卫星遥感数据的时间分辨率和空间分辨率相对较低，可能无法捕捉森林地上生物量的细微变化和误差导致的空间异质性[12]。

近年，无人机遥感技术以其低成本、高安全性、高时效性、高分辨率和低空影像获取等优势弥补了卫星遥感技术的局限[13]，在林业方面得到广泛应用。申鑫等[14]利用无人机获取的高光谱和高空间分辨率数据提取遥感因子，并结合样地实测数据，建立亚热带天然次生林林分尺度生物量估测模型，进行生物量反演；Jones 等[15]利用无人机影像生成高分辨率三维树冠模型，在林分尺度上对红树林地上生物量进行估测；Tamiminia 等[16]利用无人机多光谱影像进行纽约州柳树（Salixspp.）林分尺度上森林地上生物量的反演；李滨等[17]基于无人机遥感技术获取目标样地樟子松（Pinussylvestrisvar.mongolica）的遥感影像，并结合样地实测数据，在单株尺度上建立胸径-生物量模型，用于估测东北林业大学城市林业示范基地内目标样地的地上生物量；杨雪峰等[18]利用无人机低空遥感和高分辨率卫星影像，获取塔里木河下游胡杨（Populuseuphratica）林单株尺度的森林结构参数，估测该地区的地上生物量。目前的研究主要集中在单株和林分尺度上的森林地上生物量估测，对于最佳观测窗口和低纬度高海拔地区森林地上生物量估测，尚缺乏较深入的研究。探索最佳观测窗口，并在不同地理环境和生态系统中验证无人机遥感技术的适用性，有助于推动生物量估测方法的发展及其精确性的提高。

本研究通过地统计学变异函数确定高山松（Pinusdensat）地上生物量最佳观测窗口；采用偏最小二乘法（Partial Least Squares Regression，PLS）和随机森林（Random Forest，RF）模型对遥感影像中筛选的特征变量和地面实测生物量建模，对比分析两种模型的精度；采用筛选出的最优模型对无人机飞行区高山松地上生物量进行估测，并分析其空间分布情况，可为进一步开展高寒山区森林生态系统物质与能量循环研究及保护和改造脆弱生态环境提供参考。

1 材料与方法

1.1 研究区概况

研究区位于云南省迪庆藏族自治州香格里拉市小中甸镇下吉沙林场（99°49′E，27°24′～27°25′N），平均海拔3 300～3 500 m，总面积76.39 hm2，年均气温6.3 ℃，年均降水量649.4 mm。林场主要地貌为山地，地形呈东北高、西南低的趋势。高山松为林场从寒温带向亚热带过渡的先锋针叶树种，适应性广，再生能力强，具有喜光、耐旱和耐瘠薄等特点。具体样地分布如图1所示。

图1 研究区位置与样地分布Fig.1 Location of study area and distributions of sample sites

1.2 数据采集

1.2.1 无人机影像数据获取与预处理

采用大疆四旋翼无人机搭载的航测系统和AURedEdge 多光谱相机（含蓝光、绿光、红光、红边和近红外5 个波段，表1）采集试验区高山松多光谱影像。航拍过程中，相机面向地面90°俯拍，地面控制点采用RTK 测定，坐标系统为GCS—WGS—1984。航飞时间为2021 年11 月1 日10:00～14:30，天气晴朗，无风云干扰。通过5 次试飞和航测点调整，确定无人机飞行高度为260 m，航速为8.8 m/s。航向和旁向重叠度分别为75%和70%。地面分辨率为0.128 m，航测区面积为63.73 hm2。

表1 多光谱传感器波段参数Tab.1 Band parameters of multi-spectral sensors

获得无人机多光谱影像后，进行预处理。对无人机拍摄到的多光谱影像进行甄别，剔除姿态角不正常、成像有问题的影像。采用大疆智图软件（DJI Terra）进行影像数据预处理。将多光谱相机拍摄的原照片导入软件进行二维多光谱重建，获得不同波段反射率图像、数字表面模型（Digital Surface Model，DSM）和数字正射影像（Digital Orthophoto Map，DOM）。采用ENVI 5.3 软件进行校正，并将数据存储为栅格格式（图2）。

图2 无人机航测区影像Fig.2 Aerial survey area image of UVA

1.2.2 地面调查数据获取与预处理

根据无人机获取的多光谱影像及其尺度，采集对应时间内试验区36 个高山松角规控制检尺样地的胸径、树高和株数等信息，记录样地中心经纬度和高程信息。采用三鼎星T66系列RTK进行坐标点测量，误差控制在千分之二以内；计算每个样地的平均胸径、平均树高、总蓄积量和断面积。

参考黄金君等[19]文献，采用分层贝叶斯法拟合单木生物量模型，通过公式（1）计算不同样地地上单木平均生物量；通过公式（2）计算每公顷林木株数，并结合地上单木平均生物量，求算样地公顷生物量；将样地公顷生物量转换为与影像像元对应的样地生物量。36块样地数据描述性分析见表2。

表2 高山松样地地上生物量统计分析Tab.2 Statistical analysis on above-ground biomass of P.densata sample sites

式中，M为地上单木平均生物量（t/hm2）；D为平均胸径（cm）；H为平均树高（m）。

式中，N为每公顷林木株数（株/公顷）；Fg为断面积系数；k为样地计数木总株数；gi为样地中第i株计数木的断面积（cm2）；ki为样地中第i株计数木的株数。

1.3 研究方法

1.3.1 最佳观测窗口选择

最佳观测窗口常被称为最佳分辨率，即影像中能识别最小地物的最大空间分辨率，能清楚探索和表述这一尺度层的地理过程和现象内容，同时也能明确地分析这一尺度层地理空间信息的组合规律[20]。

半方差也称变异函数，是建立在二阶平稳假设基础上的空间变异结构分析工具，也是地统计学的主要理论基础，在地物空间变异分析方面得到广泛应用。1997 年，Atkinson[21]将变异函数法应用到遥感影像最佳观测尺度研究中，通过分析变异函数中变量因子变异的空间程度，确定最佳分辨率。传统统计学仅对像素值进行统计分析，忽视像素值间的随机性和相关性；变异函数可以揭示区域化变量在全尺度上的空间变异模式[21]，计算公式为：

式中，γ(h)为区域变量Z（w）的变异函数值；N（h）为空间间隔距离等于h的样本点对数量；Z（w）i和Z（wi+h）分别为区域变量在点wi和点wi+h位置上的像素值；h为两个像素点的分隔距离。

在变异函数的计算中，结果值为散点图。在对整个区域进行定量描述时，需对散点进行拟合。通常采用指数模型、球状模型和高斯模型3 种方式进行拟合[20，22]。

（1）指数模型表达式为：

h=3α时，γ（h）≈C0+C1，此模型的有效变程为3α。

（2）球状模型表达式为：

h=α时，γ（h）=C0+C，此模型的有效变程为α。

（3）高斯模型表达式为：

式中，C0为块金值；C1为偏基台值；C（C=C0+C1）为基台值；α为变程或相关距。模型的变程α即最佳观测窗口。实际应用时，以高山松重采样尺度大小（即重采样时像元大小）为h，以对应尺度下的像元光谱值为γ（h），求算出的变程α为最佳观测窗口。在模型优选基础上，利用块金值、基台值、变程和结构比等模型参数，对区域化变量进行空间分布特征分析，一般选用具有较大决定系数和较小残差的理论模型拟合变异函数[23]。

变异程度为块金值与基台值的比值（C0/C），能反映整个空间的变异状态，主要包括0～25%、25%～75%和75%以上3 个层次，代表空间自相关程度。空间自相关程度高（变异程度低），表明主要由结构性因素引起变异；空间自相关程度低（变异程度高），表明以随机部分为主。

1.3.2 最佳观测窗口设计方案

选取拼接后的无人机多光谱影像，筛选与实测生物量均值最接近的样方作为中心点，在ENVI 5.3软件中裁剪出大小为100 m×100 m、质量较好的一景作为观测尺度试验影像样本，利用Band Math 工具对试验影像样本进行纹理特征指数计算。

本研究中，无人机多光谱的空间分辨率为0.125 m，数据量大，冗余多，所以利用ENVI 5.3软件中的Resize Data工具对进行纹理特征指数计算后的影像样本进行重采样，将空间分辨率提升至1 m。常见的重采样内插方法有最邻近法、双线性内插法和三次卷积内插法。本研究选择最邻近法进行重采样。最邻近法将最邻近的像元值赋予新像元，输出的图像仍保持原来图像的像元值，处理速度快。

采用ArcGIS 软件将重采样后栅格数据中每个像元的纹理特征值提取至点，计算每个点的X和Y坐标；借助地统计学GS+软件，分别采用指数模型、球状模型和高斯模型作为变异函数拟合模型，对高山松地上生物量最优观测尺度进行分析，通过模型决定系数（R2）和残差平方和（Residual Sum of Squares，RSS）对模型进行评价，获得最佳观测窗口。R2越大，RSS越小，模型效果越好。

1.3.3 遥感因子提取

研究表明，在高分辨率影像中，生物量与纹理特征的相关性较强，且纹理特征不会出现植被指数的光饱和特性[24-27]。采用二阶纹理算法中的灰度共生矩阵（Gray Level Co-occurrence Matrix，GLCM）提取纹理特征，提取均值（Mean，Mea）、相异性（Dissimilarity，Dis）、角二阶矩（Second Moment，SM）、方差（Variance，Var）、对比度（Contrast，Con）、熵（Entropy，Ent）、同质性（Homogeneity，Hom）和相关性（Correlation，Cor）8种纹理特征变量，具体计算公式见文献[28-29]，纹理窗口大小为3×3、5×5、7×7、9×9和11×11。

1.3.4 生物量反演模型构建

随机抽取70%数据集采用PLS和RF方法建模，剩余30%数据集用于验证。同一种变化类型中，PLS和RF模型输入的建模因子相同。

PLS 是一种集合主成分分析、典型相关分析和多元线性回归分析的数据统计方法，可有效解决多元回归分析中变量的多重相关性和噪声，在此基础上具有预测功能并能消除参数结构不确定和模型无法辨识的影响[30]。

RF 是一种采用多棵决策树构建的增强型分类（回归器），当新的数据被带入随机森林后，全部决策树生成分类或预测结果，随机森林以这些结果的众数或均值为数据输出，可在不选取特征的情况下对高维数据进行处理。由于随机选取样本，在每次学习中，决策树均会采用不同的训练集进行训练，所以随机森林能一定程度避免过拟合现象[31]。本研究中，随机森林模型算法的实现基于R语言平台。

1.3.5 模型精度评价

选取R2、均方根误差（Root Mean Squared Error，RMSE）和估测精度（P）为模型精度评价指标。R2越大，RMSE越小，模型拟合度越好。计算公式[19，30]为：

式中，yi为样本实测值；i为模型预测值；为样本实测值均值；n为样本数。

2 结果与分析

2.1 最佳观测窗口选择

3 种模型中，球状模型R2最大（0.889），RSS 最小（0.93×10-6），说明球状模型拟合结果最佳；其最佳观测窗口为5.2 m（表3）。球状模型的变异程度为7.7%，表明高山松无人机多光谱影像空间尺度变异主要由结构性因素引起，其内在空间自相关性强。

表3 变异函数模型参数Tab.3 Parameters of variance function models

2.2 遥感因子相关性分析

以高山松地上生物量为因变量，提取的遥感因子为自变量；采用Pearson 相关系数作为评价指标，其绝对值在0～1 之间，越靠近1，相关性越强。为获取相关性强的因子，将窗口依次设置为3×3、5×5、7×7、9×9和11×11；优选出每个窗口相关系数较强的前5个因子（表4）。高山松地上生物量模型中，红边与近红外波段敏感性较强，相关系数为0.358～0.442，在提取窗口为9×9时，相关性较强。

表4 遥感因子与生物量的相关系数Tab.4 Correlation coefficients among remote sensing factors and biomass

在以上筛选的基础上，将相关系数≥0.428 设为阈值，且在0.01水平显著；5个相关性较高的遥感因子依次为RE9x9_Dis、RE9x9_Con、RE9x9_Var、NIR9x9_Dis和RE5x5_Con，相关系数分别为0.442、0.439、0.430、0.429和0.428。

2.3 飞行区高山松地上生物量模型选择

基于RF 构建的生物量模型优于PLS，R2较大（0.90），RMSE 较小（17.96 t/hm2），P较大（84.98%）（表5，图3）。选择RF 模型进行飞行区高山松地上生物量估算和反演。

表5 模型精度对比Tab.5 Comparisons on model accuracies

图3 预测值与实测值比较（a：RF；b：PLS）Fig.3 Comparisons between predicted values and measured values(a:RF;b:PLS)

2.4 飞行区高山松地上生物量反演

飞行区高山松地上生物量均值为130.48 t/hm2，总生物量为8 343.53 t（图4）。反演结果显示，飞行区高山松地上生物量空间分布呈中间高、四周低的趋势，生物量集中在123.99～147.32 t/hm2，少数处于172.48～211.13 t/hm2，与实际情况相符，反演结果较好。

图4 飞行区高山松地上生物量空间分布Fig.4 Spatial distributions of P.densata above-ground biomass in flight area

3 讨论与结论

本研究采用无人机多光谱影像为遥感数据源，结合36块样地实测数据，首先通过地统计学变异函数确定高山松地上生物量最佳观测窗口；结果显示，变异函数3 种模型中，球状模型拟合效果最佳，其最佳观测窗口为5.2 m。其次，采用ArcGIS 软件对无人机多光谱影像进行重采样至5.2 m，提取200个纹理特征参数，计算其与高山松生物量的Pearson相关系数，并选取相关系数≥0.428 且在0.01 水平显著的5 个纹理特征因子，建立PLS 和RF 地上生物量模型；结果显示，RF 模型优于PLS 模型。基于RF 模型，无人机飞行区的高山松地上生物量均值为130.48 t/hm2，总生物量为8 343.53 t。本研究结果可为高寒山地森林生物量便捷、准确估测提供技术参考。本研究中，样地生物量均值为130.48 t/hm2，超出谢福明[2]研究中香格里拉市高山松地上生物量（17.45～107.46 t/hm2），可能是因为本研究区处于林场，人工抚育强度大，保护措施较好，树木长势好。

因条件限制，本研究布设样地为1 hm2，无人机多光谱影像分辨率为厘米级，将样地公顷生物量推算至影像像元尺度生物量时，结果存在一定误差。

确定最佳观测尺度的方法主要有3 种，包括变异函数、局部方差和离散度。本研究仅探讨了基于变异函数的最佳观测尺度，对基于局部方差和离散度的最佳观测尺度有待进一步探索。

利益冲突：所有作者声明无利益冲突。

作者贡献声明：杨正道负责试验调查与设计、数据收集与分析、论文撰写和文献检索；舒清态负责论文指导与修改；黄金君负责试验调查与设计和数据收集与分析；周文武负责试验调查与设计；胥丽、王书伟负责文献检索；罗绍龙负责数据收集与分析。