邓再春,张 超,2,朱夏力,范金明,钱 慧,李成荣
(1.西南林业大学 林学院,云南 昆明 650224;2.西南林业大学 云南省山地农村生态环境演变与污染治理重点实验室,云南 昆明 650224)
森林蓄积量是一定森林面积上所有活立木材积的总和[1]。作为森林生物量和碳储量的重要评价指标,森林蓄积量能直接反映森林资源的数量与质量,是森林资源调查的重要因子之一。随着无人机遥感技术的快速发展,无人机载可见光/多光谱遥感影像在森林资源调查领域得以快速推广,为森林蓄积量调查提供了快速高效的技术手段。传统的森林蓄积量调查主要以地面调查为主,此类调查周期长,对人力、物力的需求量巨大[2]。20 世纪90 年代以来,国内外学者通过获取单一或多源遥感影像,以地面调查数据作为因变量,以各类植被指数、纹理特征、地形因子等作为自变量,采用主成分分析、偏最小二乘法、逐步回归、随机森林、k邻近模型等方法建立回归模型,估测森林蓄积量,进行了较多有益的探索[3-5]。卫星遥感影像具有长时序、大尺度、易获取等优势,但其影像易受天气影响,且难以兼顾分辨率和成本。无人机具有成本低、机动灵活、影像分辨率高等优点,作为传统遥感估测的补充手段,在森林资源调查中得到了广泛应用。通过搭载可见光、多光谱、高光谱、激光雷达等多种传感器,可获得低空地表有关森林资源的多层面数据[6]。大量研究表明:基于无人机航拍影像估测森林蓄积量具有较高的可行性。利用无人机航拍影像估测森林蓄积量主要包括2 个角度:① 基于数字正射影像(DOM)、数字表面模型(DSM)、冠层高度模型(CHM)获取林分株数、胸径、冠幅、树高等因子,从单木、林分2 个角度进行蓄积量估测[7-12];② 基于影像提取各类植被指数、纹理特征、地形因子等,建立森林蓄积量估测模型,或是建立单株材积估测模型,再进一步计算森林蓄积量[13-14]。
无人机多光谱影像较可见光影像具有更丰富的光谱信息,可计算对森林蓄积量敏感的各类植被指数,已被广泛用于植被参数信息的提取研究,在未来森林蓄积量估测研究中具有较大的潜力[14]。现有研究大多以无人机获取的可见光影像为基础,对无人机载多光谱影像的尝试相对较少。本研究基于大疆精灵4 多光谱版无人机拍摄的多光谱影像,在不进行影像分割的情景下,提取研究区单波段反射率、各类植被指数、纹理特征等因子,计算标准地范围内的均值,建立云南松Pinusyunnanensis林蓄积量估测模型,分析其在森林蓄积量估测建模中的可行性和适用性,旨在为今后森林蓄积量的遥感估算研究提供有益的方法参考。
研究区位于云南省昆明市富民县罗免乡(25°16′21″~25°25′26″N,102°20′46″~102°29′14″E)。该地区属天然云南松林分布的典型区域,地势西南、西北高,东南低。属于低纬度亚热带高原季风气候区,四季温差小,干湿季分明,年平均气温为15.8 ℃,年平均降水量为847.0 mm,研究区内以典型的天然云南松纯林为主。
于2022 年10 月在研究区内选取能充分代表林分总体特征平均水平的地块,设置大小为25 m×25 m 的标准地,共63 个。对标准地内胸径≥5.0 cm 的所有活立木进行每木定位,并测量胸径、树高、最长冠幅、最短冠幅。标准地林分因子见表1。根据实测的胸径、树高,利用二元立木材积公式计算单木材积,在此基础上计算标准地的蓄积量。云南松林二元立木材积公式:
表1 标准地林分因子汇总Table 1 Summary of stand factors in sample plots
其中:V为材积(m3);D为胸径(cm);H为树高(m)。
大疆精灵4 多光谱无人机作为多光谱影像采集平台,集成了1 个可见光相机和5 个多光谱相机[包括红光(B1)、绿光(B2)、蓝光(B3)、红边(B4)和近红外(B5)波段]。于2022 年11 月14 日利用大疆精灵4 多光谱无人机,采用DJI GS PRO 地面站软件从各标准地获取多光谱影像数据。飞行高度设置为100 m,航向和旁向重叠率均设置为85%。利用DJI Terra 软件对原始图像进行预处理,生成数字正射影像(DOM)和数字表面模型(DSM)。将5 个单波段的合成图像在ArcGIS 中合成为多光谱影像(影像分辨率为5.3 cm),计算所需的各类特征变量。
1.4.1 单波段反射率及植被指数 植被指数是指多光谱遥感数据经过线性或非线性数学运算,产生能反映植被生长状况的数值,已广泛用于森林蓄积量建模反演[15]。提取多光谱影像中红光、绿光、蓝光、红边、近红外等5 个波段的反射率(b1、b2、b3、b4、b5),计算蓄积量估测中常用的植被指数:归一化植被指数(NDVI),比值植被指数(RVI),差值植被指数(DVI),大气抗阻植被指数(ARVI),以及根据可见光波段计算的植被指数:过绿指数(EXG)[16],绿蓝比值指数(GBRI)[17],绿红比值指数(GRRI)[18],归一化绿蓝差异指数(NGBDI)[19],归一化绿红差异指数(NGRDI)[20],可见光波段差异植被指数(VDVI)[21]等10 个植被指数。
1.4.2 纹理特征 在蓄积量估测中加入纹理特征有助于提高蓄积量的估算精度[22-24]。为避免影像高频空间信息的丢失,选择较小的3×3 窗口提取纹理特征[25]。借助ENVI 5.3 的纹理提取工具,在3×3 窗口下,通过灰度共生矩阵提取纹理特征,主要包括方差(VA),均值(ME),协同性(HO),熵(EN),对比度(CO),二阶矩(SM),相异性(DI)和相关性(CC)[14]。5 个波段共40 个纹理特征。
1.4.3 均值计算与特征变量筛选 由于标准地为方形标准地,且部分样地郁闭度较低,所以以标准地中的某一点提取各特征变量值不能充分反映标准地特征。本研究以标准地边界为矢量区域,借助ArcGIS 的分区统计工具计算标准地范围内各特征变量的平均值作为自变量,建立蓄积量估测模型[26]。可用于建立蓄积量估测模型的因子会随着研究区、数据源、成像时间等的差异而不同,在建立模型之前对蓄积量与各特征变量进行Pearson 相关性分析,筛选与蓄积量相关性较高的特征变量构建模型。
根据相关性分析的结果,选择在0.01 水平与蓄积量极显著相关的特征变量为自变量,按照7∶3 比例随机划分训练集和测试集,采用多元线性回归(MLR)、支持向量机(SVR)、随机森林(RF)等3 种回归方法建立蓄积量估测模型。
利用决定系数(R2)、均方根误差(ERMS)、平均绝对误差(EMA)、平均相对误差(EMR)进行精度评价[26]。
将单波段反射率、植被指数、纹理特征与蓄积量进行相关分析(表2),5 个波段反射率中,蓄积量与b1、b2无显著相关性,与b3呈极显著负相关(P<0.01),与b4、b5呈极显著正相关(P<0.01)。
表2 蓄积量与单波段反射率、植被指数的相关性Table 2 Correlation between forest volume and single band reflectance and vegetation indexes
在植被指数中,NGRDI 与蓄积量不相关,NGBDI 与蓄积量呈显著正相关(P<0.05)外,其余植被指数均与蓄积量呈极显著正相关(P<0.01)。
如表3 所示:在8 个纹理特征中,HO、EN 与蓄积量的相关性较高,VA、CO、DI 次之,CC、SM、ME 与蓄积量的相关性较低。比较5 个波段提取的纹理特征,与蓄积量显著相关的纹理特征数由大到小依次为B5、B4、B1、B2、B3。
表3 蓄积量与纹理特征的相关性Table 3 Correlation between forest volume and texture factors
根据相关性分析结果,筛选出26 个在0.01 水平与蓄积量显著相关的因子,分别为单波段反射率b3、b4、b5,植被指数RVI、NDVI、DVI、VDVI、EXG、ARVI、GBRI、GRRI,纹理特征B1-EN、B1-HO、B1-ME、B2-EN、B2-HO、B3-EN、B4-CO、B4-DI、B4-HO、B4-VA、B5-CO、B5-DI、B5-EN、B5-HO、B5-VA。
2.2.1 多元线性回归 传统的线性回归模型易受自变量间共线性的影响,利用主成分分析可将原始的多个变量转化为少数主成分因子,每个主成分之间相互独立,克服自变量间的多重共线性的同时保留了原始变量的绝大部分信息[3-4]。根据特征根和累计方差解释率确定主成分个数(表4)。当主成分个数为4 个时,累计方差解释率达93.66%,因此确定主成分为C1、C2、C3、C4。
表4 主成分分析结果Table 4 Principal component analysis
提取的4 个主成分因子与原始变量的线性关系如表5,利用提取出的4 个主成分因子建立多元线性回归模型:M=41.569+11.949C1-2.162C2+2.020C3-2.706C4。其中:M为蓄积量(m3·hm-2);C1、C2、C3、C4为主成分分析提取的主成分因子。
表5 主成分与原始变量的线性关系Table 5 Linear relationship between principal components and original variables
2.2.2 随机森林 本研究使用Matlab 建立随机森林回归模型,决策树数目为100,最小叶子数为1,其余参数均为默认值。如图1 所示:随机森林模型中,变量重要性前5 位的均为纹理特征,说明纹理特征对于蓄积量估测模型的重要性不可忽视。纹理特征间的重要性程度差异大,不同波段间也存在较大差异;植被指数对模型的影响比较稳定,植被指数间除GRRI 外无明显差异。
图1 随机森林模型变量重要性Figure 1 Importance of random forest model variables
2.2.3 支持向量机 基于Matlab 借助LIBSVM 工具箱构建支持向量回归(SVR)模型,支持向量机(SVM)类型为e-SVR,函数选择径向基核函数(RBF),涉及惩罚系数(c)和径向基核函数的参数(γ)这2 个重要参数。利用格网化寻优,得到最佳c、γ分别为4、0.062 5。
精度评价结果:在训练集上,随机森林模型精度最高(R2=0.89,EMA=4.69 m3·hm-2,ERMS=5.45 m3·hm-2,EMR=14.5%),支持向量机次之(R2=0.74,EMA=5.27 m3·hm-2,ERMS=8.31 m3·hm-2,EMR=13.1%),多元线性回归精度最低(R2=0.35,EMA=10.12 m3·hm-2,ERMS=12.85 m3·hm-2,EMR=28.1%)。在测试集上,随机森林精度仍最高(R2=0.69,EMA=7.59 m3·hm-2,ERMS=9.05 m3·hm-2,EMR=26.4%),其次是支持向量机(R2=0.55,EMA=10.51 m3·hm-2,ERMS=12.45 m3·hm-2,EMR=38.5%),多元线性回归最低(R2=0.33,EMA=10.08 m3·hm-2,ERMS=12.21 m3·hm-2,EMR=34.9%)。3 种模型的测试集精度均有降低。
根据估测蓄积量与实测蓄积量绘制散点图。在训练集上(图2A~C),随机森林和支持向量机模型离散程度较小,散点在两侧分布较均匀,模型拟合效果较好;多元线性回归模型散点在两侧分布较均匀,但离散程度较高,模型拟合较差。在测试集上(图2D~F),3 种模型散点在两侧分布不均匀,随机森林离散程度小,其余2 种模型离散程度大。3 种模型在训练集、测试集上均存在一定的低值高估和高值低估现象。
图2 实测蓄积量与估算蓄积量之间的相关关系Figure 2 Correlation between measured forest volume and estimated forest volume
本研究基于无人机多光谱影像提取单波段反射率、各类植被指数、纹理特征,并计算其标准地均值,利用结合主成分分析的多元线性回归、随机森林、支持向量机等建立蓄积量估测模型,结果表明:随机森林模型的精度最高,支持向量机次之,多元线性回归最低。
苏迪等[10]根据平均胸径、平均树高、坡度、坡向等因子建立蓄积量估测模型(模型R2=0.738,ERMS=5.135 8 m3·hm-2),曾霞辉[11]提取平均冠幅、郁闭度、平均树高、株数密度建立蓄积量估测模型(模型对应的R2、EMR分别为0.64~0.78、18.93%~39.68%)。本研究所建立随机森林模型和支持向量机模型的误差水平与其基本一致,多元线性回归模型精度较低。
于东海[13]在对辽宁老秃顶子国家级自然保护区的蓄积量估测研究中,提取植被指数、纹理特征、地形因子,并通过主成分分析提取主成分因子建立多元线性回归模型,模型R2为0.88,拟合效果良好,测试集平均相对误差为9.96%。与该研究相比,本研究结合主成分分析法建立的多元线性回归模型R2为0.35,测试集的平均相对误差为34.9%,精度较低。对研究区进行对比发现,于东海[13]所选研究区森林覆盖率97%,郁闭度、样地单位面积蓄积量高;本研究区郁闭度差异大,标准地单位面积蓄积量低,最大值为79.49 m3·hm-2,最小值低至7.90 m3·hm-2,平均值为41.97 m3·hm-2。2 个研究区森林郁闭度存在较大差异,估测精度可能受到郁闭度的影响。
不同窗口下提取的纹理特征对蓄积量估测精度存在一定影响,提取窗口过大或过小都会造成纹理特征的错误分割进而影响估测精度[27]。基于无人机多光谱影像,在不同窗口大小下提取的纹理特征对蓄积量估测精度的影响有待进一步研究。
研究区林分较稀疏,在部分标准地中,有少量的树木位于标准地边界附近,其部分树冠在标准地边界之外。计算各变量均值时,此部分树冠所对应的特征变量没有参与均值计算,可能对估测精度产生一定的影响。
以大疆精灵4 无人机获取滇中地区典型天然云南松林多光谱影像,提取并计算标准地单波段反射率、植被指数、纹理特征均值,共55 个特征变量,经相关性分析筛选出26 个变量参与模型构建。通过对随机森林、支持向量机、多元线性回归3 种模型进行精度评价,最佳估测模型为随机森林。3 种模型在估测蓄积量时均存在一定的低值高估和高值低估现象。与同类型研究的蓄积量估测误差水平基本一致,表明在不进行单木分割的情况下,利用无人机多光谱影像提取各因子,以各因子的标准地均值建立模型估测蓄积量具有一定的可行性。