Point cloud registration method based on hierarchical grouping strategy

Abstract

The invention relates to a point cloud registration method based on a hierarchical grouping strategy. The method comprises the following steps: acquiring a corresponding point pair set; calculating to obtain a coordinate filling matrix of the (k-1) th group of point pairs; calculating a parameter correction vector of the (k-1) th group; calculating a first matrix and a second matrix of the current number of iterations; calculating a parameter correction vector of the kth group; updating the registration parameter to obtain a current registration parameter of the current iteration number; calculating error vectors of the kth group of point pairs and the (k-1) th group of point pairs; judging whether the parameter correction vector is smaller than a threshold value or not, if yes, stopping iteration, and if not, updating the number of iterations; and calculating an error in the unit weight after verification to carry out precision estimation, and carrying out precise registration on the corresponding point pair set based on the current registration parameter to obtain a point cloud registration result. Compared with the prior art, the method has the advantages of reducing point cloud registration calculation consumption and the like.

Description

Technical Field

[0001] The present invention relates to, and in particular to, a point cloud registration method based on a hierarchical grouping strategy. Background Art

[0002] Due to the high efficiency and high precision of data acquisition, point cloud data obtained by 3D laser scanning technology is widely used in many fields such as autonomous driving, deformation monitoring, industrial inspection, 3D reconstruction, etc. As an important part of point cloud data processing, the quality of point cloud registration directly determines the accuracy of the subsequent 3D model. Therefore, accurate registration is the key to the subsequent application of point cloud data. The essential work of registration is to unify the coordinate system of the point cloud source data set and the target data set. According to the degree of registration precision, it can be divided into coarse registration and fine registration. Coarse registration roughly aligns the point cloud through the registration reference points provided by the target, and then uses the fine registration method to accurately align the point cloud. The research object of this paper is the fine registration method of point cloud.

[0003] In the academic research of point cloud precision registration, the iterative closest point (ICP) method proposed by Besl and Mckay is undoubtedly a milestone classic method. This method is based on minimizing the sum of squares of distances between corresponding point pairs as the objective function, and continuously optimizes the position of sampling points in space through an iterative process to achieve accurate alignment between the point cloud to be registered and the reference point cloud. ICP has been widely used in many fields such as three-dimensional reconstruction, robot positioning, and medical imaging due to its intuitive objective function and high efficiency of convergence speed. However, in-depth exploration in the academic community has revealed some inherent limitations of the ICP method. First, the efficiency of the nearest point search has become one of the key factors restricting the performance of ICP, especially when processing large-scale or high-density point cloud data. This challenge is particularly prominent. Secondly, ICP is sensitive to the choice of initial position. Inappropriate initial position may cause the algorithm to fall into local optimality, thus affecting the accuracy and stability of registration.

[0004] In addition to the classic ICP method, the least square matching technology (LSM) is also a powerful and efficient tool that has shown excellent performance in the field of point cloud precision registration. This technology is rooted in the image stereo matching theory in photogrammetry. After continuous evolution and expansion, its application scope has far exceeded the initial two-dimensional image processing. It is particularly worth mentioning that the pioneering work of Gruen et al. (2005), who proposed the LS3D method, successfully extended the application scope of LSM technology from the two-dimensional image field to the fine processing of three-dimensional point cloud data, greatly enriching the means of three-dimensional spatial data processing. Subsequently, Akca (2010) made in-depth improvements on the basis of the LS3D method, aiming to further improve the computational efficiency and accuracy of the algorithm. This move not only optimized the processing flow, but also promoted the application and development of LSM technology in the fields of three-dimensional reconstruction and map making. However, it is worth noting that the generalized Gauss-Markov mathematical model relied on by the traditional LS3D method only considers the random error of the target point cloud when constructing it, while ignoring the error factors that may also exist in the source point cloud. This one-sided error consideration is particularly insufficient in modern applications that pursue high-precision registration, because it fails to fully capture and handle all potential error sources that affect the registration accuracy.

[0005] The EIV model can simultaneously consider the random errors in the source point cloud and the target point cloud, and is more rigorous in theory. In the field of measurement data processing, the data processing method based on the EIV model is called the overall least squares. It is worth noting that the overall least squares is not an estimation criterion but an estimation method, and it still follows the criterion of minimizing the sum of squares of errors. Although the overall least squares method based on the EIV model can obtain relatively good estimation results, due to the existence of repeated variables and fixed elements, the variance-covariance matrix of the coefficient matrix may be a singular matrix. The existing solution directly and violently inverts it, which obviously has theoretical deficiencies. With the deepening of research, it should be considered that the concept of all errors in the entire observation space has been deeply rooted in people's minds, and it should no longer be distinguished whether the error is located in the observation vector or the design matrix. At this time, the error processing model that combines the observation error with the estimated parameter is called the nonlinear Gauss-Helmert (GH) model. Wang et al. (2023, 2024) constructed a point cloud fitting method based on the GH model and obtained high-precision estimation results. In terms of registration, Ge et al. (2016) proposed the GH-LS3D method based on the GH model, which achieved better registration parameter estimation results and faster computational efficiency than ICP.

[0006] Unfortunately, the GH-LSD method indiscriminately includes all point pairs in its mathematical model when calculating the registration parameters. This approach directly leads to a sharp expansion of the dimensions of core operation matrices such as coefficient matrices, observation vectors, and weight matrices, which in turn solidifies the bottleneck of the GH-LS3D algorithm in terms of computational efficiency. This high-dimensional complexity not only increases the consumption of computing resources, but also makes the computational time a problem that cannot be ignored and needs to be overcome urgently. Summary of the invention

[0007] The purpose of the present invention is to provide a point cloud registration method based on a hierarchical grouping strategy in order to overcome the problems of high computational consumption and low efficiency of existing point cloud registration.

[0008] The purpose of the present invention can be achieved by the following technical solutions:

[0009] A point cloud registration method based on a hierarchical grouping strategy, the method comprising the following steps:

[0010] S1, obtain point cloud data and preprocess it to obtain an initial point cloud;

[0011] S2, performing rough registration of the initial point cloud to obtain rough registration point cloud data;

[0012] S3, constructing a KD-Tree to recursively divide the coarse registration point cloud data into multiple regions;

[0013] S4, select each point in the coarse registration point cloud data as a query point in turn, calculate the distance between the query point and all points in the KD-Tree, and determine the number of nearest neighbor points corresponding to the query point based on the distance;

[0014] S5. Use KNN search to calculate the distance between each pair of corresponding points, where the pair of corresponding points is a query point and a corresponding nearest neighbor point, and calculate the RMSE value based on the distance between the corresponding points;

[0015] S6, using the RMSE value as a threshold to eliminate corresponding points, and obtain a corresponding point pair set;

[0016] S7, dividing the point pairs in the corresponding point pair set into m groups, setting the k-1th group of the m groups of point pairs as the first group, and the kth group as the second group;

[0017] S8, obtaining the coordinate filling matrix of the k-1th group of point pairs;

[0018] S9, calculating the parameter correction vector of the k-1th group based on the coordinate filling matrix of the k-1th group of point pairs;

[0019] S10, calculating the error vector of the k-1th group of point pairs;

[0020] S11, obtaining the coordinate filling matrix of the kth group of point pairs;

[0021] S12, calculating the first matrix and the second matrix of the current iteration number;

[0022] S13, calculating the parameter correction vector of the kth group based on the parameter correction vector of the k-1th group, the first matrix and the second matrix;

[0023] S14, update the registration parameters, obtain the current registration parameters of the current iteration number, update k to k+1, and return to S11 until k reaches the maximum threshold;

[0024] S15, judging whether the kth group parameter correction vector is less than the threshold value, if so, stopping the iteration and executing S16, if not, updating the coefficient matrix A based on the error vector, updating the number of iterations, and returning to S8;

[0025] S16. Calculate the post-verification unit weighted mean error to perform accuracy estimation, and accurately align the corresponding point pair set based on the current alignment parameters to obtain a point cloud alignment result.

[0026] Furthermore, the parameter correction vector of the k-1th group is:

[0027]

[0028] in, represents the parameter correction vector of the k-1th group at the jth iteration, represents the cofactor matrix of parameter estimates for the k-1th group at the jth iteration, represents the equivalent cofactor matrix of the k-1th set of coordinate values ​​for the jth iteration, represents the coefficient matrix of the k-1th set of error equations for the jth iteration, The constant vector representing the k-1th set of error equations for the jth iteration.

[0029] Furthermore, the equivalent cofactor matrix is:

[0030]

[0031] Among them, Q k-1 The prior cofactor matrix J represents the k-1th group of coordinate values, which consists of the unit matrix and the linearized error coefficient matrix B. k-1 =[I3 B k-1 ].

[0032] Furthermore, the covariance matrix of the parameter estimation of the k-1th group is:

[0033]

[0034] Furthermore, the first matrix is:

[0035]

[0036] Among them, G k-1,k represents the first matrix, represents the cofactor matrix of parameter estimates after the k-1th set of iterations converge, The equivalent cofactor matrix A represents the parameter estimates of the k-1th set of point pairs. k Represents the coefficient matrix of the error equation for the kth set of point pairs.

[0037] Furthermore, the second matrix is:

[0038]

[0039] in, represents the second matrix, L k Represents the observation vector of the kth group of point pairs.

[0040] Furthermore, the parameter correction vector of the kth group is:

[0041]

[0042] in, represents the parameter correction vector of the kth group at the jth iteration, represents the first matrix of the j-th iteration.

[0043] Furthermore, the current registration parameter of the current iteration number is:

[0044]

[0045] in, is the current registration parameter of the j-1th iteration, is the current registration parameter for the current iteration number j, is the current registration correction parameter, which is the parameter correction vector of the kth group when k reaches the maximum threshold.

[0046] Furthermore, the error vector is:

[0047]

[0048] Q k-1 Represents the prior cofactor matrix of the k-1th group of coordinate values, L k represents the observation vector of the kth point pair, A k Represents the coefficient matrix of the error equation for the kth set of point pairs.

[0049] Furthermore, when the error vector prediction value updates the coefficient matrix A, the error vector replaces the random error in the coefficient matrix A.

[0050] Compared with the prior art, the present invention has the following beneficial effects:

[0051] The present invention not only considers the error terms in the A array but also considers the error terms in L. The GM model only considers the error terms in L. Therefore, the solution accuracy is significantly better than the estimation result based on the nonlinear GM model, and is consistent with the overall solution of the nonlinear GH model, which is due to the absolute stability of its mathematical model and objective function. At the same time, due to the grouping strategy, the processing starts from the first subset, and only the processing result of the subset and the data of the next subset are used for processing to obtain a new parameter estimation result, and then a new parameter estimation result is obtained again together with the third subset until the last subset data is processed. Due to the large amount of point cloud data, this grouping processing can avoid the inversion calculation of high-dimensional matrices, and its computational efficiency is greatly improved. Compared with the nonlinear GH model, it shows higher practicality and timeliness, especially when processing multi-period observation data or large-scale unrelated data sets, it shows excellent decomposition and processing capabilities. The iterative solution of the present invention in the parameter solution process divides the observation values ​​into several groups, and performs adjustment group by group, which overcomes the problem of low efficiency when processing massive observation data. BRIEF DESCRIPTION OF THE DRAWINGS

[0052] Figure 1 is a flow chart of the present invention;

[0053] Figure 2 is the Stanford point cloud dataset of the present invention, where Figure 2 (a) is the rabbit point cloud dataset, and (b) is the dragon point cloud dataset. Figure 2 (c) is the monkey point cloud dataset;

[0054] Figure 3 is the point cloud registration result based on the overall solution, where Figure 3 (a) is the rabbit point cloud registration result of the overall solution, and (b) in 3 is the dragon point cloud registration result of the overall solution. Figure 3 (c) is the monkey point cloud registration result of the overall solution;

[0055] Figure 4 is the point cloud registration result based on the sequential solution, where Figure 4 (a) is the rabbit point cloud registration result of the sequential solution, and (b) is the dragon point cloud registration result of the sequential solution. Figure 4 (c) is the monkey point cloud registration result of the sequential solution;

[0056] Figure 5 Comparison of the efficiency of the overall solution and sequential solution for Rabbit point cloud data;

[0057] Figure 6 Comparison of the efficiency of the overall solution and sequential solution for dragon point cloud data;

[0058] Figure 7 Comparison of the efficiency of overall solution and sequential solution for monkey point cloud data;

[0059] Figure 8 is the registration deviation under different random model conditions, where Figure 8 (a) is the registration deviation under equal weight, Figure 8 (b) is the nominal accuracy registration deviation;

[0060] Fig. 9 Comparison of the efficiency of sequential solution and overall solution for measured data. DETAILED DESCRIPTION

[0061] The present invention is described in detail below in conjunction with the accompanying drawings and specific embodiments. This embodiment is implemented based on the technical solution of the present invention, and provides a detailed implementation method and specific operation process, but the protection scope of the present invention is not limited to the following embodiments.

[0062] The present invention proposes a point cloud registration method based on a hierarchical grouping strategy. The flowchart of the method is as follows: Figure 1 As shown, the method comprises the following steps:

[0063] S1, obtain point cloud data and preprocess it to obtain an initial point cloud;

[0064] S2, performing rough registration of the initial point cloud to obtain rough registration point cloud data;

[0065] S3, constructing a KD-Tree to recursively divide the coarse registration point cloud data into multiple regions;

[0066] S4, select each point in the coarse registration point cloud data as a query point in turn, calculate the distance between the query point and all points in the KD-Tree, and determine the number of nearest neighbor points corresponding to the query point based on the distance;

[0067] S5. Use KNN search to calculate the distance between each pair of corresponding points, where the pair of corresponding points is a query point and a corresponding nearest neighbor point, and calculate the RMSE value based on the distance between the corresponding points;

[0068] S6, using the RMSE value as a threshold to eliminate corresponding points, and obtain a corresponding point pair set;

[0069] S7, dividing the point pairs in the corresponding point pair set into m groups, setting the k-1th group of the m groups of point pairs as the first group, and the kth group as the second group;

[0070] S8, obtaining the coordinate filling matrix of the k-1th group of point pairs;

[0071] S9, calculating the parameter correction vector of the k-1th group based on the coordinate filling matrix of the k-1th group of point pairs;

[0072] S10, calculating the error vector of the k-1th group of point pairs;

[0073] S11, obtaining the coordinate filling matrix of the kth group of point pairs;

[0074] S12, calculating the first matrix and the second matrix of the current iteration number;

[0075] S13, calculating the parameter correction vector of the kth group based on the parameter correction vector of the k-1th group, the first matrix and the second matrix;

[0076] S14, update the registration parameters, obtain the current registration parameters of the current iteration number, update k to k+1, and return to S11 until k reaches the maximum threshold;

[0077] S15, judging whether the kth group parameter correction vector is less than the threshold value, if so, stopping the iteration and executing S16, if not, updating the coefficient matrix A based on the error vector, updating the number of iterations, and returning to S8;

[0078] S16. Calculate the post-verification unit weighted mean error to perform accuracy estimation, and accurately align the corresponding point pair set based on the current alignment parameters to obtain a point cloud alignment result. Figure 1 The error vector update matrix A in is only performed when the iteration number is updated.

[0079] Firstly, the GH model is linearized and its parameters are estimated.

[0080] If both surfaces are analytically continuous without any deterministic or random differences, assume that each surface element f on the target surface has an exact corresponding surface element on the search surface g, and consider that the random errors on the target surface and the search surface caused by factors such as sensors, environmental conditions or measurement methods are e f =[e X ,e Y ,e Z ] T and e g =[e x ,e y ,e z ] T , according to GH-LS3D description, the function model of the i-th surface element is:

[0081]

[0082] Where T is the translation parameter vector, T = [ΔxΔyΔz] T; R is the rotation matrix, and there is

[0083]

[0084] In order to perform least squares estimation, equation (1) must be linearized through Taylor expansion:

[0085]

[0086] In the formula,

[0087]

[0088] make, J i =[I3B i ]. At this time, formula (2) can be rewritten as:

[0089] l i -J i e i =A i dξ (3)

[0090] For n surface elements, equation (3) can be transformed into:

[0091] l-Je=Adξ (4)

[0092] In the formula, l=[l1 l2…l n ] T , represents the Kronecker inner product.

[0093] For the error vector e:

[0094]

[0095] Where Q, P and D represent the prior cofactor matrix, weight matrix and variance-covariance matrix respectively; The difference represents the unit weight variance.

[0096] According to the least squares criterion, the solution of the GH-LS3D model is:

[0097]

[0098] In the formula, Q J =JQJ T .

[0099] From formula (6), we can see that the dimension of the coefficient matrix A is 3n×6. Due to the high-dimensional characteristics of point cloud data, 3n may be equal to tens of millions or even hundreds of millions, which will undoubtedly increase the time consumption of parameter calculation. For this purpose, Section 2.2 will introduce our proposed efficient point cloud registration method based on grouping strategy, which we call HiGoReg. It can obtain the same estimation results as GH-LS3D and significantly improve the computational efficiency.

[0100] Furthermore, the point cloud registration sequential solution HiGoReg based on grouping strategy is analyzed.

[0101] We first divide the data into m groups of observation vectors, and the grouping of observation vectors and random models is as follows:

[0102]

[0103] In the formula, v, u and q represent numerical values ​​and satisfy 1<v<u<q<n.

[0104] In order to derive the HiGoReg method more vividly, we take two adjacent groups as an example, assuming they are group k-1 and group k. According to equations (4) and (7), the function equations of the two groups are:

[0105]

[0106] The corresponding random model is:

[0107]

[0108] According to the least squares principle, the overall solution of the two sets of data (the calculation method refers to Section 2.1) can be expressed as:

[0109]

[0110] In the formula, Therefore, it is assumed that there is a unique solution in the k-1th group. According to the least squares principle, the parameter estimate of the observation equation for the k-1th period is:

[0111]

[0112] According to the least squares principle, the cofactor matrix of the overall solution of the two sets of data is:

[0113]

[0114] Accordingly, formula (10) can be rewritten as:

[0115]

[0116] Reconstruct equations (11) and (12) as follows:

[0117]

[0118] Replace the formula (15) with the formula (14) Substituting into formula (13), we can obtain the preliminary form of HiGoReg:

[0119]

[0120] From the above formula, it can be found that the overall solution of the two sets of data can be obtained by using the calculation results of the current set of data and the previous set, but there is an unknown matrix in formula (16): Therefore, in order to obtain the recursive relationship between the parameter cofactor matrices of adjacent groups, multiply both ends of the formula by And right multiply Then the recursive relationship between the cofactor matrices of the parameter vectors of the k-1th and kth groups satisfies:

[0121]

[0122] Using the Woodbury matrix transformation formula for equation (12), we can get:

[0123]

[0124] Comparing equations (17) and (18), we can obtain:

[0125]

[0126] make and And rewriting formula (16) can obtain the final form of HiGoReg:

[0127]

[0128] Accordingly, the error vector predicts:

[0129]

[0130] Post-test unit weighted mean error:

[0131]

[0132] In the formula, r is the degree of freedom.

[0133] According to the above derivation process, the calculation steps of HiGoReg are shown in Table 1.

[0134] Table 1 HiGoReg calculation process

[0135]

[0136] When updating the matrix A, the error vector and replace In parameter.

[0137] Based on the above detailed theoretical derivation and clear comparison flow chart, the control process of the present invention and the characteristics of each step can be summarized as follows:

[0138] (1) Import point cloud data and perform preprocessing (denoising, sampling, etc.)

[0139] The point cloud preprocessing stage of the present invention integrates efficient denoising and sampling technology. Through the statistical filtering method, the distance between each point and the points in the neighborhood is calculated, and the noise points are identified and eliminated according to the set threshold range. Parameters such as the neighborhood radius and the standard deviation multiple are adjustable to flexibly control the denoising intensity. Subsequently, the voxel grid sampling strategy is adopted to divide the point cloud space into uniform small cubes (voxels), and representative points (such as center points or mean points) are selected in each voxel to achieve sparse point cloud data while retaining the main shape features and improving the efficiency of subsequent processing.

[0140] (2) Point cloud coarse registration

[0141] To achieve preliminary alignment between the source point cloud and the target point cloud, the present invention uses CloudCompare software to manually select multiple pairs of corresponding points and automatically calculates the transformation matrix through its "Align" function. This matrix contains the rotation matrix and translation vector, which accurately describes the transformation operation required to adjust the source point cloud to roughly align with the target point cloud, laying a solid foundation for subsequent precise registration.

[0142] (3) KD-Tree Construction

[0143] For the processed point cloud data, the present invention constructs a KD-Tree to optimize the spatial search efficiency. KD-Tree recursively divides the point cloud data into multiple regions through multi-dimensional space segmentation technology, and creates index nodes to achieve the function of quickly locating the neighboring areas of the query point. By adjusting the minimum number of subset points or tree depth and other parameters, the structure of KD-Tree can be flexibly controlled, and finally an efficient spatial index structure is obtained to support efficient data access and retrieval.

[0144] (4) Set the initial value of the fine registration parameter to 0, and KNN search for the nearest point pair

[0145] After rough registration, the source point cloud and the target point cloud are roughly aligned. Since HiGoReg needs to be solved iteratively, the initial value of the registration parameter is first set to 0. In the fine registration process, KNN search (K value is set to 1) is used to find the nearest neighbor point in the KD-Tree for each point, using the Euclidean distance as the metric. This step provides an accurate set of corresponding point pairs for subsequent registration optimization.

[0146] (5) RMSE evaluation

[0147] By calculating the Euclidean distance between corresponding point pairs and obtaining the root mean square error (RMSE) based on the distance, the present invention realizes a quantitative evaluation of the registration effect. The RMSE value directly reflects the accuracy of the registration and provides a key indicator for optimizing the iterative process.

[0148] (6) Point pair elimination

[0149] In order to further improve the registration accuracy, the present invention implements a point pair elimination strategy. Based on the mean error of the point pair distance calculation, 2 times the mean error is used as the threshold to eliminate point pairs with large deviations, thereby obtaining a purer set of corresponding point pairs and reducing the interference of noise and errors.

[0150] (7) Divide the retained point pairs into m groups

[0151] In order to avoid the inherent high-dimensional matrix operations of GH-LS3D, the retained point pairs are grouped, each group containing several point pair data. The HiGoReg algorithm proposed in the present invention is used to calculate and update the registration parameters group by group, which can effectively avoid the inherent efficiency problems of high-dimensional matrix calculations.

[0152] (8) Fill in matrices A, J, L according to the coordinates of the k-1th set of point pairs

[0153] Since the point pairs are divided into several groups, two adjacent groups are taken as an example, where the k-1th group is the former group and the kth group is the latter group. According to the coordinates of the k-1th group of point pairs, the matrices A, J, and L are filled respectively to calculate the registration parameter correction number of the group.

[0154] (9) Calculate the parameter correction vector of the k-1th group

[0155] According to the filled matrix, the parameter correction number of the k-1th group is calculated based on the least squares principle and the parameter estimation equation (Equation (11)), which is used to calculate the kth group of registration parameter correction numbers.

[0156] (10) Fill in matrices A, J, L according to the coordinates of the kth set of point pairs

[0157] The matrices A, J, and L are filled according to the coordinates of the kth group of point pairs. It should be noted that if the number of observations in each group cannot be evenly divided, the matrix dimensions here may be different from those of the k-1th group. Similarly, the filled matrices are used to calculate the kth group of registration parameter corrections.

[0158] (11) Calculate matrices G and L

[0159] The Woodbury matrix transformation formula can be used to construct the internal connection between the k-1th group and the kth group, and the matrices G and L are calculated based on the parameter correction number estimation results of the k-1th group, the filling matrix of the kth group, and the cofactor matrix of the corresponding parameter correction number. The matrix obtained in this step still serves the calculation of the kth group registration parameter correction number.

[0160] (12) Estimate the parameter correction vector of the kth group of observations

[0161] According to the k-1th group of registration parameter corrections, matrices G and L, the parameter correction vector of the kth group of observations is calculated. This process realizes the continuous optimization and adjustment of parameters, laying a solid foundation for the accuracy of the overall registration process. Through such iterations and improvements, the real registration parameters can be more efficiently approached, thereby significantly improving the quality and efficiency of data processing.

[0162] (13) Calculate or update the parameter vector

[0163] With the help of the pre-set initial parameter values ​​and the estimated correction values ​​of the kth group of registration parameters obtained through precise calculation, the registration parameters in the current iteration are systematically calculated and updated in detail until all the data in the group are calculated. This process not only makes full use of the prior knowledge of the initial parameters, but also integrates the correction value information based on the latest data feedback, ensuring that the registration parameters can gradually approach the optimal solution.

[0164] (14) Calculate the error vector used to update the matrix A

[0165] Since HiGoReg itself is a nonlinear model, parameter estimation needs to be performed iteratively. Therefore, the error vector of each group needs to be estimated simultaneously during the calculation process for updating the matrix A in the next iteration.

[0166] (15) Convergence

[0167] Determine whether the parameter correction number is less than the preset threshold. If so, stop the iteration and perform accuracy estimation. Otherwise, jump to step (8) to continue execution.

[0168] (16) Accuracy estimation

[0169] The posterior unit weighted mean error is calculated using formula (22) to estimate the accuracy.

[0170] The present invention proposes a sequential solution for point cloud registration based on a hierarchical grouping strategy, where hierarchical refers to iterative solution in the parameter solution process, and each iteration represents a layer; grouping refers to dividing the observation values ​​into several groups and performing adjustment group by group. The new method overcomes the problem of low efficiency in processing massive observation data, and is theoretically derived and experimentally verified by taking point cloud registration as an example.

[0171] The present invention has the following advantages:

[0172] (1) Sequential modeling innovation that fully considers random errors and information: When constructing the model, the new method fully incorporates the random errors and random information existing in all observations, making the model structure more rigorous and closer to the actual data characteristics. The introduction of the sequential solution method not only extends the traditional GH model in theoretical depth, but also greatly enhances the adaptability and accuracy of the model in complex and changing data environments.

[0173] (2) Efficient and accurate solution strategy: The new method is significantly better than the estimation results based on the nonlinear GM model in terms of solution accuracy, and is consistent with the overall solution of the nonlinear GH model, thanks to the absolute stability of its mathematical model and objective function. At the same time, its computational efficiency is greatly improved, showing higher practicality and timeliness than the nonlinear GH model, especially when dealing with multi-period observation data or large-scale unrelated data sets, showing excellent decomposition and processing capabilities.

[0174] Wide applicability and flexible configuration: The new method is presented in the form of an abstract function, which makes it possible to apply it to problems in a wider range of fields. By adjusting the analytical model according to the specific problem, such as choosing a sequential solution based on the GM model to deal with scenarios without random error effects, it shows a high degree of flexibility and customization potential. However, it is also pointed out that the performance of the method may be limited in the case of specific data anomalies (such as gross error contamination), ill-conditioned or rank-deficient method matrices, suggesting that future research needs to further explore optimization strategies to enhance its robustness and universality.

[0175] To verify the effectiveness of the new method, the experimental data comes from the point cloud dataset publicly available at Stanford University, such as Figure 2 As shown in Figure 1, red represents the source point cloud and green represents the target point cloud. In the experiment, the unit matrix or nominal accuracy is still used to weight the point cloud observations.

[0176] The present invention needs to find all the nearest neighbor points of the source point cloud in the target point cloud, and then use these point pairs to calculate the parameters. Here, the source point cloud is first downsampled using a voxel grid filter, and then the ISS key points of the point cloud are used as calculation objects to achieve the registration of the point cloud data. After processing, the number of coordinate points in the three point clouds of Bunny, Dragon and Monkey is reduced from 35947, 437645, and 125952 to 330, 936, and 233 respectively.

[0177] The experimental design is as follows:

[0178] Solution 1: GH-LS3D;

[0179] Scheme 2: EGH-LS3D.

[0180] Since the above datasets do not provide information about the instruments, all observations are treated equally in the calculation. The three point cloud data are processed using Scheme 1 and Scheme 2 respectively; in Scheme 2, the point cloud data are registered when the number of groups is 2, 3, 4, 5, and 6, and the data in each group is distributed as evenly as possible (for example, when 330 is divided into 4 groups, the number of each group is 82, 82, 82, and 84). During the registration process, the minimum and maximum number of iterations are set to 5 and 20, respectively, and the iteration threshold is the same as Ge and Wunderlich (2016).

[0181] The registration result is as follows Figures 3-4 As shown (blue is the source point cloud data after registration); according to the number of groups, each solution is run 10 times, and the efficiency statistics are as follows Figures 5 to 7 The RMSE and conversion parameter estimation values ​​after registration are shown in Table 3.

[0182] Table 3. Root mean square error, posterior unit weighted mean error and conversion parameter estimates for point cloud registration

[0183]

[0184] In addition to the simulation experiment, the second example is based on the measured point cloud data of a statue. The Faro Focus3DX330 scanner is used to obtain data from two different angles, with distance and angle measurement accuracies of 0.4mm and 0.009° respectively. For each scan, the scanner was placed 6 to 8m away from the statue. The source surface and target surface contain 32546 and 27894 data points respectively. In the experiment, the initial approximation of the registration is calculated by manually selecting three points in CloudCompare. The function model selected for the experiment is still GH-LS3D, and the random model is divided into two types: equal weight and nominal accuracy weighting. After registration, the 3D Compare module of Geomagic Studio 2012 software is still used to compare the overlapping areas between the two point clouds, such as Figure 8 shown.

[0185] exist Figures 5 to 7 and Fig. 9 In the legend, BS represents the overall solution of GH-LS3D; S2, S3, S4, S5, S6 represent the statistical results of the sequential solution when the number of groups is 2, 3, 4, 5, 6 respectively; the suffixes E and N represent the equal weighting and nominal accuracy weighting methods respectively. Figures 3 to 7 It can be found that the registration accuracy and parameter estimation values ​​of the global solution and sequential solution of point cloud registration based on GH-LS3D are the same, which is consistent with the conclusion of the previous section, proving the equivalence of the sequential solution with the global solution in terms of parameter estimation. In addition, the RMSE of the three point cloud data is extremely small and almost negligible. This is because the target point cloud is obtained by performing a slight rotation and translation on the source point cloud in a manual interactive manner, that is, the points in the source point cloud and the target point cloud after registration are one-to-one corresponding, and there are no redundant points (such as noise points, drift points, etc.) that may affect feature matching. In terms of efficiency, the computational efficiency of the sequential solution is better than that of the global solution under given grouping conditions; but as the number of groups increases, the efficiency is not always improved. On the contrary, too many groups will lead to reduced efficiency. Figure 6 In the above data, the computational efficiency has been improving. This is because the amount of Dragon point cloud data is higher than the other two data sets, so more groups are needed to show the above trend. Continuing to group, it can be found that the Dragon data set can achieve the maximum computational efficiency when the number of groups is 9. This shows that there should not be too many groups in sequential adjustment, and there should be an optimal number of groups.

[0186] In the measured data, Figure 8 As shown in Table 4, the registration results based on equal weights and nominal accuracy weighting are almost the same (RMSE has a difference of about 1%). This may be because the TLS is close to the target and the point cloud data itself has extremely high accuracy. Therefore, the contribution of weighting to improving the registration accuracy is not significant. Fig. 9 , the calculation under equal weights is obviously faster, which depends on the fact that the unit matrix is ​​faster when inverted, while the cofactor matrix determined based on nominal accuracy is more time-consuming when inverted. In addition, the two weighting methods achieve the highest computational efficiency when the number of groups is 3 and 2 respectively. The optimal grouping under equal weights can improve the efficiency by about 15%. In the case of nominal accuracy weighting, the use of the optimal grouping can improve the efficiency by about 30% compared with the overall solution, indicating that the weighting method has a certain influence on the determination of the optimal number of groups. According to the existing experimental results, the optimal grouping is when the number of observations in each group is about 100, but a more detailed analysis should also consider factors such as time complexity, space complexity, and computer configuration, which urgently needs further research.

[0187] Table 4 Registration accuracy of different random models

[0188] Random Model RMSE(m) Average distance (m) Iterations Equal weight <![CDATA[2.7992×10 -4 ]]> 0.003 5 Nominal accuracy <![CDATA[2.7722×10 -4 ]]> 0.003 5

[0189] The preferred specific embodiments of the present invention are described in detail above. It should be understood that a person skilled in the art can make many modifications and changes based on the concept of the present invention without creative work. Therefore, any technical solution that can be obtained by a person skilled in the art through logical analysis, reasoning or limited experiments based on the concept of the present invention on the basis of the prior art should be within the scope of protection determined by the claims.

Claims

1. A point cloud registration method based on a hierarchical grouping strategy, characterized in that: The method comprises the following steps: S1, obtain point cloud data and preprocess it to obtain an initial point cloud; S2, performing rough registration of the initial point cloud to obtain rough registration point cloud data; S3, constructing a KD-Tree to recursively divide the coarse registration point cloud data into multiple regions; S4, select each point in the coarse registration point cloud data as a query point in turn, calculate the distance between the query point and all points in the KD-Tree, and determine the number of nearest neighbor points corresponding to the query point based on the distance; S5. Use KNN search to calculate the distance between each pair of corresponding points, where the pair of corresponding points is a query point and a corresponding nearest neighbor point, and calculate the RMSE value based on the distance between the corresponding points; S6, using the RMSE value as a threshold to eliminate corresponding points, and obtain a corresponding point pair set; S7, dividing the point pairs in the corresponding point pair set into m groups, setting the k-1th group of the m groups of point pairs as the first group, and the kth group as the second group; S8, obtaining the coordinate filling matrix of the k-1th group of point pairs; S9, calculating the parameter correction vector of the k-1th group based on the coordinate filling matrix of the k-1th group of point pairs; S10, calculating the error vector of the k-1th group of point pairs; S11, obtaining the coordinate filling matrix of the kth group of point pairs; S12, calculating the first matrix and the second matrix of the current iteration number; S13, calculating the parameter correction vector of the kth group based on the parameter correction vector of the k-1th group, the first matrix and the second matrix; S14, update the registration parameters, obtain the current registration parameters of the current iteration number, update k to k+1, and return to S11 until k reaches the threshold; S15, judging whether the kth group parameter correction vector is less than the threshold value, if so, stopping the iteration and executing S16, if not, updating the coefficient matrix A based on the error vector, updating the number of iterations, and returning to S8; S16. Calculate the post-verification unit weighted mean error to perform accuracy estimation, and accurately align the corresponding point pair set based on the current alignment parameters to obtain a point cloud alignment result.

2. A point cloud registration method based on a hierarchical grouping strategy according to claim 1, characterized in that: The parameter correction vector of the k-1th group is: in, represents the parameter correction vector of the k-1th group at the jth iteration, represents the cofactor matrix of parameter estimates for the k-1th group at the jth iteration, represents the equivalent cofactor matrix of the k-1th set of coordinate values ​​for the jth iteration, represents the coefficient matrix of the k-1th set of error equations for the jth iteration, The constant vector representing the k-1th set of error equations for the jth iteration.

3. A point cloud registration method based on a hierarchical grouping strategy according to claim 2, characterized in that: The equivalent cofactor matrix is: Among them, Q k-1 The prior cofactor matrix J represents the k-1th group of coordinate values, which consists of the unit matrix and the linearized error coefficient matrix B. k-1 =[I3 B k-1 ].

4. A point cloud registration method based on a hierarchical grouping strategy according to claim 3, characterized in that: The cofactor matrix of the parameter estimation of the k-1th group is:

5. A point cloud registration method based on a hierarchical grouping strategy according to claim 4, characterized in that: The first matrix is: Among them, G k-1,k represents the first matrix, represents the cofactor matrix of parameter estimates after the k-1th set of iterations converge, The equivalent cofactor matrix A represents the parameter estimates of the k-1th set of point pairs. k Represents the coefficient matrix of the error equation for the kth set of point pairs.

6. A point cloud registration method based on a hierarchical grouping strategy according to claim 5, characterized in that: The second matrix is: in, represents the second matrix, L k Represents the observation vector of the kth group of point pairs.

7. A point cloud registration method based on a hierarchical grouping strategy according to claim 6, characterized in that: The parameter correction vector of the kth group is: in, represents the parameter correction vector of the kth group at the jth iteration, represents the first matrix of the j-th iteration.

8. The point cloud registration method based on hierarchical grouping strategy according to claim 1, characterized in that: The current registration parameters for the current number of iterations are: in, is the current registration parameter of the j-1th iteration, is the current registration parameter for the current iteration number j, is the current registration correction parameter, which is the parameter correction vector of the kth group when k reaches the maximum threshold.

9. The point cloud registration method based on hierarchical grouping strategy according to claim 1, characterized in that: The error vector is: Q k-1 Represents the prior cofactor matrix of the k-1th group of coordinate values, L k represents the observation vector of the kth point pair, A k Represents the coefficient matrix of the error equation for the kth set of point pairs.

10. The point cloud registration method based on hierarchical grouping strategy according to claim 1, characterized in that: When updating the coefficient matrix A based on the error vector prediction value, the error vector replaces the random error in the coefficient matrix A.

Images