基于树冠体素特征的森林机载-地基激光点云无控配准

摘要

高精度森林机载-地基点云配准能够极大提高点云数据完整性，为森林生态参数反演提供良好的数据支撑。针对森林环境特征稀疏下的机载-地基点云配准问题，提出了一种基于树冠体素特征的机载-地基激光点云数据无控配准方法。通过包围盒位置方向初始化、树高-频数的树冠点云滤波、树冠点云体素化、基于滑动模板匹配的特征提取以及旋转矩阵求解等步骤实现机载-地基点云配准。通过不同森林结构复杂度的4组数据对所提出的方法进行了实验验证，获得的平均距离差异分别为0.245m、0.238m、0.184m以及0.020m。与最近邻迭代法、人工配准及二者组合方法这3种方法进行对比实验，结果表明，本文提出的无控配准方法在精度和稳定性上均优于这3种方法。

关键词

森林; 激光雷达遥感; 点云无控配准; 体素化; 滑动模板匹配

激光雷达技术（light detection and ranging， LiDAR）作为一种主动遥感技术已广泛应用于各类林业资源调查任务中^［1］。其在信息更新频率、采集成本、数据精度、量化分析等方面比传统人工测量方式更具备优势^［2］。LiDAR遥感通过激光脉冲信号，以点云数据形式感知森林三维场景，通过分析处理激光点云数据，可计算得到树高^［3］、胸径^［4］、冠幅^［5］等森林结构参数，进而为蓄木量^［6］、固碳量^［7］、生物量^［8］等宏观森林生态参数信息反演提供数据基础。激光雷达遥感成果在森林动态监测、林业资源精细化管理等任务中日趋重要^［9］。

地基激光雷达（terrestrial laser scanning， TLS）遥感和机载激光雷达（airborne laser scanning， ALS）遥感是目前该领域任务中主流的森林三维结构感知方式^［10］。TLS以地面基站为平台，常用于单木、样地尺度测量^［11］。TLS基站通过水平和竖直面内的二轴旋转来获取关于目标对象千万级以上点云数据，能够准确捕捉到树干等林下结构^［12］。ALS以有人驾驶飞机或无人机为载体，扫描覆盖面广，在区域尺度森林资源调查等任务中较为常见^［13］。ALS激光脉冲穿透树冠间隙，并收集多次回波信息，以此探测不同冠层和森林垂直结构^［14-15］。但是，由于数据采集过程中激光传感器视角限制、激光脉冲能量衰减、树冠遮挡等因素，这2种采集方式均存在单一源数据缺失、数据质量不佳等问题^［16］。因此，如何高效、准确配准机载-地基点云，对于多源点云融合、提高数据完整性以及观测精度具有重要意义。

针对机载和地基激光点云匹配问题，最直接的方法是通过手工方式提取实地勘测过程中所布设的地面控制点^［17］或者从树顶^［18］和地形^［19］中提取特征，实现匹配。但匹配特征选取易受主观影响，其结果具有一定偏向性，推广性有限。近年来的研究更多集中于自动化的无控匹配，主要应用在针叶林等较为简单森林结构场景中。Liu等^［20］以冠层高度模型为基础，运用滑动窗口搜索方式提取具备相似灰度值的点对作为特征对，以此求解旋转矩阵，并运用最近邻点迭代（ICP）方法精化配准结果，配准结果距离差异为0.10~0.43m。但该方法在生成冠层高度模型中造成的信息损失会影响特征点对提取与最终匹配结果。Fekry等^［21］和Polewski等^［22］以ALS和TLS分割后的单木点云为基础，提取单木平面质心、树顶等特征对，进而求解旋转矩阵，配准精度为0.3~0.7m。这类方法的问题在于需预先进行点云分割，无疑大幅增加了计算成本。此外，Dai等^［23］通过提取树冠点云分布模型，运用相关点漂移方法提取特征，实现粗配准，再利用迭代最近点（iterative closest point，ICP）算法精化配准结果，精度约为0.3m。

针对上述问题，提出一种精度更高且可应用于多种森林复杂度结构的自动化空地点云融合匹配方法。该方法首先利用ALS和TLS轴对齐包围盒中心，与主方向向量进行初步定向，初始化相对位置，在此基础上分别提取不同源树冠层点云，并分别进行体素化，将数据空间转换为更加有序的数据形式，然后以TLS点云为模板窗口进行滑动，在ALS空间中寻找具备最大化匹配得分的位置，并返回匹配特征对，最后将问题转换为最小化ALS和TLS点云空间差异，通过奇异值分解解算旋转矩阵，并计算平移矩阵，最终实现机载与地基激光点云配准。

1 机载-地基点云无控配准方法

1.1　无控配准流程

如图1，所提出的机载-地基激光点云无控配准方法主要包括：点云包围盒定向、树冠点云滤波、树冠点云体素化、基于体素模板的特征对提取以及旋转平移矩阵计算。

图 1  机载-地基无控匹配流程

Fig. 1  Flowchart to show mark-free ALS and TLS point clouds registration

1.2　基于点云包围盒的位置初始化

该步骤目的是消除较大的平移旋转差异，先对ALS和TLS点云构建轴对齐包围盒，分别确定各自包围盒x、y、z主方向（图2）。在此基础上，对齐ALS和TLS点云包围盒中心，减小平移差异，按式（1）计算ALS点云和TLS点云各主方向向量夹角。

$$A=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{{\mathit{V}}_i^{\mathrm{a}\mathrm{l}\mathrm{s}}·{\mathit{V}}_i^{\mathrm{t}\mathrm{l}\mathrm{s}}}{\left|{\mathit{V}}_i^{\mathrm{a}\mathrm{l}\mathrm{s}}\right|\left|{\mathit{V}}_i^{\mathrm{t}\mathrm{l}\mathrm{s}}\right|}\right),$$$$A\in \{\gamma \mathrm{、}\beta \mathrm{、}\alpha ]$$

（1）

式中：A为单一分量的旋转角度，其取值范围$\mathrm{为}\gamma$、$\beta$、$\alpha ,$分别为相对于x、y、z轴的旋转角度；${\mathit{V}}_i^{\mathrm{a}\mathrm{l}\mathrm{s}}$和${\mathit{V}}_i^{\mathrm{t}\mathrm{l}\mathrm{s}}$分别为ALS点云和TLS点云单一主方向分量向量。

图 2  包围盒定向示意

Fig. 2  ALS and TLS position initialization with bounding box orientation

当角度计算完成后，分别按照式（2）、（3）、（4）分别计算相对于x、y、z轴的旋转矩阵${\mathit{R}}_x(\gamma )$、${\mathit{R}}_y(\beta )$、${\mathit{R}}_z\left(\alpha \right)$，最后依次将各主方向旋转矩阵应用于原始TLS点云上，减小其与目标ALS点云的差异，实现位置初始化。

$${\mathit{R}}_x(\gamma )=\left[\begin{array}{ccc}\mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{c}\mathrm{o}\mathrm{s}\gamma & -\mathrm{s}\mathrm{i}\mathrm{n}\gamma \\ \mathrm{0}& \mathrm{s}\mathrm{i}\mathrm{n}\gamma & \mathrm{c}\mathrm{o}\mathrm{s}\gamma \end{array}\right]$$

（2）

$${\mathit{R}}_y(\beta )=\left[\begin{array}{ccc}\mathrm{c}\mathrm{o}\mathrm{s}\beta & \mathrm{0}& \mathrm{s}\mathrm{i}\mathrm{n}\beta \\ \mathrm{0}& \mathrm{1}& \mathrm{0}\\ -\mathrm{s}\mathrm{i}\mathrm{n}\beta & \mathrm{0}& \mathrm{c}\mathrm{o}\mathrm{s}\beta \end{array}\right]$$

（3）

$${\mathit{R}}_z\left(\alpha \right)=\left[\begin{array}{ccc}\mathrm{c}\mathrm{o}\mathrm{s}\alpha & -\mathrm{s}\mathrm{i}\mathrm{n}\alpha & \mathrm{0}\\ \mathrm{s}\mathrm{i}\mathrm{n}\alpha & \mathrm{c}\mathrm{o}\mathrm{s}\alpha & \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]$$

（4）

1.3　基于高度与密度分布的树冠点云滤波

森林包括不同层级树冠。在不同高度，点云频数也呈现出一定规律性，同时高大树冠比低矮树冠形态更为清晰和稳定，因此从中提取的特征更为清晰与稳定。由图3a树高与点云频数分布折线图可知，示例森林样区（图3b）主要包括2级树冠层，分别对应2个波峰。通过区间定位波谷位置（局部极小值点），得到滤波阈值点，进而提取树冠上层点云，作为后续特征提取的数据基础。

图 3  基于树高与频数分布的树冠点云滤波

Fig. 3  Canopy point clouds filtering with height and density

1.4　树冠点云体素化

分别对ALS和TLS中所提取的冠层点云进行体素化，划分体素格网。通过格网表达，将无序点云转换为规则化的数据结构（图4），为后续模板匹配提供数据基础，同时也有利于降低数据量，提升算法整体效率。

图 4  树冠点云体素化示意

Fig. 4  Canopy point clouds voxelization

1.5　基于体素模板的特征对提取

针对森林点云数据难以有效提取不变特征的问题，提出一种区域模板匹配的方法。以最大化匹配得分（S_total）为目标函数，寻找最佳搜索位置，从而获取特征对。该模块主要包括体素过滤、模板搜索、最大化得分以及特征对提取四部分（图5）。

图 5  模板匹配流程

Fig. 5  Workflow for template matching process

在ALS和TLS体素质心集中，相同平面位置可能存在多个不同z值冗余体素中心，因此进一步通过限定最大高度值提取树冠上表层体素集。

在模板匹配过程中，以整体面积更小的TLS体素空间为模板，并将其坐标通过式（5）与ALS体素空间对齐，并通过x和y方向上的增量实现模板滑动。

$$\left\{\begin{array}{c}{P}_x={P_x^{\mathrm{t}\mathrm{l}\mathrm{s}}}_{}-P_{x\mathrm{m}\mathrm{i}\mathrm{n}}^{\mathrm{t}\mathrm{l}\mathrm{s}}+P_{x\mathrm{m}\mathrm{i}\mathrm{n}}^{\mathrm{a}\mathrm{l}\mathrm{s}}\\ P_y=P_y^{\mathrm{t}\mathrm{l}\mathrm{s}}-P_{y\mathrm{m}\mathrm{i}\mathrm{n}}^{\mathrm{t}\mathrm{l}\mathrm{s}}+P_{y\mathrm{m}\mathrm{i}\mathrm{n}}^{\mathrm{a}\mathrm{l}\mathrm{s}}\end{array}\right.$$

（5）

式中：$P_x$和$P_y$分别表示体素空间对齐后的x和y方向的坐标值。

在每个搜索位置，对ALS和TLS体素中心集整体建立KDTree^［24］，对KDTree中每一个TLS体素质心点检索其最近邻的ALS体素质心点，若存在最近邻点，且两者间Mahalanobis距离^［25］小于一个体素长度（$\epsilon$），则认定其为潜在匹配点对，增加一个匹配得分，见式（6）：

$$S_i=$$

$$\left\{\begin{array}{c}\mathrm{1},\mathrm{若}\sqrt[]{{\sum }_{i\in \{x,y,z\}}(P_i^{\mathrm{t}\mathrm{l}\mathrm{s}}-P_i^{\mathrm{a}\mathrm{l}\mathrm{s}}{)}^{\mathrm{2}}}<\epsilon \\ \mathrm{0},\mathrm{其他}\end{array}\right.$$

（6）

式中：$S_i$为单一质心点的匹配得分；$P_i^{\mathrm{t}\mathrm{l}\mathrm{s}}$、$P_i^{\mathrm{a}\mathrm{l}\mathrm{s}}$分别表示TLS和ALS质心集坐标分量值。因此，单一搜索位置匹配总得分则为潜在点对的匹配得分总和，见式（7）：

$$S_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}=\sum S_i$$

（7）

式中：$S_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}$为当前位置匹配总得分。

完成所有位置搜索后返回具有最大匹配得分位置的TLS和ALS特征点对，作为旋转平移矩阵求解步骤中的输入。

1.6　旋转平移矩阵求解

在森林场景下，ALS和TLS点云匹配为刚体变换，不考虑形变和缩放问题，因此，该问题转换为最小化TLS和ALS点云差异，见式（8）：

$$(\mathit{R},\mathit{t})=\begin{array}{c}\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}\\ \mathit{R},\mathit{t}\end{array}{‖(\mathit{R}\times {\mathrm{P}}_i^{\mathrm{t}\mathrm{l}\mathrm{s}}+\mathit{t})-P_i^{\mathrm{a}\mathrm{l}\mathrm{s}}‖}^{\mathrm{2}}$$

（8）

式中：R、t分别为旋转矩阵和平移矩阵。

$$\mathit{R}={\mathit{R}}_z(\alpha ){\mathit{R}}_y(\beta ){\mathit{R}}_x(\gamma )$$

（9）

分别对ALS和TLS点云去中心化，并构建协方差矩阵，见式（10）：

$$C\left({\mathit{P}}^{\mathrm{a}\mathrm{l}\mathrm{s}},{\mathit{P}}^{\mathrm{t}\mathrm{l}\mathrm{s}}\right)=$$

$$({\mathit{P}}_{i\in \{x,y,z\}}^{\mathrm{a}\mathrm{l}\mathrm{s}}-{\mathit{\mu }}_{\mathrm{a}\mathrm{l}\mathrm{s}})({\mathit{P}}_{i\in \{x,y,z\}}^{\mathrm{t}\mathrm{l}\mathrm{s}}-{\mathit{\mu }}_{\mathrm{t}\mathrm{l}\mathrm{s}}{)}^{\mathrm{T}}$$

（10）

式中：${\mathit{\mu }}_{\mathrm{a}\mathrm{l}\mathrm{s}}$和${\mathit{\mu }}_{\mathrm{t}\mathrm{l}\mathrm{s}}$分别为ALS与TLS点云质心坐标向量，在此基础上通过SVD奇异值分解方法获取U、S、V矩阵（式（11）），通过式（12）解算旋转矩阵R，最后将R矩阵回代式（8），利用式（13）计算得到平移矩阵t。

$$\left[\mathit{U},\mathit{S},\mathit{V}\right]=C(P^{\mathrm{a}\mathrm{l}\mathrm{s}},P^{\mathrm{t}\mathrm{l}\mathrm{s}})$$

（11）

$$\mathit{R}=\mathit{V}{\mathit{U}}^{\mathrm{T}}$$

（12）

$$\mathit{t}=-\mathit{R}{\mu }_{\mathrm{a}\mathrm{l}\mathrm{s}}+{\mu }_{\mathrm{t}\mathrm{l}\mathrm{s}}$$

（13）

最终将计算得到的旋转矩阵R和平移矩阵t应用于TLS点云，实现机载-地基点云配准。

2 实验与分析

2.1　实验数据

采用4种不同森林结构复杂度的空地点云数据，对所提出的无控配准方法进行验证与比较（表1）。

表 1  4种不同森林结构的机载激光雷达（ALS）与地基激光雷达（TLS）点云数据

Tab. 1  Four sets of ALS and TLS point clouds with different forest structures

点云数据类型	样地1	样地2	样地3	样地4
机载点云
地基点云

样地1为马来西亚婆罗洲某自然保护区热带雨林，机载点云区域大小约125m×125m，地基点云范围约为70m×70m，区域中主要包括石南树、龙果树等树种，最大树高约达158m。样地2为德国卡尔斯鲁厄地区阔叶、针叶混合林，点云数据来源于文献［19］，机载点云区域大小约为145m×145m，地基点云范围约为200m×200m，其中主要包括樟子松、橡木等树种，最大树高约为62m。样地3和样地4为分别为荷兰Veluwe地区的阔叶林和针叶林，点云数据均来源于文献［17］，机载点云区域大小分别约为95m×95m和60m×60m，地基点云范围分别约为102m×102m和51m×51m，主要包括桦木、云杉等树种，最大树高分别约45m和40m。4组数据按森林结构由复杂到简单排序为：样地1、样地2、样地3、样地4。4组样地ALS点云量均为百万至千万级，其密度为139~3 865个·m^-2，TLS点云数据量为千万至亿级，其密度为2 375~20 381个·m^-2。其他参数见表2。

表 2  机载激光雷达（ALS）点云与地基激光雷达（TLS）点云数据指标详情

Tab. 2  Details for ALS and TLS point clouds

指标	ALS测区尺寸/（m×m）	TLS测区尺寸/（m×m）	机载点云总量/个	地基点云总量/个	机载点云密度/（个·m-2）	地基点云密度/（个·m-2）
样地1	125×125	70×70	2 191 898	27 223 843	139	5 400
样地2	145×145	200×200	2 638 216	98 018 685	125	2 375
样地3	95×95	102×102	35 442 634	181 012 500	3 865	17 149
样地4	60×60	51×51	9 025 398	54 251 479	2 508	20 381

2.2　机载-地基点云配准精度评价指标

为避免倒树等可能造成无法匹配最近点的情况，通过人工采样ALS和TLS点云中树冠与树干位置的平均距离差异来描述机载-地基点云匹配成果质量。即构建KDtree建立空间索引，以转换后点云为基准，按照式（14）计算得到其到目标点云中最近点的平均距离差异。

$$D=\frac{{\sum }_{i=\mathrm{1}}^N‖{\widehat{\mathit{P}}}_i^{\mathrm{t}\mathrm{l}\mathrm{s}}-{\mathit{P}}_i^{\mathrm{a}\mathrm{l}\mathrm{s}}‖}{N}$$

（14）

式中：D为转换后的TLS点云与目标ALS点云的平均距离差异；${\widehat{\mathit{P}}}_i^{\mathrm{t}\mathrm{l}\mathrm{s}}$和${\mathit{P}}_i^{\mathrm{a}\mathrm{l}\mathrm{s}}$分别代表转换后的TLS点云坐标向量和目标ALS点云坐标向量。

2.3　实验参数选择

基于Python3.7.12软件实现所提出的方法，分别运用Open3D构建KDTree与检索最近邻点，以Numpy库实现矩阵运算与SVD分解，以multiprocess库并行加速。实验硬件环境为64G内存、英特尔2.6GHz CPU。

在基于树高-频数的树冠点云滤波环节中，通过区间对比可知，检测到的4块林区上层树冠点云高度分别约为120m、45m、32m、34m。此外，所选择的体素格网尺寸均定义为1。

对比实验选择了人工配准、最近邻迭代法^［26］（iterative closest point， ICP）以及上述二者组合的方法。人工配准方法是通过手动选取测区周边和中心区域9组几何特征显著点对作为匹配特征，进而求解旋转平移矩阵，实现点云匹配。在样地1中，特征点为树主干与枝干交叉点。在样地2、3和4中，树顶点作为特征点。单独使用ICP方法时，由于其对初始位置要求较高，因此手动将TLS点云移动至与对应ALS点云有重叠区域的位置，并且迭代参数选择为1 500，该参数在组合方法中选择为1 000。

2.4　机载-地基点云配准结果与分析

对4组方法在4组样地中机载、地基点云匹配结果按式（15）计算平均距离差异。

表3表明，本文所提出的方法在4组不同森林复杂度场景中均表现出较低的平均距离差异，分别为0.245m、0.238m、0.184m和0.020m。同时，人工方法也具有较高的精度，最大距离差异低于0.37m。相比下，基于位置初始化的单一ICP方法在这4组结果中均出现较为显著的不一致性（图6），其最小数值超过0.48m，最大差异达18m。并且在组合方法中，ICP算法在多数情况下对最终结果造成负面影响。除了样地2数据外，组合方法匹配在平均距离差异上相较于单一人工配准方法分别上升了0.378m、0.082m和0.325m。

表 3  机载点云与地基点云匹配结果平均距离差异

Tab. 3  Average spatial difference of ALS-TLS registration result ( m )

配准方法	样地1	样地2	样地3	样地4
ICP	18.013	23.044	0.486	5.245
人工配准	0.219	0.317	0.361	0.096
人工配准+ICP	0.597	0.289	0.443	0.421
本文	0.245	0.238	0.184	0.020

图 6  机载-地基点云配准结果细节

Fig. 6  Details of ALS and TLS point clouds registration results

此外，本文所提出方法在应对不同森林场景时也表现出优良的稳定性。其中，在结构最简单的针叶林中效果最好，其平均距离差异仅为0.020m。在其他更为复杂的混合林（样地3）、阔叶林（样地2）以及热带雨林（样地1）中，精度浮动在0.061m内。而人工配准与组合方法精度出现较大波动，上述3个样地中浮动分别在0.142m和0.308m。

对于ICP算法而言，尽管在执行算法前通过人工方式将ALS和TLS点云初始化为相似位置，但森林环境数据更加杂乱，特征对应关系更加模糊，若无稳定不变性的输入特征，很可能导致算法失效。这也是单一ICP算法失效以及其组合方法在大多数样地实验中出现精度下降的主要原因。

人工配准是人为选取树顶和树干交叉作为特征点，具备较好的稳定性，因而配准精度更高。然而树顶特征是不稳定的，通过图6可看到，样地3和4的匹配结果在冠层以下部分依然出现明显偏差，而树干交叉特征相对稳定，其配准结果与本文所提出方法基本一致。

相比而言，本文方法更优的原因主要有两方面：第一，树冠点云体素化解决了原数据空间的点云无序性与散乱性问题，体素分布更加规律，有助于后续特征提取。第二，在特征对选择中通过定义体素模板并通过网格滑动方式寻找在最优匹配位置下的特征点对，该过程本质是基于局部采样点集提取具备最佳树冠层体素拓扑结构相似度的配准特征对，而非简单地基于点特征或基于面特征的配准，因此算法精度和稳定性更优。

3 结论

针对森林场景，提出了一种基于树冠体素特征的机载-地基激光点云无控配准方法。首先初始化ALS和TLS点云包围盒方位来消除显著方位差异，然后以树冠层点云体素质心为基础，通过体素模板搜索过程计算最大匹配得分位置提取匹配特征对，并以此计算旋转平移矩阵，最终实现机载-地基点云配准。基于针叶林、混合林、阔叶林和热带雨林4组不同结构复杂度场景的机载-地基点云数据，与ICP、人工配准、人工配准与ICP结合这3种方法进行对比实验，表明所本文方法在精确度与鲁棒性方面均最优，平均距离差异分别为0.245m、0.238m、0.184m、0.020m。配准结果极大提高了森林场景点云数据完整性，有助于多源点云融合、森林结构生态参数反演等工作。

作者贡献声明

林 怡：研究方向确定、算法指导、论文修改。

曹宇杰：算法实现、数据采集、数据分析、论文撰写与修改。

程效军 ：研究方向确定。