联合UAV-LiDAR和HMLS技术的森林样地点云数据融合

王楚虹，刘浩然，钟 浩，林文树

（东北林业大学 工程技术学院，黑龙江 哈尔滨 150040）

激光雷达作为一种高精度的主动遥感技术，在林业应用领域具有传统光学遥感无可替代的优势[1-3]。按照平台的不同，激光雷达主要分为星载激光雷达、机载激光雷达、车载激光雷达及地基激光雷达等[4]。单一激光雷达技术平台测量森林结构参数时具有一定局限性，如地面移动激光雷达技术可快速准确地获得树干点云数据信息，但对于树冠点云信息的获得存在一定缺失，而无人机激光雷达技术虽然可以有效获取林分冠层信息，但对于树干点云信息却无法做到精确获取。因此，要利用激光雷达技术同时获取森林的水平和垂直结构信息，就需要进行多源遥感平台数据的融合来弥补单一激光雷达点云数据在林业中应用的技术瓶颈。

目前已经有相关学者将遥感影像作为激光雷达数据的补充应用于森林参数提取及三维模型建立中，其中部分研究是将遥感影像与激光雷达点云数据配准后进行联合应用。张吴明等[5]基于人为设置控制点的方法将无人机影像与地基激光雷达数据进行配准，并转换为在同一坐标系下的点云，以此提取树高。曹明兰等[6]将无人机倾斜影像生成密集点云后，使用迭代最近点（ICP,Iterative Closest Point）算法将密集点云与背包式三维激光扫描系统获取的点云配准，并构建三维模型，使2 种技术实现了优势互补。李丹[7]采用联合地基激光雷达与机载成像激光雷达的方法对森林参数进行提取，并以地基提取的 DEM 为基准，对机载影像解算的 DSM 进行地理坐标配准。激光雷达数据联合无需配准的同参考坐标系下的平台影像数据也可以有效提高森林参数获取的精度。骆社周等[8]处理星载激光雷达GLAS 数据建立反演森林叶面积指数模型，并利用神经网络融合星载激光雷达GLAS 与TM 光学遥感数据，以此实行区域连续高精度森林叶面积指数的反演，为生态环境研究精准输入参数提供了新思路。云增鑫等[9]利用航空激光雷达数据对高光谱数据进行林下植被信息的剔除，研究发现处理后的高光谱遥感影像能够有效改善森林冠层有效叶面积指数的估算精度。

不同平台的激光雷达点云数据也可以相互补充，因而有学者针对无标识不同平台激光雷达点云数据的融合进行了研究。2017年，Paris 等[10]利用基于点云数据栅格化的图像配准方法实现了地基和机载点云数据融合，但该方法适用于树木形状不规则的且较为开阔的林区，算法复杂度高。Polewski 等[11]使用模拟退火算法对二维树的位置点之间相对距离进行比较，完成地面移动激光雷达点云数据与空中激光雷达点云数据的配准，精度较好。Guan 等[12]提出了一种多平台激光雷达数据自动配准框架，即建立不规则三角网并通过三角形相似度投票策略来完成粗配准，为无人机激光雷达与地面移动激光雷达森林点云数据的融合提出了有效解决方法。Dai 等[13]通过树冠密度分析来融合机载激光扫描和地面激光扫描的森林点云，利用均值偏移法提取基于模型的关键点，并利用最大似然估计对关键点进行排序，该方法避免了对森林点云的几何描述。

随着以地面手持式、背包式及无人机为平台的轻小型激光雷达扫描系统的快速发展[14]，其在林业中的应用也越来越广泛。目前，基于地面移动与无人机平台的林地点云数据配准研究尚在起步阶段，大多数地空平台的激光雷达点云数据融合研究多应用于城市建筑方面[15-17]，而在林业应用方面的研究较少。另外，林业激光雷达点云数据配准中，部分研究通过利用周围建筑物或设立明显标志物完成[18]，选择的森林样地中林分也比较稀疏。因此面对复杂的林分环境，在不易设立地面标志物的情况下，如何将地面移动激光雷达数据与无人机激光雷达数据的精准与高效融合仍有待进一步研究。基于此，本研究提出一种基于Delaunay 三角网的手持移动激光雷达数据与无人机激光雷达数据的无标识融合方法，旨在以较高的效率及精准度完成2 种平台的配准，为提取完整的森林结构参数奠定基础，从而推动激光雷达技术在林业上的广泛应用。

1 研究样地与激光雷达数据获取

1.1 样地概况

研究样地位于哈尔滨城市林业示范基地，地理坐标为45°43′10″N，126°37′15″E。在樟子松人工林与蒙古栎人工林中分别选取2 块15 m×15 m的区域作为研究样地，其中，樟子松样地树木分布较松散，林中灌木较多；蒙古栎样地树木分布较密集，林木生长繁茂，郁闭度较高。选取的样地与基地内道路相邻，便于在道路上设置用于评价配准精度的标志物。为了使无人机激光雷达能够清晰地扫描到标志物，利用放置于三脚架上的白色标靶球作为标志物，并将三脚架放置在无树冠遮挡的位置，确保标志物被手持移动激光雷达与无人机激光雷达2 个平台完整扫描。

1.2 激光雷达数据获取

2020年10月，分别使用手持移动激光雷达系统和无人机激光雷达系统进行数据获取。手持移动激光雷达扫描系统选取GeoSLAM 手持移动SLAM 激光雷达ZEB-HORIZON，其测程为100 m，采集速度为每秒30 万点，获取点云的坐标系为相对坐标系。设置间隔3 m 左右的蛇形的行走路线，采集过程平稳缓慢，林木密集复杂处以及转弯处进行停留采集数据，确保树木数据的完整性。数据采集完成后，使用GEOSLAM-Hub 对数据进行SLAM 拼接，生成LAS 格式文件。

无人机激光雷达系统使用蜂鸟Genius 微型无人机LiDAR 系统，最大测距优于250 m，最大测量速度64 万点/s，实际作业飞行密度点大于200点/m2。实际飞行高度100 m，选取东北林业大学哈尔滨城市林业示范基地300 m × 160 m 的区域进行无人机激光雷达数据获取，确保研究样地包含在内。通过POS 数据解算，完成对无人机激光雷达点云的整合，点云坐标系为WGS-84 大地坐标系。

2 HMLS 和UAV-LiDAR 数据融合方法

在融合前，需要对手持移动激光雷达数据和无人机激光雷达数据进行预处理。首先，将2 种平台得到的点云数据进行滤波去噪。然后将无人机激光雷达点云数据根据分离出的地面点归一化处理；将手持移动激光雷达点云数据进行均匀抽稀处理，保证其在不丢失树木自身特征的前提下与无人机激光雷达点云的密度基本一致。最后将处理后的手持移动激光雷达点云数据和无人机激光雷达点云数据进行单木分割，分别提取树干的中心点与树冠的几何中心点作为树木的位置点。

融合手持移动激光雷达数据与无人机激光雷达数据的方法分为粗配准和精配准2 个步骤。将预处理后的2 个不同平台点云数据基于树木位置点构建Delaunay 三角网进行粗配准并利用模拟退火算法优化配准过程。选取一定高度区域，将粗配准的手持移动激光雷达点云数据基于无人机激光雷达点云数据运用ICP 算法进一步精配准调整。选取樟子松样地和蒙古栎样地对融合方法进行验证，下文中的点云融合过程以樟子松样地为例，点云融合平台基于Python 编译。具体流程如图1所示。

图1 点云配准流程Fig.1 Point cloud registration flow chart

2.1 Delaunay 三角网的构建

构建Delaunay 三角网是粗配准的第一步，也是提高粗配准效率的关键步骤。Delaunay 三角网是一系列相连的、不重叠的三角形的集合，该集合内任意一个三角形的外接圆中不包含平面内其他点[19]。对于同一个点集而言，从任意一点构建Delaunay 三角网，所得结果都是一致的。本研究分别将预处理提取的手持移动激光雷达和无人机激光雷达数据树木位置的三维坐标投影于平面，形成2组二维平面位置点。由于树木位置是固定的，又因相同离散点构建的Delaunay 三角网具有唯一性，故使用2 组树木的二维位置点集合分别构建的2 个Delaunay 三角网形状应基本相同。

本研究选择效率高、使用最为广泛的逐点插入算法构建Delaunay 三角网[20]。使用逐点插入算法构建基于树木位置点的Delaunay 三角网的基本步骤如下：1）根据树木位置点的最大分布来求得随机一个超级三角形，并将其加入临时三角形集合中；2）依次插入树木位置点，找出临时三角形集合中外接圆包含插入点的三角形T，删除三角形T 的公共边，并将插入的树木位置点与T 的顶点连接，完成该点的插入，将形成的三角形放入临时三角形集合；3）重复2）步骤，直到所有树木位置点插入完成。

2.2 对应点的确定

理想情况下，2 组树木位置点构建的Delaunay三角网形状应完全相同，但是树木郁闭度高等因素会影响无人机激光雷达点云数据单木分割的精度，导致例如树木位置点未被记录或记录错误等偏差出现。另外，手持移动激光雷达点云单木位置是树干中心，而无人机激光雷达点云单木位置是树冠中心，在树干倾斜的情况下2 点在同一平面的投影点不能重合。因此在2 个由树木位置点构造的三角网仍可能会出现较大差别的情况下，需通过设定限制确保2 组树木位置点之间建立准确的对应关系。

以Delaunay 三角网为基础，通过以下限制搜索2 种平台所得到的树木位置点。1）搜索所有线段，设定合理阈值，记录差值在阈值内的线段及其对应线段，两两构成线段组。2）将线段组整合，记录手持移动激光雷达点云中同一树木位置点出发的线段所对应的无人机激光雷达点云线段。3）寻找对应线段相交的公共点，该点即为手持移动激光雷达点云中树木位置点的对应点。在搜索对应点的过程中，设置线段差阈值为可调节的参数。当线段长度差的阈值设置极小时，搜索过程中该方法对于线段与其对应线段的长度偏差的容忍度极低。由于手持移动激光雷达点云树木位置点与无人机激光雷达点云树木位置点并非全部垂直对应，当线段差阈值设置极小时，容易出现无法搜索到足够对应点的情况。反之，将线段长度差的阈值设置极大，表明搜索过程中容忍度极高，这会大大增加搜索进程时间，使得搜索效率降低，同时有可能出现树木位置点对应错误的问题，导致配准失败。

为了提高粗配准的效率，本研究利用模拟退火算法优化搜索无人机激光雷达树木位置点对应的手持移动激光雷达树木位置点这一过程。模拟退火算法是一个搜索全局最优解的算法，是寻找温度最低时概率分布中具有最大概率的状态，并在搜索过程中设置了一个变化且最终趋于零的概率突跳性，避免搜索陷入局部最优解[21]。基于模拟退火优化算法的具体配准过程如下：1）设定较高的初始温度T，随机选取可能的无人机激光雷达点云树木位置点与其手持移动激光雷达点云对应点作为初始模型A0，计算其相应的目标函数值S(A0)，目标函数见式（1）；2）对当前模型进行扰动产生一个新的可能的无人机激光雷达点云树木位置点与其手持移动激光雷达点云对应点作为模型A，计算相应的目标函数值S(A)，得到ΔS=S(A)-S(A0)；3）若ΔS＜0，则新模型A被接受；若ΔS＞0，则新模型A按概率P=exp(-ΔS/T)进行接受。当模型被接受时，A0=A；4）温度T下，多次进行扰动和接受；5）逐渐降低温度T；6）重复步骤，至结果收敛或者到达迭代次数，迭代次数这一参数可在配准过程中根据点云数据的数据量进行相应调节。模拟退火算法的目标函数设为下式。

式（1）中，(x,y)表示无人机激光雷达点云树木位置目标点在平面的投影点，(xi,yi)为Delaunay 三角网中(x,y)的相邻点，(xi,yi)与(x,y)两点构成目标线段。(X,Y)为疑似手持移动激光雷达点云树木位置对应点在平面的投影点，(X,Y)与其相邻点(Xi,Yi)的距离为与目标线段长度相似的对应线段。通过计算，确定(X,Y)是否为(x,y)真正的对应点。寻找所有目标点的对应点，使用模拟退火算法找到最优对应点集合，舍弃一些无法确定的点对。

2.3 旋转平移矩阵

利用模拟退火算法搜索得到手持移动激光雷达树木位置点与无人机激光雷达树木位置点对应的点集合，利用奇异值分解法根据对应点对的集合求得旋转矩阵R与平移矩阵t[22]。

2.4 迭代最近点算法（ICP）

使用ICP 算法对粗配准后的手持移动激光雷达点云数据与归一化的无人机激光雷达数据进行精配准。ICP 算法是基于最小二乘法的最优匹配，通过将点与点不断进行匹配后，对点云旋转平移，直到找到最优变换矩阵，使目标点云与待配准点云重合。ICP 算法需要2 组点云形成一致或者包含的关系，才能确保配准正常运行[23]。手持移动激光雷达点云数据包括树干以及部分树冠信息，而无人机激光雷达点云数据包括树冠、部分树干及少量地面点，由此可知2 组点云数据中，重合的点云大部分属于树干上部及树冠中下部区域。为保证ICP 算法能够获取理想配准结果，选取一定高度区域的手持移动激光雷达点云与无人机激光雷达点云进行精配准。通过ICP 算法，得到最终精配准结果，使手持移动激光雷达点云转换到无人机激光雷达点云的同一坐标系下。

2.5 融合精度评价

本研究从2 方面对融合结果进行精度评价：

1）根据配准后2 个平台的树木位置投影点坐标之间的偏移距离评价。提取手持移动激光雷达点云树干中心在平面上的投影点坐标(XH,YH)及无人机激光雷达点云树冠中心在同一平面上的投影点坐标(XU,YU)，计算无人机激光雷达点云树冠中心点坐标相较于手持移动激光雷达点云树干中心点坐标的偏移距离D，计算公式见式2 以检验2个样地的融合效果。

2）根据配准后2 个平台的标志物点云数据中心点坐标之间的均方误差（mean-square error,MSE）评价。设定标志物的无人机激光雷达数据为真值，手持移动激光雷达数据为估计值，由此来计算均方误差。均方误差的计算公式见式3，其中(xi,yi,zi)为标志物的手持移动激光雷达点云中心点坐标，(x,y,z)为标志物的无人机激光雷达点云中心点坐标。

3 结果与分析

3.1 树木位置提取

经过预处理后，手持移动激光雷达点云与无人机激光雷达点云的树木位置点如图2所示,与样地中树木实际位置相比较可得分割结果较为精准。将图2a 顺时针旋转约90°后，图2a 与图2b 形状相似，即图2a 中树木位置点的分布与图2b 中树木位置点的分布大致相同。

图2 树木位置点Fig.2 Tree locations

3.2 点云数据粗配准

如图3所示，使用由2 个平台激光雷达点云提取的树木位置点分别构造2个Delaunay三角网。图3a 三角网中大部分三角形可以在图3b 三角网中找到对应的三角形。Polewski 等通过比较2 种点云数据中的单个位置点与其他位置点间的所有距离来实现对应点的确定[11]，而本研究提出的建立Delaunay 三角网后再比较位置点距离的方法来搜索对应点的效率更高。

图3 树木位置点构造Delaunay 三角网Fig.3 Delaunay triangulation of tree locations

经过反复试验可得，粗配准中线段长度差阈值为0.8 m 时，树木匹配数量最多且结果准确。模拟退火算法中迭代次数的设置对于配准过程的影响较小，如果迭代次数太小可能会使得结果出现偏差，同样经过试验得到设置迭代次数为15 时配准结果较准确且运行效率较高。

经过粗配准后，不同平台数据得到的树木位置点的对应关系如图4所示，图4中同种颜色圈出的点为一组对应点，由于树木倾斜及无人机激光雷达点云数据树木位置点分割不准确等原因，部分树木位置点不能很好地找到其对应点，如图4中用“×”符号标出的树木位置点。只要求得3组以上对应点即可计算旋转平移矩阵，因此存在未搜索到对应点的目标点不会影响粗配准的进行。

图4 树木位置点的对应关系Fig.4 The corresponding relationship between tree locations

手持移动激光雷达点云去除高度值后的二维坐标利用旋转平移矩阵进行二维旋转平移。将旋转平移后得到的粗配准结果点云数据匹配点云数据原有的高度值，从而获取粗配准后的三维点。选取樟子松样地中心位置的任意一棵树木作为示例，其粗配准效果如图5所示，图中黄色代表手持移动激光雷达点云，蓝色代表无人机激光雷达点云。由图5可得，树干与树冠部分均有所偏移。

图5 粗配准结果Fig.5 Rough registration results

3.3 点云数据精配准

针对有偏移、误差较大的粗配准结果，对点云进一步进行精配准。选取一定高度区域的点云数据，即手持移动激光雷达点云与无人机激光雷达点云均较密集的区域，运用ICP 算法配准，实现2 种平台点云数据的精准融合。ICP 算法选取的点云区域与获取激光雷达点云数据的仪器性能、树木平均高度以及林区郁闭度均有很大相关性。由于手持移动激光雷达点云数据与无人机激光雷达点云数据重合度高的区域为树冠区域，因此本研究通过对2 块样地中不同树木点云数据高度区域范围进行选择，当处于平均树高4～6 m 的范围时可以获取较好的精配准结果。

针对樟子松样地的点云数据，选取14～15 m高度区域的点云数据进行精配准。精配准计算得到旋转矩阵：平移矩阵：

将粗配准后的手持移动激光雷达点云利用旋转平移矩阵进行坐标转换，从而实现点云数据的精配准。精配准后，与图5粗配准结果相比，图6中显示的同一棵树的树冠和树干部分融合效果得到了改进，树干与树冠部分基本没有偏移。

图6 精配准结果Fig.6 Precise registration results

3.4 融合效果与精度

3.4.1 融合效果

经过2 次配准后，基于手持移动激光雷达与无人机激光雷达扫描的2 块样地的点云数据融合效果较好，结果如图7所示，图中蓝色点云为无人机激光雷达点云。由图7可以看出，无人机激光雷达点云与手持移动激光雷达点云的树冠部分基本重合，2 个平台点云相互补充。

图7显示了样地整体的配准效果，然而为了进一步检验样地内部树木的融合效果，随机选取样地内部区域的树木进行检验，选取的区域见图8中的方框区域。

图7 样地配准结果Fig.7 Registration results of plots

图8 随机区域选取Fig.8 Selection of random regions

图9和图10分别是樟子松样地与蒙古栎样地选取区域内部树木的激光雷达点云图，其中图9a与图10a 为无人机激光雷达点云图。由方框区域可见，无人机激光雷达的穿透性使其能够获取样地中树干的部分点云数据。由于樟子松样地的郁闭度更低且树木分布更稀疏，与图10所示蒙古栎样地无人机激光点云树干处点云相比，图9所示樟子松样地的树干处点云更加清晰可见。由无人机激光雷达点云数据可以大致确定树木的生长方向，该方向与图9b、图10b 中图所示手持移动激光雷达点云树干生长方向基本一致。图9c、图10c中黄色点云为手持移动激光雷达点云，蓝色点云为无人机激光雷达点云，可以看到随机选择的样地内部区域中，2 个平台获取的数据在树干和树冠区域均得到有效配准。

图9 樟子松样地Fig.9 Pinus sylvestris plot

图10 蒙古栎样地Fig.10 Quercus mongolica plot

3.4.2 融合精度

1）根据2 个平台的部分树木位置坐标计算偏移距离。针对随机选择的区域，经过提取区域内树木的手持移动激光雷达点云树干中心及相应的无人机激光雷达点云树冠中心，计算两个中心投影点的偏移量可得表1～2，其中XH与YH分别表示手持移动激光雷达点云树干中心点投影点的平面X 坐标与Y 坐标，XU与YU分别表示无人机激光雷达点云树冠中心点投影点的平面X 坐标与Y坐标。通过计算可得，樟子松的投影点坐标偏移距离D 最大和最小值分别为0.30 和0.09 m，而蒙古栎的投影点坐标偏移距离最大和最小值分别为0.59 和0.06 m。

表1 樟子松样地树木位置点坐标偏移量Table 1 Offset of tree locations in Pinus sylvestris plot

表2 蒙古栎样地树木位置点坐标偏移量Table 2 Offset of tree locations in Quercus mongolica plot

2）根据同一标志物在2 个平台中的坐标计算均方误差。通过计算样地上放置的标志物在2 种平台中的坐标之间的偏差来得到点云配准后的均方误差，结果表明，樟子松和蒙古栎样地融合后均方误差分别为0.051 2 和0.080 2，与蒙古栎样地相比，樟子松样地树木融合精度更高。这可能是因为樟子松样地为针叶林，无人机激光雷达点云中树冠更易被准确分割，因而由此提取的树木位置点更加精准，同时樟子松样地的郁闭度较低、树木分布较离散，所以樟子松样地的融合精度更高。

4 讨 论

与普通测绘相比，林业测绘具有地形复杂，树木树冠不规则、非刚体，森林点云特征较普通测绘中建筑物、道路点云特征边线更难提取等特点，因此特征的选取是森林点云配准的关键。本研究通过树木位置构成的Delaunay 三角网作为配准特征，有效解决了森林特征难提取的问题，达成手持移动激光雷达点云与无人机激光雷达点云精确配准的目标。

手持移动激光雷达和无人机激光雷达扫描所得数据尺度均为树木实际尺度，但由于不同平台设备采集数据的角度不同，所以不同平台点云数据的特点也不同，点云的精准融合需要考虑不同平台点云数据的特点。手持移动激光雷达自下而上采集数据，点云密度自下而上递减。无人机激光雷达自上而下采集数据，树冠处点云较密集，而森林下层数据采集数量较少。2 种设备的不同特点导致采集到的点云数据整体重叠度较低，但是针对树木中上层点云数据，手持移动激光雷达数据与无人机激光雷达数据重叠度较高，较高的重叠度是使用ICP 算法进行精配准的前提。

对于不同的树种，融合精度有一定差异。针对阔叶林，如本研究中的蒙古栎，由于郁闭度较高，以及阔叶树树冠形状无规则的特征，无人机激光雷达点云数据的单木位置提取效果较差，导致后期配准结果容易出现误差。而针对冠幅较小、树干生长垂直于地面的针叶林，如本研究中的樟子松，无人机激光雷达点云数据单木分割效果较好，故配准结果也更加精准。

本研究结果为多源遥感数据在林业中的广泛应用提供了可能，有利于克服单一平台无法获取树木完整点云的局限性，对于森林结构参数的精确提取和树木三维模型的精准构建有重要意义，但仍存在以下不足：1）本研究区域为人工林，地形平缓，树木种类较单一，针对树种及地形更加复杂的天然林，本算法的适用性仍有待验证。2）虽然模拟退火算法使得融合方法运行效率有所提升，但当本方法应用于大区域尺度时，仍存在花费时间较长的问题。同时随着区域的增大，树木位置对应点的搜索结果也更易出现误差。3）本研究虽然将手持移动激光雷达数据与无人机激光雷达数据精准融合，但是并未对融合后的数据进一步精简处理，数据仍存在冗余。

在未来研究中将尝试选取地形起伏较大、林分结构更复杂的天然林作为研究对象，进一步验证算法的适用性。另外，对算法进行优化以提高算法的自动化程度及运行效率，减少时间成本。针对融合后数据冗余问题，基于森林结构特点提出点云压缩算法，为后续提取高精度森林结构参数、构建三维可视化模型奠定基础。

5 结 论

本研究利用树木位置点，提出了一种基于Delaunay 三角网和ICP 算法的不同平台点云数据的融合方法，该方法属于无标识配准方法，并通过2 块复杂程度不同的样地获取得到的点云数据进行了验证。结果表明，在不同的参数设置下，树木匹配的精度可能不同。粗配准中设置的线段差阈值和精配准中设置的高度区域范围这2 个参数直接影响配准结果的精准度。而粗配准中设置的迭代次数的变化对精准度结果的影响较小，但是迭代次数的合理选择可以提高方法运行效率。

本研究方法可以较为准确地将手持移动激光雷达点云和无人机激光雷达点云进行融合，融合效果较好。融合后2 块样地随机选取区域中的2种平台数据树木位置投影点的偏移距离均值分别为0.19 和0.25 m。根据配准后标志物在2 种平台中的位置坐标，计算可得融合后2 块样地的均方误差均小于0.1，融合精度较高。本研究结果可为单木结构参数精确提取、树木地上生物量估算以及三维模型精准构建提供基础，为多源数据在林业中的应用提供参考。