0.   引言

地面激光扫描是一种有效的林业调查工具，可以用于植被点云提取^[1]单木位置^[2]、胸径^[3]、叶面积指数^[4]、树种^[5]等森林参数的提取和估算。由于前视遮挡、扫描范围有限等原因，单测站激光扫描往往难以覆盖整个场景，因而需要架设多个测站并将不同测站扫描数据的坐标转换至同一坐标系下，该过程为点云配准。森林激光点云配准可以为林业调查提供完整的基础数据，是森林激光扫描数据处理的重要环节。

传统森林激光点云配准方法有人工配准^[6]、布设标靶^[7]等，但森林场景中树干、树冠等主要地物相似性较高，且前视遮挡较为严重，人工配准和布设标靶方法在大范围森林场景下效率低。因此，无标靶的森林点云配准方法对提高森林环境下的激光扫描数据采集和使用效率有重要的意义。

森林环境下的无标靶激光点云配准通常使用树干或树冠提取得到的配准基元。LIANG等人使用单木位置点构建水平向量作为配准基元，通过长度、旋转角等特征进行筛选并计算水平转换参数，但由于复杂度较高，该方法仅计算了水平转换参数^[8]。在此基础上上，LIU等人提取树干不同高度的中心坐标和半径，用以计算竖直方向的平移参数^[9]。KELBE等人基于树干位置点集构建三元组，并使用坐标、胸径和协方差矩阵特征值作为特征描述，查找可能的同名树干位置点实现森林点云配准^[10]。由于需要遍历所有三元组对，该方法复杂度较高。针对这个问题，TREMBLAY等人将多线程运算应用于三元组匹配过程，并使用三元组边长替代协方差矩阵特征值，通过调整存储方案使算法更适用于多线程策略，以提高配准的效率^[11]。与TREMBLAY等人思路相近，JIN等人将三元组的边长作为特征值，通过K维(K-dimensional, KD)树索引提高三元组匹配的效率^[12]。DAI等人通过均值漂移在树冠处提取特征点实现粗配准，然后筛选重叠区域的单木位置点进行精配准，与直接使用所有树干位置点的配准结果相比，该方法在保证配准精度的同时提高了配准效率^[13]。在面向单木的配准基元检测思路之外，GUAN等人提出从树干阴影区域检测特征基元的方法^[14]。直接从激光点云中提取配准基元往往复杂度较高，在复杂森林存在稳定性不足的缺点，部分研究人员通过使用额外仪器获取辅助信息来简化配准过程，例如，通过测斜器获取激光点云绕x轴和y轴的旋转参数^[15]，或者通过反射片后视定向获取测站间平移参数^[16]以及通过实时动态(real-time kinematic，RTK)设备获取测站位置^[17]。

已有的研究表明，激光扫描辅助信息可以提供部分转换参数初值^[18]，简化森林激光点云配准过程，但需携带额外设备或对扫描仪功能提出更多要求，会增加林下作业负担或限制方法的适用范围。此外，基于树干位置点的配准方法大多使用所有位置点进行配准。尽管有很多研究显示，由于地面点或树干点密度的下降甚至缺失，远距离树干位置点精度会随之降低，进而影响配准精度^{[11, 13]}，但对树干位置点进行筛选的研究仍然不足。针对上述问题，本文作者提出一种以智能手机辅助的森林激光点云配准方法，使用智能手机获取测站定位和定向数据等辅助信息实现粗配准，然后构建基于2级目标函数的精配准方法，同时在该过程中基于测站间距对单木位置点进行筛选以提高配准精度。

1.   方法

1.1   智能手机数据采集与点云预处理

本文中使用智能手机的内置传感器辅助激光扫描和粗配准过程，由于智能手机定位定向数据仅用于粗配准，因此在每站扫描前将手机置于扫描仪上，无需精确校准二者相对位置关系。通过手机陀螺仪实现扫描仪粗略置平，并测量每站起始扫描方向的方位角，即起始扫描方向与北向的顺时针夹角。由于地面激光扫描仪通常以起始扫描方向作为局部坐标系x轴，因此可得测站局部坐标系x轴的方位角，实现定向数据获取。通过手机定位设备获取设站点的大地坐标(即设站点的经纬度)，可实现定位数据获取。智能手机的定位和定向功能可通过软件开发工具包(software development kit，SDK)调用或手机地图应用获取。

对每站激光点云数据，使用树干位置点集作为输入，即一系列用于描述每颗树木位置的3维坐标集合。本文作者使用参考文献[19]中的方法提取单木位置点，通过随机森林算法提取树干点，并分段拟合树干的圆柱模型，计算每棵树对应的最低圆柱轴线与地面模型的交点作为单木位置点。

1.2   使用定位定向数据的粗配准

本文中采用多测站逐对配准的方式进行转换参数计算，基本思路是在全局投影坐标系下计算测站间的平面转换参数。首先基于智能手机定向信息将两站点云转换至x轴为北向时的扫描结果，如图 1a所示；然后基于参考站和待配准站定位数据计算测站间的平移参数，如图 1b所示。

图  1  使用定位和定向数据表示测站间的转换关系

Figure  1.  Relation between two scans based on the positioning and orientation data

下载: 全尺寸图片 幻灯片

转换参数计算如下:

式中：$\boldsymbol{p}_{\mathrm{r}}=\left[x_{\mathrm{r}}, y_{\mathrm{r}}, z_{\mathrm{r}}\right]^{\mathrm{T}}$和$\boldsymbol{p}_{\mathrm{s}}=\left[x_{\mathrm{s}}, y_{\mathrm{s}}, z_{\mathrm{s}}\right]^{\mathrm{T}}$分别表示两测站间的对应点; α_{r, x}和α_{s, x}分别为p_r和p_s对应的x轴方位角，地面激光扫描通常将起始扫描方向作为x轴，因此α_{r, x}和α_{s, x}可通过智能手机定向数据直接获取; R (0, 0, α)为旋转矩阵，表示绕z轴逆旋转α角度; M为平移矩阵，可通过智能手机定位数据得到:

式中：(X_r, Y_r)和(X_s, Y_s)分别为p_r和p_s的全局投影坐标，可由智能手机定位得到的大地坐标转换得到，由于投影坐标系和测站坐标系分别为左手系和右手系，二者对应的y轴的平移方向相反；D_{r, s}为设站位置p_r和p_s间的水平距离; β_{r, s}为p_r和p_s连线与y轴正方向之间的顺时针夹角。

在此基础上，通过变换式(1)可得粗配准转换参数，即旋转矩阵R和平移矩阵t。

1.3   使用两级目标函数的精配准

由于智能手机定位和定向数据可能存在较大误差，由式(4)得到的转换参数通常存在较大误差。在精配准中，本文作者提出基于树干位置点的两级评价函数，以参考站和待配准站的树干位置点集M_r和M_s作为输入，构建一系列转换参数并进行评估。

1.3.1   第1级目标函数

对式(4)，若将β_{r, s}作为常数，则转换参数由智能手机定位定向所得数据α_{r, x}、α_{s, x}和D_{r, s}计算得到。以α_{r, x}、α_{s, x}和D_{r, s}作为初始值，可构建相应搜索空间:

式中: α_r和α_s分别为参考站和待配准站x轴方位角; D为测站位置间距; 设角度α_r和α_s的搜索间隔设为δ_α; D的搜索间隔为δ_D。在搜索区间内，共构建N_t组参数，对每组(α_r, α_s, D)使用式(4)构建转换参数，对树干位置点集M_s进行坐标转换，并使用互为最邻近点的方式搜索其在M_r内的水平邻近点，得到点对集合$M_{x, y}=\left\{\left(\boldsymbol{p}_i \in M_{\mathrm{s}}, \boldsymbol{q}_i \in M_{\mathrm{r}}\right)\right\}$, 且i=1, 2, …, n。使用水平距离均方根误差构建第1级评价函数:

式中: (p_i, q_i)表示M_{x, y}中第i对点；向量n =(1, 1, 0)^T；符号“ ”表示在搜索间隔δ_α、δ_D下得到每一组(α_r, α_s, D)；$\|\cdot\|^2$表示求范数，即点至原点的距离。选取函数值最小的N_r组参数并记录其对应转换参数R和t，然后使用迭代最近点(iterative closest point, ICP)算法^[20]精化相应的转换参数，得到转换参数集合{(R_j, t_j)}, 且j=1, 2, …, N_r。由于点云密度随距离增加而减小，同时遮挡随扫描距离而增多，导致树干位置点的精度通常随扫描距离增加而降低，因此本文中仅使用到原点水平距离不超过D的树干位置点参与ICP迭代。

1.3.2   第2级目标函数

使用每组转换参数(R_j, t_j)对M_s进行坐标转换，然后基于互为最邻近点的方式搜索其在M_r内的3维邻近点，得到M_{x, y, z}={(p_k∈M_s, q_k∈M_r)}, 且k=1, 2, …, m。在此基础上，以对应点数目构建第2级评价函数：

式中: (p_k, q_k)表示M_{x, y, z}中第k对点; m为点对数目; ε为距离阈值，表示对应位置点的最大允许偏差。本文中将ε设为0.2 m。当点对间距超过ε时，通过函数f(γ)将其标记为无效点对。

式中: $\gamma=\left|\boldsymbol{R} \boldsymbol{p}_k+\boldsymbol{t}-\boldsymbol{q}_k\right| <\varepsilon ; N\left(\boldsymbol{R}_j, \boldsymbol{t}_j \mid \varepsilon\right)$表示满足条件的邻近点对数目。最终取最大评价函数值对应的转换参数R ^*和t ^*作为最后的精配准结果。

若多组转换参数得到的临近点对数目相同，则取均方根误差最小的转换参数作为配准结果。

1.4   配准结果的精度评定

通过计算上述转换参数与参考值之间的差异进行精度评估。对配准结果，使用旋转矩阵R ^*可得x、y和z轴分别对应的旋转角，记为φ_e、ω_e和κ_e; 将平移矩阵t ^*对应的分量记为t_{x, e}、t_{y, e}和t_{z, e}，即沿x、y和z轴的平移分量。将旋转角和平移量的参考值分别记为φ_r、ω_r、κ_r、t_{x, r}、t_{y, r}和t_{z, r}，则本文中使用旋转角误差(rotation angle error, RAE)、水平平移误差(translation error, HTE)和竖直平移误差(vertical translation error, VTE)评估配准结果，分别用e_RAE、e_HTE和e_VTE表示:

2.   结果与分析

2.1   实验数据介绍

实验数据采集于湖北省武汉市某树林区域(30.53°N，114.35°E)，采集时间为6月，在场景内部和边缘共架设5个测站，设站位置如图 2所示。场景树种以梧桐树和栾树为主，为每公顷370棵树的密度，并有少量灌木，地面以草类植物为主，并有少量步道、石桌、垃圾桶等人工地物，场景边缘有少量建筑物。扫描仪为Reigl-VZ400，角分辨率为0.04°，扫描视角为360°×100°，通过华为P20手机(集成单频定位模块)测量设站位置大地坐标和起始扫描方向的方位角，测站信息如表 1所示，定向误差通过与罗盘仪对比得到。

图  2  试验区域与设站位置

Figure  2.  Test area and scanner positions

下载: 全尺寸图片 幻灯片

表  1  实验数据信息

Table  1.  Information of test data

station	number of points	trunk position count	canopy thickness/m	irientation error/(°)
s1	22671334	214	3.3	11°
s2	21586721	227	2	17°
s3	24659841	186	7.3	5°
s4	24208730	166	11.5	12°
s5	24783758	259	9.5	54°

下载: 导出CSV 

| 显示表格

2.2   配准结果对比

将本文中的方法与参考文献[11]和参考文献[12]中的方法进行对比，结果如表 2所示。选取的6对测站间距分布在19 m~40 m，选取测站间重叠区域内的部分人造地物和裸露地面点，通过人工选取同名点并进行ICP精配准的方式得到转换参数参考值，用于精度评定。本文中方法的主要参数δ_α，δ_D和N_r分别设置为4°、1 m和1%×N_t(即式(5)所得参数组合总数的1%)。

表  2  配准结果对比

Table  2.  Comparison of registration results

		s1-s5	s2-s5	s3-s5	s4-s5	s4-s1	s2-s4	mean
station spacing	reference/m	23	21	21	19	37	40	26.8
	error/m	8.789	6.982	1.952	8.75	0.44	12.378	6.5
our method	RAE/(′)	2.4	2	8.3	1.1	7.2	4.6	4.3
	HTE/mm	3	9	9	11	11	7	8.3
	VTE/mm	29	18	44	3	66	50	35
	time/s	69	74	79	68	71	75	73
reference [11]	RAE/(′)	18	22.5	8.1	6.4	18.2	29.4	17.1
	HTE/mm	20	17	13	10	28	27	19.2
	VTE/mm	151	36	59	40	116	158	93.3
	time/s	18974	57070	12459	6380	827	4231	16657
reference [12]	RAE/(′)	18.2	21.6	8.1	6.7	12.8	31.6	16.5
	HTE/mm	20	19	9	11	24	19	17
	VTE/mm	145	24	57	34	173	42	79
	time/s	58	75	53	47	27	39	50

下载: 导出CSV 

| 显示表格

从配准精度看，参考文献[11]中方法与参考文献[12]中方法配准精度相近，本文中方法的整体配准精度高于三角匹配法，平均旋转角误差、水平和竖直平移误差分别为4.3′、8.3 mm和35 mm, 尤其在测站间距超过35 m时(s₄-s₁和s₂-s₄)，本文中方法的精度仍然保持稳定，而三角匹配法的精度有明显降低。主要原因在于本文中的方法在粗配准的基础上使用测站间距对树干位置点进行了筛选，剔除了测站间距之外的树干位置点。考虑到已有研究的测站间距通常不高于25 m^[6-10]，本文作者通过智能手机辅助森林点云配准的方法提供了一种增大设站间隔、提高林下作业效率的可行思路。此外，本实验中测站间距误差为0.440 m~ 12.378 m，与参考文献[21]和参考文献[22]中得到的消费级定位设备的林下定位误差相符(最大定位误差为7 m~8 m，即对应最大测距误差为14 m~16 m)，表明本文中方法可以在消费级定位设备的误差下为森林点云配准提供有效辅助信息。

从配准效率看，本文中方法的运行时间与参考文献[12]中方法相近且较为稳定，而参考文献[11]中方法的运行时间受单木数量影响较大，主要原因在于本文中的方法通过智能手机辅助信息构建搜索空间，计算量不受单木数量影响，而参考文献[11]中方法需要遍历所有单木以查找最佳匹配三元组，运行效率受单木数量影响较大。

考虑到树干位置点提取存在漏检测和误检测，将一定比例的原始位置点替换为随机误差点，分别使用本文中方法、参考文献[11]中方法和参考文献[12]中方法进行配准，以对比两种方法的稳定性。对每个测站，使用6组替换比例值(5%~30%，间隔为5%)分别随机生成50组树干位置点集，然后选取4对测站s₁-s₅、s₄-s₁、s₂-s₄和s₃-s₅进行配准实验，每种方法共6×50×4=1200组配准实验，结果如图 3所示。实验结果与表 2中的结论相近，在不同的替换比例下, 本文中方法配准精度整体高于参考文献[11]中方法和参考文献[12]中方法(见图 3a~图 3c)。与参考文献[11]中方法相比，在添加随机随机误差后，参考文献[12]中方法的配准精度整体减小(见图 3a、图 3b)，说明参考文献[12]中方法的配准结果受树干位置的漏检或误检影响较大。在图 3d中，本文中将水平平移误差、竖直平移误差和旋转角误差均不高于0.02 m、0.20 m和20′的结果记为正确配准，可得不同替换比例下的配准成功率。由结果可知，本文中方法的成功率高于三角匹配法，且不低于80%。

图  3  不同替换比例下的配准结果

Figure  3.  Registration results in terms of different replacement ratio of stem position points

下载: 全尺寸图片 幻灯片

为研究本文中方法在不同精度激光雷达传感器数据条件下的精度，本文作者通过往树干点中添加不同标准差的高斯噪声来模拟更低精度的测量数据进行实验, 结果如图 4所示。可以看到，配准误差随噪声增加而呈现整体增加的趋势，其中水平误差受噪声影响最大；此外，对多数测站对，在标准差不超过0.06 m时，噪声对配准精度影响较小，超过0.06 m后，配准误差增加显著。证明对于一般精度测量型激光扫描仪所获取的数据，本文中方法仍能达到较高的配准精度。

图  4  不同精度下的配准结果

Figure  4.  Registration results under different precision data

下载: 全尺寸图片 幻灯片

图 5以s₁和s₅测站为例，展示了本文中方法不同阶段的单木位置点拼合结果。可以看到, 位置点在配准后可以准确贴合。图 6展示了本文中方法得到的部分树干点轮廓的拼合效果。图中, 绿色为参考点云，蓝色为待配准点云。可以看到，本文中方法可以准确拼合重叠部分较少或无重叠部分的树干点，在测站(s₂-s₄)间距为40 m时可以取得和测站(s₁-s₅)间距为20 m相同的拼合效果。

图  5  不同阶段的位置点拼合效果

a—粗配准b—第1级目标函数筛选后结果c—最终配准结果

Figure  5.  Alignment of stem positions at different stages

a—coarse registration b—result after evaluation of the first-level objective function c—final result

下载: 全尺寸图片 幻灯片

图  6  树干轮廓拼合结果

Figure  6.  Result of alignment of the stem contours

下载: 全尺寸图片 幻灯片

2.3   参数实验

本文中方法共有3个参数需要设置：方位角搜索间距δ_α，测站搜索间距δ_D和参数总数目N_t(由式(5)得到)的固定比例N_r。图 7为不同参数设置对配准成功率的影响。选取4对测站s₁-s₅、s₄-s₁、s₂-s₄和s₃-s₅进行配准实验，每个测站使用第2.2节中替换比例30%的50组单木位置点集，共进行50×4=200组配准试验。由图 7可知，当δ_α和δ_D分别高于4°和1 m，且N_r取值不超过N_t的1%时，配准成功率变化较大，并随着δ_α和δ_D的增大而减小，随N_r的增大而增大。这是因为在越大的搜索间隔(δ_α、δ_D)下，得到定位定向真值的概率越小，而更大比例(N_r)的参数使误差更小，因此配准成功率越大。当δ_α和δ_D分别低于4°和1 m，且N_r取值超过N_t的1%时，3个参数对配准成功率的影响趋于平稳，这是由于配准精度达到一定限值。在第2.2节中，本文作者将3个参数分别设置为这3个临界值，以验证本文中方法的效率。同时可以看出，这3个值的确定是一个相对随机的过程，对于一般的智能手机所获取的定位定向数据(误差在同一数量级)，本文中参数同样具有较大的参考价值。

图  7  不同参数对配准的影响

Figure  7.  Effect of parameters on registration result

下载: 全尺寸图片 幻灯片

3.   结论

本文中使用智能手机的定位定向数据和树干位置点实现森林点云的配准，通过手机定位定向数据得到粗配准参数, 并构建转换参数搜索空间，在搜索空间内对不同转换参数使用两级目标函数进行筛选, 得到最终配准结果。实验表明，智能手机在冠层厚度不大于11.5 m时，可以为森林点云配准提供有效辅助信息，同时本文中方法在测站间距由20 m增加至40 m时仍可以实现树干轮廓的准确拼合。因此，该方法可为森林数据采集时减少额外设备携带、扩展测站间距，进而为提高林下作业效率提供了一种可行思路。此外，对远距离(本文中为40 m)单木位置点的剔除可以有效提高配准精度。在后续工作中，智能手机定位定向数据在更厚冠层下的有效性仍需验证；同时本文中方法只使用了水平间距对树干位置点进行筛选，未来将探索使用多种因子对树干位置点进行筛选的方法。