基于无人机影像的无瓣海桑单木提取与地上生物量估算

余楚滢，龚 辉，曹晶晶，刘燕君，刘 凯

[ 1. 中山大学 地理科学与规划学院//广东省公共安全与灾害工程技术研究中心//广东省城市化与地理环境空间模拟重点实验室，广州510006；2. 南方海洋科学与工程广东省实验室（珠海），广东 珠海 519000；3. 广州市城市规划勘测设计研究院，广州 510060 ]

红树林是热带亚热带地区海岸潮间带的木本植物群落，具有显著的防风消浪、促淤造陆和净化海洋等效果（陈玉军 等，2012；卢昌义 等，2019），其在保护海岸带生态安全中发挥着重要作用。红树林拥有强大的固碳能力，准确的红树林生物量估算是定量研究其生态系统生产力和生态系统碳循环的重要内容（田义超 等，2019；王子予 等，2020）。无瓣海桑（Sonneratia apetala）是海桑科海桑属乔木，其具有良好的环境适应性，且生物量积累快、固碳速率高，生产力水平处于中国红树林群落中的较高位置，是人工红树林营造应用中重要的优质红树树种之一（朱可峰 等，2011；周元满 等，2012）。相比于传统的红树林生物量估算方式，遥感技术已成为红树林生物量估算研究的重要手段和数据源。然而，基于卫星遥感影像估算精度通常受限于较低的影像分辨率，且大多数的红树林生物量估算研究对所有树种进行整体建模，较少考虑不同树种的生物量差异。由于红树林群落冠层密集、结构复杂，目前，树种尺度的红树林生物量估算研究仍较少。

单木尺度的森林结构信息有助于提高红树林生物量估算精度。基于遥感技术的单木提取包括单株立木探测和单木树冠提取两部分（董新宇 等，2019）。单株立木探测指对单木树冠的所在位置进行探测，找出每个树冠的中心点；单木树冠提取指将探测到的树冠中心点作为参照，找到树木树冠的边界点，对树冠轮廓进行描绘（刘晓双 等，2010）。目前，国内外单木提取研究主要集中在针叶林（Dalponte et al., 2014; Yang et al., 2016），这是因为针叶林的树冠呈尖塔形，树冠高点明显，而红树林由于树冠比较密集，相邻树木之间高度差较小，在遥感数据中难以识别树冠边界，红树林的单木提取精度受限。近年来，广受关注的低空无人机（Unmanned aerial vehicle, UAV）遥感技术为从高分辨率影像中提取单木树冠信息提供新的机会，且已在红树林单木提取中得到应用。如Yin 等（2019）基于无人机载激光雷达（Light Detection And Ranging, LiDAR）数据实现了红树林的树冠分割和参数估计，并分析了树冠聚集度和点云密度对单木提取精度的影响；Navarro 等（2020）使用无人机影像生成高分辨率的冠层高度模型（Canopy Height Model, CHM），基于可变窗口滤波（Variable Window Filter, VWF）和控制分水岭分割方法（Controlled Watershed Segmentation Method, CWSM）实现红树林单木提取。

以往研究主要利用地面样方数据，采用皆伐法、平均木法和异速生长法等方法（薛立 等，2004）构建红树林地上生物量模型。其中，异速生长方程是一种通过部分单木参数得到单木生物量的方法，相比于皆伐法、平均木法，该法较为简便，破坏性小，便于推广（Castedo-Dorado et al.,2012； 金川 等，2012）。基于遥感技术的红树林地上生物量估算方法可分为两类：一是提取光谱、指数、纹理等遥感特征，建立生物量反演模型（Lei et al., 2022）；二是通过高分辨率影像、LiDAR数据等获取树高、冠幅等林分因子，基于异速生长方程构建生物量估算模型（Pham et al., 2019）。前者结合多个特征变量，可达到较高的精度，但反演模型可推广性较差。后者基于异速生长法理论，具有较强解释能力。但卫星遥感影像通常受限于较低的影像分辨率（Zhu et al., 2017），难以实现单木尺度的红树林生物量估算。近年来，无人机遥感在精细尺度的植被生物量监测中备受关注，部分学者也开始探讨其在红树林生物量估算中的应用潜力。如Qiu 等（2019）结合WorldView-2 高分辨率影像和无人机载LiDAR 数据进行单木和格网尺度的红树林地上生物量估算；Navarro 等（2020）利用无人机影像进行红树林单木检测，并结合异速生长方程估算生物量。然而，无人机遥感在红树林生物量监测中的应用潜力仍有待进一步探究，且当前基于无人机遥感的单木尺度无瓣海桑生物量估算研究仍较少。

因此，本文以广东省珠海市淇澳岛红树林自然保护区为研究区，基于无人机影像和种子区域生长（Seed Region Growing, SRG）算法进行无瓣海桑单木提取，通过确定研究区无瓣海桑树高和胸径之间的回归关系，构建基于树高的无瓣海桑异速生长方程以准确估算其地上生物量。以期为红树林生态恢复、保护与管理提供决策支持。

1 研究区与数据

1.1 研究区介绍

广东省珠海市淇澳岛（22°23′40″-22°27′38″N，113°36′40″-113°23′40″ E）位于珠江口内伶仃洋西侧，东望深圳、香港，北邻中山、东莞，海岛总面积约24 km2（王树功 等，2005；Cao et al.,2021）。淇澳岛属热带海洋性季风气候，雨量充沛、光照充足、雨热同期（黄建荣 等，2011）。淇澳岛红树林保护区内有秋茄（Kandelia candel）、桐花树（Aegiceras corniculatum）、银叶树（Heritiera littoralis Dryand）、老鼠簕（Acanthus ilicifolius）和无瓣海桑（Sonneratia apetala）等主要红树植物。本研究区位于淇澳岛红树林保护区的中、高潮区域（图1），面积为2.89 hm2。研究区的植被群落中上层木主要为无瓣海桑，中层木主要为桐花树和老鼠簕，其中无瓣海桑为研究区的优势种。研究区的无瓣海桑为2001-2002年左右种植（Yu et al., 2020），2019年采集本文使用的无人机影像时无瓣海桑的树龄约为17～18 a。

图 1 广东珠海淇澳岛地理位置（a）与研究区域（b、c）及野外调查样方点位置（c）Fig.1 Geographical location of Qi'ao Island, Zhuhai, Guangdong (a);Location of the study area (b, c); Quadrat positions in the study area (c)

1.2 无人机影像获取与预处理

采用大疆DJI精灵4 PRO V2.0四旋翼消费级可见光无人机设备进行研究区数据采集。该无人机质量轻、体积小，搭载CMOS 传感器，镜头焦距为8.8 mm，无人机平台与RGB 传感器集成一体，整机重量1 375 g。共有红、绿、蓝3个波段，传感器分辨率为5 472×3 648（大疆创新，2021）。无人机影像的采集时间为2019-11-10，天气状况良好，飞行时段选择正午T 10:00-14:00，风速约为1～4 m/s。使用Altizure生成无人机倾斜摄影航线，航向重叠率设置为80%，旁向重叠率设置为75%，航高设置为100 m，飞行速度约为9 m/s，倾斜航线的俯仰角度为45°。共获取覆盖研究区的无人机倾斜摄影影像566幅。

基于倾斜摄影建模软件ContextCapture 生成研究区无人机正射影像图（Digital Orthophoto Map, DOM） 和数字表面模型（Digital Surface Model, DSM）以及密度点云（图2）。其中DOM空间分辨率约为0.04 m，DSM 空间分辨率为0.25 m，密度点云面积为0.343 km2，空间分辨率为0.03 m。

图2 基于无人机航片生成的研究区正射影像图DOM（a）、数字表面模型DSM（b）和密度点云（c）Fig.2 Digital Orthographic Model DOM (a), Digital Surface Model DSM (b) and point cloud (c) of the study area based on UAV Imagery

1.3 地面数据采集与预处理

1.3.1 实测单株无瓣海桑树高和胸径数据 使用Trimble SX10三维激光扫描仪、测高仪和卷尺对无瓣海桑进行调查，获取无瓣海桑的树高，并测量树干1.3 m处的直径（胸径），以构建无瓣海桑树高与胸径的回归关系。由于红树林生长在潮间带且生境复杂，研究区内难以架设激光扫描仪设备，因而，采样地点主要选取研究区边缘的无瓣海桑群落。采集时间为2018-06-28，共采集了48棵无瓣海桑的树高和胸径数据，用于构建两者的回归关系。

1.3.2 样方数据 针对研究区开展无瓣海桑地面样方调查，数据采集时间为2020-11-03，在研究区内构建5 个10 m×10 m 的样方（见图1）。借助测高仪和近距离目视估测无瓣海桑的树高；使用软尺测量树干1.3 m处的直径（胸径）。获取的树高和胸径数据主要用于代入异速生长方程以获取样方地上生物量，验证无人机获取的结果。

1.4 技术路线

基于无人机倾斜摄影影像的无瓣海桑单木提取与地上生物量估算技术路线如图3所示。首先，基于无人机影像生成冠层高度模型；然后，采用局部极大值算法和种子区域生长算法进行无瓣海桑单木提取；同时，利用地面调查数据确定树高和胸径之间的关系，构建基于树高的无瓣海桑异速生长方程；最后，估算研究区的无瓣海桑地上生物量，并结合地面调查数据进行验证分析。

图3 无瓣海桑单木树冠提取与地上生物量估算流程Fig.3 Flow chart of tree crown delineation and aboveground biomass estimation of Sonneratia apetala

2 基于无人机点云的无瓣海桑单木树冠提取

2.1 点云滤波与冠层高度模型生成

基于无人机航片生成的研究区原始点云数据可分为植被点云和地面点云两部分。为了精确获取植被的冠层高度模型（Canopy Height Model, CHM），通常需先将地面点与非地面点分离开。已有研究证明，布料模拟滤波（Cloth Simulation Filtering,CSF）算法效果较好且所需参数较少（Zhang et al.,2016； 张昌赛 等，2018）。在CloudCompare 中采用CSF算法进行点云滤波，进而分离地面点与植被点。最终得到的点云滤波结果如图4-a 所示，可以看出，地面点云主要包括林分外围滩涂点云和林窗间隙内的地面点云。对CSF滤波处理得到的地面点云数据进行处理，采用最邻近算法进行插值，并在ArcGIS中使用低通滤波进行平滑处理，生成数字高程模型（Digital Elevation Model, DEM），将DSM和DEM 进行差值运算得到研究区冠层高度模型（Canopy Height Model, CHM），重采样后空间分辨率为0.25 m（图4-b）。运用CHM 模型描述研究区林木从地表至冠层的高度分布。

图4 研究区地面点云（a）和CHM（b）Fig.4 Ground point cloud (a) and CHM (b) of the study area

2.2 单木位置探测与树冠分割

单木提取通常包括2个步骤：首先，进行单木位置探测，即将搜索到的树冠顶点作为单木的空间位置；再以提取的树冠顶点为标记点，进行单木的树冠分割。采用经典的局部极大值算法（Yin et al.,2019）进行单木顶点提取。该算法将树冠最高点作为单木位置，通过移动窗口逐步对栅格数据进行搜索，并判断搜索窗口的中心点是否为局部极大值，若为局部极大值，则将此像素标记为单木顶点。单木探测结果的准确性通常与搜索窗口大小、影像分辨率和冠幅大小等因素有关。采用可变窗口滤波（Variable Window Filter, VWF） 算法（Popescu et al., 2004）确定搜索窗口大小。VWF算法基于树高与冠幅的相关关系假设，即树木越高则冠幅越大。该算法根据窗口中心的像素值确定移动窗口的大小，其中可变窗口函数的选择非常关键。根据前期地面调查数据和多次试验，进而确定函数y=0.06x为窗口变换函数（其中y为窗口大小，单位为m；x为树高，单位为m）。以上步骤通过R 语言Forest Tools工具包实现，最终得到研究区红树林单木顶点的局部探测结果（图5）。

图5 无瓣海桑单木顶点探测结果局部（a. 局部位置示意图；b. 局部区域正射影像；c. 局部区域单木顶点检测结果）Fig.5 Partial result of individual tree detection of Sonneratia apetala (a. Location of the partial area; b. DOM of the partial area;c. Individual tree detection result of the partial area)

结合研究区CHM 影像，在SAGA GIS 中采用种子区域生长（Seed Region Growing,SRG）算法（Yin et al., 2019）进行单木树冠分割。SRG 算法将图像种子点标记为生长区域，进一步判断相邻像素与种子点的相似度，若相似度在阈值以内则并入生长区域。不断循环上述操作直到满足停止条件或遍历结束。该算法利用特征阈值和距离阈值控制分割对象的生长。使用树冠顶点作为种子点进行区域生长，以进行树冠分割，得到的研究区域无瓣海桑群落分割结果及局部分割效果（图6）。

图6 无瓣海桑单木分割结果（a. 研究区总体单木分割结果与局部区域位置示意图；b. 局部区域正射影像；c. 局部区域单木分割结果）Fig.6 The individual tree crown delineation result of Sonneratia apetala (a. Individual tree crown delineation result and the location of the partial area; b. DOM of the partial area; c. Individual tree crown delineation of the partial area)

2.3 单木提取精度

利用研究区的DOM影像、CHM影像，采用目视解译，同时结合前期实地调研，人工随机选取100 个无瓣海桑样点，勾绘样点的树冠边界，进而生成单木提取结果的验证样本，验证效果如图7所示。考虑到树冠分割结果与样本之间存在多种空间关系，参考已有研究，将分割结果与验证样本之间存在重合且重合部分同时占分割结果与验证样本树冠50%以上（Zhen et al., 2013）定义为正确分割。采用探测精度（Detection Accuracy, DA）指标进行树冠分割精度评价（Yin et al., 2019），公式为：

图7 研究区局部无瓣海桑单木树冠提取结果与验证样本（a. 局部位置示意图；b. 局部区域正射影像；c. 局部区域单木分割结果与验证样本）Fig.7 Individual tree crown delineation result and reference samples of Sonneratia apetala (a. Location of the partial area;b. DOM of the partial area; c. Individual tree crown delineation result and reference crowns of the partial area)

式中：T为正确分割的树冠数量（株）；R为所有样本的数量；DA 为正确分割的树冠数量T与所有样本数量R的比例。

由于研究区内无瓣海桑冠层十分密集、群落内部冠层相互交叠、树冠间缝隙较小，且树冠之间的高度变化较为平滑，使得单木尺度的无瓣海桑提取工作具有挑战性。使用SRG算法得到的无瓣海桑单木树冠探测精度DA 为67%，基本能较为准确地将研究区的无瓣海桑单木树冠分割出来。

3 单木尺度的无瓣海桑地上生物量估算

3.1 无瓣海桑生物量估算模型构建

异速生长方程是林木生物量估算的经典方法，其通过建立树高、冠幅、胸径等林木参数与生物量之间的经验方程估算林木生物量（武高洁 等，2014）。已有研究表明，异速生长方程具有较高的生物量估算精度，且相比传统的皆伐法和平均木法，其能够大幅降低对森林的破坏性采样，适用于自然保护区等特殊政策环境下的受保护树种（Navarro et al., 2020）。同种树种的异速生长关系往往有较高的相似性，区域差异较小。因此，对不同地域的同一树种可以使用同一异速生长方程（Komiyama et al., 2008）。

参考Ren 等（2010）利用实地调查数据建立的无瓣海桑异速生长方程（表1）构建无瓣海桑单木地上质量估算模型。Ren 等的实验区域为广东省雷州湾红树林，该区域与广东省珠海市淇澳岛处于相近的纬度，且具有相似的红树林生境与树种组成。同时，Zhu等（2021）也将该异速生长方程应用于深圳福田红树林保护区17～18年树龄的无瓣海桑单木地上质量估算。基于表1中地上部分的干、枝、皮和叶生物量求和，得到无瓣海桑单木地上质量。

表1 基于胸径和树高建立无瓣海桑异速生长方程Table 1 The allometric growth equation of Sonneratia apetala based on DBH and tree height

已有研究表明，无瓣海桑的树高和胸径之间有较强的相关性（朱可峰 等，2011），因此，通过建立树高和胸径的回归方程，利用树高推算胸径，进而代入异速生长方程以估算无瓣海桑单木地上质量。利用地面采集工作获取的48株无瓣海桑单木树高和胸径，并采用最小二乘法进行回归拟合，得到研究区无瓣海桑的树高和胸径的回归方程，其公式为：

式中：D为胸径，单位为cm；H为树高，单位为m。

基于地面调查数据构建的无瓣海桑的树高和胸径线性回归趋势线（图8），可以看出，树高和胸径具有较高的相关性，相关系数R2为0.871 3。

图8 基于地面调查数据得到的无瓣海桑树高-胸径回归线Fig.8 Regression line of height-DBH of Sonneratia apetala based on ground survey data

利用上述树高-胸径回归方程，将异速生长方程（见表1）的胸径参数用树高参数表示，从而得到优化后的异速生长方程（表2）。

结合无瓣海桑单木提取结果，搜索每个树冠多边形内的最大值，将其作为无瓣海桑的单木树高，代入优化后的异速生长方程（见表2），得到无瓣海桑单木地上质量。进而计算无瓣海桑地上生物量，公式为：

表2 优化后的无瓣海桑异速生长方程Table 2 The optimized allometric growth equation of Sonneratia apetala

式中：AGB是单位面积内无瓣海桑地上生物量，单位为t/hm2；W为单位面积内无瓣海桑单木地上质量总和，单位是kg；A为单位面积，单位为hm2。

3.2 地上生物量估算结果与精度评价

从研究区无瓣海桑单木地上质量分布（图9）可以看出，无瓣海桑单木地上质量的范围为29.60～388.44 kg，平均值为145.72 kg，地上质量总和为368.97 t。

总体上，研究区内无瓣海桑单木地上质量分布存在部分空间集聚现象。采用全局Moran'sI指数（Moran, 1950）进一步评估研究区内无瓣海桑单木地上质量分布的空间集聚程度，得到其Moran'sI指数值为0.594（显著性检验结果p值为0，Z检验结果为50.934），表现为较强的空间聚集性。可以看出，林分边缘和林窗部分的无瓣海桑单木地上质量较小，如图10所示，主要可能有2个原因，一是部分研究区内的无瓣海桑正常死亡后留下林窗，其周围的幼苗仍处于生物量累积阶段，二是可能是由于林分边缘的无瓣海桑易受自然和人为因素的干扰，其单木地上质量常常受到损失。

图10 林分内部（a）与林分边缘和林窗处（b）的无瓣海桑单木地上质量对比Fig.10 Comparison of individual tree mass of Sonneratia apetala inside the forest (a) with that at the edge of the forest and at the gap (b)

通过对比基于无人机影像估算的无瓣海桑单木地上生物量、他人估算的研究区地上生物量与实地样方调查获得的地上生物量，进而验证利用无人机影像估算地上生物量的适用性和准确性。参考地面调查样方大小，将无人机估算得到的研究区无瓣海桑地上生物量分布分割为399个10 m×10 m的格网。图9-b 为按研究区10 m×10 m 样方尺度统计的无瓣海桑地上生物量分布，计算得到每个样方内无瓣海桑平均单木数为6.3棵，样方内最多单木数为17棵，最少为1棵；得到的样方平均地上生物量为92.14 t/hm2，最小值为2.99 t/hm2，最大值为247.24 t/hm2。根据前期实地样方数据显示，研究区内5 个10 m×10 m 无瓣海桑样方的平均单木数为5.3 棵，平均地上生物量为130.89 t/hm2，最小值为52.69 t/hm2，最大值为269.09 t/hm2。可以看出，基于无人机影像估算得到的研究区无瓣海桑生物量总体上低于地面实测生物量数据。Zhu 等（2020）利用WorldView-2高分辨率影像和无人机DSM 数据，采用随机森林算法对广东珠海淇澳岛红树林保护区地上生物量进行估算。首先将淇澳岛内植被主要分为秋茄、无瓣海桑与其他植物3种类型，将遥感影像的波段特征和植被指数作为输入变量，实地考察并通过异速生长方程得到的地上生物量作为输出变量对模型进行训练，进而得到整个保护区的地上生物量，RMSEr为30.48%。其得到的2016 年研究区平均地上生物量为92.53 t/hm2，最小值为39.24 t/hm2，最大值为249.99 t/hm2，总体上与本文估算的地上生物量结果一致。

图9 研究区无瓣海桑单木地上质量分布（a）和样方尺度地上生物量分布（b）Fig.9 Aboveground mass distribution (a) and quadrat scale biomass distribution (b) of Sonneratia apetala in the study area

4 结论

针对广东珠海淇澳岛红树林保护区的无瓣海桑群落，基于无人机影像和种子区域生长算法实现单木树冠提取，构建了无瓣海桑树高-胸径回归模型以优化生物量异速生长方程，进而实现了研究区无瓣海桑地上生物量的精确估算。得到的主要结论为：

1）基于无人机影像能够实现无瓣海桑的单木分割，且单木树冠探测精度为67%，并验证其在单木尺度无瓣海桑地上生物量估算中的可用性，可为后续基于无人机遥感的红树林生物量监测提供技术参考。

2）利用地面生物量调查数据分析了研究区无瓣海桑的树高和胸径间较高的相关性，相关系数R2为0.871 3，进而构建了基于树高估算无瓣海桑地上生物量异速生长方程。

3）基于无人机遥感估算得到研究区无瓣海桑的平均地上生物量为92.14 t/hm2，最低为2.99 t/hm2，最高为247.24 t/hm2，整体上略低于实地调查的地上生物量。

此外，本研究存在以下几点不足：1）由于无人机未携带RTK 系统，可能会造成垂直定高上约1 m的误差，同时无瓣海桑的树种特点为同一植株树冠可能存在多个树顶点，这些因素会造成单木分割的误差从而影响地上生物量计算；2）受到局地气候条件与潮汐等因素的影响，研究区内有部分无瓣海桑处于倒伏状态，树冠被其他未倒伏无瓣海桑树冠遮盖，而这些情况在高空视角往往难以获取，该部分无瓣海桑地上生物量在基于遥感的研究中通常不被考虑，但在研究区林下实地调查时却可以被统计，因而，基于无人机影像与高分辨率卫星影像估算的无瓣海桑地上生物量相比地面实测数据总体上偏低；3）由于研究区的无瓣海桑属于人工种植，相比滩涂区自然生长的无瓣海桑群落密度较低，因而本文单木分割算法中阈值的适用性在群落密度较高的区域仍有待探究。