A multi-period point cloud registration method and system under a ground station scene

Abstract

The application discloses a multi-period point cloud registration method and system under a ground station scene, and the method comprises the following steps: collecting point cloud data at different time points by using a ground three-dimensional laser scanner, and cutting out monitoring and non-monitoring areas; calculating local attributes of point clouds in the non-monitoring area, searching for generalized feature points, calculating advanced feature descriptors of the generalized feature points, and preliminarily registering based on the same-named points; comparing features of the same-named points in the preliminarily registered point clouds, constructing a normalized weight function, determining point position changes according to function values, sorting weight values, selecting the first N points with the smallest changes in the two-period point clouds as global invariable or slightly changed points, and taking the global invariable or slightly changed points as registration key points; based on the key points, performing IGGIII method weighted iterative reweighted least squares constraint ICP registration, and removing points with large differences to realize accurate registration; and based on the registered point clouds, building an elevation change contour line map, rendering and displaying deformation, and analyzing trends according to the density of the contour lines.

Description

Technical Field

[0001] This invention belongs to the field of three-dimensional laser scanning technology, and more specifically, relates to a multi-phase point cloud registration method and system in a ground station scenario. Background Technology

[0002] In the field of deformation monitoring technology based on 3D laser scanning point cloud data, with the increasing demands for accuracy and reliability in deformation monitoring from infrastructure construction, geological disaster prevention and other projects, traditional monitoring methods have revealed many problems that urgently need to be solved in practical applications.

[0003] Most existing deformation monitoring technologies are built upon idealized scenarios. Their core idea typically involves pre-assuming a reference surface that closely matches the target's ideal shape, and then calculating the distance between the actual point cloud and this ideal reference surface to determine the target's deformation. However, in actual engineering practice, the spatial geometry of monitored objects exhibits significant complexity and diversity. For example, large ancient buildings and complex geological slopes do not have regular geometric shapes, making it difficult to establish accurate spatial geometric equations using simple mathematical models. This results in the inability to accurately determine the ideal reference surface in real-world scenarios. In this context, the applicability of traditional monitoring methods relying on ideal reference surfaces is greatly reduced, making deformation monitoring based on the comparison of point cloud data collected in two phases an inevitable choice.

[0004] Meanwhile, the actual environmental conditions faced by deformation monitoring are extremely complex and challenging. In actual monitoring scenarios, the monitoring environment generally exhibits varying degrees of vibration and continuous deformation. For example, when monitoring building deformation in busy traffic areas, the vibrations of passing vehicles and disturbances from surrounding construction activities all affect the monitoring environment. When monitoring mountain deformation in areas prone to geological disasters, factors such as minute crustal movements and soil loosening caused by rainfall also keep the monitoring environment in a state of dynamic change. These environmental vibration factors directly cause the location of points to shift at different observation periods. Even with fixed facilities such as observation piers, the piers themselves are difficult to withstand the effects of environmental vibrations and may even vibrate or deform with the environment. This complex environmental condition seriously affects the accuracy and consistency of 3D laser scanning point cloud data, making data registration for deformation monitoring based on multi-period point clouds extremely difficult. Because the spatial location of point cloud data changes at different times, traditional methods are unable to effectively extract stable deformation features and cannot accurately assess key information such as the degree of deformation and development trend of the deformed area.

[0005] Therefore, how to efficiently and reliably find globally invariant points or points with minute changes of the same name from 3D laser scanning point cloud data in complex and ever-changing real-world environments, and use these points as reference points for point cloud registration, thereby achieving accurate registration of multi-period point cloud data; at the same time, after completing point cloud registration, how to establish a scientific and reasonable evaluation system to accurately analyze the deformation state of deformed areas, including the magnitude of deformation, deformation direction, and deformation development trend, has become a key bottleneck in the development of deformation monitoring technology based on multi-period point clouds.

[0006] Existing technologies have shortcomings in areas such as strategies for selecting globally invariant or minutely changing points, methods for accurate point cloud registration, and deformation state assessment systems, and therefore cannot effectively solve the aforementioned problems. There is an urgent need for an innovative technical solution to overcome the limitations of existing technologies and achieve high-precision, high-reliability deformation monitoring based on 3D laser scanning point cloud data, so as to meet the growing monitoring needs in practical engineering applications. Summary of the Invention

[0007] This invention aims to solve the challenges of existing deformation monitoring technologies based on 3D laser scanning point cloud data. Addressing issues such as the difficulty in determining an ideal reference surface in real-world scenarios, the complexity of the monitoring environment leading to difficulties in point cloud registration, and inaccurate deformation region assessment, it provides a multi-phase point cloud registration method and system for ground station scenarios. By selecting globally invariant or minimally changing points from non-monitoring targets for registration, and generating contour maps based on the registered point clouds to assess deformation, high-precision and reliable deformation monitoring is achieved.

[0008] To address the aforementioned deficiencies or improvement needs of existing technologies, as a first aspect of this invention, the present invention provides a multi-phase point cloud registration method in a ground station scenario, comprising:

[0009] S1. Use a terrestrial 3D laser scanner to collect point cloud data at different times, and crop out the monitoring area, leaving the rest as the non-monitoring area; calculate the local attributes of the point cloud in the non-monitoring area and search for generalized feature points; for the searched generalized feature points, calculate their high-level feature descriptors; obtain the preliminary registration point cloud based on the corresponding points with the same descriptor.

[0010] S2. Compare the features of the corresponding points in the preliminary registration point cloud, construct the normalized weight function of the overall feature descriptor, determine the point position change based on the function value, sort the weight values, and select the N points with the smallest changes in the two point clouds as global invariant points or points with small changes, which are used as key points for point cloud registration.

[0011] S3. Based on the selected key points, the ICP registration method is constrained by the iteratively reweighted least squares method with weights determined by the IGGIII method. The weight function factor is determined by calculating the Euclidean distance between corresponding points of the registered point clouds. Points with large distance differences are removed during the iteration process, and the truly invariant points are retained. Furthermore, the spatial transformation relationship between the two-phase point clouds is calculated to achieve precise registration.

[0012] S4. Based on the registered point clouds, an elevation change contour map of the deformation area is established, and different distances represented by the contour map are rendered to display the deformation situation of the monitoring area. The deformation trend is analyzed according to the density of the contour lines.

[0013] Furthermore, in S1, the specific method for calculating the local attributes of the point cloud in the non-monitoring area and searching for generalized feature points is as follows:

[0014] The local attributes of the point cloud include information such as the local point cloud normal vector V, the local point cloud principal curvature K, the local point cloud density M, and the local point cloud intensity I.

[0015] The calculation method of the local normal vector V is as follows:

[0016] For any point P in the point cloud of the non-monitoring area, its neighborhood point set can be expressed as: N(P) = {P i |‖P - P i ‖ < r}, where r is the neighborhood radius and i is the serial number of the points in the neighborhood point set; then the centroid of the corresponding neighborhood point set of this point is:

[0017]

[0018] All points are translated to the centroid to obtain the translated point set N ′ (P) = {P i ′ |P i ∈ N(P)}, where

[0019] The covariance matrix of the translated point set is:

[0020]

[0021] The covariance matrix is subjected to eigenvalue decomposition:

[0022] Cv i = λ i v i

[0023] Among them, the dimension of the covariance matrix C is 3×3, λ1, λ2, and λ3 are the eigenvalues of the covariance matrix, and λ1 ≤ λ2 ≤ λ3; v1, v2, and v3 are the corresponding eigenvectors;

[0024] Then the local normal vector V = v1;

[0025] The calculation method of the local point cloud principal curvature K is as follows:

[0026]

[0027] The calculation methods of the local point cloud density M and the local point cloud intensity I are as follows:

[0028] For any point P in the non-monitoring area point cloud and its neighborhood point set N(P) = {P i |‖P - P i ‖ < r}, there are m points, then the calculation expression of its local point cloud density M is:

[0029]

[0030] In the formula, π is the pi;

[0031] The calculation expression of its local point cloud intensity I is:

[0032]

[0033] Through the above calculation results, search for corner points, inflection points, wall points, high-reflection intensity points, and high-density points in the non-monitoring area as the generalized feature points.

[0034] Furthermore, for the generalized feature points searched in S1, the specific method for calculating their high-level feature descriptors is as follows:

[0035] Calculate the FPFH feature and SHOT feature descriptors of the generalized feature points;

[0036] The FPFH feature descriptor is a fast point feature histogram, which describes the angle and distance distribution of points and points in the neighborhood. Its statistical principle is:

[0037] For any point P and its neighborhood point Q in the generalized feature area point cloud, their respective normal vectors are V P and V Q , then there are three angle and one distance relationships between the two points;

[0038] The vector between the two points is Normalize d:

[0039]

[0040] The included angle between the normal vector V P and d is:

[0041]

[0042] In V Q The projection angle on is:

[0043]

[0044] V P and V Q The included angle is:

[0045] θ = arccos(V) P ·V Q )

[0046] These relationships are statistically represented as a simplified histogram (SPFH). Then, the SPFH values ​​for each neighboring point are weighted and summed to obtain the FPFH features.

[0047]

[0048] In the formula, N k It is the number of SPFHs in the neighborhood; W PQ It is the weight between point P and its neighboring point Q, reflecting the contribution of the SPFH feature of neighboring point Q to the calculation of the FPFH feature of point P.

[0049] Furthermore, the SHOT feature descriptor specifically includes:

[0050] The SHOT feature descriptor is a local feature descriptor used for 3D point cloud data, capable of capturing the local geometric and texture features of the point cloud surface. The SHOT descriptor generates stable features by analyzing the geometric and color information within the local neighborhood of a point; its statistical principle is as follows:

[0051] A local reference frame is constructed for each keypoint, making the features invariant to rotation; the x, y, and z axes of the local reference frame are determined by calculating the covariance matrix around the keypoint and performing eigenvalue decomposition.

[0052] The neighborhood is divided into multiple three-dimensional cubes, and the normal vector direction and color information within each cube are statistically analyzed to construct a descriptor histogram.

[0053] In each cube, the normal vector direction and color information are statistically analyzed to construct a descriptor histogram; the robustness to point cloud density and noise is enhanced by normalizing the histogram.

[0054] Furthermore, the specific method for constructing the normalized weight function of the overall feature descriptor in S2 is as follows:

[0055] For any point (x) with index j j ,y j ,z j The local descriptive features in the point cloud measured in the k-th period are respectively Where t is the number of features, the following feature vectors are constructed. This represents the comprehensive local descriptive features of point j in the point cloud measured in period k:

[0056]

[0057] Then, in the two measurements, the residual vectors for period k and period (k+1) are expressed as:

[0058]

[0059] Each term in the residual vector at the current point is weighted and normalized to obtain the normalized weight function of the overall feature descriptor as follows:

[0060]

[0061] In the formula, This represents the normalized weight function value of the overall feature descriptor obtained after weighting and normalizing each term of the residual vector at the residual points; a u It is the weight when the u-th term in the residual vector is weighted, and t represents the number of features in the residual vector.

[0062] Furthermore, the iterative reweighted least squares method constraining the ICP registration method of the IGGIII legal weights in S3 is specifically as follows:

[0063] For the two phases of point cloud P s and P t In the t point cloud registration key points, P t It is the target point cloud, P s It is a source point cloud, and we need to register the source point cloud to the target point cloud; therefore, we take point cloud P as an example. t For reference, calculate the point cloud P. s With reference point cloud P t The Euclidean distance between the corresponding registration keypoints is given by the formula:

[0064]

[0065] In the formula, (X s ,Y s Z s (Source point cloud P) s Key point coordinates in (X) t ,Y t Z t ) is the target point cloud P t The corresponding key point coordinates in the image;

[0066] After calculating the distances between j keypoint pairs, the distances are normalized to obtain the normalized distance residual value V.j The calculation method is shown in the following formula:

[0067]

[0068] In the formula, mind E and maxd E These are the maximum and minimum values ​​of the registration keypoint spacing, respectively; ε is a small amount that prevents the normalized distance residual from being zero.

[0069] Using the above normalized distance residual value V j As the basis for determining the weights in the IGGⅢ method, the formula for calculating the weight function factor is:

[0070]

[0071] In the formula, w j It is the weighting factor, used to determine the weights of corresponding data points in the iterative reweighted least squares method; V j σ is the normalized distance residual between corresponding points in the above-mentioned registered point cloud, reflecting the degree of spatial difference between corresponding points during the point cloud registration process; σ is the mean square error of the residual, reflecting the degree of dispersion of the residual; k0, k1 are the set threshold parameters.

[0072] During the process of finding corresponding points using the ICP algorithm, the reweighted iterative least squares constraint can gradually eliminate registration points with large distance differences, while retaining the truly invariant points in the two point clouds.

[0073] Finally, the spatial transformation relationship between the two point clouds is calculated:

[0074]

[0075] In the formula, ρ is the scaling parameter; R is the rotation matrix; by calculating the rotation matrix and scaling coefficient through the correspondence, the precise registration of the two phases of point cloud data can be achieved.

[0076] Furthermore, the reweighted iterative least squares constraint in S3 is specifically as follows:

[0077] The weighted total least squares regression model is as follows:

[0078]

[0079] The parameter estimation criteria are:

[0080]

[0081] In the formula, Y represents the observation vector; A represents the coefficient matrix, and E represents the coefficient matrix. A The error matrix represents the coefficient matrix A; X represents the parameter vector to be estimated; vecE A This indicates that matrix EA The vector obtained by straightening the columns; Q,Q A Represents the weight matrix; P, P0, P x It is a matrix related to observation accuracy, point cloud characteristics, etc.; P A The weight matrix representing the error of the coefficient matrix, and Q A They are inverse matrices;

[0082] By minimizing To estimate the parameter X, the observation error and coefficient matrix error are minimized in a weighted sense.

[0083] Furthermore, the specific method for establishing the contour map of the elevation change of the deformed region in S4 is as follows:

[0084] Based on the registered point cloud of the monitoring area, taking one phase as the benchmark, the elevation difference of the corresponding points in the point cloud of the other phase is projected onto the image, and the pixel value is normalized to the elevation difference range of 0 to 255; the pixel position of the image represents the planar position of the point cloud, and the pixel value on the image represents the deformation of the point cloud at the corresponding position.

[0085] As a second aspect of the present invention, a multi-phase point cloud registration system for a ground station scenario is provided, comprising:

[0086] The data preprocessing and preliminary registration unit is used to collect point cloud data at different times using a terrestrial 3D laser scanner, and to crop out the monitoring area, with the remainder as the non-monitoring area; to calculate the local attributes of the point cloud in the non-monitoring area and search for generalized feature points; to calculate the high-level feature descriptors of the searched generalized feature points; and to obtain the preliminary registration point cloud based on the corresponding points with the same descriptors.

[0087] The key point screening unit is used to compare the features of the corresponding points in the preliminary registration point cloud, construct the normalized weight function of the overall feature descriptor, determine the point position change based on the function value, sort the weight values, and select the N points with the smallest change in the two point clouds as global invariant points or points with small changes, which serve as key points for point cloud registration.

[0088] The point cloud fine registration unit is used to constrain the ICP registration method based on the selected key points using the iterative reweighted least squares method with IGGIII weights. The weight function factor is determined by calculating the Euclidean distance between corresponding points in the registered point cloud. During the iteration process, points with large distance differences are eliminated, and truly invariant points are retained. Then, the spatial transformation relationship between the two point clouds is calculated to achieve fine registration.

[0089] The deformation monitoring and analysis unit is used to establish a contour map of elevation changes in the deformation area based on the registered point cloud, render different distances represented by the contour map, display the deformation of the monitoring area, and analyze the deformation trend based on the density of the contour lines.

[0090] As a third aspect of the present invention, a computer-readable storage medium is also provided, on which a computer program is stored, characterized in that the computer program is processed by any step of the multi-phase point cloud registration method in the ground station scenario described herein.

[0091] In summary, compared with the prior art, the above-described technical solutions conceived by this invention can achieve the following beneficial effects:

[0092] 1. The multi-phase point cloud registration method for ground station scenarios of the present invention divides the area into monitored and non-monitored regions, calculates local attributes and searches for generalized feature points in the non-monitored regions, and then calculates high-level feature descriptors for preliminary registration. This effectively solves the problem of directly determining the correspondence between point clouds in complex real-world scenarios. Preliminary registration lays the foundation for subsequent precise registration, enabling point clouds from different phases to match in approximate locations. This overcomes the disorder of the initial state of point cloud data and provides a feasible data basis for subsequent selection of stable points and fine registration, improving the efficiency and accuracy of point cloud registration.

[0093] 2. The multi-phase point cloud registration method for ground station scenarios of the present invention, by comparing the features of corresponding points in the initially registered point cloud, constructs a normalized weight function to screen globally invariant points or points with minor changes, thus addressing the point offset problem caused by factors such as vibration in the monitoring environment. These screened key points are less affected by environmental interference and can serve as reliable registration benchmarks, enabling stable correspondences to be found even in complex environments. This greatly enhances the reliability of point cloud registration, effectively avoids the accumulation of registration errors caused by environmental factors, and provides strong support for accurate deformation assessment.

[0094] 3. The multi-stage point cloud registration method for ground station scenarios of this invention achieves precise registration of two-stage point clouds by employing an iterative reweighted least squares method constrained by IGGIII weighting for ICP registration, and generates a deformation monitoring contour map based on the registered point cloud. Precise registration ensures accurate spatial correspondence of the point clouds, while the deformation monitoring contour map visually displays the deformation situation, facilitating the analysis of deformation trends. This entire process achieves complete and accurate monitoring from point cloud registration to deformation analysis, meeting the needs of practical engineering for high-precision and high-reliability deformation monitoring. Attached Figure Description

[0095] Figure 1 This is a flowchart of a multi-phase point cloud registration method in a ground station scenario according to an embodiment of the present invention;

[0096] Figure 2 This is a detailed flowchart of an embodiment of the present invention;

[0097] Figure 3 This is a system unit diagram of an embodiment of the present invention. Detailed Implementation

[0098] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention. Furthermore, the technical features involved in the various embodiments of this invention described below can be combined with each other as long as they do not conflict with each other.

[0099] Example 1

[0100] Please refer to Figure 1 This embodiment 1 provides a multi-phase point cloud registration method in a ground station scenario, including:

[0101] S1. Use a terrestrial 3D laser scanner to collect point cloud data at different times, and crop out the monitoring area, leaving the rest as the non-monitoring area; calculate the local attributes of the point cloud in the non-monitoring area and search for generalized feature points; for the searched generalized feature points, calculate their high-level feature descriptors; obtain the preliminary registration point cloud based on the corresponding points with the same descriptor.

[0102] S2. Compare the features of the corresponding points in the preliminary registration point cloud, construct the normalized weight function of the overall feature descriptor, determine the point position change based on the function value, sort the weight values, and select the N points with the smallest changes in the two point clouds as global invariant points or points with small changes, which are used as key points for point cloud registration.

[0103] S3. Based on the selected key points, the ICP registration method is constrained by the iterative reweighted least squares method with IGGIII weights. The weighting factor is determined by calculating the Euclidean distance between corresponding points in the registered point cloud. During the iteration process, points with large distance differences are eliminated, and truly invariant points are retained. Then, the spatial transformation relationship between the two point clouds is calculated to achieve fine registration.

[0104] S4. Based on the registered point cloud, establish a contour map of elevation changes in the deformed area, render the different distances represented by the contour map, display the deformation of the monitored area, and analyze the deformation trend based on the density of the contour lines.

[0105] Please refer to Figure 2 This embodiment 1 further elaborates on the above steps:

[0106] (1) Point cloud data preprocessing and preliminary registration

[0107] 1.1 Delineation of monitoring and non-monitoring areas

[0108] Collect point cloud data at different times using a terrestrial three-dimensional laser scanner. Process the original point cloud data using point cloud data processing software such as CloudCompare, cut out the deformation monitoring area, and use the remaining part as the non-monitoring area; the monitoring area is used for deformation analysis after registration in step S4, and the non-monitoring area is used for subsequent steps of S1.

[0109] The original point cloud data refers to the three-dimensional point cloud collected by a terrestrial three-dimensional laser scanner. Its general characteristics are that the point density is higher for points closer to the scanner and lower for points farther from the scanner. In addition, the original point cloud data can be the original point cloud data collected at a single station or the data after stitching the original point cloud data collected at multiple stations within the same time.

[0110] 1.2 Search for generalized feature points

[0111] Within the non-monitoring area, calculate the local properties of the point cloud and search for generalized feature points. The local properties of the point cloud include information such as the local point cloud normal vector V, the local point cloud principal curvature K, the local point cloud density M, and the local point cloud intensity I.

[0112] The calculation method of the local normal vector V is as follows:

[0113] For any point P in the point cloud of the non-monitoring area, its neighborhood point set can be expressed as: N(P) = {P i |‖P - P i ‖ < r}, where r is the neighborhood radius; then the centroid of the corresponding neighborhood point set of this point is:

[0114]

[0115] Translate all points to the centroid to obtain the translated point set N ′ (P) = {P i ′ |P i ∈ N(P)}, where

[0116] The covariance matrix of the translated point set is:

[0117]

[0118] Perform eigenvalue decomposition on the covariance matrix:

[0119] Cv i = λ i v i

[0120] Among them, the covariance matrix C has a dimension of 3×3, λ1, λ2, and λ3 are the eigenvalues of the covariance matrix, and λ1 ≤ λ2 ≤ λ3; v1, v2, and v3 are the corresponding eigenvectors;

[0121] Then the local normal vector V = v1;

[0122] The calculation method of the local point cloud principal curvature K is as follows:

[0123]

[0124] The calculation methods of the local point cloud density M and the local point cloud intensity I are as follows:

[0125] For any point P in the non-monitoring area point cloud and its neighborhood point set N(P) = {P i |‖P - P i ‖ < r}, there are m points in total, then the calculation expression of its local point cloud density M is:

[0126]

[0127] In the formula, π is the pi;

[0128] The calculation expression of its local point cloud intensity I is:

[0129]

[0130] Through the above calculation results, search for corner points, inflection points, wall points, high-reflection intensity points, and high-density points in the non-monitoring area as the generalized feature points.

[0131] 1.3 Calculate feature descriptors and match

[0132] Calculate the FPFH feature and the SHOT feature descriptor of the generalized feature points;

[0133] The FPFH feature descriptor is a fast point feature histogram, which describes the angle and distance distribution of points and points in the neighborhood. Its statistical principle is:

[0134] For any point P in the generalized feature area point cloud and its neighborhood point Q, their respective normal vectors are V P and V Q , then there are three angle and one distance relationships between the two points;

[0135] The vector between the two points is Normalize d:

[0136]

[0137] The included angle between the normal vector V P and d is:

[0138]

[0139] In V Q The projection angle on is:

[0140]

[0141] V P and V Q The included angle is:

[0142] θ = arccos(V) P ·V Q )

[0143] These relationships are statistically represented as a simplified histogram (SPFH). Then, the SPFH values ​​for each neighboring point are weighted and summed to obtain the FPFH features.

[0144]

[0145] In the formula, N k It is the number of SPFHs in the neighborhood; W PQ It is the weight between point P and its neighboring point Q, reflecting the contribution of the SPFH feature of neighboring point Q to the calculation of the FPFH feature of point P.

[0146] The SHOT (Signature of Histograms of Orientations) feature descriptor is a local feature descriptor for 3D point cloud data, capable of capturing local geometric and textural features of the point cloud surface. The SHOT descriptor generates stable features by analyzing geometric and color information within the local neighborhood of a point. Its statistical principle is as follows:

[0147] A local reference frame is constructed for each keypoint, making the features rotation-invariant. The x, y, and z axes of the local reference frame are determined by calculating the covariance matrix around the keypoint and performing eigenvalue decomposition.

[0148] The neighborhood is divided into multiple 3D angle cubes (bins), and the normal vector direction and color information within each cube are statistically analyzed to construct a descriptor histogram;

[0149] Within each cube, the normal vector direction and color information are statistically analyzed to construct a descriptor histogram. Normalization of the histogram enhances robustness to point cloud density and noise.

[0150] Due to the local reference frame constructed by the SHOT features, the spherical support region is divided into 32 small regions: 2 in the radial direction, 2 in the latitudinal direction, and 8 in the longitudinal direction. A local histogram is constructed for each small region. Voting is performed based on the cosine value of the angle between the normal vector of the feature point and the Z-axis of the local reference frame. Each local histogram is divided into 11 cubes. Therefore, the SHOT descriptor is 32×11=352-dimensional.

[0151] After calculating the FPFH and SHOT features of points within the generalized feature region, the feature descriptors of each region are matched to establish a preliminary registration relationship between the two point clouds.

[0152] (2) Filter points that remain unchanged globally or have minor changes

[0153] By comparing the features of the corresponding points in the preliminary registered point cloud, a normalized weight function of the overall feature descriptor is constructed. The changes in the points are determined based on the value of the normalized weight function. The weight values ​​are sorted, and the top N points with the smallest changes in the two point clouds are selected as globally unchanged points or points with minor changes, and used as key points for point cloud registration.

[0154] For any point (x) with index j j ,y j ,z j The local descriptive features in the point cloud measured in the k-th period are respectively Where t is the number of features, the following feature vectors are constructed. This represents the comprehensive local descriptive features of point j in the point cloud measured in period k:

[0155]

[0156] Then, in the two measurements, the residual vectors for period k and period (k+1) are expressed as:

[0157]

[0158] Each term in the residual vector at the current point is weighted and normalized to obtain the normalized weight function of the overall feature descriptor as follows:

[0159]

[0160] In the formula, This represents the normalized weight function value of the overall feature descriptor obtained after weighting and normalizing each term of the residual vector at the residual points; a u It is the weight when the u-th term in the residual vector is weighted, and t represents the number of features in the residual vector.

[0161] (3) Fine registration of point clouds in two phases

[0162] For the two phases of point cloud Ps and P t After the operations S1-S2, corresponding globally invariant points or points of slight change are obtained, but slight differences still exist between the corresponding feature points. To improve the robustness of registration, an iterative reweighted least squares method based on IGGIII weights is used to constrain the ICP registration algorithm. The weighted global least squares regression model is as follows:

[0163]

[0164] The parameter estimation criteria are:

[0165]

[0166] In the formula, Y represents the observation vector; A represents the coefficient matrix, and E represents the coefficient matrix. A The error matrix represents the coefficient matrix A; X represents the parameter vector to be estimated; vecE A This indicates that matrix E A The vector obtained by straightening the columns; Q,Q A Represents the weight matrix; P, P0, P x It is a matrix related to observation accuracy, point cloud characteristics, etc.; P A The weight matrix representing the error of the coefficient matrix, and Q A They are inverse matrices;

[0167] By minimizing To estimate the parameter X, the observation error and coefficient matrix error are minimized in a weighted sense.

[0168] For the two phases of point cloud P s and P t In the t point cloud registration key points, P t It is the target point cloud, P s It is a source point cloud, and we need to register the source point cloud to the target point cloud; therefore, we take point cloud P as an example. t For reference, calculate the point cloud P. s With reference point cloud P t The Euclidean distance between the corresponding registration keypoints is given by the formula:

[0169]

[0170] In the formula, (X s ,Y s Z s (Source point cloud P) s Key point coordinates in (X) t ,Y t Z t ) is the target point cloud P t The corresponding key point coordinates in the image;

[0171] After calculating the distances between j keypoint pairs, the distances are normalized to obtain the normalized distance residual value V. j The calculation method is shown in the following formula:

[0172]

[0173] In the formula, min d E and max d E These are the maximum and minimum values ​​of the registration keypoint spacing, respectively; ε is a small amount that prevents the normalized distance residual from being zero.

[0174] Using the normalized distance described above as the basis for weighting in the IGGⅢ method, the formula for calculating the weight factor is as follows:

[0175]

[0176] In the formula, w j It is the weighting factor, used to determine the weights of corresponding data points in the iterative reweighted least squares method; V j σ is the normalized distance residual between corresponding points in the above-mentioned registered point cloud, reflecting the degree of spatial difference between corresponding points during the point cloud registration process; σ is the mean square error of the residual, reflecting the degree of dispersion of the residual; k0, k1 are the set threshold parameters.

[0177] During the process of finding corresponding points using the ICP algorithm, the reweighted iterative least squares constraint can gradually eliminate registration points with large distance differences, while retaining the truly invariant points in the two point clouds.

[0178] Finally, the spatial transformation relationship between the two point clouds is calculated:

[0179]

[0180] In the formula, ρ is the scaling parameter; R is the rotation matrix; by calculating the rotation matrix and scaling coefficient through the correspondence, the precise registration of the two phases of point cloud data can be achieved.

[0181] (4) Deformation monitoring and analysis

[0182] Contour map generation for deformed areas: Based on the registered point cloud of the monitoring area, using one phase as a baseline, the elevation difference of corresponding points in the point cloud of the other phase is projected onto the image. The pixel values ​​are normalized to the elevation difference range of 0-255. The pixel position in the image represents the planar position of the point cloud, and the pixel value on the image represents the deformation of the point cloud at the corresponding position.

[0183] Based on the generated contour map, the deformation of the monitored area can be determined. Specifically, normalized pixel values ​​are converted into color gradients using a specified color level for intuitive rendering. Areas with dense elevation changes are analyzed; these areas correspond to dense contour lines, indicating significant deformation in the two point cloud phases and requiring close monitoring.

[0184] Example 2

[0185] Please refer to Figure 3 This embodiment 2 provides a multi-phase point cloud registration system for a ground station scenario, including:

[0186] The data preprocessing and preliminary registration unit is used to collect point cloud data at different times using a terrestrial 3D laser scanner, and to crop out the monitoring area, with the remainder as the non-monitoring area; to calculate the local attributes of the point cloud in the non-monitoring area and search for generalized feature points; to calculate the high-level feature descriptors of the searched generalized feature points; and to obtain the preliminary registration point cloud based on the corresponding points with the same descriptors.

[0187] The key point screening unit is used to compare the features of the corresponding points in the preliminary registration point cloud, construct the normalized weight function of the overall feature descriptor, determine the point position change based on the function value, sort the weight values, and select the N points with the smallest change in the two point clouds as global invariant points or points with small changes, which serve as key points for point cloud registration.

[0188] The point cloud fine registration unit is used to constrain the ICP registration method based on the selected key points using the iterative reweighted least squares method with IGGIII weights. The weight function factor is determined by calculating the Euclidean distance between corresponding points in the registered point cloud. During the iteration process, points with large distance differences are eliminated, and truly invariant points are retained. Then, the spatial transformation relationship between the two point clouds is calculated to achieve fine registration.

[0189] The deformation monitoring and analysis unit is used to establish a contour map of elevation changes in the deformation area based on the registered point cloud, render different distances represented by the contour map, display the deformation of the monitoring area, and analyze the deformation trend based on the density of the contour lines.

[0190] Example 3

[0191] This embodiment 3 also provides a computer-readable storage medium storing a computer program. When the computer program is executed by a processor, it can implement any step of a multi-phase point cloud registration method in a ground station scenario.

[0192] The computer-readable storage medium may include various media capable of storing program code, such as USB flash drives, portable hard drives, read-only memory (ROM), random access memory (RAM), magnetic disks, or optical disks.

[0193] For a description of the computer-readable storage medium provided in this application, please refer to the above method embodiments; further details will not be repeated here.

[0194] Those skilled in the art will readily understand that the above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A multi-phase point cloud registration method in a ground station scenario, characterized in that, include: S1. Use a terrestrial 3D laser scanner to collect point cloud data at different times, and crop out the monitoring area, leaving the rest as the non-monitoring area; calculate the local attributes of the point cloud in the non-monitoring area and search for generalized feature points; for the searched generalized feature points, calculate their high-level feature descriptors; obtain the preliminary registration point cloud based on the corresponding points with the same descriptor. S2. Compare the features of corresponding points in the preliminary registered point cloud, construct a normalized weight function for the overall feature descriptor, determine the point location changes based on the function values, sort the weight values, and select the points with the smallest changes between the two point cloud periods. These points, either globally invariant or with minor variations, are used as key points for point cloud registration. S3. Based on the selected key points, the ICP registration method is constrained by the iterative reweighted least squares method with IGGIII weights. The weighting factor is determined by calculating the Euclidean distance between corresponding points in the registered point cloud. During the iteration process, points with large distance differences are eliminated, and truly invariant points are retained. Then, the spatial transformation relationship between the two point clouds is calculated to achieve fine registration. S4. Based on the registered point cloud, establish a contour map of elevation changes in the deformed area, render the different distances represented by the contour map, display the deformation of the monitored area, and analyze the deformation trend based on the density of the contour lines.

2. The multi-phase point cloud registration method in a ground station scenario according to claim 1, characterized in that, The specific method for calculating the local attributes of the point cloud in the non-monitored area and searching for generalized feature points in S1 is as follows: The local attributes of the point cloud include the local point cloud normal vector. Local point cloud principal curvature Local point cloud density and local point cloud intensity Information including; The local point cloud normal vector The calculation method is as follows: For any point in the point cloud of the non-monitored area Its neighborhood point set is represented as: ,in It is the neighborhood radius. If the index of a point is a point within its neighborhood set, then the centroid of that point's corresponding neighborhood set is: ， Translate all points to the centroid to obtain the translated point set. ,in ; The covariance matrix of the translated point set is: ， Perform eigenvalue decomposition on the covariance matrix: ， Wherein, the covariance matrix Its dimensions are 3×3. , , These are the eigenvalues ​​of the covariance matrix, and ; , , It is the corresponding feature vector; Then the local point cloud normal vector ; The principal curvature of the local point cloud The calculation method is as follows: ， The local point cloud density and local point cloud intensity The calculation method is as follows: For any point in the point cloud of the non-monitored area and its neighborhood point set There are a total of If there are 1 point, then its local point cloud density is... The calculation expression is: ， In the formula, Pi; Its local point cloud intensity The calculation expression is: ， Based on the above calculation results, corner points, inflection points, wall points, high reflectivity points, and high density points in the non-monitored area are searched as the generalized feature points.

3. The multi-phase point cloud registration method in a ground station scenario according to claim 1, characterized in that, The specific method for calculating the high-level feature descriptor of the searched generalized feature points in S1 is as follows: Calculate the FPFH and SHOT feature descriptors for the generalized feature points; The FPFH feature descriptor is a fast point feature histogram, which describes the angle and distance distribution between a point and its neighbors. Its statistical principle is as follows: For any point in the generalized feature region point cloud and its neighboring points Their respective normal vectors are and If so, then there are three angles and one distance relationship between the two points; The vector between two points is ,Will Normalization: ， normal vector and The included angle is: ， exist The projection angle on is: ， and The included angle is: ， These relationships are statistically represented as a simplified histogram (SPFH). Then, the SPFH values ​​for each neighboring point are weighted and summed to obtain the FPFH features. ， In the formula, It is the number of SPFHs in the neighborhood; It is a point With neighboring points The weights between them reflect the neighborhood points of Features at calculation points of The degree of contribution when considering features.

4. The multi-phase point cloud registration method in a ground station scenario according to claim 3, characterized in that, The SHOT feature descriptor is specifically: The SHOT feature descriptor is a local feature descriptor used for 3D point cloud data, capable of capturing the local geometric and texture features of the point cloud surface. The SHOT descriptor generates stable features by analyzing the geometric and color information within the local neighborhood of a point; its statistical principle is as follows: A local reference frame is constructed for each keypoint, making the features invariant to rotation; the x, y, and z axes of the local reference frame are determined by calculating the covariance matrix around the keypoint and performing eigenvalue decomposition. The neighborhood is divided into multiple three-dimensional cubes, and the normal vector direction and color information within each cube are statistically analyzed to construct a descriptor histogram. In each cube, the normal vector direction and color information are statistically analyzed to construct a descriptor histogram; the robustness to point cloud density and noise is enhanced by normalizing the histogram.

5. The multi-phase point cloud registration method in a ground station scenario according to claim 1, characterized in that, The specific method for constructing the normalized weight function of the overall feature descriptor in S2 is as follows: For the serial number is any point In the The local descriptive features in the period measurement point cloud are respectively ,in If the number of features is [number], then construct the following feature vectors. , indicating that the serial number is The point at the In the measured point cloud, the comprehensive situation of local descriptive features at this point: ， Then, in the two measurements, the residual vectors for period k and period (k+1) are expressed as: ， Each term in the residual vector at the current point is weighted and normalized to obtain the normalized weight function of the overall feature descriptor as follows: ， In the formula, This represents the normalized weight function value of the overall feature descriptor obtained after weighting and normalizing each term of the residual vector at the residual points; For the th in the residual vector The weights when the items are weighted. This represents the number of features in the residual vector.

6. The multi-phase point cloud registration method in a ground station scenario according to claim 1, characterized in that, The iterative reweighted least squares method constraining the ICP registration method of the IGGIII legal weights in S3 is specifically as follows: For the two phases of point cloud and In Register key points in a point cloud, among which... It is a target point cloud. It is a source point cloud, and we need to register the source point cloud to the target point cloud; therefore, we use the point cloud as the reference. For reference, calculate the point cloud. With reference point cloud The Euclidean distance between the corresponding registration keypoints is given by the formula: ， In the formula, Source Point Cloud Key point coordinates in For target point cloud The corresponding key point coordinates in the image; Seeking After calculating the distance between each keypoint pair, the distance is normalized to obtain the normalized distance residual value. The calculation method is shown in the following formula: ， In the formula, and These are the maximum and minimum values ​​of the registration keypoint spacing, respectively; It is a tiny quantity that prevents the normalized distance residual from being zero; Using the above normalized distance residual value As the basis for determining the weights in the IGGⅢ method, the formula for calculating the weight function factor is: ， In the formula, It is a weighting factor used to determine the weight of the corresponding data point in the iterative reweighted least squares method; It is the normalized distance residual value between corresponding points in the above-mentioned registered point cloud, which reflects the degree of spatial difference between corresponding points in the point cloud registration process; It is the standard error of the residuals, reflecting the degree of dispersion of the residuals; It is the set threshold parameter; During the process of finding corresponding points using the ICP algorithm, the reweighted iterative least squares constraint can gradually eliminate registration points with large distance differences, while retaining the truly invariant points in the two point clouds. Finally, the spatial transformation relationship between the two point clouds is calculated: ， In the formula, For scaling parameters; The rotation matrix is ​​used; by calculating the rotation matrix and scaling factor through the correspondence, the precise registration of the two phases of point cloud data can be achieved.

7. The multi-phase point cloud registration method in a ground station scenario according to claim 6, characterized in that, The reweighted iterative least squares constraint in S3 is specifically as follows: The weighted total least squares regression model is as follows: ， The parameter estimation criteria are: ， In the formula, Represents the observation vector; Represents the coefficient matrix. Represents the coefficient matrix The error matrix; Represents the vector of parameters to be estimated; Indicates the matrix The vector obtained by straightening the columns; Represents the weight matrix; It is a matrix related to observation accuracy and point cloud characteristics; The weight matrix representing the error of the coefficient matrix, and They are inverse matrices; By minimizing To estimate parameters This minimizes the observation error and coefficient matrix error in a weighted sense.

8. The multi-phase point cloud registration method in a ground station scenario according to claim 1, characterized in that, The specific method for establishing the contour map of the elevation change of the deformation region in S4 is as follows: Based on the registered monitoring area point cloud, taking one phase as the benchmark, the elevation difference of the corresponding points in the other phase point cloud is projected onto the image, and the pixel value is the elevation difference normalized to the range of 0~255. The pixel position in the image represents the planar position of the point cloud, and the pixel value in the image represents the deformation of the point cloud at the corresponding position.

9. A multi-phase point cloud registration system for ground station scenarios, characterized in that, include: The data preprocessing and preliminary registration unit is used to collect point cloud data at different times using a terrestrial 3D laser scanner, and to crop out the monitoring area, with the remainder as the non-monitoring area; to calculate the local attributes of the point cloud in the non-monitoring area and search for generalized feature points; to calculate the high-level feature descriptors of the searched generalized feature points; and to obtain the preliminary registration point cloud based on the corresponding points with the same descriptors. The key point filtering unit is used to compare the features of corresponding points in the preliminary registration point cloud, construct a normalized weight function for the overall feature descriptor, determine the point location changes based on the function values, sort the weight values, and select the points with the smallest changes between the two point clouds. Each point serves as a globally invariant point or a point with minute changes, acting as a key point for point cloud registration. The point cloud fine registration unit is used to constrain the ICP registration method based on the selected key points using the iterative reweighted least squares method with IGGIII weights. The weight function factor is determined by calculating the Euclidean distance between corresponding points in the registered point cloud. During the iteration process, points with large distance differences are eliminated, and truly invariant points are retained. Then, the spatial transformation relationship between the two point clouds is calculated to achieve fine registration. The deformation monitoring and analysis unit is used to establish a contour map of elevation changes in the deformation area based on the registered point cloud, render different distances represented by the contour map, display the deformation of the monitoring area, and analyze the deformation trend based on the density of the contour lines.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that, The computer program is executed by a processor using the multi-phase point cloud registration method for ground station scenarios as described in any one of claims 1-8.

Images