CN117455963A - Natural forest region foundation airborne laser point cloud registration method - Google Patents

The invention discloses a natural forest region foundation airborne laser point cloud registration method, which comprises the following steps: s1, cutting, denoising and filtering ULS point cloud and TLS point cloud acquired in a measuring area range respectively; s2, converting the ULS coordinate system into a TLS coordinate system; s3, rasterizing and tree high-point detection are respectively carried out on the ULS point cloud and the TLS point cloud; s4, extracting homonymous point pairs of the ULS point cloud and the TLS point cloud through distance sorting and similar distance searching; s5, completing coarse registration of ULS point cloud and TLS point cloud by adopting a singular value decomposition method based on homonymy point pairs; s6, combining the nearest point fine registration method to finish final registration of TLS and ULS forest region point clouds; according to the method, the characteristic that the positions among single woods in a natural forest area are unchanged is utilized, the high points of the tree are detected by utilizing point cloud data of different platforms to approximately replace the positions of the single woods, and the registration accuracy can be effectively improved by extracting homonymous point pairs of a foundation point cloud and an airborne point cloud and completing coarse registration and fine registration methods of ULS point clouds and TLS point clouds based on the homonymous point pairs.

一种天然林区地基机载激光点云配准方法A method for airborne laser point cloud registration in natural forest areas

技术领域Technical field

本发明涉及计算机图形学领域，具体涉及一种天然林区地基机载激光点云配准方法。The invention relates to the field of computer graphics, and in particular to a method for airborne laser point cloud registration in natural forest areas.

背景技术Background technique

激光雷达可以获取目标地物不同尺度的三维结构信息，通过进一步的数据处理可以提取单木结构参数，在森林资源管理方面有无穷潜力。但是单一平台的激光雷达数据存在点云数据缺陷，一定程度上影响林木参数的反演。因此，多平台点云数据的融合显得尤为重要。Lidar can obtain three-dimensional structural information of target features at different scales. Through further data processing, single wood structure parameters can be extracted, which has unlimited potential in forest resource management. However, the lidar data from a single platform has point cloud data defects, which affects the inversion of forest parameters to a certain extent. Therefore, the fusion of multi-platform point cloud data is particularly important.

点云配准是指将两个或者多个坐标系中的大容量三维空间数据点通过旋转、平移转换到同一坐标系下的数学计算过程，其关键是同名特征的获取和坐标转换参数的稳定计。点云配准算法发展以来已经出现了许多算法和变体算法，根据侧重点不同，基本可以将其划分为基于局部特征描述的算法、基于全局搜索策略的算法和基于统计概率的算法。基于局部特征描述算法是通过邻域信息对特征点进行特征描述，并根据局部特征搜索点云对应点完成配准。基于全局搜索策略的算法以迭代最近点算法(Iterative Closest Point,ICP)算法最为经典，ICP算法通过在全局范围内搜索最近点建立对应点对关系，迭代计算直到残差平方最小，得到变换参数，这种方法计算时间成本较高，而且要求点云具有较好的初始位。基于统计概率的算法中主要有正态变换分布算法(Normal DistributionTransform,NDT)、一致性点漂移算法(Coherent Point Drift,CPD)。NDT算法构建参照数据多维变量的正态分布，求解待变换点在参考坐标系中的概率密度之和，计算其变换参数实现配准。CPD算法将配准问题看作高斯混合模型概率密度函数估计问题，通过期望最大值算法求解最大似然估计。这些方法在许多场景比如道路、建筑都得到了广泛应用，但是天然森林样地点云分布不规则，不同平台点云密度差异大，扫描视场不同，这对配准方法提出了新要求。Point cloud registration refers to the mathematical calculation process of converting large-capacity three-dimensional space data points in two or more coordinate systems to the same coordinate system through rotation and translation. The key is the acquisition of features with the same name and the stability of coordinate conversion parameters. count. Since the development of point cloud registration algorithms, many algorithms and variant algorithms have emerged. According to different focuses, they can be basically divided into algorithms based on local feature description, algorithms based on global search strategies, and algorithms based on statistical probability. The algorithm based on local feature description describes feature points through neighborhood information, and searches for corresponding points in the point cloud based on local features to complete registration. The algorithm based on the global search strategy is the most classic algorithm, the Iterative Closest Point (ICP) algorithm. The ICP algorithm establishes corresponding point pair relationships by searching for the closest points in the global scope, and iteratively calculates until the square of the residual is minimized to obtain the transformation parameters. This method has a high computational time cost and requires a good initial position of the point cloud. Algorithms based on statistical probability mainly include Normal Distribution Transform (NDT) and Coherent Point Drift (CPD). The NDT algorithm constructs the normal distribution of the multi-dimensional variables of the reference data, solves the sum of probability densities of the points to be transformed in the reference coordinate system, and calculates its transformation parameters to achieve registration. The CPD algorithm regards the registration problem as a Gaussian mixture model probability density function estimation problem, and solves the maximum likelihood estimation through the expected maximum algorithm. These methods have been widely used in many scenes such as roads and buildings. However, the distribution of point clouds in natural forest samples is irregular, the density of point clouds on different platforms varies greatly, and the scanning fields of view are different, which puts forward new requirements for registration methods.

为解决天然林地的点云配准问题，许多学者针对天然林样地的特征提出了一些配准方法。Paris等将ALS点云分割林木冠层，将冠层边界在TLS中进行分割，根据冠层位置进行配准，但是当林分密度较高时，冠层分割不确定性增强，配准方法适用性降低。Dai等通过均值漂移算法实现机载平台和地基平台点云数据的关键点提取，结合ICP算法和CPD算法完成配准，平均精度在0.05m-0.08m。Polewski通过机载点云和地基点云分别检测单木位置，将单木垂直和水平距离作为约束条件，通过影像匹配完成配准，但是在地基点云中利用胸径中心作为单木位置具有一定的偏差，部分特殊树种并不适用。Liu等通过栅格化后的CHM树高点，以树高值和单木距样地中心距离作为约束条件寻找同名点对，结合最邻近迭代算法完成点云配准，配准精度在0.18m-0.69m之间。以上研究方法证明，在林区机载点云和地基点云配准中，单木位置是值得利用的刚性特征点，但是由于实际天然林条件的不同和具体树种的区别，各种方法的适用性仍然值得讨论，配准精度也仍然有提升空间。In order to solve the point cloud registration problem of natural forest land, many scholars have proposed some registration methods based on the characteristics of natural forest plots. Paris et al. segmented the forest canopy from ALS point cloud, segmented the canopy boundary in TLS, and performed registration according to the canopy position. However, when the forest stand density is high, the uncertainty of canopy segmentation increases, and the registration method is suitable. Sexuality is reduced. Dai et al. used the mean drift algorithm to extract key points from the point cloud data of the airborne platform and the ground-based platform, and combined the ICP algorithm and the CPD algorithm to complete the registration, with an average accuracy of 0.05m-0.08m. Polewski detects the position of single trees through airborne point clouds and ground point clouds respectively, uses the vertical and horizontal distances of single trees as constraints, and completes registration through image matching. However, using the center of the diameter at breast height as the position of single trees in the ground point cloud has certain limitations. Deviation, some special tree species are not applicable. Liu et al. used the rasterized CHM tree height points to find point pairs with the same name using the tree height value and the distance between a single tree and the plot center as constraints, and combined the nearest neighbor iterative algorithm to complete point cloud registration, with a registration accuracy of 0.18m. -0.69m. The above research methods prove that in the registration of airborne point clouds and ground point clouds in forest areas, the position of a single tree is a rigid feature point worthy of use. However, due to the differences in actual natural forest conditions and differences in specific tree species, the applicability of various methods The accuracy is still worthy of discussion, and there is still room for improvement in registration accuracy.

发明内容Contents of the invention

本发明提供了一种天然林区地基机载激光点云配准方法，该天然林区地基机载激光点云配准方法利用天然林区单木间位置不变的特点，利用不同平台的点云数据检测树高点近似代替单木位置，通过提取地基点云和机载点云的同名点对，并基于同名点对完成ULS点云和TLS点云的粗配准和精配准方法，可有效提高配准精度。The invention provides an airborne laser point cloud registration method for a natural forest area foundation. The airborne laser point cloud registration method for a natural forest area foundation utilizes the characteristics of the constant position between individual trees in the natural forest area and uses points from different platforms. Cloud data detects tree height points to approximately replace the position of individual trees. By extracting point pairs of the same name from the ground point cloud and the airborne point cloud, and completing the coarse registration and fine registration methods of ULS point cloud and TLS point cloud based on the same name point pair, It can effectively improve the registration accuracy.

一种天然林区地基机载激光点云配准方法，包括以下步骤：A method for airborne laser point cloud registration in natural forest areas, including the following steps:

S1、分别对测区范围获取的ULS点云和TLS点云进行裁剪、去噪和滤波；S1. Cut, denoise and filter the ULS point cloud and TLS point cloud obtained within the measurement area;

S2、将ULS坐标系转换到TLS坐标系下；S2. Convert the ULS coordinate system to the TLS coordinate system;

S3、分别对ULS点云和TLS点云进行栅格化与树高点检测；S3. Perform rasterization and tree height detection on ULS point cloud and TLS point cloud respectively;

S4、通过距离排序和相似距离搜索提取ULS点云和TLS点云的同名点对；S4. Extract the same-name point pairs of ULS point cloud and TLS point cloud through distance sorting and similar distance search;

S5、基于同名点对采用奇异值分解方法完成ULS点云和TLS点云的粗配准；S5. Use singular value decomposition method to complete the rough registration of ULS point cloud and TLS point cloud based on the same-name point pairs;

S6、结合最邻近点精配准方法，完成TLS和ULS林区点云的最终配准。S6. Combined with the nearest neighbor point precise registration method, complete the final registration of TLS and ULS forest area point clouds.

优选地，步骤S1中裁剪过程为：根据样地中心坐标将ULS点云数据裁剪为半径25m的圆形样地，将TLS点云数据裁剪为半径20m的圆形样地。Preferably, the clipping process in step S1 is: clip the ULS point cloud data into a circular sample plot with a radius of 25m according to the sample plot center coordinates, and clip the TLS point cloud data into a circular sample plot with a radius of 20m.

优选地，步骤S1中去噪过程为：分别对ULS点云数据集和TLS点云数据集中的每个点进行K邻域统计分析，计算每个点到其K个邻近点的平均距离，假设结果服从高斯分布，将平均距离大于阈值之外的点视作噪声点，保留非噪声点。Preferably, the denoising process in step S1 is: perform K neighborhood statistical analysis on each point in the ULS point cloud data set and TLS point cloud data set, and calculate the average distance from each point to its K neighboring points. Assume The results obey Gaussian distribution, and points whose average distance is greater than the threshold are regarded as noise points, and non-noise points are retained.

优选地，步骤S1中滤波过程为：通过地面点滤波分别从ULS点云数据集和TLS点云数据集中识别出地面点和非地面点，利用地面点将点云数据进行高度归一化，其中，地面点滤波的具体过程为：Preferably, the filtering process in step S1 is: identifying ground points and non-ground points from the ULS point cloud data set and TLS point cloud data set respectively through ground point filtering, and using the ground points to highly normalize the point cloud data, where , the specific process of ground point filtering is:

S11、将点云倒置；S11. Invert the point cloud;

S12、设置模拟布料，设置布料网格分辨率，确定模拟粒子数，布料的位置设置在点云最高点以上；S12. Set the simulated cloth, set the cloth grid resolution, determine the number of simulated particles, and set the position of the cloth above the highest point of the point cloud;

S13、将布料模拟点和点云投影到水平面，为每个布料模拟点找到最相近的点云的高度值；S13. Project the cloth simulation points and point clouds to the horizontal plane, and find the height value of the closest point cloud for each cloth simulation point;

S14、将布料粒子设置为可移动，粒子首先受到重力作用，当粒子高度小于最邻近点云高度时，将邻近点云高度赋予粒子高度，粒子设置为不可移动；S14. Set the cloth particles to be movable. The particles are first affected by gravity. When the particle height is less than the height of the nearest point cloud, the height of the neighboring point cloud is assigned to the particle height, and the particles are set to be immovable;

S15、计算布料粒子之间的内力作用，根据设置的布料刚性参数，调整布料粒子之间的相对位置；S15. Calculate the internal force between cloth particles, and adjust the relative positions between cloth particles according to the set cloth rigid parameters;

S16、设置迭代次数，重复步骤S14和步骤S15的计算；S16. Set the number of iterations and repeat the calculations of steps S14 and S15;

S17、计算激光点云和对应布料模拟点的距离，小于阈值则标记为地面点，大于阈值则标记为非地面点。S17. Calculate the distance between the laser point cloud and the corresponding cloth simulation point. If it is less than the threshold, it will be marked as a ground point, and if it is greater than the threshold, it will be marked as a non-ground point.

优选地，步骤S2的具体过程为：通过步骤S17中提取出地面点云，计算地面点云的中心值，将非地面点云减去地面点云中心值，更新转换后的机载点云坐标，完成ULS点云数据坐标转换，且两个平台点云数据范围具有重叠度，具体转换公式如下：Preferably, the specific process of step S2 is: extract the ground point cloud in step S17, calculate the center value of the ground point cloud, subtract the center value of the ground point cloud from the non-ground point cloud, and update the converted airborne point cloud coordinates. , complete the ULS point cloud data coordinate conversion, and the point cloud data ranges of the two platforms overlap. The specific conversion formula is as follows:

ULS(x,y,z)_i＝ULS_nonground(x,y,z)_i-mean(ULS_ground(x,y,z)),i＝1,2,…,n (1)ULS(x,y,z) _i =ULS _nonground (x,y,z) _i -mean(ULS _ground (x,y,z)),i＝1,2,…,n (1)

其中，ULS(x,y,z)_i表示坐标转换后机载点云非地面点，ULS_nonground(x,y,z)_i表示坐标转换前第i个机载点云非地面点，mean(ULS_ground(x,y,z))表示坐标转换前机载点云地面中心点。Among them, ULS(x,y,z) _i represents the non-ground point of the airborne point cloud after coordinate conversion, ULS _nonground (x,y,z) _i represents the i-th non-ground point of the airborne point cloud before coordinate conversion, mean( ULS _ground (x,y,z)) represents the ground center point of the airborne point cloud before coordinate conversion.

优选地，步骤S3中栅格化与树高点检测的具体过程为：Preferably, the specific process of rasterization and tree height detection in step S3 is:

S31、选定栅格影像分辨率；S31. Select the raster image resolution;

S32、划分格网，建立格网与点云的索引关系；S32. Divide the grid and establish the index relationship between the grid and the point cloud;

S33、确定插值方式，向格网中填充属性值；S33. Determine the interpolation method and fill the attribute values into the grid;

S34、翻转格网；S34, flip the grid;

S35、为每个格网交叉点补充坐标，生成分辨率为0.1m的栅格影像；S35. Supplement coordinates for each grid intersection and generate a raster image with a resolution of 0.1m;

S36、通过3×3大小的滑动窗口对栅格影像进行最大值检测，所检测到的点即为树高点；S36. Perform maximum value detection on the raster image through a 3×3 sliding window, and the detected point is the tree height point;

其中，ULS点云经栅格化与树高点检测得到树高点TH_ULS，TLS点云经栅格化与树高点检测得到树高点TH_TLS，树高点TH_ULS和树高点TH_TLS均为单木位置点。Among them, the ULS point cloud is rasterized and the tree height point is detected to obtain the tree height point TH _ULS , the TLS point cloud is rasterized and the tree height point is detected to obtain the tree height point TH _TLS , the tree height point TH _ULS and the tree height point TH _TLS are all single wood position points.

优选地，步骤S4的具体过程为：Preferably, the specific process of step S4 is:

S41、分别计算树高点TH_ULS和树高点TH_TLS至样地中心距离Dcenter，分别排序成为DcenterULS列表和DcenterTLS列表；S41. Calculate the distance Dcenter from the tree height point TH _ULS and the tree height point TH _TLS to the plot center respectively, and sort them into a DcenterULS list and a DcenterTLS list respectively;

S42、分别从列表中抽取距离样地中心最近树高点，比较距离，遍历所有树高点，距离相似点作为待选点对；S42. Extract the nearest tree height points from the center of the sample plot from the list, compare the distances, traverse all tree height points, and use points with similar distances as candidate point pairs;

S43、从待选点对中分别搜索最邻近树高点，判断最邻近距离相似性；S43. Search the nearest neighbor tree height points from the pairs of points to be selected, and determine the nearest neighbor distance similarity;

S44、不相似，跳回步骤S43；相似，计算最邻近树高点至样地中心距离，并判断该距离相似性；S44. Not similar, jump back to step S43; similar, calculate the distance from the nearest tree height point to the center of the sample plot, and determine the similarity of the distance;

S45、不相似，跳回步骤S43；相似，保留该树高点对，记为同名点对，存储在matchpoint集合中，直至遍历完所有树高点，输出match point集合。S45. If they are not similar, jump back to step S43; if they are similar, keep the tree high point pair, record it as a point pair with the same name, and store it in the matchpoint set until all tree high points are traversed and the match point set is output.

优选地，步骤S5通过步骤S4提取出的同名点对集合，利用奇异值分解得到旋转矩阵R和平移向量t，实现林区ULS点云和TLS点云的粗配准融合，粗配准的具体过程为：假设Q＝{q1,q2,…,qn}和P＝{p1,p2,…,pn}是两组d维空间中的对应点集，Q表示TLS树高点集，P表示ULS树高点集，根据TLS树高点集和ULS树高点集计算出它们之间的刚性转换信息，利用最小二乘法求解最优解，计算公式为：Preferably, step S5 uses the set of point pairs with the same name extracted in step S4, and uses singular value decomposition to obtain the rotation matrix R and the translation vector t, so as to realize the coarse registration fusion of the ULS point cloud and the TLS point cloud in the forest area. The details of the coarse registration are: The process is: Suppose Q={q1,q2,…,qn} and P={p1,p2,…,pn} are corresponding point sets in two sets of d-dimensional spaces, Q represents the TLS tree high point set, and P represents ULS Tree height point set, calculate the rigid transformation information between them based on the TLS tree high point set and ULS tree high point set, and use the least squares method to find the optimal solution. The calculation formula is:

其中，R表示旋转矩阵，t表示平移向量，w_i是点集中第i个点的权重，w_i>0，p_i表示机载点云中提取的第i个树高点，q_i表示地基点云中提取的第i个树高点；Among them, R represents the rotation matrix, t represents the translation vector, w _i is the weight of the i-th point in the point set, w _i >0, p _i represents the i-th tree height point extracted from the airborne point cloud, and q _i represents the ground point. The i-th tree height point extracted from the base point cloud;

求解最优旋转矩阵和平移向量的具体过程为：The specific process of solving the optimal rotation matrix and translation vector is:

S51、计算TLS树高点集的加权质心和ULS树高点集的加权质心/> S51. Calculate the weighted centroid of the TLS tree high point set and the weighted centroid of the ULS tree high point set/>

S52、计算中心向量：S52. Calculate the center vector:

S53、计算d×d协方差矩阵：S53. Calculate the d×d covariance matrix:

S＝XWY^T (5)S＝XWY ^T (5)

其中，X和Y分别是具有x_i列和y_i列的d×n矩阵，where, X and Y are d×n matrices with columns xi _{and yi} respectively,

W＝diag(w₁,w₂,...,w_n)；W=diag(w ₁ ,w ₂ ,...,w _n );

S54、计算奇异值分解S＝UΣV^T，得到旋转矩阵R：S54. Calculate the singular value decomposition S=UΣV ^T and obtain the rotation matrix R:

S55、计算最优平移向量t：S55. Calculate the optimal translation vector t:

计算得到的旋转矩阵R和平移向量t代入待配准点集即ULS点云数据集中，计算得到配准后的机载点云数据集ULS_tr：The calculated rotation matrix R and translation vector t are substituted into the point set to be registered, that is, the ULS point cloud data set, and the registered airborne point cloud data set ULS _tr is calculated:

ULS_tr＝p_i×R+t (8)。ULS _tr = _pi ×R+t (8).

优选地，步骤S6的具体过程为：Preferably, the specific process of step S6 is:

S61、搜索粗配准后的机载点云数据集ULS_tr上的每一点在地基点云数据集中的空间最邻近点；S61. Search for the spatial nearest point in the ground-based point cloud data set of each point on the coarsely registered airborne point cloud data set ULS _tr ;

S62、对搜索到的最邻近点对进行奇异值分解，求解旋转矩阵和平移向量，并应用到待精配准的机载点云中；S62. Perform singular value decomposition on the searched nearest point pairs, solve the rotation matrix and translation vector, and apply them to the airborne point cloud to be precisely registered;

S63、进行误差迭代计算，直至求解得到使最近点之间的距离均方差最小的旋转矩阵和平移向量；S63. Perform error iterative calculation until the rotation matrix and translation vector that minimize the mean square error of the distance between the closest points are obtained;

S64、将经过最邻近点精配准后的旋转矩阵R和平移向量t，代入到公式(4)中，得到最终配准后的机载点云，并将经过精配准的机载点云和地基点云进行合并，得到完整的融合TLS和ULS的融合点云。S64. Substitute the rotation matrix R and the translation vector t after precise registration of the nearest neighbor points into formula (4) to obtain the final registered airborne point cloud, and use the precisely registered airborne point cloud to Merge with the ground point cloud to obtain a complete fusion point cloud that combines TLS and ULS.

采用上述技术方案后，本发明具有如下有益效果：本发明采用树高点相似距离搜索的林区点云配准算法，基于单木位置刚性不变的特征，利用林区点云数据提取树高点替代单木位置，通过单木相互位置及单木与样地中心距离比对搜索同名点对完成粗配准，利用最邻近迭代完成精配准，本申请所提出的方法对于四个树种的天然林样地均具有适用性，相比部分已有方法精度有所提升，其配准流程简洁明了。After adopting the above technical solution, the present invention has the following beneficial effects: The present invention adopts a forest area point cloud registration algorithm based on similar distance search of tree height points. Based on the characteristics of single tree position rigidity, the invention uses forest area point cloud data to extract tree height. Points replace the position of a single tree. By comparing the mutual positions of single trees and the distance between the single tree and the center of the sample plot, the same-name point pairs are searched to complete the rough registration. The nearest neighbor iteration is used to complete the fine registration. The method proposed in this application is suitable for four tree species. It is applicable to all natural forest sample plots, has improved accuracy compared with some existing methods, and its registration process is simple and clear.

附图说明Description of the drawings

图1为本发明的ULS点云和TLS点云配准流程图；Figure 1 is a flow chart of ULS point cloud and TLS point cloud registration according to the present invention;

图2为本发明的树高点检测示意图；(a)TLS树高点检测，(b)ULS树高点检测；Figure 2 is a schematic diagram of tree high point detection according to the present invention; (a) TLS tree high point detection, (b) ULS tree high point detection;

图3为本发明的粗配准效果图；(a)样地1粗配准前，(b)样地1粗配准后，(c)样地2粗配准前，(d)样地2粗配准后，(e)样地3粗配准前，(f)样地3粗配准后，(g)样地4粗配准前，(h)样地4粗配准后；Figure 3 is a rough registration effect diagram of the present invention; (a) before rough registration of sample plot 1, (b) after rough registration of sample plot 1, (c) before rough registration of sample plot 2, (d) sample plot 2 after rough registration, (e) before rough registration of sample plot 3, (f) after rough registration of sample plot 3, (g) before rough registration of sample plot 4, (h) after rough registration of sample plot 4;

图4为本发明的样地1高山松精配准全局效果图；(a)、(d):X-Z视角；(b)、(e):Y-Z视角；(c)、(f):X-Y视角；Figure 4 is a global rendering of alpine pine precision registration in sample plot 1 of the present invention; (a), (d): X-Z perspective; (b), (e): Y-Z perspective; (c), (f): X-Y perspective ;

图5为本发明的样地2云杉精配准全局效果图；(a)、(d):X-Z视角；(b)、(e):Y-Z视角；(c)、(f):X-Y视角；Figure 5 is a global rendering of spruce precision registration in sample plot 2 of the present invention; (a), (d): X-Z perspective; (b), (e): Y-Z perspective; (c), (f): X-Y perspective ;

图6为本发明的样地3云南松精配准全局效果图；(a)、(d):X-Z视角；(b)、(e):Y-Z视角；(c)、(f):X-Y视角；Figure 6 is a global rendering of Yunnan pine fine registration in sample plot 3 of the present invention; (a), (d): X-Z perspective; (b), (e): Y-Z perspective; (c), (f): X-Y perspective ;

图7为本发明的样地4冷杉精配准全局效果图；(a)、(d):X-Z视角；(b)、(e):Y-Z视角；(c)、(f):X-Y视角；Figure 7 is a global effect diagram of precise registration of fir trees in sample plot 4 of the present invention; (a), (d): X-Z perspective; (b), (e): Y-Z perspective; (c), (f): X-Y perspective;

图8为本发明的样地1高山松精配准效果局部图；(a)：TLS点云；(b):ULS点云；(c)：配准融合后点云X-Y-Z视角；(d):配准融合后点云X-Z视角；(e)：配准融合后点云Y-Z视角；Figure 8 is a partial view of the alpine pine precision registration effect in sample plot 1 of the present invention; (a): TLS point cloud; (b): ULS point cloud; (c): X-Y-Z perspective of point cloud after registration and fusion; (d) : X-Z perspective of the point cloud after registration and fusion; (e): Y-Z perspective of the point cloud after registration and fusion;

图9为本发明的样地2云杉精配准效果局部图；(a)：TLS点云；(b):ULS点云；(c)：配准融合后点云X-Y-Z视角；(d):配准融合后点云X-Z视角；(e)：配准融合后点云Y-Z视角；Figure 9 is a partial view of the precise registration effect of spruce in sample plot 2 of the present invention; (a): TLS point cloud; (b): ULS point cloud; (c): X-Y-Z perspective of point cloud after registration and fusion; (d) : X-Z perspective of the point cloud after registration and fusion; (e): Y-Z perspective of the point cloud after registration and fusion;

图10为本发明的样地3云南松精配准效果局部图；(a)：TLS点云；(b):ULS点云；(c)：配准融合后点云X-Y-Z视角；(d):配准融合后点云X-Z视角；(e)：配准融合后点云Y-Z视角；Figure 10 is a partial view of the registration effect of Yunnan pines in sample plot 3 of the present invention; (a): TLS point cloud; (b): ULS point cloud; (c): X-Y-Z perspective of point cloud after registration and fusion; (d) : X-Z perspective of the point cloud after registration and fusion; (e): Y-Z perspective of the point cloud after registration and fusion;

图11为本发明的样地4冷杉精配准效果局部图；(a)：TLS点云；(b):ULS点云；(c)：配准融合后点云X-Y-Z视角；(d):配准融合后点云X-Z视角；(e)：配准融合后点云Y-Z视角；Figure 11 is a partial view of the precise registration effect of fir trees in sample plot 4 of the present invention; (a): TLS point cloud; (b): ULS point cloud; (c): X-Y-Z perspective of point cloud after registration and fusion; (d): X-Z perspective of the point cloud after registration and fusion; (e): Y-Z perspective of the point cloud after registration and fusion;

图12为本发明的样地1ULS-TLS最邻近点距离统计图；Figure 12 is a statistical diagram of the nearest neighbor point distance of sample plot 1ULS-TLS of the present invention;

图13为本发明的样地2ULS-TLS最邻近点距离统计图；Figure 13 is a statistical diagram of the nearest neighbor point distance of sample plot 2ULS-TLS of the present invention;

图14为本发明的样地3ULS-TLS最邻近点距离统计图；Figure 14 is a statistical diagram of the nearest neighbor point distance of sample plot 3ULS-TLS of the present invention;

图15为本发明的样地4ULS-TLS最邻近点距离统计图。Figure 15 is a statistical diagram of the nearest neighbor point distance of the sample plot 4ULS-TLS of the present invention.

具体实施方式Detailed ways

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。In order to make the purpose, technical solutions and advantages of the present invention more clear, the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described here are only used to explain the present invention and are not intended to limit the present invention.

参见图1至图15；一种天然林区地基机载激光点云配准方法，包括以下步骤：See Figures 1 to 15; a method for airborne laser point cloud registration in natural forest areas, including the following steps:

S1、分别对测区范围获取的ULS点云和TLS点云进行裁剪、去噪和滤波；S1. Cut, denoise and filter the ULS point cloud and TLS point cloud obtained within the measurement area;

步骤S1中裁剪过程为：根据样地中心坐标将ULS点云数据裁剪为半径25m的圆形样地，将TLS点云数据裁剪为半径20m的圆形样地；The clipping process in step S1 is: clip the ULS point cloud data into a circular sample plot with a radius of 25m according to the sample plot center coordinates, and clip the TLS point cloud data into a circular sample plot with a radius of 20m;

步骤S1中去噪过程为：分别对ULS点云数据集和TLS点云数据集中的每个点进行K邻域统计分析，计算每个点到其K个邻近点的平均距离，假设结果服从高斯分布，将平均距离大于阈值之外的点视作噪声点，保留非噪声点；The denoising process in step S1 is: perform K neighborhood statistical analysis on each point in the ULS point cloud data set and TLS point cloud data set, calculate the average distance from each point to its K neighboring points, assuming that the result obeys Gaussian Distribution, points whose average distance is greater than the threshold are regarded as noise points, and non-noise points are retained;

步骤S1中滤波过程为：通过地面点滤波分别从ULS点云数据集和TLS点云数据集中识别出地面点和非地面点，利用地面点将点云数据进行高度归一化，其中，地面点滤波的具体过程为：The filtering process in step S1 is: identify ground points and non-ground points from the ULS point cloud data set and TLS point cloud data set respectively through ground point filtering, and use the ground points to highly normalize the point cloud data. Among them, the ground points The specific process of filtering is:

S11、将点云倒置；S11. Invert the point cloud;

S12、设置模拟布料，设置布料网格分辨率，确定模拟粒子数，布料的位置设置在点云最高点以上；S12. Set the simulated cloth, set the cloth grid resolution, determine the number of simulated particles, and set the position of the cloth above the highest point of the point cloud;

S13、将布料模拟点和点云投影到水平面，为每个布料模拟点找到最相近的点云的高度值；S13. Project the cloth simulation points and point clouds to the horizontal plane, and find the height value of the closest point cloud for each cloth simulation point;

S14、将布料粒子设置为可移动，粒子首先受到重力作用，当粒子高度小于最邻近点云高度时，将邻近点云高度赋予粒子高度，粒子设置为不可移动；S14. Set the cloth particles to be movable. The particles are first affected by gravity. When the particle height is less than the height of the nearest point cloud, the height of the neighboring point cloud is assigned to the particle height, and the particles are set to be immovable;

S15、计算布料粒子之间的内力作用，根据设置的布料刚性参数，调整布料粒子之间的相对位置；S15. Calculate the internal force between cloth particles, and adjust the relative positions between cloth particles according to the set cloth rigid parameters;

S16、设置迭代次数，重复步骤S14和步骤S15的计算；S16. Set the number of iterations and repeat the calculations of steps S14 and S15;

S17、计算激光点云和对应布料模拟点的距离，小于阈值则标记为地面点，大于阈值则标记为非地面点；S17. Calculate the distance between the laser point cloud and the corresponding cloth simulation point. If it is less than the threshold, it will be marked as a ground point, and if it is greater than the threshold, it will be marked as a non-ground point;

S2、将ULS坐标系转换到TLS坐标系下；S2. Convert the ULS coordinate system to the TLS coordinate system;

步骤S2的具体过程为：通过步骤S17中提取出地面点云，计算地面点云的中心值，将非地面点云减去地面点云中心值，更新转换后的机载点云坐标，完成ULS点云数据坐标转换，且两个平台点云数据范围具有重叠度，具体转换公式如下：The specific process of step S2 is: extract the ground point cloud in step S17, calculate the center value of the ground point cloud, subtract the center value of the ground point cloud from the non-ground point cloud, update the converted airborne point cloud coordinates, and complete ULS. Point cloud data coordinate conversion, and the point cloud data ranges of the two platforms overlap. The specific conversion formula is as follows:

ULS(x,y,z)_i＝ULS_nonground(x,y,z)_i-mean(ULS_ground(x,y,z)),i＝1,2,…,n(1)ULS(x,y,z) _i =ULS _nonground (x,y,z) _i -mean(ULS _ground (x,y,z)),i＝1,2,…,n(1)

其中，ULS(x,y,z)_i表示坐标转换后机载点云非地面点，ULS_nonground(x,y,z)_i表示坐标转换前第i个机载点云非地面点，mean(ULS_ground(x,y,z))表示坐标转换前机载点云地面中心点；S3、分别对ULS点云和TLS点云进行栅格化与树高点检测；Among them, ULS(x,y,z) _i represents the non-ground point of the airborne point cloud after coordinate conversion, ULS _nonground (x,y,z) _i represents the i-th non-ground point of the airborne point cloud before coordinate conversion, mean( ULS _ground (x, y, z)) represents the ground center point of the airborne point cloud before coordinate conversion; S3, perform rasterization and tree height point detection on the ULS point cloud and TLS point cloud respectively;

步骤S3中栅格化与树高点检测的具体过程为：The specific process of rasterization and tree height detection in step S3 is:

S31、选定栅格影像分辨率；S31. Select the raster image resolution;

S32、划分格网，建立格网与点云的索引关系；S32. Divide the grid and establish the index relationship between the grid and the point cloud;

S33、确定插值方式，向格网中填充属性值；S33. Determine the interpolation method and fill the attribute values into the grid;

S34、翻转格网；S34, flip the grid;

S35、为每个格网交叉点补充坐标，生成分辨率为0.1m的栅格影像；S35. Supplement coordinates for each grid intersection and generate a raster image with a resolution of 0.1m;

S36、通过3×3大小的滑动窗口对栅格影像进行最大值检测，所检测到的点即为树高点；S36. Perform maximum value detection on the raster image through a 3×3 sliding window, and the detected point is the tree height point;

其中，ULS点云经栅格化与树高点检测得到树高点TH_ULS，TLS点云经栅格化与树高点检测得到树高点TH_TLS，树高点TH_ULS和树高点TH_TLS均为单木位置点；Among them, the ULS point cloud is rasterized and the tree height point is detected to obtain the tree height point TH _ULS , the TLS point cloud is rasterized and the tree height point is detected to obtain the tree height point TH _TLS , the tree height point TH _ULS and the tree height point TH _TLS are all single wood position points;

S4、通过距离排序和相似距离搜索提取ULS点云和TLS点云的同名点对；步骤S4的具体过程为：S4. Extract the same-name point pairs of ULS point cloud and TLS point cloud through distance sorting and similar distance search; the specific process of step S4 is:

S41、分别计算树高点TH_ULS和树高点TH_TLS至样地中心距离Dcenter，分别排序成为DcenterULS列表和DcenterTLS列表；S41. Calculate the distance Dcenter from the tree height point TH _ULS and the tree height point TH _TLS to the plot center respectively, and sort them into a DcenterULS list and a DcenterTLS list respectively;

S42、分别从列表中抽取距离样地中心最近树高点，比较距离，遍历所有树高点，距离相似点作为待选点对；S42. Extract the nearest tree height points from the center of the sample plot from the list, compare the distances, traverse all tree height points, and use points with similar distances as candidate point pairs;

S43、从待选点对中分别搜索最邻近树高点，判断最邻近距离相似性；S43. Search the nearest neighbor tree height points from the pairs of points to be selected, and determine the nearest neighbor distance similarity;

S44、不相似，跳回步骤S43；相似，计算最邻近树高点至样地中心距离，并判断该距离相似性；S44. Not similar, jump back to step S43; similar, calculate the distance from the nearest tree height point to the center of the sample plot, and determine the similarity of the distance;

S45、不相似，跳回步骤S43；相似，保留该树高点对，记为同名点对，存储在matchpoint集合中，直至遍历完所有树高点，输出match point集合；S45. If they are not similar, jump back to step S43; if they are similar, retain the tree high point pair, record it as a point pair with the same name, and store it in the matchpoint set until all tree high points are traversed and the match point set is output;

S5、基于同名点对采用奇异值分解方法完成ULS点云和TLS点云的粗配准；步骤S5通过步骤S4提取出的同名点对集合，利用奇异值分解得到旋转矩阵R和平移向量t，实现林区ULS点云和TLS点云的粗配准融合，粗配准的具体过程为：假设Q＝{q1,q2,…,qn}和P＝{p1,p2,…,pn}是两组d维空间中的对应点集，Q表示TLS树高点集，P表示ULS树高点集，根据TLS树高点集和ULS树高点集计算出它们之间的刚性转换信息，利用最小二乘法求解最优解，计算公式为：S5. Use the singular value decomposition method to complete the rough registration of the ULS point cloud and the TLS point cloud based on the point pairs with the same name. Step S5 uses the set of point pairs with the same name extracted in step S4 to obtain the rotation matrix R and the translation vector t using singular value decomposition. Realize the coarse registration and fusion of ULS point cloud and TLS point cloud in the forest area. The specific process of coarse registration is: assuming that Q = {q1, q2,..., qn} and P = {p1, p2,..., pn} are two Set the corresponding point set in the d-dimensional space, Q represents the TLS tree high point set, P represents the ULS tree high point set, calculate the rigid transformation information between them according to the TLS tree high point set and the ULS tree high point set, and use the minimum The square method is used to find the optimal solution, and the calculation formula is:

其中，R表示旋转矩阵，t表示平移向量，w_i是点集中第i个点的权重，w_i>0，p_i表示机载点云中提取的第i个树高点，q_i表示地基点云中提取的第i个树高点；Among them, R represents the rotation matrix, t represents the translation vector, w _i is the weight of the i-th point in the point set, w _i >0, p _i represents the i-th tree height point extracted from the airborne point cloud, and q _i represents the ground point. The i-th tree height point extracted from the base point cloud;

求解最优旋转矩阵和平移向量的具体过程为：The specific process of solving the optimal rotation matrix and translation vector is:

S51、计算TLS树高点集的加权质心和ULS树高点集的加权质心/> S51. Calculate the weighted centroid of the TLS tree high point set and the weighted centroid of the ULS tree high point set/>

S52、计算中心向量：S52. Calculate the center vector:

S53、计算d×d协方差矩阵：S53. Calculate the d×d covariance matrix:

S＝XWY^T (5)S＝XWY ^T (5)

其中，X和Y分别是具有x_i列和y_i列的d×n矩阵，where, X and Y are d×n matrices with columns xi _{and yi} respectively,

W＝diag(w₁,w₂,...,w_n)；W=diag(w ₁ ,w ₂ ,...,w _n );

S54、计算奇异值分解S＝UΣV^T，得到旋转矩阵R：S54. Calculate the singular value decomposition S=UΣV ^T and obtain the rotation matrix R:

S55、计算最优平移向量t：S55. Calculate the optimal translation vector t:

计算得到的旋转矩阵R和平移向量t代入待配准点集即ULS点云数据集中，计算得到配准后的机载点云数据集ULS_tr：The calculated rotation matrix R and translation vector t are substituted into the point set to be registered, that is, the ULS point cloud data set, and the registered airborne point cloud data set ULS _tr is calculated:

ULS_tr＝p_i×R+t (8)；ULS _tr =p _i ×R+t (8);

S6、结合最邻近点精配准方法，完成TLS和ULS林区点云的最终配准；步骤S6的具体过程为：S6. Combine the nearest neighbor point precise registration method to complete the final registration of TLS and ULS forest area point clouds; the specific process of step S6 is:

S61、搜索粗配准后的机载点云数据集ULS_tr上的每一点在地基点云数据集中的空间最邻近点；S61. Search for the spatial nearest point in the ground-based point cloud data set of each point on the coarsely registered airborne point cloud data set ULS _tr ;

S62、对搜索到的最邻近点对进行奇异值分解，求解旋转矩阵和平移向量，并应用到待精配准的机载点云中；S62. Perform singular value decomposition on the searched nearest neighbor point pairs, solve the rotation matrix and translation vector, and apply them to the airborne point cloud to be precisely registered;

S63、进行误差迭代计算，直至求解得到使最近点之间的距离均方差最小的旋转矩阵和平移向量；S63. Perform error iterative calculation until the rotation matrix and translation vector that minimize the mean square error of the distance between the closest points are obtained;

S64、将经过最邻近点精配准后的旋转矩阵R和平移向量t，代入到公式(4)中，得到最终配准后的机载点云，并将经过精配准的机载点云和地基点云进行合并，得到完整的融合TLS和ULS的融合点云。S64. Substitute the rotation matrix R and the translation vector t after precise registration of the nearest neighbor points into formula (4) to obtain the final registered airborne point cloud, and use the precisely registered airborne point cloud to Merge with the ground point cloud to obtain a complete fusion point cloud that combines TLS and ULS.

配准效果对比：Registration effect comparison:

本发明通过检测树高点代替单木位置，进而搜索相似距离寻找同名特征点的方式进行ULS点云与TLS点云的初始配准即点云粗配准。通过图3(深色点代表ULS树高点，浅色点代表TLS树高点)可以看出，在进行粗配准前的树高点是杂乱无章的，通过肉眼无法判断同名单木点，经过粗配准后，所选的四块样地中的机载点云树高点与地基点云树高点，出现了相似的空间分布形态，能够较为清晰地辨别出两个平台数据的同名树高点。事实上，并不是所有的ULS树高点都有对应的TLS树高点与之相匹配，这是因为机载激光雷达自上而下的扫描方式和地基自下而上的扫描方式，导致树高检测时并不能完全检测出对应树高。机载点云范围要大于地基点云范围，因此，机载点云的边缘树高点，并未有与之相匹配的地基点云。This invention performs initial registration of ULS point clouds and TLS point clouds, that is, point cloud coarse registration, by detecting tree height points instead of single tree positions, and then searching for similar distances to find feature points with the same name. It can be seen from Figure 3 (the dark dots represent the high points of the ULS tree, the light dots represent the high points of the TLS tree) that the tree heights before coarse registration are messy, and it is impossible to judge the same single tree point with the naked eye. After rough registration, the high points of the airborne point cloud trees and the high points of the ground point cloud trees in the four selected plots have similar spatial distribution patterns, and the trees with the same name in the two platform data can be clearly distinguished. High Point. In fact, not all ULS tree heights have corresponding TLS tree heights that match them. This is because of the top-down scanning method of airborne lidar and the bottom-up scanning method of the ground. The corresponding tree height cannot be completely detected during high detection. The range of the airborne point cloud is larger than the range of the ground point cloud. Therefore, the edge tree heights of the airborne point cloud do not have matching ground point clouds.

通过近似距离搜索检测特征同名点对，通过奇异值分解实现ULS点云和TLS点云的初始配准，为进一步提高配准精度，又利用最邻近点迭代实现精配准，完成ULS点云和TLS点云的最终配准。最终配准效果如图4-图7所示，其中(a)(b)(c)表示配准前示意图，(d)(e)(f)表示配准后示意图。Detect feature point pairs with the same name through approximate distance search, and achieve initial registration of ULS point cloud and TLS point cloud through singular value decomposition. In order to further improve the registration accuracy, the nearest neighbor point iteration is used to achieve precise registration, and the ULS point cloud and TLS point cloud are completed. Final registration of TLS point clouds. The final registration effect is shown in Figures 4 to 7, where (a) (b) (c) represents the schematic diagram before registration, and (d) (e) (f) represents the schematic diagram after registration.

从图4-图7中可以看出，经过粗配准和精配准后，ULS点云和TLS点云实现了一定程度的重叠融合。为进一步查看配准效果，研究从每个样地中的ULS点云和TLS点云中裁剪提取了单木点云，以查看局部配准效果，如图8-图11所示，从局部单木图可以看出，将ULS点云与TLS点云配准融合后，一定程度上弥补了两个平台数据各自的缺陷。相比较单一TLS点云而言，融合后的点云获得了更为完整的冠层信息，尤其是冠层高度信息，相比较单一ULS点云而言，融合后的点云弥补了缺失的林下信息和单木主干信息。As can be seen from Figures 4 to 7, after coarse registration and fine registration, the ULS point cloud and TLS point cloud have achieved a certain degree of overlap and fusion. In order to further check the registration effect, the study extracted single tree point clouds from the ULS point cloud and TLS point cloud in each plot to check the local registration effect. As shown in Figures 8 to 11, from the local single point cloud It can be seen from the map that after integrating ULS point cloud and TLS point cloud registration, the shortcomings of the data of the two platforms are compensated to a certain extent. Compared with a single TLS point cloud, the fused point cloud obtains more complete canopy information, especially canopy height information. Compared with a single ULS point cloud, the fused point cloud makes up for the missing forest information. Lower information and single tree trunk information.

本发明选择利用最邻近点距离作为配准精度验证的指标，即搜索配准后的ULS点云中与TLS点云的最邻近点，计算其距离之和的平均值p，即为最终的配准精度。其计算公式如公式(9)所示。研究计算并统计了粗配准后和精配准后，ULS点云在TLS点云中的最邻近点距离分布情况，如图12-图15所示，其中(a)表示粗配准后距离分布情况，(b)表示精配准后距离分布情况。可以发现，经过精配准后，最邻近距离的分布趋向于更小距离，精配准前较远距离的点得到改善。This invention chooses to use the nearest neighbor point distance as an indicator for registration accuracy verification, that is, to search for the nearest neighbor point between the registered ULS point cloud and the TLS point cloud, and calculate the average p of the sum of their distances, which is the final registration. Accuracy. Its calculation formula is shown in formula (9). The study calculated and counted the distance distribution of the nearest neighbor points of the ULS point cloud in the TLS point cloud after coarse registration and fine registration, as shown in Figure 12-Figure 15, where (a) represents the distance after coarse registration. Distribution, (b) represents the distance distribution after precise registration. It can be found that after fine registration, the distribution of nearest neighbor distances tends to be smaller, and points farther away are improved before fine registration.

其中，n表示点总个数，d_knn(i)表示第i个点的最邻近距离。Among them, n represents the total number of points, and d _knn(i) represents the nearest neighbor distance of the i-th point.

选择了两种配准方法进行了配准精度比较，以证明本申请所提出配准方法的有效性和可靠性。配准精度表如表1所示。通过与其他方法的对比发现，本发明所提出的搜索相似距离法进行粗配准就取得了不错的效果，在样地1、3、4中，粗配准效果比ICP算法精配准还要高。而每块样地的精配准效果相比TR算法的精配准精度都有所提升。证明本发明提出的配准方法具有较好的有效性和可靠性。Two registration methods were selected and the registration accuracy was compared to prove the effectiveness and reliability of the registration method proposed in this application. The registration accuracy table is shown in Table 1. Through comparison with other methods, it was found that the search similarity distance method proposed by the present invention achieved good results in rough registration. In sample plots 1, 3, and 4, the rough registration effect was better than the fine registration using the ICP algorithm. high. The fine registration effect of each plot is improved compared to the fine registration accuracy of the TR algorithm. It is proved that the registration method proposed by the present invention has better effectiveness and reliability.

表1各配准方法精度对比Table 1 Comparison of accuracy of various registration methods

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求的保护范围为准。The above are only preferred specific embodiments of the present invention, but the protection scope of the present invention is not limited thereto. Any person familiar with the technical field can easily think of changes or modifications within the technical scope disclosed in the present invention. All substitutions are within the scope of the present invention. Therefore, the protection scope of the present invention should be subject to the protection scope of the claims.