基于多时相无人机影像的高郁闭度森林采伐生物量估算

周小成 王 佩 谭芳林 陈崇成 黄洪宇 林 宇

(1.福州大学空间数据挖掘与信息共享教育部重点实验室, 福州 350108;2.福建省林业科学研究院, 福州 350012; 3.福建省闽侯白沙国有林场, 福州 350102)

0 引言

森林作为陆地生态系统中最重要的组成部分之一,在促进生物多样性保护以及全球碳循环过程中发挥着至关重要的作用。森林地上生物量(Above ground biomass,AGB)是反映森林碳储量和森林碳汇能力的重要因子,准确估算森林地上生物量的动态变化有助于为政府及有关部门进行碳排放调控和碳汇交易提供科学决策依据[1]。传统的森林地上生物量直接调查方法人力成本高、劳动强度大、有效覆盖面积较小,还会对森林造成一定的损害[2]。当前采用遥感技术估算森林地上生物量已经成为重要的研究方向。按照数据源与研究方法的不同,大体可以分成以下两类:数据源上主要包括光学遥感影像、SAR、LiDAR点云及多源遥感数据的结合[3-6];研究方法上主要包括参数回归(例如多元回归、异速生长模型等)和非参数机器学习模型(例如随机森林(Random forest,RF)、支持向量机(Support vector machine,SVM)、最邻近回归(K-nearest neighbor,KNN)等)[7-8]。从数据源获取途径来说,航天卫星遥感影像分辨率一般较低,精度有限;机载LiDAR在获取林分参数时精度较高,但存在成本高、数据处理复杂的限制;而可见光无人机遥感具有飞行时效强、成本低和可获取厘米级空间分辨率影像等优势,现已被广泛应用在林业资源调查和精准农业[9]等多个领域。

近年来,许多学者利用可见光无人机影像,使用了不同的方法来对地上生物量进行估算。但由于估算难度大,缺少关于高郁闭度森林采伐生物量的估算研究[10]。由于高郁闭度的森林内部及冠层结构的复杂性,仅使用单一时相的可见光无人机影像,难以穿透林冠,获取准确的地面信息,使得树高信息提取难以达到精度要求[11],进而使得估算高郁闭度地上生物量比较困难。因此,当前的地上生物量研究区域主要集中在中低郁闭度样地。如LIN等[3]在稀疏的亚高山针叶林地样地,对可见光无人机影像采用冠层高度模型(Canopy height model,CHM)并叠加多尺度分割的树冠数据获取树高,建立生物量-树高的异速生长模型,根据实测单株树木生物量验证,结果表明单株生物量估算的相对均方根误差(RMSE)为54.9 kg。何游云等[12]以岷江冷杉为研究对象,因分布间隙较大,对可见光无人机影像采用人机交互方式提取树冠面积,建立生物量-树冠面积模型估算生物量,与实测的胸径估算生物量进行对比,结果表明偏差为5.15%。LATIFI[13]在稀疏的树丛中,采用消费级可见光无人机影像提取树高和树冠面积估算胸径,采用异速生长模型估算地上生物量,结果表明对地上生物量估算结果较好。

而对于森林采伐,结合采伐后无人机影像则可得到较为准确的数字高程模型(DEM)和伐桩信息。综上所述,本文根据多时相可见光无人机影像对高郁闭度森林采伐生物量进行估算。首先是采用多时相可见光无人机影像数据,生成冠层高度模型(CHM),运用动态窗口局部最大值法提取采伐株数和单木树高信息,再结合采伐后高分辨可见光影像,运用YOLO v5(You only look once versions 5)方法检测并提取伐桩直径信息,建立胸径-伐桩直径模型,利用二元生物量公式估算森林采伐生物量,为实现高郁闭度森林采伐生物量估算提供思路与方法。

1 研究区与数据源

1.1 研究区概况

研究区位于福建省闽侯白沙国有林场一个松林采伐小班(26°15′13″N,119°04′48″E),地处闽侯县西北部、闽江之畔,多为中、低山地貌,海拔200 m以上,试验区坡度40°,小班测区面积27 533 m2,其中采伐带面积6 333 m2,主要树种为马尾松,占比达到90%以上,间杂零星杉木,小班郁闭度0.8,林龄33 a,采伐时间为2021年4月7日,实际调查采伐株数192株,采伐类型为条带带状采伐,根据多时相无人机影像(图1)可看到森林采伐前的树冠与森林采伐后的地形地貌。

图1 研究区位置Fig.1 Location of test region

1.2 数据获取

本研究使用的多时相可见光影像数据采集设备为大疆精灵4 RTK(Real-time kinematic)型无人机(Unmanned aerial vehicle,UAV),该款无人机具备厘米级导航定位系统及2 000万像素CMOS传感器,便携易用。本试验区的森林采伐前无人机影像获取时间为2021年3月27日,航线设计为航向重叠率80%,旁向重叠率70%,相对飞行高度300 m,影像正射空间分辨率9.7 cm。对于森林采伐后无人机影像数据的采集策略是:考虑到该测区地形起伏较大且该型无人机具备仿地形飞行的能力,为保证地面分辨率不至于变化过大导致拼接效果差的情况,因此对测区采取仿地飞行,拍摄时间为2021年5月6日,最终航线设计为航向重叠率80%,旁向重叠率70%,仿地飞行高度37 m,影像正射空间分辨率1.8 cm。

采用胸径尺及图帕斯激光测高仪在白沙林场现场收集146株马尾松伐桩直径(离地5 cm)、胸径、冠幅和树高信息,其中66株有树高、冠幅对应信息。为提高回归模型的泛化能力,同时获得福建省将乐县国有林场收集的179株马尾松幼林地的胸径、树高信息和102株树高、冠幅信息,现场实测数据统计如表1所示。

表1 实测数据统计Tab.1 Statistical results of measured data

2 研究方法

根据多时相可见光无人机影像对高郁闭度森林伐区进行采伐生物量估算,整体思路为:基于采伐前后可见光影像DSM特征差值生成CHM,运用动态窗口局部最大值法提取采伐株数与单木树高参数,再进一步运用YOLO v5算法从采伐后影像中提取伐桩直径,以实测数据建立胸径-伐桩直径模型,代入二元生物量公式估算森林采伐生物量。为开展多时相与单时相估算精度的有效对比,基于采伐前单一时相可见光影像分割单木冠幅,以实测数据建立树高-冠幅面积模型,结合胸径-树高模型,代入二元生物量公式估算森林采伐生物量,将本研究设计的一种多时相估算生物量新方法与采伐前单期影像估算结果进行对比分析。总体技术路线如图2所示。

图2 技术路线图Fig.2 Technical flow chart

2.1 无人机数据预处理

无人机数据使用Pix4D软件进行全自动空中三角测量计算,得到研究区的数字正射影像图(Digital orthophoto map,DOM)、数字表面模型(Digital surface model,DSM)和UAV影像点云产品。获取的可见光正射影像色彩自然、协调、无明显拼接痕迹。

由于树高提取需要高精度的DSM数据,因此须对无人机获取的多时相DSM数据进行精度验证。本文数据采集设备为大疆精灵4 RTK型无人机,该型无人机自带RTK系统,可将航测图像的POS信息精确到厘米级,无需设置像控点。随机从伐区中选取14个采伐前后两期DSM的同名地面点,用于进行高程精度验证。

2.2 单木冠幅分割

由于高郁闭度(小班郁闭度为0.8)森林冠层结构复杂,仅使用单一时相的可见光无人机影像,难以穿透林冠,获取准确的地面信息,使得树高信息提取难以达到精度要求。而可见光无人机影像能有效获取地物顶部的纹理特征,本研究则转换思路,首先采用分水岭分割算法对单木冠幅进行提取,再利用实测数据建立树高-冠幅面积回归模型估算树高信息。分水岭分割算法[14]原理是将待分割的影像视为测地学的拓扑地貌,采用影像中的每一点像元值表示该点的海拔,那么像元的变化即可描述成地形图的山脊、谷底和集水盆,每个局域极小值及其影响的区域构成集水盆,而集水盆的边界范围则形成分水岭。

2.3 采伐株数与树高估算

根据采伐带的矢量边界,首先对多时相影像的DSM数据进行掩膜分析,再通过采伐前后DSM差值运算生成树冠层高度模型(CHM),此时获得的林冠高度最为准确,可直接用来获得高精度的单木树高。

国内外检测树冠顶点的研究方法中局部最大值法应用的最多[15],该方法窗口类型一般分为固定窗口与动态窗口。由于利用固定窗口的局部最大值法在高郁闭度森林中检测林冠顶点时经常会出现漏提或多提现象[16],因此,选择合适的策略来检测树冠顶点十分重要。

本研究通过IDL语言编写动态窗口局部最大值法来检测树冠顶点位置。首先是对生成的CHM采用一个较小的固定窗口检测林区潜在的树冠顶点,然后运用自适应动态窗口对提取的顶点位置进行判别,假如当前的顶点为对应窗口范围中的最大值点则保存,否则剔除。动态窗口尺寸根据计算潜在顶点的8个不同剖面方向的经验半方差变化值确定[17],其计算公式为

(1)

式中γ(h)——经验半方差

xi——图像中像元点位置

h——两个像元间的分割距离

Z(xi)——对应影像xi像元值

N(h)——在分隔距离为h时像元对数量

由于研究区的主要树种为马尾松,属于针叶林,而针叶林的最高点一般为树冠中心点位置,若动态窗口的中心为最大值,则判断为树冠顶点,而其对应像元属性值则为树高。

2.4 伐桩目标检测与伐径提取

对于本研究的伐区场景下的伐桩检测与其他目标检测相比具有一定的特殊性,伐桩目标尺寸一般较小,背景信息量大(有类圆形石块、造林坑穴等影响),且许多伐桩被残枝和杂草遮挡,为伐桩检测增加了难度。由于YOLO v5具有目标检测速度快、背景误检率较低且模型泛化能力较强等优点[18],因此本研究拟采用YOLO v5方法对伐桩目标进行检测。

YOLO v5是由UItralytics团队于2020年发布的一种能直接得到预测目标的精准边界框和类别的端到端网络模型。它的基本框架[19]主要分为输入端(Input)、主干网络(Backbone)、颈部模块(Neck)和输出预测端(Prediction),其结构原理如图3所示。

图3 YOLO v5网络结构图Fig.3 YOLO v5 network structure diagram

2.4.1伐桩目标检测

本研究的样本集来自于福建省闽候白沙国有林场的可见光无人机原始图像。本研究对获取的伐区图像按一定重叠的滑动裁切,筛选出208幅高纯度的伐桩图像。为了增加样本的多样性和提高模型的泛化能力,对伐桩图像采用随机旋转、随机改变亮度、镜像变换的离线数据增广方法,增加到440幅伐桩图像,将训练集和测试集按照9∶1进行随机抽取划分,使用Labelme工具对训练集进行标注,完成样本数据集的划分。将正射影像按一定规则格网进行滑动裁切,然后利用YOLO v5方法对上寨工区的伐桩进行检测,再根据空间坐标映射方法进行镶嵌处理。

为了定量评估YOLO v5方法的性能以及对伐桩检测的有效性,本文将目视解译的伐桩数量作为真实值,引入精确率P(Precision)、召回率R(Recall)和总体精度F1值作为评价指标。

2.4.2伐桩直径提取

由于YOLO v5方法是一种能直接得到预测目标精准边界框的端到端网络模型,将检测到的伐桩目标矩形框栅格数据转为矢量数据,在ArcGIS可自动计算得到检测目标矩形框的面积及周长。由于伐桩直径实地测量时是量取东西、南北两方向长度并计算平均值获取。伐桩一般属于类椭圆形目标,根据矩形内切椭圆的原理(图4),假设以矩形的中心为椭圆的中心,矩形的长为a、宽为b。因此,椭圆的长轴长为a、短轴长为b,可得出矩形周长的1/4即为伐桩直径。

图4 矩形内切椭圆原理图Fig.4 Schematic of rectangular inscribed ellipse

2.5 树高和胸径回归模型构建

2.5.1树高-冠幅面积回归模型

对于采伐前单期无人机影像树高的估算,本研究通过在福建省将乐县林场及福州市白沙林场野外实测的168株马尾松的单木树高以及冠幅为基础数据,使用SPSS分析常用的5种模型:一元线性、对数函数、幂函数、S曲线、增长函数,建立适用于福建省闽侯县白沙林场的树高-冠幅面积回归模型。随机选取134株成活立木的树高(H)、冠幅面积(W)为试验数据,建立5种树高-冠幅面积回归模型,以剩余的34株树高、冠幅面积作为验证数据。

2.5.2胸径-树高回归模型

通过在福建省将乐县林场及福州市白沙林场野外实测的325株马尾松的单木树高(H)以及胸径(D)为基础数据,根据已有的4个树高-胸径经验模型(线性回归、对数模型、异速生长模型、Bates and Watts模型)[20],反推出胸径-树高模型,以选取最优的马尾松胸径-树高回归模型。

随机选取260株马尾松的树高、胸径作为试验数据。运用Origin 2018编译出4种候选胸径-树高模型公式,可得到胸径-树高回归模型结果,以剩余的65株胸径、树高作为验证数据。

2.5.3胸径-伐桩直径回归模型

由于马尾松树种的胸径-伐桩直径经验模型较少,本研究采用在白沙国有林场收集的146株马尾松的伐桩直径(离地5 cm)、胸径的野外实测值为基础数据,使用SPSS分析常用的5种模型(一元线性、对数函数、幂函数、S曲线、增长函数),建立适用于福建省闽侯县白沙国有林场的胸径-伐桩直径回归模型。随机选取116株马尾松胸径(D)、伐桩直径(D0)试验数据,建立5种胸径-伐桩直径回归模型,以剩余的30株胸径、伐桩直径作为验证数据。

2.6 生物量估算

计算生物量时,先计算单株立木的生物量,再计算采伐区域的总生物量。胸径获取的2种途径:①利用采伐前无人机影像所得的单木冠幅面积,根据树高-冠幅面积模型估算出树高,再利用胸径-树高模型估算胸径。代入国家林业局2014年颁布的中华人民共和国林业标准《立木生物量模型及碳计量参数—马尾松》的二元立木生物量公式[21]中,计算单株树木生物量和总生物量。②结合采伐后YOLO v5方法提取的伐桩直径进行回归估算,再将多时相影像提取的采伐株数和平均树高代入二元立木生物量公式中,总生物量计算式为

(2)

式中Dj——单木胸径

Hj——单木树高

Mj——单木二元生物量

2.7 精度评价

将决定系数(R2)和均方根误差(RMSE)作为回归模型精度评价指标,确定最优模型。

为验证生物量估算精度,以野外实地测量得到的平均树高、平均胸径计算获得的生物量为参考,对无人机遥感采伐生物量估算精度进行验证。

3 结果与分析

3.1 无人机数据预处理验证

通过对伐区的采伐前后两期DSM的同名地面点高程精度验证结果表明,采伐前后两期DSM中的同名点高程表现出极强的相关性,两者的R2达到0.999 6,其中RMSE为0.066 78 m(图5),根据林业生产要求,该精度能满足伐区树高估算的要求。

图5 采伐前后DSM精度验证Fig.5 Verification of DSM accuracy before and after deforestation

3.2 单木冠幅分割

根据分水岭分割算法对高郁闭度的伐区进行单木冠幅分割,其中分割结果如图6所示。由于高郁闭度森林的林木冠幅存在相连、遮挡等情况,无法获取准确的单木冠幅与面积,得到分割株数结果为162株,野外采伐株数为192株,株数估算精度为84%,林木冠幅平均面积为23.91 m2,野外实测该地区的林木冠幅平均面积为16.15 m2,估算精度为67.55%。

图6 分水岭分割冠幅结果Fig.6 Crown segmentation results based on watershed segmentation

3.3 采伐株数与树高估算

采用多时相DSM差值特征(CHM)的动态窗口局部最大值法估算出的株数为185株,实际采伐株数为192株,株数估算精度为96.35%,提取平均树高为19.07 m,野外实测树高为18.88 m,估算精度为99.01%,图7为冠层高度模型提取的单木树高分级显示。

图7 冠层高度模型的树高分级显示图Fig.7 Tree height classification diagram of canopy height model

3.4 伐桩目标检测与伐径提取

采用YOLO v5方法对上寨工区的伐桩进行检测,伐区内伐桩检测的结果及局部放大区域如图8所示。

图8 YOLO v5检测结果Fig.8 YOLO v5 detection results

由于伐桩容易受到影像空间分辨率、杂草和泥块遮挡等因素影响,仅在伐区目视解译到54个伐桩。YOLO v5方法在上寨工区从采伐带中正确检测出47个,正确检测数量较高,精确率为69.12%,召回率为87.04%,伐桩总体F1值为77.05%。YOLO v5方法提取出47个伐桩直径平均为31.41 cm,结合人工目视修正伐桩直径结果为29.61 cm。

3.5 树高和胸径估算回归模型结果

3.5.1树高-冠幅面积回归模型

使用SPSS分析常用的5类模型结果(表2)可知,马尾松的树高和冠幅面积表现出较强相关性。根据R2最大、RMSE最小的标准,本研究则选取幂函数为马尾松的树高-冠幅面积回归模型。

表2 选用的5种树高-冠幅面积回归模型结果Tab.2 Results of five selected tree height-crown area regression models

选定马尾松的最优树高-冠幅面积回归模型后,将剩余的34株马尾松野外实测的单木树高作为模型的验证数据,验证模型精度。由图9可知,马尾松的实测树高和估算树高的决定系数为0.760 1,表明建立的树高-冠幅面积模型可用于白沙国有林场。

图9 树高-冠幅面积建模和验证Fig.9 Tree height-crown area modeling and verification

3.5.2胸径-树高回归模型

使用Origin 2018编译的4个胸径-树高模型结果(表3)表明,马尾松的胸径和树高表现出较强的相关性,其中R2均在0.6以上。根据R2最大、RMSE最小的标准,本研究选取异速生长模型,作为马尾松的胸径-树高回归模型。

表3 马尾松的胸径-树高模型Tab.3 DBH-tree height model of Masson pine

选定马尾松的最优胸径-树高模型后,将剩余的65株马尾松野外实测胸径作为模型的验证数据,验证模型精度。由图10可知,马尾松的实测胸径和估算胸径的决定系数为0.628 1,在一定程度上表明建立的马尾松胸径-树高模型可用于白沙国有林场。

图10 胸径-树高建模和验证Fig.10 DBH-tree height modeling and verification

3.5.3胸径-伐桩直径回归模型

使用SPSS分析常用的5类模型结果(表4)可以看出,马尾松的胸径和伐桩直径表现出极强相关性,其中R2均在0.9以上。根据R2最大、RMSE最小的标准,本研究则选取幂函数模型,为马尾松的胸径-伐桩直径回归模型。

表4 选用的5种胸径-伐桩直径回归模型结果Tab.4 Results of five selected DBH-pile diameter regression models

选定马尾松的最优胸径-伐桩直径模型后,将剩余的30株马尾松野外实测伐桩直径作为模型的验证数据,验证模型精度。由图11可知,马尾松的实测胸径和估算胸径的决定系数为0.949 6,表明建立的胸径-伐桩直径模型可用于白沙国有林场。

图11 胸径-伐桩直径建模和验证Fig.11 DBH-stump diameter modeling and verification

3.6 两种生物量估算方法结果分析

表5为前面提到的几种方法获得的采伐生物量结果。其中,仅利用采伐前单一可见光无人机影像估算的高郁闭度森林采伐生物量精度仅为53.13%,这是由于高郁闭度森林的林木冠幅存在相连、遮挡等情况,无法获取准确的单木冠幅面积,导致采伐生物量提取精度较低。根据多时相影像获取的平均树高与采伐株数,再结合YOLO v5方法检测伐桩并提取直径信息,估算的采伐生物量精度为83.08%,表明结合伐桩信息来估算森林采伐生物量这一思路与方法,具有较强的可行性。

表5 生物量精度验证Tab.5 Validation of biomass accuracy

4 讨论

胸径和树高是估算森林地上生物量的重要因子,传统的森林采伐生物量直接调查方法需要对树木砍伐后进行干燥称量,这样费时费力,且不适用于大规模森林的采伐生物量估算。由于仅利用单一时相可见光无人机影像对高郁闭度森林采伐生物量估算时,冠幅存在相连、遮挡等情况,提取单木冠幅困难,后续可尝试将影像转为灰度梯度图像,对树冠边缘特征进行增强后再进行分割[22]。本研究则提出根据多时相无人机影像估算森林采伐生物量的途径,该方法估算精度结果较好,且与传统卫星影像和机载LiDAR数据相比具有飞行时效强、作业高效和成本低等优势[23]。

本研究中的伐桩检测与其他目标检测相比具有一定的特殊性,伐桩目标尺寸较小,伐区背景信息量大,且许多伐桩被杂草和泥块遮挡,为伐桩准确检测增加了一定难度。目前采用无人机遥感影像进行伐桩检测并推算森林采伐生物量的研究非常有限。现有的少量文献对伐桩进行检测,主要可分为2方面:①利用传统的计算机视觉及机器学习方法对伐桩检测,如SATHISHKUMAR等[24]使用可见光无人机影像通过模板匹配、霍夫变换和相位编码方法对火炬松伐桩的检测精度分别为41%、77.3%和72%,其中霍夫变换、相位编码方法提取伐桩直径RMSE分别为7.2 cm和4.3 cm。STEFANO等[25]采用随机森林分类对空间分辨率为1.7 cm的无人机影像进行伐桩检测,根据11个57 m×57 m的地块研究表明,伐桩总体检测精度为68%～80%,伐桩直径的RMSE为6.4～7.5 cm。KAMARULZAMAN[26]使用分辨率为0.06 cm的无人机影像采用支持向量机分类对伐桩进行检测,根据4个地块研究表明,检测伐桩总体精度为66.67%～71.43%。这些研究主要是选取无人机影像中伐桩的纹理特征参与检测分类,分类过程繁琐。②利用深度学习方法进行的伐桩检测研究,如WINDRIM等[27]使用分辨率0.2～0.4 cm的无人机影像通过深度学习卷积神经网络(Convolutional neural networks,CNN)方法,对9个10 m×10 m地块的松树伐桩检测的精度为66.6%～97.7%。由于经典卷积神经网络,缺乏自适应瞄框,需要进行大量调参,会导致目标检测边界框难以准确描绘。

与文献[24-27]的方法相比,本文采用的YOLO v5方法具有自适应瞄框计算,是一种能直接得到预测目标的精准边界框的端到端模型。本研究通过森林采伐后高分辨率无人机影像,根据YOLO v5方法检测到伐桩数并提取伐桩直径信息,通过建立胸径-伐桩直径模型,估算的平均胸径为 26.38 cm,相较于实测的平均胸径23.78 cm,偏高 2.6 cm。森林采伐生物量结果偏高的主要原因可能是由于本研究的伐桩样本数量较少,从而对胸径估算精度产生较大的波动影响,导致估算的采伐生物量偏高。后续的研究将扩大研究区域范围,增加伐桩样本数量,以提升森林采伐生物量估算精度。

由于时间和试验条件有限,本文仅对主要树种为马尾松的样地进行试验,后期的研究将采用混交林、阔叶林等不同类型的样地进行试验,并建立优化的胸径、树高和冠幅之间的模型,以进一步提升该技术方法估算精度。

5 结论

(1)根据采伐前后多时相DSM数据,获取采伐带冠层高度模型(CHM),运用动态窗口局部最大值法检测林冠顶点,得到估算的株数为185株,采伐株数为192株,精度达到96.35%,估算的平均树高为19.07 m,野外实测的平均树高为18.88 m,精度达到99.01%。

(2)根据森林采伐后的高分辨无人机可见光影像,运用YOLO v5方法检测伐桩的总体检测精度为77.05%,通过提取伐桩直径,建立胸径-伐桩直径回归模型,估算的平均胸径为26.38 cm,野外实测胸径为23.78 cm,精度为90.14%。

(3)根据采伐前后多时相影像获取的采伐株数和平均树高信息,结合YOLO v5方法从采伐后的可见光无人机影像中检测伐桩直径,建立胸径-伐桩回归模型,再利用二元立木生物量模型,估算出采伐样地森林采伐生物量为50 458.09 kg,野外实测数据生物量41 922.76 kg,精度为83.08%,估算结果精度较好。结果表明,采伐前后多时相无人机可见光影像获取采伐株数与树高,再结合伐桩直径信息,估算森林采伐生物量的思路,具备一定有效性,需进一步改进伐桩直径的估算精度,从而提高生物量估算精度。