Multi-level feature fusion adaptive point cloud registration method

Abstract

The invention discloses a multilevel feature fusion adaptive point cloud registration method, and relates to the field of point cloud processing technology and industrial visual inspection technology, and the method comprises the steps: carrying out the preprocessing of a point cloud pair; multi-scale feature extraction is carried out on the preprocessed point cloud pair, a feature pyramid is constructed through the extracted multi-scale features, and a top-down guide mechanism is adopted between layers of the feature pyramid; performing preliminary feature fusion on feature points of the point cloud pair based on the extracted multi-scale features; performing multi-level registration on the point cloud pair; the multiple layers comprise first-layer registration, second-layer registration and third-layer registration, a transformation matrix of each layer is obtained after each layer of registration is completed, and registration is carried out on the to-be-registered point cloud through the transformation matrix of each layer; according to the method, a multi-scale feature extraction strategy is adopted, three feature descriptors with different scales, namely ESF, SHOT and FPFH feature descriptors, are selected for the degradation problem in point cloud registration, and feature extraction and parameter adaptive adjustment from macroscopic to microscopic are achieved.

Description

Technical Field

[0001] The present invention relates to the technical fields of point cloud processing technology and industrial vision detection technology, and particularly relates to a multi-level feature fusion adaptive point cloud registration method. Background Art

[0002] Point cloud registration technology is a key link in 3D data processing and has important applications in fields such as industrial silo monitoring and volume measurement. For example, the feeding and discharging of the silo is a continuous process. If we want to measure the change amount between two point clouds in the granary before and after, due to the influence of different external environments (the influence of temperature, humidity, dust, etc. on radar scanning and the small offset caused by physical vibration to the hardware device on the point cloud angle), it is very difficult to keep the current scanning state consistent with the previous scanning state. It is necessary to select appropriate features to align the point clouds, so as to calculate the difference by combining the historical point cloud and the current point cloud, and then calculate the change amount.

[0003] With the rapid development of Light Detection and Ranging (LiDAR) technology, continuous monitoring of the shape of the material pile can be achieved through a scanning device fixedly installed on the top of the silo. However, in the actual industrial environment, the acquisition and processing of point cloud data still face many challenges. In practical applications, the quality of point cloud data is often degraded due to various factors. First, the non-rigid transformation caused by the jitter and temperature drift of the device itself makes the deviation between continuously collected point clouds; second, the noise interference caused by factors such as dust and water vapor in the industrial environment. In addition, due to the irregularity of the material pile surface and the dust generated when the material falls, the point cloud data often has quality problems such as uneven density and local point loss. Existing point cloud registration methods have obvious deficiencies when dealing with industrial silo scenarios. Although the traditional Iterative Closest Point (ICP) algorithm is widely used, it is easy to fall into local optimum when dealing with large-scale deformations; methods based on features such as Fast Point Feature Histogram (FPFH) are unstable in the presence of a large amount of noise and deformations. Most existing methods adopt a fixed parameter strategy and lack an adaptive mechanism for point cloud quality, making it difficult to cope with the complex and variable data quality problems in the industrial environment. These problems seriously restrict the application effect of point cloud registration technology in industrial silo monitoring. Therefore, it is of great engineering practical significance to develop a registration method that can effectively process point cloud degradation, adapt to data quality, and has high accuracy. Summary of the Invention

[0004] To solve the technical problems proposed in the background art, the present invention proposes a multi-level feature fusion adaptive point cloud registration method, including:

[0005] S1. Preprocess the point cloud pair;

[0006] S2. Extract multi-scale features from the preprocessed point cloud pairs, construct a feature pyramid through the extracted multi-scale features, and adopt a top-down guidance mechanism between the levels of the feature pyramid;

[0007] S3. Perform preliminary feature fusion on the feature points of the point cloud pairs based on the extracted multi-scale features;

[0008] S4. Perform multi-level registration on the point cloud pairs; the multi-level includes the first-level registration, the second-level registration, and the third-level registration. After each level of registration is completed, the transformation matrix of that level is obtained, and the point cloud to be registered is registered through the transformation matrix of each level.

[0009] Specifically, the S1 includes the following steps:

[0010] S1.1 Perform noise reduction processing on the point cloud, identify and remove abnormal points by analyzing the neighborhood distribution characteristics of the points;

[0011] S1.2 Downsample the point cloud using the voxel grid filtering algorithm;

[0012] S1.3 Calculate the average point spacing d of the downsampled point cloud, and use it as the benchmark parameter for feature extraction;

[0013] S1.4 Construct an octree-based spatial index structure for neighborhood search.

[0014] Specifically, the S2 includes the following steps:

[0015] S2.1. Extract the first-scale features using the ESF features to obtain the ESF feature descriptor ESF(P). The first-scale features serve as the top layer of the pyramid, perform silo structure recognition and region segmentation on the preprocessed point cloud, obtain the cylindrical section, conical section, and transition section, and verify the rationality of the region segmentation through the ESF features;

[0016] S2.2. Extract the second-scale features using the SHOT features to obtain the SHOT feature descriptor SHOT(p); the second-scale features serve as the middle layer of the pyramid, are used to extract the local geometric structure, and transfer the information of the middle layer of the feature pyramid to the bottom layer to provide guidance for the bottom layer;

[0017] S2.3. Extract the third-scale features using the FPFH features to obtain the FPFH feature descriptor FPFH(p), and the third-scale features serve as the bottom layer of the pyramid for extracting the fine geometric structure;

[0018] For each point cloud, its set of feature points is its set of FPFH feature points, and the feature points in the subsequent S3 and S4 are the feature points of this set of FPFH feature points.

[0019] Specifically, the S2.1 includes the following steps:

[0020] S2.1.1 Perform cylindrical segment detection on the point cloud, extract the cylinder radius r and the axial direction vector v, and use the axial direction as the z-axis direction;

[0021] S2.1.2 Extract the ESF feature descriptor ESF(P) of the point cloud, and perform reliability verification on the cylindrical segment through the ESF feature. If the verification passes, record the cylinder parameters, the cylinder radius r and the axial direction vector v, and mark the points in the point cloud that conform to the cylinder model characteristics as cylindrical segments, and enter S2.1.4; if the verification fails, enter S2.1.3; the reliability verification criteria include two aspects: (1) whether the deviation of the D2 distribution of the ESF feature from the theoretical cylinder established based on the cylinder radius r extracted by RANSAC is within the range of ε1, and (2) whether the deviation of the A3 component of the ESF feature from the roundness of the cross-section of the theoretical cylinder is within the range of ε2. Both aspects need to be satisfied simultaneously for the verification to pass.

[0022] S2.1.3 Crop the point cloud in the z-axis direction, remove the point cloud data within δh×H near the highest point and the lowest point, where H is the total height of the point cloud and δh is the cropping ratio; after cropping, return to step S2.1.1 and accumulate the cropping times; when the cropping times reach the maximum number of attempts N, if the verification still fails, the quality of the point cloud is insufficient to extract reliable cylinder features, and the registration is exited;

[0023] S2.1.4 Further divide the point cloud into cylindrical segments, transition segments and conical segments based on the segmented cylindrical segments, where the transition segment is the area within the boundary (-tr, tr) of the cylindrical segment, and the conical segment is the remaining area;

[0024] Specifically, S2.2 includes the following steps:

[0025] S2.2.1 Define the structure indication function I(p),

[0026]

[0027] where p represents the points in the point cloud that have been divided, β is the sampling correction coefficient for the conical segment, μ is the sampling density coefficient for the transition segment, β∈(1, 2), μ∈(1, 2);

[0028] S2.2.2 Calculate the basic sampling interval,

[0029] ds_base(p) = d / ρ(p)

[0030] ρ(p) = α·I(p)

[0031] where d is the average point spacing of the point cloud, α is the basic sampling coefficient, β is the sampling correction coefficient for the conical segment, and ρ(p) is the sampling density function;

[0032] S2.2.3. Obtain the local quality score Q(p),

[0033] Q(p) = (Q_d(p) + Q_n(p) + Q_r(p)) / 3

[0034] Q_d(p)= min(1, |N(p)| / N_exp)

[0035] Q_n(p) = 1 - (1 / |N(p)|)∑(arccos(|n_i·n_p|) / π)

[0036] Q_r(p) = exp(-σ / d)

[0037] Among them, Q_d(p) is the local point density score of point p, N(p) is the local neighborhood point set of point p. For any point p in the point cloud, its neighborhood N(p) is defined as the set of all points within a spherical space centered at p with a radius of r,

[0038] N(p) = {q | ||q - p|| ≤ r}

[0039] where r is the neighborhood radius, ||·|| represents the Euclidean distance, N_exp is the expected number of neighborhood points, Q_n(p) is the normal vector consistency score of point p, n_i and n_p are the normal vectors of the neighborhood point and the center point respectively, Q_r(p) is the local plane fitting residual score of point p, and σ is the root mean square error of local plane fitting based on the principal component analysis method;

[0040] S2.2.4. Obtain the adaptive sampling interval ds(p) based on the local quality score,

[0041] ds(p) = ds_base(p)·(1 + γ(1-Q(p)))

[0042] where γ is the sampling adjustment coefficient, Q(p)∈[0,1]. When Q(p)=1, the sampling interval is maintained. When Q(p) decreases, the sampling interval gradually increases through 1+γ(1-Q(p));

[0043] S2.2.5. Construct the SHOT local reference frame;

[0044] For the cylindrical section and the part of the transition section on the cylindrical structure, use the identified cylindrical axis as the main direction of the reference frame, and then determine the secondary direction in combination with the normal vectors of the local neighborhood point sets of each point on the cylindrical section point cloud;

[0045] For the conical section and the transition section on the conical structure, the principal directions of the local surface are calculated, and the complete reference frame is determined using the principal component analysis method. Specifically, first, the local neighborhood point set of each point on the point cloud of the conical section is obtained, and the principal component analysis is performed on this point set. The eigenvector corresponding to the largest eigenvalue is used as the principal direction of the reference frame, the eigenvector corresponding to the second-largest eigenvalue is used as the secondary direction, and the eigenvector corresponding to the smallest eigenvalue is used as the third direction.

[0046] S2.2.6. Extract SHOT feature points from the point cloud. For each point p in the point cloud, the points in the point cloud are sorted based on its local quality score Q(p). Traverse the sorted point list, and for the neighborhood N1(p) of the current point p, N1(p) = {q | ||q - p|| ≤ ds(p)}

[0047] where ds(p) is the adaptive sampling distance; if there is no SHOT feature point already added in the neighborhood N1(p) of the current point p, then p is added to the SHOT feature point set.

[0048] S2.2.7. Calculate the SHOT feature descriptor SHOT(p), and calculate the feature response value S(p) for each point p in the SHOT feature point set of the point cloud, including the following steps:

[0049] (1) Based on the local reference frame established in S2.2.5, divide the local spherical support region of point p into 32 spatial grids, calculate the normal vector distribution histogram of 11 bins in each grid, and obtain a 352-dimensional SHOT feature vector. This SHOT feature vector is the SHOT feature descriptor SHOT(p);

[0050] (2) Reorganize the feature vector into a 32×11 matrix M with every 11 elements as a row, which is used to calculate the feature response value S(p).

[0051] S(p)= mean(diff(i,j))=mean(|M[i] - M[j]|)

[0052] where M[i] is the 11-dimensional histogram vector of the i-th spatial grid, diff(i,j) is the histogram difference between adjacent grids, and mean(·) is the averaging operation.

[0053] Specifically, the 2.3 includes the following steps:

[0054] S2.3.1. Sort the points in the SHOT feature point set according to the feature response value S(p). Traverse the sorted point list, and check its neighborhood N2(p) for the current point, which is defined as

[0055] N2(p) = {q | ||q - p|| ≤ df(p)}

[0056] df(p) is the sampling spacing, df(p) = d, that is, the average point spacing d of the point cloud downsampling. If there are no selected FPFH feature points in the neighborhood, then add p to the FPFH feature point set of this point cloud;

[0057] S2.3.2. Based on the feature response value S(p) of the SHOT feature descriptor, for each point p in the FPFH feature point set, construct an adaptive support radius,

[0058] R(p) = d·[γ·H(S(p)-θ) + λ·(1-H(S(p)-θ))]

[0059] where γ and λ are both radius adjustment coefficients and satisfy λ > γ > 0;

[0060] Obtain the support radius neighborhood of point p for use in subsequent FPFH feature descriptor calculations,

[0061] N3(p) = {q | ||q - p|| ≤ R(p)}

[0062] S2.3.3. Based on the feature response value S(p) of the SHOT feature descriptor, for each point p in the FPFH feature point set, construct an adaptive feature weight mapping function,

[0063] Wf(p) = η·S(p)

[0064] where η is the weight adjustment coefficient;

[0065] S2.3.4. Calculate the FPFH feature descriptor FPFH(p) for each point p in the FPFH feature point set,

[0066] FPFH(p) = SPFH(p) + (1 / |N3(p)|) * ∑(Wf(pk) / dk * SPFH(pk))

[0067] where FPFH(p) is the FPFH feature descriptor of point p, SPFH(p) is the local geometric feature histogram of point p for describing local geometric features, N3(p) is the support radius neighborhood centered at point p, |N3(p)| represents the number of points in the support radius neighborhood, dk is the Euclidean distance from the neighborhood point pk to the center point p, and Wf(pk) is the feature weight of the neighborhood point pk.

[0068] Specifically, the said S3 includes the following steps:

[0069] S3.1. Calculate the feature distances of the point cloud pairs, including the overall point cloud distance dist_ESF based on the global ESF feature, the distance dist_SHOT(p,q) between the corresponding feature points of the local features based on the SHOT feature, and the distance dist_FPFH(p,q) between the corresponding feature points of the local features based on the FPFH feature.

[0070] dist_ESF = ||ESF(P) - ESF(Q)||2

[0071] dist_SHOT(p,q) = ||SHOT(p) - SHOT(q)||2

[0072] dist_FPFH(p,q) = ||FPFH(p) - FPFH(q)||2

[0073] Among them, dist_ESF is the overall point cloud distance based on the global ESF feature, ESF(P) and ESF(Q) are the ESF feature descriptors of the source point cloud P and the target point cloud Q respectively, || ||2 is the Euclidean distance calculation, p and q are the corresponding feature points in the source point cloud and the target point cloud respectively, SHOT(p) and FPFH(p) are the SHOT feature descriptor and the FPFH feature descriptor of the feature point p, and dist_SHOT(p,q) and dist_FPFH(p,q) are the SHOT and FPFH feature distances of the feature point pair (p,q) respectively.

[0074] S3.1. Normalize the feature distances of the point cloud pairs.

[0075] Among them, represents the maximum value of the SHOT feature distances of all feature point pairs. represents the maximum value of the FPFH feature distances of all feature point pairs; dist_ESF_norm is the normalized ESF feature distance, and dist_SHOT_norm(p,q) and dist_FPFH_norm(p,q) are the normalized SHOT and FPFH feature distances of the feature point pair (p,q) respectively.

[0076] Specifically, the S4 includes the following steps:

[0077] S4.1. Obtain the first-layer transformation matrix T_coarse, perform feature matching in the global range, and achieve the first-layer registration of the point cloud.

[0078] S4.2. Obtain the second-layer transformation matrix T_optimize and achieve the second-layer registration of the point cloud.

[0079] S4.3. Obtain the third - layer transformation matrix \(T_{fine}\) to achieve the third - layer registration of the point cloud;

[0080] S4.4. Obtain the registered point cloud \(Q_p\), and calculate \(Q_p = T_{fine}\times T_{optimize}\times T_{coarse}\times Q\).

[0081] Specifically, in the above - mentioned S4, the registration process for each layer is as follows:

[0082] Determine the search space for the registration of this layer;

[0083] Determine the weights of the ESF feature distance, SHOT feature distance, and FPFH feature distance for the point - cloud pairs in the registration of this layer;

[0084] Calculate the fused feature distance in a weighted manner in the search space based on the three feature distances and their respective weights;

[0085] Determine the matching point pairs based on the fused feature distance and the region - adaptive matching strategy;

[0086] Conduct geometric consistency verification on the results of the matching point pairs;

[0087] Calculate the transformation matrix according to the set of matching point pairs that pass the geometric consistency verification, and register the point cloud to be registered through the transformation matrix of each layer.

[0088] Specifically, the above - mentioned S4.1 includes the following specific steps:

[0089] S4.1.1. Determine the search space; within the global range, the search space for any feature point \(p\) is,

[0090] \(\Omega_1(p)=\{q\in Q|\vert\vert q - p\vert\vert\leq R_{max}\}\)

[0091] where \(\Omega_1(p)\) represents the search space for the feature point \(p\) within the global range \(R_{max}\), and \(R_{max}\) is the maximum search radius;

[0092] S4.1.2. Determine the feature - weight vector,

[0093]

[0094] where \(w_1\), \(w_2\), and \(w_3\) respectively correspond to the initial weights of the ESF feature, SHOT feature, and FPFH feature, \(0.2\leq w_i\leq0.5\), \(w_1 + w_2 + w_3 = 1\), and \(w_1>w_2\), \(w_1>w_3\),

[0095] S4.1.3. Calculate the fused feature distance \(D_{final}(p,q)\) of the point pair:

[0096]

[0097] S4.1.4 Determine the matching point pairs based on the fused feature distance D_final(p,q) and the region - adaptive matching strategy; including:

[0098] (1) Calculate the fused distance D_final(p,q) between all feature points q in the point cloud Q to be registered within the search space Ω(p) for each feature point p in the target point cloud P;

[0099] (2) Find two points q1 and q2 with the smallest distances to the feature point p in the point cloud Q to be registered. The corresponding distances are D1 and D2, and different matching criteria are adopted according to the region where p is located;

[0100] For the feature points in the transition section, if D1 / D2 < λ1 and the cross - validation rule is required to be satisfied, then the matching point of this feature point p in P in Q is q1; Similarly, for the feature point q1 in Q, if its nearest neighbor in P is p, it is considered that the matching point of this feature point q1 in Q in P is p; If both of the above are satisfied, the matching point pair (p,q1) is retained;

[0101] For the feature points in the cylindrical section, if D1 / D2 < λ2, for the feature point q1 in Q, if its nearest neighbor in P is p, it is considered that the matching point of this feature point q1 in Q in P is p;

[0102] Where λ1 < λ2 < 1, both are the ratio thresholds for feature point applications;

[0103] S4.1.5. Conduct geometric consistency verification on the matching point pair results; Calculate the geometric consistency index for any two randomly selected pairs of matching point pairs (pi,qi) and (pj, qj),

[0104] ε_ij = |||pi - pj|| / ||qi - qj|| - 1|

[0105] If ε_ij < εd, this matching point pair passes the verification; otherwise, it fails the geometric consistency verification; where, (pi,pj) is a point pair in the source point cloud, (qi, qj) is a point pair in the target point cloud, ε_ij is the distance ratio deviation, and εd is the geometric consistency threshold;

[0106] Define Q1 = Nc / Nt, where Nc is the number of matching point pairs that pass the geometric consistency constraint, Nt is the total number of matching point pairs, and Q1 is the proportion of matching point pairs that satisfy the geometric constraint;

[0107] When Q1 < 0.6, it indicates insufficient geometric consistency. Perform the following update, return to step 4.1.3 to recalculate D_final, and accumulate the update times k;

[0108]

[0109] Among them, σpace is the weight adjustment step size, σpace = 0.05. If Q1 >= 0.6 or the cumulative update times k > 5, then the weight is not updated;

[0110] S4.1.6 Based on the set of filtered matching point pairs, use the RANSAC - SVD method to calculate the transformation matrix T_coarse;

[0111] The said S4.2 includes the following steps:

[0112] S4.2.1. Determine the search range based on the transformation matrix T_coarse: For any feature point p, take its position p = T_coarse·p after being transformed by the transformation matrix T_coarse as the search center, and its local search space Ω2(p) is defined as,

[0113] Ω2(p) = {q ∈ Q | ||q - p|| ≤ Rmiddle}

[0114] Rmiddle = βRmax

[0115] Among them, Rmiddle is the local search radius, and β is the search radius reduction coefficient;

[0116] S4.2.2. Perform weight configuration based on the weights defined in S4.1.2, w1 = w3 = 0.25, w2 = 0.5;

[0117] S4.2.3. Based on the local search space Ω2(p), weights, and the formula defined in S4.1.3, calculate the fusion feature distance D_final in the medium registration stage;

[0118] S4.2.4. Calculate new matching point pairs based on the local search space Ω2(p) and the definition in S4.1.4;

[0119] S4.2.5. Calculate Q1 based on the new matching point pairs and the definition in S4.1.5. The specific judgment criteria, weight update logic, and transformation matrix calculation results are the same as in step 4.1.5;

[0120] S4.3 includes the following steps:

[0121] S4.3.1. Determine the search range based on the optimized transformation matrix T_optimize: Rmin = γRmax, γ ≈ 0.1, where Rmin is the fine search radius. For any feature point p, the position p = T_optimize · p after being transformed by the transformation matrix T_optimize is used as the search center, and its local search space Ω3(p) is defined as,

[0122] Ω3(p) = {q ∈ Q | ||q - p|| ≤ Rmin}

[0123] S4.3.2. Perform weight configuration based on the weights defined in step 4.1.2, w1 = 0.2, w2 = 0.3, w3 = 0.5;

[0124] S4.3.3. Calculate D_final in the fine registration stage based on the above local search space Ω3(p), weights, and the definition in step 4.1.3;

[0125] S4.3.4. Calculate new matching point pairs based on the above local search space Ω3(p) and the definition in step 4.1.4;

[0126] S4.3.5. Define local shape similarity and calculate the local shape similarity index of the matching point pairs obtained in step 4.3.4,

[0127] s_i = cos(SHOT(pi), SHOT(qi))

[0128] where SHOT(pi) and SHOT(qi) are the SHOT descriptors of the matching point pi and the point qi respectively;

[0129] Define Q2 = (1 / Nt) * ∑s_i as the overall local shape similarity score, and Nt is the total number of matching point pairs;

[0130] When Q2 < 0.7, it indicates insufficient local similarity. Perform the following updates, return to step 4.3.3 to recalculate D_final, and accumulate the update times k,

[0131]

[0132] If Q2 >= 0.7 or the cumulative update times k > 5, then no weight update is performed. Based on the set of selected matching point pairs, use the RANSAC + SVD method to calculate the fine optimization transformation matrix T_fine.

[0133] After adopting the above scheme, the beneficial effects of the present invention are as follows:

[0134] The weight adjustment strategy based on geometric consistency and local shape similarity realizes the adaptive optimization of the feature fusion process by dynamically balancing the contributions of global features and local features. When the geometric consistency is poor, the weight of the ESF feature is increased to enhance the global structure matching; when the local similarity is insufficient, the weights of the SHOT and FPFH features are increased to improve the local matching accuracy. This dual verification mechanism provides a more comprehensive and reliable quality evaluation standard for feature fusion. Based on this verification mechanism, the problems in the matching can be identified: when the geometric consistency is poor, it indicates that there is a deviation in the global structure matching; when the local similarity is low, it means that the local feature matching is not accurate enough. This problem positioning ability provides clear direction guidance for subsequent parameter adjustment. Description of the Drawings

[0135] Figure 1 This is the main flowchart of the present invention. Detailed Implementation Manner

[0136] The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative efforts fall within the protection scope of the present invention.

[0137] In the description of the present invention, the terms "first" and "second" are only used for descriptive purposes, and cannot be construed as indicating or implying relative importance or implicitly specifying the quantity of the indicated technical features. Thus, the features defined with "first" and "second" may explicitly or implicitly include one or more of the described features. In the description of the present invention, "a plurality of" means two or more, unless otherwise specifically defined.

[0138] In the description of the present invention, the term "for example" is used to mean "serving as an example, illustration, or explanation". Any embodiment described as "for example" in the present invention is not necessarily to be construed as more preferred or more advantageous than other embodiments. In order for any person skilled in the art to implement and use the present invention, the following description is given. In the following description, details are set forth for the purpose of explanation. It should be understood that those skilled in the art can realize that the present invention can be implemented without using these specific details. In other instances, well-known structures and processes are not described in detail to avoid unnecessary details from obscuring the description of the present invention. Therefore, the present invention is not intended to be limited to the embodiments shown, but is to be accorded the widest scope consistent with the principles and features disclosed herein.

[0139] This specific implementation will explain the method of the present invention more clearly and completely, such asFigure 1 As shown in Figure 1 , a multi-level feature fusion adaptive point cloud registration method of the present invention includes the following steps:

[0140] S1. Preprocess the point cloud pair; before feature extraction, it is necessary to preprocess the original point cloud data to improve the reliability of subsequent processing. The S1 includes the following steps:

[0141] S1.1 Use the statistical outlier removal algorithm to denoise the point cloud, and identify and remove outliers by analyzing the neighborhood distribution characteristics of points.

[0142] S1.2 Use the voxel grid filtering algorithm to downsample the point cloud; achieve the uniformization of point cloud density, and at the same time reduce the data volume to improve the processing efficiency.

[0143] S1.3 Calculate the average point spacing d of the downsampled point cloud, and use it as the reference parameter for feature extraction. Specifically, this parameter will be used to determine the basic sampling spacing ds_base(p) for SHOT feature extraction, the sampling spacing df(p) for FPFH feature, and its support radius R(p).

[0144] S1.4 Construct an octree-based spatial index structure for neighborhood search; this structure can improve the efficiency of neighborhood search in the subsequent feature calculation process. In the following text, the neighborhood is based on the octree to determine the adjacent relationship and distance limit of points, and the set of adjacent points within a certain distance.

[0145] In the above steps, the point cloud pair refers to the source point cloud and the target point cloud, or the first point cloud to be registered and the second point cloud to be registered. In the subsequent steps, if there is no special limitation, this interpretation will be used; the preprocessing lays the foundation for the subsequent feature extraction and registration process.

[0146] S2. Extract multi-scale features from the preprocessed point cloud pair, and construct a feature pyramid through the extracted multi-scale features; a top-down guidance mechanism is adopted between the levels of the feature pyramid.

[0147] This solution adopts a multi-scale feature extraction strategy. Aiming at the degradation problem in point cloud registration, three mature feature descriptors at different scales are selected: the ESF feature descriptor at the macro scale is used to capture global structure information, the SHOT feature descriptor at the meso scale describes local geometric structures, and the FPFH feature descriptor at the micro scale provides fine geometric descriptions; and a feature pyramid is constructed through the ESF, SHOT, and FPFH feature descriptors to achieve feature extraction and parameter adaptive adjustment from macro to micro. The hierarchical structure of the feature pyramid from top to bottom is as follows:

[0148] 1) ESF feature (top layer): Capture global structure information and provide a basis for structure classification;

[0149] 2) SHOT features (middle layer): Describe local geometric structures and are guided by ESF;

[0150] 3) FPFH features (bottom layer): Provide fine geometric descriptions and are guided by SHOT;

[0151] In the present invention, macro, meso, and micro are only relative concepts, aiming to distinguish the scales of the three feature descriptors and not representing a specific size. In the following description, these three concepts are written as the first scale, the second scale, and the third scale respectively.

[0152] S2.1. Extract the first-scale features using ESF features to obtain the ESF feature descriptor ESF(P). The first-scale features serve as the top layer of the pyramid. Perform silo structure recognition and region segmentation on the preprocessed point cloud to obtain the cylindrical section, conical section, and transition section, and verify the rationality of the region segmentation through ESF features. Through its global structure description ability, ESF features provide structural-level guidance for SHOT features. This top-down guidance mechanism ensures that the middle-layer features can obtain the support of global structure information while maintaining local description ability. Here, structure recognition and region segmentation specifically involve segmenting the point cloud into the cylindrical section, conical section, and transition section. The principle of lidar is the reflection of optical signals, so the transition section and conical section are greatly affected by the inclination angle, and there is a high possibility of point cloud defects (i.e., expansion, contraction, irregular deformation, missing, etc.) in actual operation. To solve this problem, the structure is divided into 3 sections, and in subsequent steps, separate processing will be carried out according to the point cloud distribution.

[0153] S2.1.1 Detect the cylindrical section of the point cloud through the RANSAC technique, extract the cylindrical radius r and the axial direction vector v, and use the axial direction as the z-axis direction.

[0154] S2.1.2 Extract the ESF feature descriptor ESF(P) of the point cloud, and verify the reliability of the cylindrical section through ESF features. If the verification passes, record the cylindrical parameters, the cylindrical radius r and the axial direction vector v, and mark the points in the point cloud that conform to the cylindrical model features as the cylindrical section and enter S2.1.4; if the verification fails, enter S2.1.3; The reliability verification criteria include two aspects: (1) whether the deviation of the D2 distribution of the ESF feature from the theoretical cylinder established based on the cylindrical radius r extracted by RANSAC is within the range of ε1, and (2) whether the deviation of the A3 component of the ESF feature from the roundness of the cross-section of the theoretical cylinder is within the range of ε2. Both aspects need to be satisfied simultaneously for the verification to pass. In this specific implementation, ε1 is taken as 0.1 and ε2 is taken as 0.05, which means that the D2 distribution deviation is allowed to be within 10% and the roundness deviation is allowed to be within 5%.

[0155] S2.1.3 Crop the point cloud in the z-axis direction, removing the point cloud data within δh×H near the highest and lowest points, where H is the total height of the point cloud and δh is the cropping ratio; after cropping, return to step S2.1.1 and accumulate the cropping times; when the cropping times reach the maximum attempt times N, if the verification still fails, the point cloud quality is insufficient to extract reliable cylindrical features, and the registration is exited.

[0156] S2.1.4 Further divide the point cloud into cylindrical segments, transition segments, and conical segments based on the segmented cylindrical segments, where the transition segment is the area within the range (-tr, tr) of the cylindrical segment boundary, and the conical segment is the remaining area. In this specific implementation, tr is taken as 0.2m.

[0157] S2.2. Extract the second-scale features using SHOT features to obtain the SHOT feature descriptor SHOT(p); the second-scale features are used as the middle layer of the pyramid to achieve local geometric structure extraction; and the information of the middle layer of the feature pyramid is transmitted to the bottom layer to provide guidance for the bottom layer; the SHOT features use their stable description ability of local geometric structures to provide local geometric significance guidance for the FPFH features, which enables the bottom-layer features to obtain more reliable local structure support while maintaining fine description ability; the S2.2 includes:

[0158] S2.2.1. Define the structure indication function I(p),

[0159]

[0160] where p represents the points in the point cloud that have been divided, β is the conical segment sampling correction coefficient, μ is the transition segment sampling density coefficient, β ∈ (1, 2), and μ ∈ (1, 2).

[0161] S2.2.2. Calculate the basic sampling interval:

[0162] ds_base(p) = d / ρ(p)

[0163] ρ(p) = α·I(p)

[0164] where d is the average point spacing of the point cloud, α is the basic sampling coefficient, β is the conical segment sampling correction coefficient, and ρ(p) is the sampling density function.

[0165] S2.2.3. To improve the reliability of feature extraction, design a local quality assessment. Obtain the local quality score Q(p) and comprehensively evaluate three key features,

[0166] Q(p) = (Q_d(p) + Q_n(p) + Q_r(p)) / 3

[0167] Q_d(p)= min(1, |N(p)| / N_exp)

[0168] Q_n(p) = 1 - (1 / |N(p)|)∑(arccos(|n_i·n_p|) / π)

[0169] Q_r(p) = exp(-σ / d)

[0170] Among them, Q_d(p) is the local point density score of point p, N(p) is the set of local neighborhood points of point p. For any point p in the point cloud, its neighborhood N(p) is defined as the set of all points within a spherical space centered at p with a radius of r.

[0171] N(p) = {q | ||q - p|| ≤ r}

[0172] Among them, r is the neighborhood radius, ||·|| represents the Euclidean distance, and all subsequent neighborhood definitions are based on this. N_exp is the expected number of neighborhood points, Q_n(p) is the normal vector consistency score of point p, n_i and n_p are the normal vectors of the neighborhood point and the center point respectively, Q_r(p) is the local plane fitting residual score of point p, and σ is the root mean square error of local plane fitting based on the principal component analysis method.

[0173] S2.2.4. Based on the local quality score, construct an adaptive mapping to achieve dynamic adjustment of the sampling interval, and obtain the adaptive sampling interval ds(p):

[0174] ds(p) = ds_base(p)·(1 + γ(1 - Q(p)))

[0175] Among them, γ is the sampling adjustment coefficient, Q(p) ∈ [0,1]. When Q(p) = 1, the sampling interval is maintained. When Q(p) decreases, the sampling interval gradually increases through 1 + γ(1 - Q(p)) to cope with unstable regions.

[0176] This adaptive sampling strategy realizes the dynamic optimization of sampling parameters through the organic combination of structural feature classification and local quality assessment, which not only ensures the reliability of feature extraction but also improves the computational efficiency.

[0177] S2.2.5. Construct a SHOT local reference frame; the construction of the local reference frame also adopts different strategies according to the region type. The local reference frame refers to the local coordinate system established for each point when calculating the SHOT feature. This local coordinate system is used to transform the local region of the point cloud into a standardized reference space, making the feature description rotation-invariant, ensuring consistent feature descriptions for the same geometric structure under different perspectives, and improving the accuracy of feature matching.

[0178] Specifically, for the parts of the cylindrical section and the transition section on the cylindrical structure, the identified axial direction of the cylinder is used as the main direction (z-axis) of the reference frame, and the secondary directions (x-axis and y-axis) are determined by combining the normal vectors of the local neighborhood point sets of each point on the point cloud of the cylindrical section, thereby establishing a stable local coordinate system.

[0179] For the parts of the conical section and the transition section on the conical structure, the main direction of the local surface is calculated, and the complete reference frame (the three directions of x, y, and z) is determined using the principal component analysis method; specifically, first, the local neighborhood point sets of each point on the point cloud of the conical section are obtained (in this article, two expressions, N(p) and local neighborhood point set, are specifically used, and their concepts are exactly the same. When not emphasizing the point p, the local neighborhood point set is used, and when emphasizing the point p, N(p) is used). The principal component analysis is performed on this point set, and the eigenvector corresponding to the largest eigenvalue is used as the main direction (z-axis) of the reference frame, the eigenvector corresponding to the second largest eigenvalue is used as the secondary direction (x-axis or y-axis), and the eigenvector corresponding to the smallest eigenvalue is used as the third direction (y-axis or x-axis), thereby establishing a complete reference frame; this method can effectively capture the main change directions of the local geometric structure.

[0180] S2.2.6. Extract SHOT feature points from the point cloud. For each point p in the point cloud, the points in the point cloud are sorted based on its local quality score Q(p). Traverse the sorted point list. For the neighborhood N1(p) of the current point p, N1(p) = {q | ||q - p|| ≤ ds(p)}

[0181] where ds(p) is the adaptive sampling spacing; if there is no SHOT feature point already added in the neighborhood N1(p) of the current point p, then p is added to the SHOT feature point set.

[0182] S2.2.7. Calculate the SHOT feature descriptor SHOT(p), and calculate the feature response value S(p) for each point p in the SHOT feature point set of the point cloud, including the following steps:

[0183] (1) Based on the local reference frame established in S2.2.5, the local spherical support region of the point p is divided into 32 spatial grids, and the normal vector distribution histograms of 11 bins are calculated in each grid to obtain a 352-dimensional SHOT feature vector, and this SHOT feature vector is the SHOT feature descriptor SHOT(p).

[0184] (2) The feature vector is reorganized into a 32×11 matrix M with every 11 elements as a row for calculating the feature response value S(p).

[0185] S(p)= mean(diff(i,j))=mean(|M[i] - M[j]|)

[0186] Among them, M[i] is the 11-dimensional histogram vector of the i-th spatial grid, diff(i, j) is the histogram difference between adjacent grids, and mean(·) is the averaging operation; the response value S(p) reflects the significance of the local geometric features at point p: when S(p) is large, it indicates that the point is located in an area with obvious geometric features, such as edges or corners; when S(p) is small, it indicates that the geometric change in the area where the point is located is relatively gentle, such as a planar area. This calculation method of the feature response value provides an important basis for the subsequent adaptive adjustment of feature extraction parameters.

[0187] S2.3. Use FPFH features to extract the third-scale features and obtain the FPFH feature descriptor FPFH(p). The third-scale features are used as the lower layer of the pyramid to achieve fine geometric structure extraction.

[0188] The specific content of 2.3 includes:

[0189] S2.3.1. Take dense features, and use a larger spacing value νd at locations where the features are not obvious to improve the calculation efficiency;

[0190] Sort the points in the SHOT feature point set according to the feature response value S(p), traverse the sorted point list, check the neighborhood N2(p) of the current point, which is defined as

[0191] N2(p) = {q | ||q - p|| ≤ df(p)}

[0192] df(p) is the sampling spacing, df(p)=d, that is, the average point spacing d of the point cloud downsampling. If there is no selected FPFH feature point in the neighborhood, add p to the FPFH feature point set of this point cloud.

[0193] S2.3.2. Based on the feature response value S(p) of the SHOT feature descriptor, for each point p in the FPFH feature point set, construct an adaptive support radius

[0194] R(p) = d·[γ·H(S(p)-θ) + λ·(1-H(S(p)-θ))]

[0195] Among them, both γ and λ are radius adjustment coefficients and satisfy λ>γ>0; this enables the use of a smaller support radius γd to retain details in areas with significant geometric features (S(p)≥θ), and a larger support radius λd to improve stability in areas with insignificant features (S(p)<θ);

[0196] Obtain the support radius neighborhood of point p for use in subsequent FPFH feature descriptor calculations.

[0197] N3(p) = {q | ||q - p|| ≤ R(p)}

[0198] S2.3.3. Based on the feature response value S(p) of the SHOT feature descriptor, for each point p in the FPFH feature point set, construct an adaptive feature weight mapping function.

[0199] Wf(p) = η·S(p)

[0200] where η is the weight adjustment coefficient; the weight is proportional to the local geometric significance, ensuring that the geometric feature significant region has a greater contribution in feature calculation.

[0201] S2.3.4. Calculate the FPFH feature descriptor FPFH(p) for each point p in the FPFH feature point set.

[0202] FPFH(p) = SPFH(p) + (1 / |N3(p)|) * ∑(Wf(pk) / dk * SPFH(pk))

[0203] where FPFH(p) is the FPFH feature descriptor of point p, SPFH(p) is the local geometric feature histogram of point p, used to describe local geometric features, N3(p) is the support radius neighborhood centered at point p, |N3(p)| represents the number of points in the support radius neighborhood, dk is the Euclidean distance from neighborhood point pk to the center point p, and Wf(pk) is the feature weight of neighborhood point pk. The specific calculation process is an existing mature theory and will not be elaborated here.

[0204] Finally, for each point cloud, its feature point set is its FPFH feature point set. The feature points in S3 and S4 later are the feature points of this FPFH feature point set. These points are only used to select registration points. The specific registration features also need to use the ESF features and SHOT features calculated previously.

[0205] S3. Based on the extracted multi-scale features, perform preliminary feature fusion on the feature points of the point cloud pair; specifically, it includes the following steps:

[0206] S3.1. Calculate the feature distances of the point cloud pair, including the overall point cloud distance dist_ESF based on the global ESF feature, the distance dist_SHOT(p,q) between the corresponding local feature points based on the SHOT feature, and the distance dist_FPFH(p,q) between the corresponding local feature points based on the FPFH feature.

[0207] dist_ESF = ||ESF(P) - ESF(Q)||2

[0208] dist_SHOT(p,q) = ||SHOT(p) - SHOT(q)||2

[0209] dist_FPFH(p,q) = ||FPFH(p) - FPFH(q)||2

[0210] Among them, dist_ESF is the overall point cloud distance based on the global ESF feature, which is used to measure the global structural similarity. ESF(P) and ESF(Q) are the ESF feature descriptors of the source point cloud P and the target point cloud Q respectively. ||||2 is the Euclidean distance calculation. p and q are the corresponding feature points in the source point cloud and the target point cloud respectively. SHOT(p) and FPFH(p) are the SHOT feature descriptor and FPFH feature descriptor of the feature point p. dist_SHOT(p,q) and dist_FPFH(p,q) are the SHOT and FPFH feature distances of the feature point pair (p,q) respectively.

[0211] S3.1. Normalize the feature distances of the point cloud pair,

[0212] Among them, represents the maximum value of the SHOT feature distances of all feature point pairs, represents the maximum value of the FPFH feature distances of all feature point pairs; dist_ESF_norm is the normalized ESF feature distance, and dist_SHOT_norm(p,q) and dist_FPFH_norm(p,q) are the normalized SHOT and FPFH feature distances of the feature point pair (p,q) respectively. By dividing by their respective maximum values, all feature distances are mapped to the unified [0,1] interval, providing a comparable measurement standard for subsequent feature fusion.

[0213] S4. Perform multi-level registration on the point cloud pair; the multi-level registration proposed by the present invention realizes a progressive registration process from coarse to fine through the dynamic adjustment of the feature-dominant weights and the information transfer mechanism between levels; the multi-level includes the first-level registration, the second-level registration, and the third-level registration. After each level of registration is completed, the transformation matrix of that level is obtained, and the point cloud to be registered is registered through the transformation matrix of each level; these three levels of registration can be regarded as coarse registration, medium registration, and fine registration. The specific process of each level of registration is as follows:

[0214] (1) Determine the search space for the registration of that level.

[0215] (2) Determine the respective weights of the ESF feature distance, SHOT feature distance, and FPFH feature distance of the point cloud pair in the registration of that level.

[0216] (3) Calculate the fused feature distance in a weighted manner in the search space based on the three feature distances and their respective weights.

[0217] (4) Determine the matching point pairs based on the matching strategy that combines feature distance and region - adaptive segmentation.

[0218] (5) Conduct geometric consistency verification on the results of the matching point pairs.

[0219] (6) Calculate the transformation matrix according to the set of matching point pairs that pass the geometric consistency verification, and register the point cloud to be registered through the transformation matrix of each layer.

[0220] It includes the following steps:

[0221] S4.1. The first step of multi - level registration is coarse registration. The main goal of this stage is to obtain the initial transformation matrix, that is, the first - layer transformation matrix T_coarse, perform feature matching in the global range, and achieve the first - layer registration of the point cloud. In terms of feature configuration, the ESF global feature is the dominant one. The registration process includes the following steps:

[0222] S4.1.1. Determine the search space; in the global range, the search space for any feature point p is,

[0223] Ω1(p) = {q ∈ Q | ||q - p|| ≤ Rmax}

[0224] where Ω1(p) represents the search space for feature point p within the global range Rmax, that is, the set of points in the point cloud Q to be registered that may match p. Rmax is the maximum search radius, which is a threshold for restricting the search range;

[0225] S4.1.2. Determine the feature weight vector,

[0226]

[0227] where w1, w2, and w3 correspond to the initial weights of the ESF feature, SHOT feature, and FPFH feature respectively. 0.2 ≤ wi ≤ 0.5 to ensure the basic contributions of each feature. w1 + w2 + w3 = 1 to ensure weight normalization, and w1 > w2, w1 > w3.

[0228] In the coarse registration stage of this specific implementation, the weight configuration satisfies: w1 = 0.5, w2 = w3 = 0.25, to ensure the reliability of the overall pose estimation by increasing the weight of the global feature.

[0229] S4.1.3. Calculate the fused feature distance D_final(p,q) of the point pairs,

[0230]

[0231] Through this weighted fusion method, an adaptive combination of multi-feature distances is achieved, which not only ensures that the contribution degrees of various features do not deviate excessively but also maintains a certain adjustment space.

[0232] S4.1.4 Determine the matching point pairs based on the fused feature distance D_final(p,q) and the region-adaptive matching strategy; including,

[0233] (1) Calculate the fused distance D_final(p,q) between all feature points in the point cloud Q to be registered within the search space Ω(p) for each feature point p in the target point cloud P.

[0234] (2) Find the two points q1 and q2 with the smallest distances to the feature point p in the point cloud Q to be registered. The corresponding distances are D1 and D2, and different matching criteria are adopted according to the region where p is located.

[0235] For the feature points in the transition section, if D1 / D2 < λ1 and the cross-validation rule is required to be satisfied (the nearest neighbor of q1 in P is p), then the matching point of this feature point p in P in Q is q1; similarly, for the feature point q1 in Q, if its nearest neighbor in P is p, it is considered that the matching point of this feature point q1 in Q in P is p; if both of the above are satisfied, the matching point pair (p,q1) is retained.

[0236] For the feature points in the cylindrical section, if D1 / D2 < λ2, it is considered that the matching point of the feature point p in P in Q is q1; for the feature point q1 in Q, if its nearest neighbor in P is p, it is considered that the matching point of this feature point q1 in Q in P is p. If both of the above are satisfied, the matching point pair (p,q1) is retained.

[0237] Where λ1 < λ2 < 1, both are ratio thresholds for feature point applications; the smaller the values of λ1 and λ2, the more stringent the ratio test, that is, the higher the matching degree of the point pairs passing the test. However, too small a value will result in too few matching points, and it needs to be adjusted according to actual experience.

[0238] S4.1.5. Conduct geometric consistency verification on the matching point pair results to screen out; used to examine the distance ratio relationship between matching point pairs; calculate the geometric consistency index of any two randomly selected matching point pairs (pi, qi) and (pj, qj),

[0239] ε_ij = |||pi - pj|| / ||qi - qj|| - 1|

[0240] If ε_ij < εd, the matching point pair passes the verification; otherwise, it fails the geometric consistency verification; where, (pi,pj) is a point pair in the source point cloud, (qi, qj) is a point pair in the target point cloud, ε_ij is the distance ratio deviation, and εd is the geometric consistency threshold.

[0241] Define Q1 = Nc / Nt, where Nc is the number of matched point pairs passing through geometric consistency constraints, Nt is the total number of matched point pairs, and Q1 is the proportion of matched point pairs satisfying geometric constraints;

[0242] When Q1 < 0.6, it indicates insufficient geometric consistency. Perform the following updates and return to step 4.1.3 to recalculate D_final while accumulating the update count k;

[0243]

[0244] Among them, σpace is the weight adjustment step size, σpace = 0.05. If Q1 >= 0.6 or the cumulative update count k > 5, then no weight update is performed.

[0245] S4.1.6 Based on the filtered set of matched point pairs, use the RANSAC - SVD method to calculate the transformation matrix T_coarse (the specific implementation method is the classic transformation matrix calculation method, which will not be elaborated here).

[0246] S4.2. The second step of multi - level registration is medium registration. The main goal of this stage is to obtain the optimized transformation matrix, that is, the second transformation matrix T_optimize, to achieve the second - layer registration of the point cloud and improve the alignment accuracy of structural features. In terms of feature configuration, with SHOT features as the dominant, the registration process includes the following steps:

[0247] S4.2.1. Determine the search range based on the transformation matrix T_coarse: For any feature point p, take its position p = T_coarse·p after being transformed by the transformation matrix T_coarse as the search center, and its local search space Ω2(p) is defined as,

[0248] Ω2(p) = {q ∈ Q | ||q - p|| ≤ Rmiddle}

[0249] Rmiddle = βRmax

[0250] Among them, Rmiddle is the local search radius, β is the search radius reduction coefficient, and in this specific implementation, β takes 0.3, indicating that the search radius is reduced to 30% of the global search radius through the coefficient β.

[0251] S4.2.2. Perform weight configuration based on the weights defined in step 4.1.2, w1 = w3 = 0.25, w2 = 0.5, highlighting the dominant role of structural features.

[0252] S4.2.3. Calculate the fusion feature distance D_final in the medium registration stage based on the above local search space Ω2(p), weights, and the definitions in step 4.1.3.

[0253] S4.2.4. Calculate new matching point pairs based on the above local search space Ω2(p) and the definitions in step 4.1.4; the new matching point pairs are based on a smaller search range and weights that pay more attention to structural features, ensuring a more accurate feature point matching result at the structural feature level compared to the matching point pairs in the rough registration stage of step 4.1.4.

[0254] S4.2.5. Calculate Q1 based on the above new matching point pairs and the definitions in step 4.1.5, i.e., also require geometric consistency to be ensured in the medium registration stage. The specific determination criteria, weight update logic, and transformation matrix calculation results are the same as those in step 4.1.5.

[0255] S4.3. The third step of the multi-level registration is the fine registration. The main goal of this stage is to obtain a refined optimized transformation matrix, i.e., the third-layer transformation matrix T_fine, to achieve the third-layer registration of the point cloud, i.e., fine alignment; in terms of feature configuration, FPFH features are the mainstay, and the registration process includes the following steps:

[0256] S4.3.1. Determine the search range based on the optimized transformation matrix T_optimize:

[0257] For any feature point p, take its position p = T_optimize·p after being transformed by the transformation matrix T_optimize as the search center, and its local search space Ω3(p) is defined as,

[0258] Ω3(p) = {q ∈ Q | ||q - p|| ≤ Rmin}

[0259] Rmin = γRmax

[0260] where Rmin is the fine search radius, γ ≈ 0.1, and in this specific implementation, γ is taken as 0.1, further reducing the search range to 10% of the global range through the coefficient γ.

[0261] S4.3.2. In this fine registration stage, perform weight configuration based on the weights defined in step 4.1.2, w1 = 0.2, w2 = 0.3, w3 = 0.5; highlighting the role of local fine features.

[0262] S4.3.3. Calculate D_final in the fine registration stage based on the above local search space Ω3(p), weights, and the definitions in step 4.1.3.

[0263] S4.3.4. Define new matching point pairs based on the above local search space Ω3(p) and the definition in step 4.1.4; the new matching point pairs are based on a smaller search range and weights that are more concerned with local fine features, ensuring more accurate feature point matching results at the local fine feature level compared to the matching point pairs in the intermediate registration stage of step 4.2.4.

[0264] S4.3.5. Define local shape similarity and calculate the local shape similarity index of the matching point pairs obtained in step 4.3.4.

[0265] s_i = cos(SHOT(pi), SHOT(qi))

[0266] where SHOT(pi) and SHOT(qi) are the SHOT descriptors of the matching point pi and the point qi respectively;

[0267] Define Q2 = (1 / Nt) * ∑s_i as the overall local shape similarity score, which is obtained by calculating their cosine similarity and averaging over all matching point pairs, and Nt is the total number of matching point pairs;

[0268] When Q2 < 0.7, it indicates insufficient local similarity. Perform the following updates and return to step 4.3.3 to recalculate D_final, and at the same time accumulate the update count k.

[0269]

[0270] If Q2 >= 0.7 or the cumulative update count k > 5, then no weight update is performed. Based on the set of selected matching point pairs, use the RANSAC+SVD method to calculate the fine optimization transformation matrix T_fine.

[0271] S4.4. Obtain the registered point cloud Qp, which is calculated as Qp = T_fine × T_optimize × T_coarse × Q.

[0272] Although the preferred embodiments of the present invention have been described, those skilled in the art can make additional changes and modifications to these embodiments once they learn the basic creative concept. Therefore, the appended claims are intended to be construed as including the preferred embodiments as well as all changes and modifications that fall within the scope of the present invention.

[0273] Obviously, those skilled in the art can make various changes and variations to the present invention without departing from the spirit and scope of the present invention. Thus, if these modifications and variations of the present invention fall within the scope of the claims of the present invention and their equivalent technologies, the present invention is also intended to include these modifications and variations.

Claims

1. A multi-level feature fusion adaptive point cloud registration method, characterized in that Including: S1. Preprocess the point cloud pair; S2. Extract multi-scale features from the preprocessed point cloud pair, construct a feature pyramid through the extracted multi-scale features, and adopt a top-down guidance mechanism among the levels of the feature pyramid; S3. Perform preliminary feature fusion on the feature points of the point cloud pair based on the extracted multi-scale features; S4. Perform multi-level registration on the point cloud pair; the multi-level includes first-level registration, second-level registration, and third-level registration. After each level of registration is completed, the transformation matrix of that level is obtained, and the point cloud to be registered is registered through the transformation matrix of each level.

2. The multi-level feature fusion adaptive point cloud registration method according to claim 1, wherein The S1 includes the following steps: S1.1 Perform noise reduction processing on the point cloud, identify and remove abnormal points by analyzing the neighborhood distribution characteristics of points; S1.2 Downsample the point cloud using the voxel grid filtering algorithm; S1.3 Calculate the average point spacing d of the downsampled point cloud and use it as the reference parameter for feature extraction; S1.4 Construct a spatial index structure based on an octree for neighborhood search.

3. A multi-level feature fusion adaptive point cloud registration method according to claim 1, characterized in that, The S2 includes the following steps: S2.1 Extract the first-scale features using the ESF feature to obtain the ESF feature descriptor ESF(P). The first-scale features serve as the top layer of the pyramid. Perform silo structure recognition and region segmentation on the preprocessed point cloud to obtain cylindrical sections, conical sections, and transition sections, and verify the rationality of the region segmentation through the ESF feature; S2.2 Extract the second-scale features using the SHOT feature to obtain the SHOT feature descriptor SHOT(p); the second-scale features serve as the middle layer of the pyramid, which is used to extract local geometric structures and transfer the information of the middle layer of the feature pyramid to the bottom layer to provide guidance for the bottom layer; S2.3 Extract the third-scale features using the FPFH feature to obtain the FPFH feature descriptor FPFH(p), and the third-scale features serve as the bottom layer of the pyramid, which is used to extract fine geometric structures; For each point cloud, its set of feature points is its set of FPFH feature points, and the feature points in subsequent S3 and S4 are the feature points of this set of FPFH feature points.

4. The multi-level feature fusion adaptive point cloud registration method according to claim 3, wherein The S2.1 includes the following steps: S2.1.1 Detect the cylindrical section of the point cloud, extract the cylindrical radius r and the axial direction vector v, and use the axial direction as the z-axis direction; S2.1.2 Extract the ESF feature descriptor ESF(P) of the point cloud, and verify the reliability of the cylindrical section through the ESF feature. If the verification passes, record the cylindrical parameters, the cylindrical radius r and the axial direction vector v, and mark the points in the point cloud that conform to the characteristics of the cylindrical model as the cylindrical section and enter S2.1.4; if the verification fails, enter S2.1.3; the reliability verification criteria include two aspects: (1) whether the deviation between the D2 distribution of the ESF feature and the theoretical cylinder established based on the cylindrical radius r extracted by RANSAC is within the range of ε1, and (2) whether the deviation of the A3 component of the ESF feature from the roundness of the cross-section of the theoretical cylinder is within the range of ε2. Both aspects need to be satisfied simultaneously for the verification to pass; S2.1.3 Crop the point cloud in the z-axis direction, removing the point cloud data within δh×H near the highest and lowest points, where H is the total height of the point cloud and δh is the cropping ratio; after cropping, return to step S2.1.1 and accumulate the cropping times; when the cropping times reach the maximum number of attempts N, if the verification still fails, the point cloud quality is insufficient to extract reliable cylindrical features, and registration is exited; S2.1.

4. Further divide the point cloud into cylindrical segments, transition segments, and conical segments based on the segmented cylindrical segments, where the transition segment is the area within (-tr, tr) of the cylindrical segment boundary, and the conical segment is the remaining area.

5. The multi-level feature fusion adaptive point cloud registration method according to claim 3, wherein The above S2.2 includes the following steps: S2.2.

1. Define the structure indication function I(p), ； where p represents the points in the point cloud that have been divided, β is the sampling correction coefficient for the conical segment, μ is the sampling density coefficient for the transition segment, β∈(1, 2), μ∈(1, 2); S2.2.

2. Calculate the basic sampling spacing, ds_base(p) = d / ρ(p); ρ(p) = α·I(p); where d is the average point spacing of the point cloud, α is the basic sampling coefficient, β is the sampling correction coefficient for the conical segment, and ρ(p) is the sampling density function; S2.2.

3. Obtain the local quality score Q(p), Q(p) = (Q_d(p) + Q_n(p) + Q_r(p)) / 3; Q_d(p)= min(1, |N(p)| / N_exp); Q_n(p) = 1 - (1 / |N(p)|)∑(arccos(|n_i·n_p|) / π); Q_r(p) = exp(-σ / d); where Q_d(p) is the local point density score of point p, N(p) is the local neighborhood point set of point p, for any point p in the point cloud, its neighborhood N(p) is defined as the set of all points within a spherical space centered at p with a radius of r, N(p) = {q | ||q - p|| ≤ r}; where r is the neighborhood radius, ||·|| represents the Euclidean distance, N_exp is the expected number of neighborhood points, Q_n(p) is the normal vector consistency score of point p, n_i and n_p are the normal vectors of the neighborhood point and the center point respectively, Q_r(p) is the local plane fitting residual score of point p, and σ is the root mean square error of local plane fitting based on the principal component analysis method; S2.2.

4. Obtain the adaptive sampling spacing ds(p) based on the local quality score, ds(p) = ds_base(p)·(1 + γ(1 - Q(p))); where γ is the sampling adjustment coefficient, Q(p)∈[0, 1], when Q(p)=1, keep the sampling spacing, when Q(p) decreases, gradually increase the sampling spacing through 1 + γ(1 - Q(p)); S2.2.

5. Construct the SHOT local reference frame; For the cylindrical segment and the part of the transition segment on the cylindrical structure, use the identified cylindrical axis as the main direction of the reference frame, and then determine the secondary direction in combination with the normal vectors of the local neighborhood point sets of each point on the cylindrical segment point cloud; For the parts of the conical section and the transition section on the conical structure, the principal direction of the local surface is calculated, and the complete reference frame is determined using the principal component analysis method. Specifically, first, the local neighborhood point set of each point on the point cloud of the conical section is obtained, and the principal component analysis is performed on this point set. The eigenvector corresponding to the largest eigenvalue is used as the principal direction of the reference frame, the eigenvector corresponding to the second largest eigenvalue is used as the secondary direction, and the eigenvector corresponding to the smallest eigenvalue is used as the third direction. S2.2.

6. Extract SHOT feature points from the point cloud. For each point p in the point cloud, the points in the point cloud are sorted based on its local quality score Q(p). Traverse the sorted point list. For the neighborhood N1(p) of the current point p, N1(p) = {q | ||q - p|| ≤ ds(p)}. where ds(p) is the adaptive sampling spacing. If there are no previously added SHOT feature points in the neighborhood N1(p) of the current point p, then p is added to the SHOT feature point set. S2.2.

7. Calculate the SHOT feature descriptor SHOT(p), and calculate the feature response value S(p) for each point p in the SHOT feature point set of the point cloud, including the following steps: (1) Based on the local reference frame established in S2.2.5, the local spherical support region of point p is divided into 32 spatial grids. The normal vector distribution histogram of 11 bins is calculated in each grid, and a 352-dimensional SHOT feature vector is obtained. This SHOT feature vector is the SHOT feature descriptor SHOT(p). (2) The feature vector is reorganized into a 32×11 matrix M with every 11 elements as a row, which is used to calculate the feature response value S(p). S(p)= mean(diff(i,j))=mean(|M[i] - M[j]|); where M[i] is the 11-dimensional histogram vector of the i-th spatial grid, diff(i,j) is the histogram difference between adjacent grids, and mean(·) is the averaging operation.

6. The multi-level feature fusion adaptive point cloud registration method according to claim 3, wherein The above 2.3 includes the following steps: S2.3.

1. Sort the points in the SHOT feature point set according to the feature response value S(p). Traverse the sorted point list. Check the neighborhood N2(p) of the current point, which is defined as N2(p) = {q | ||q - p|| ≤ df(p)}; df(p) is the sampling spacing, df(p)=d, that is, the average point spacing d of the downsampled point cloud. If there are no previously selected FPFH feature points in the neighborhood, then p is added to the FPFH feature point set of this point cloud. S2.3.

2. Based on the feature response value S(p) of the SHOT feature descriptor, for each point p in the FPFH feature point set, construct an adaptive support radius R(p) = d·[γ·H(S(p)-θ) + λ·(1-H(S(p)-θ))]; where γ and λ are both radius adjustment coefficients and satisfy λ>γ>0; Obtain the support radius neighborhood of point p, which is used for the subsequent calculation of the FPFH feature descriptor. N3(p) = {q | ||q - p|| ≤ R(p)}; S2.3.

3. Based on the feature response value S(p) of the SHOT feature descriptor, for each point p in the FPFH feature point set, construct an adaptive feature weight mapping function, Wf(p) = η·S(p); where η is the weight adjustment coefficient; S2.3.

4. Calculate the FPFH feature descriptor FPFH(p) for each point p in the FPFH feature point set, FPFH(p) = SPFH(p) + (1 / |N3(p)|) * ∑(Wf(pk) / dk * SPFH(pk)); where FPFH(p) is the FPFH feature descriptor of point p, SPFH(p) is the local geometric feature histogram of point p, used to describe local geometric features, N3(p) is the support radius neighborhood centered on point p, |N3(p)| represents the number of support radius neighborhood points, dk is the Euclidean distance from neighborhood point pk to the center point p, and Wf(pk) is the feature weight of neighborhood point pk.

7. A multi-level feature fusion adaptive point cloud registration method according to claim 1, characterized in that The said S3 includes the following steps: S3.

1. Calculate the feature distances of the point cloud pair, including the overall point cloud distance dist_ESF based on the global ESF feature, the distance dist_SHOT(p,q) between the corresponding feature points of the local features based on the SHOT feature, and the distance dist_FPFH(p,q) between the corresponding feature points of the local features based on the FPFH feature, dist_ESF = ||ESF(P) - ESF(Q)|| 2； dist_SHOT(p,q) = ||SHOT(p) - SHOT(q)|| 2； dist_FPFH(p,q) = ||FPFH(p) - FPFH(q)|| 2； where dist_ESF is the overall point cloud distance based on the global ESF feature, ESF(P) and ESF(Q) are the ESF feature descriptors of the source point cloud P and the target point cloud Q respectively, ||||2 is the Euclidean distance calculation, p and q are the corresponding feature points in the source point cloud and the target point cloud respectively, SHOT(p) and FPFH(p) are the SHOT feature descriptor and the FPFH feature descriptor of feature point p respectively, and dist_SHOT(p,q) and dist_FPFH(p,q) are the SHOT and FPFH feature distances of the feature point pair (p,q); S3.

1. Normalize the feature distances of the point cloud pair, ； Among them, represents the maximum value of the SHOT feature distance for all pairs of feature points, represents the maximum value of the FPFH feature distance for all pairs of feature points; dist_ESF_norm is the normalized ESF feature distance, and dist_SHOT_norm(p,q) and dist_FPFH_norm(p,q) are the normalized SHOT and FPFH feature distances of the pair of feature points (p,q) respectively.

8. A multi-level feature fusion adaptive point cloud registration method according to claim 1, characterized in that The said S4 includes the following steps: S4.

1. Obtain the first-layer transformation matrix T_coarse, perform feature matching globally, and achieve the first-layer registration of the point cloud; S4.

2. Obtain the second-layer transformation matrix T_optimize and achieve the second-layer registration of the point cloud; S4.

3. Obtain the third-layer transformation matrix T_fine and achieve the third-layer registration of the point cloud; S4.

4. Obtain the registered point cloud Qp, and calculate Qp = T_fine×T_optimize×T_coarse×Q.

9. The multi-level feature fusion adaptive point cloud registration method according to claim 8, wherein In the said S4, the process of each layer of registration is: Determine the search space for this layer of registration; Determine the weights of the ESF feature distance, the SHOT feature distance, and the FPFH feature distance of the point cloud pair in this layer of registration respectively; Based on the three feature distances and their respective weights, calculate the fused feature distance in a weighted manner in the search space; Determine the matching point pairs based on the fused feature distance and the regionally adaptive matching strategy; Perform geometric consistency verification on the matching point pair results; Calculate the transformation matrix based on the set of matching point pairs verified by geometric consistency, and register the point cloud to be registered through the transformation matrix of each layer.

10. A multi-level feature fusion adaptive point cloud registration method according to claim 9, characterized in that, S4.1 includes the following specific steps: S4.1.

1. Determine the search space; globally, the search space for any feature point p is Ω1(p) = {q ∈ Q | ||q - p|| ≤ Rmax}; where Ω1(p) represents the search space for feature point p within the global range Rmax, and Rmax is the maximum search radius; S4.1.

2. Determine the feature weight vector, ； where w1, w2, and w3 correspond to the initial weights of the ESF feature, SHOT feature, and FPFH feature respectively, 0.2 ≤ wi ≤ 0.5, w1 + w2 + w3 = 1, and w1 > w2, w1 > w3; S4.1.

3. Calculate the fusion feature distance D_final(p,q) of the point pair: D_final(p,q) = w1×dist_ESF_norm + w2×dist_SHOT_norm(p,q) + w3×dist_FPFH_norm(p,q); S4.1.

4. Determine the matching point pairs based on the fusion feature distance D_final(p,q) and the region - adaptive matching strategy; including: (1) Calculate the fusion distance D_final(p,q) of all feature points in the point cloud Q to be registered within the search space Ω(p) for each feature point p in the target point cloud P; (2) Find the two points q1 and q2 with the smallest distances to the feature point p in the point cloud Q to be registered, and the corresponding distances are D1 and D2, and adopt different matching criteria according to the region where p is located; For the feature points in the transition section, if D1 / D2 < λ1 and the cross - validation rule is required to be satisfied, then the matching point of this feature point p in P in Q is q1; similarly, for the feature point q1 in Q, if its nearest neighbor in P is p, it is considered that the matching point of this feature point q1 in Q in P is p; if both of the above are satisfied, the matching point pair (p,q1) is retained; For the feature points in the cylindrical section, if D1 / D2 < λ2, for the feature point q1 in Q, if its nearest neighbor in P is p, it is considered that the matching point of this feature point q1 in Q in P is p; where λ1 < λ2 < 1, both are the ratio thresholds for feature point applications; S4.1.

5. Conduct geometric consistency verification on the matching point pair results; calculate the geometric consistency index of any two randomly selected pairs of matching point pairs (pi, qi) and (pj, qj), ε_ij = |||pi - pj|| / ||qi - qj|| - 1|; If ε_ij < εd, the matching point pair passes the verification, otherwise, it fails the geometric consistency verification; where, (pi, pj) is a point pair in the source point cloud, (qi, qj) is a point pair in the target point cloud, ε_ij is the distance ratio deviation, and εd is the geometric consistency threshold; Define Q1 = Nc / Nt, where Nc is the number of matched point pairs that pass the geometric consistency constraint, Nt is the total number of matched point pairs, and Q1 is the proportion of matched point pairs that satisfy the geometric constraint; When Q1 < 0.6, it indicates insufficient geometric consistency. Perform the following updates and return to step 4.1.3 to recalculate D_final, and at the same time accumulate the update count k. ； Among them, σpace is the weight adjustment step size, σpace = 0.

05. If Q1 >= 0.6 or the cumulative update count k > 5, then no weight update is performed; S4.1.6 Based on the set of selected matched point pairs, use the RANSAC-SVD method to calculate the transformation matrix T_coarse; The said S4.2 includes the following steps: S4.2.1 Determine the search range based on the transformation matrix T_coarse: For any feature point p, take its position p = T_coarse·p after being transformed by the transformation matrix T_coarse as the search center, and its local search space Ω2(p) is defined as Ω2(p) = {q ∈ Q | ||q - p|| ≤ Rmiddle}; Rmiddle = βRmax; Among them, Rmiddle is the local search radius, and β is the search radius reduction coefficient; S4.2.2 Perform weight configuration based on the weights defined in S4.1.2, w1 = w3 = 0.25, w2 = 0.5; S4.2.3 Calculate the fusion feature distance D_final in the middle registration stage based on the local search space Ω2(p), weights, and the formula defined in S4.1.3; S4.2.4 Calculate new matched point pairs based on the local search space Ω2(p) and the definition in S4.1.4; S4.2.5 Calculate Q1 based on the new matched point pairs and the definition in S4.1.

5. The specific judgment criteria, weight update logic, and transformation matrix calculation results are the same as in step 4.1.5; S4.3 includes the following steps: S4.3.1 Determine the search range based on the optimized transformation matrix T_optimize: Rmin = γRmax, γ ≈ 0.

1. Among them, Rmin is the fine search radius. For any feature point p, take its position p = T_optimize·p after being transformed by the transformation matrix T_optimize as the search center, and its local search space Ω3(p) is defined as Ω3(p) = {q ∈ Q | ||q - p|| ≤ Rmin}; S4.3.2 Perform weight configuration based on the weights defined in step 4.1.2, w1 = 0.2, w2 = 0.3, w3 = 0.5; S4.3.3 Calculate D_final in the fine registration stage based on the above local search space Ω3(p), weights, and the definition in step 4.1.3; S4.3.4 Calculate new matched point pairs based on the above local search space Ω3(p) and the definition in step 4.1.4; S4.3.

5. Define the local shape similarity and calculate the local shape similarity index of the matching point pairs obtained in step 4.3.

4. s_i = cos(SHOT(pi), SHOT(qi)); where SHOT(pi) and SHOT(qi) are the SHOT descriptors of the matching point pi and the point qi respectively; Define Q2 = (1 / Nt) * ∑s_i as the overall local shape similarity score, and Nt is the total number of matching point pairs; When Q2 < 0.7, it indicates insufficient local similarity. Perform the following updates and return to step 4.3.3 to recalculate D_final, and at the same time accumulate the update times k. ； If Q2 >= 0.7 or the cumulative update times k > 5, then no weight update is performed. Based on the set of selected matching point pairs, use the RANSAC+SVD method to calculate the fine optimization transformation matrix T_fine.

Images