Multi-period point cloud registration method and system in ground station scene

Abstract

The invention discloses a multi-period point cloud registration method and system in a ground station scene, and the method comprises the steps: collecting point cloud data at different times through a ground type three-dimensional laser scanner, and cutting out a monitoring region and a non-monitoring region; calculating point cloud local attributes in a non-monitoring area, searching generalized feature points, calculating advanced feature descriptors of the generalized feature points, and performing preliminary registration based on homonymy points; comparing the homonymy point features of the preliminarily registered point clouds, constructing a normalized weighting function, determining point position changes according to function values, sorting weighted values, and selecting the first N points with the minimum changes in the point clouds in the two periods as global invariant or tiny change points as registration key points. And on the basis of the key points, iterative reweighted least square constraint ICP registration with the weight fixed by an IGGIII method is used, and points with large differences are eliminated to realize fine registration. And establishing an elevation change contour map based on the registered point cloud, rendering and displaying deformation, and analyzing the trend according to the density of contour lines. Registration error accumulation caused by environmental factors is avoided, and a guarantee is provided for accurate deformation evaluation.

Description

Technical Field

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

[0002] In the field of deformation monitoring technology based on three-dimensional laser scanning point cloud data, with the continuous improvement of the requirements for deformation monitoring accuracy and reliability in infrastructure construction, geological disaster prevention and control projects, traditional monitoring methods have exposed many problems that need to be solved in practical applications.

[0003] Most existing deformation monitoring technologies are built on the basis of idealized scenarios. Their core idea is usually to presuppose the existence of a reference surface that is highly consistent with the ideal model of the target. By calculating the distance between the actual point cloud and the ideal reference surface, this is used as the basis for judging the deformation of the target. However, in actual engineering practice, the spatial geometric form of the monitored object presents significant complexity and diversity. Taking monitoring targets such as large ancient buildings and complex geological slopes as an example, their structural forms are not regular geometric shapes, and it is difficult to establish accurate spatial geometric equations through simple mathematical models. This makes it impossible to accurately determine the ideal reference surface in actual scenarios. In this case, the applicability of traditional monitoring methods that rely 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] At the same time, the actual environmental conditions faced by deformation monitoring are extremely complex and challenging. In real-world monitoring scenarios, the monitoring environment is often subject to varying degrees of vibration and continuous deformation. For example, when monitoring building deformation in busy traffic areas, the vibrations from passing vehicles and disturbances from surrounding construction activities can impact the monitoring environment. When monitoring mountain deformation in areas prone to geological disasters, factors such as minute movements in the earth's crust and soil loosening caused by rainfall can also cause the monitoring environment to undergo dynamic changes. These environmental vibrations directly cause point positions to shift between different observation periods. Even with fixed facilities such as observation piers, the piers themselves are vulnerable to environmental vibrations and may even vibrate and deform with the environment. This complex environmental condition severely impacts the accuracy and consistency of 3D laser scanning point cloud data, making data registration for deformation monitoring based on multi-period point clouds a significant challenge. Due to the shifting spatial positions of point cloud data over time, traditional methods struggle to effectively extract stable deformation features, making it impossible to accurately assess key information such as the extent and trend of deformation within the deformed area.

[0005] Therefore, how to efficiently and reliably find globally invariant points or slight change points with the same name from three-dimensional laser scanning point cloud data in a complex and changeable actual environment, and use them as reference points for point cloud registration, so as to achieve accurate registration of multi-period point cloud data; at the same time, after completing the point cloud registration, how to establish a scientific and reasonable evaluation system to accurately analyze the deformation state of the deformed area, including the deformation amount, deformation direction and deformation development trend, has become the key bottleneck in the development of deformation monitoring technology based on multi-period point cloud.

[0006] Due to deficiencies in the existing technology in terms of strategies for selecting globally invariant points or slightly changing points, precise point cloud registration methods, and deformation state assessment systems, it is still unable to effectively solve the above problems. An innovative technical solution is urgently needed to break through the limitations of existing technology and achieve high-precision and high-reliability deformation monitoring based on three-dimensional laser scanning point cloud data to meet the growing monitoring needs in actual engineering applications. Summary of the Invention

[0007] This invention aims to address the current challenges of deformation monitoring technology based on 3D laser scanning point cloud data. To address the difficulties of determining an ideal reference surface in real-world scenarios, the complex monitoring environment that makes point cloud registration difficult, and the inaccurate assessment of deformation areas, a multi-period point cloud registration method and system for ground station scenarios are provided. By selecting globally unchanged or slightly changing points in non-monitored 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] In view of the above defects or improvement needs of the prior art, as a first aspect of the present invention, the present invention provides a multi-period 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 the monitored area, leaving the remaining area as the non-monitored area. Calculate local properties of the point cloud in the non-monitored area and search for generalized feature points. Calculate high-level feature descriptors for the searched generalized feature points. Preliminary registration of point clouds is performed based on the corresponding points with the same descriptors.

[0010] S2. Compare the features of the same-name points in the preliminary registration point clouds, construct a normalized weight function for the overall feature descriptor, determine the point position change based on the function value, sort the weight values, and select the top N points with the smallest change between the two point clouds as the global invariant points or the points with slight change, which are used as the 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 conditions of the monitoring area. The deformation trend is analyzed according to the density of the contour lines.

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

[0014] 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.

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

[0016] For any point P in the non-monitoring area point cloud, 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 [[ID=3s]]

[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 principal curvature K of the local point cloud 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 found 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 with 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 angle is:

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

[0046] These relationships are counted into a simplified histogram SPFH, and then the SPFH of each neighborhood point is weighted and accumulated to obtain the FPFH feature:

[0047]

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

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

[0050] The SHOT feature descriptor is a local feature descriptor for 3D point cloud data that can capture 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 the point. Its statistical principle is:

[0051] Construct a local reference frame for each key point so that the feature is invariant to rotation; determine the x, y, and z axes of the local reference frame by calculating the covariance matrix around the key point and performing eigenvalue decomposition;

[0052] Divide the neighborhood into multiple solid angle cubes, count the normal vector direction and color information in each cube to construct a descriptor histogram;

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

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

[0055] For any point (x j ,y j ,z j ), the local description features in the k-th measurement point cloud are Where t is the number of features, the following feature vector is constructed Indicates the comprehensive situation of the local description features of the point with serial number j in the kth measurement point cloud:

[0056]

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

[0058]

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

[0060]

[0061] Where, Represents the normalized weight function value of the overall feature descriptor obtained after weighting and normalizing the residual vector of the residual point; a u is the weight when weighting the u-th item in the residual vector, and t represents the number of features in the residual vector.

[0062] Furthermore, the IGGIII-weighted iterative reweighted least squares method constrained ICP registration method in S3 is specifically:

[0063] For the two-phase point cloud P s and P t There are t point cloud registration key points in P, where P t is the target point cloud, P s It is the source point cloud, and the source point cloud needs to be registered to the target point cloud; then the point cloud P t As a reference, calculate the point cloud P s With the reference point cloud P t The Euclidean distance corresponding to the registration key points is:

[0064]

[0065] Where, (X s ,Y s ,Z s ) is the source point cloud P s The key point coordinates in (X t ,Y t ,Z t ) is the target point cloud P t The corresponding key point coordinates in ;

[0066] After obtaining the distance of j key point pairs, the distance is normalized to obtain the normalized distance residual value Vj , the calculation method is as follows:

[0067]

[0068] In the formula, mind E and maxd E are the maximum and minimum values ​​of the distance between the registration key points 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 weights in the IGG III method, the weight function factor calculation formula is:

[0070]

[0071] Where w j is the weight function factor, which is used to determine the weight of the corresponding data point in the iterative reweighted least squares method; V j is the normalized distance residual value between the corresponding points of the above-mentioned registration point cloud, which reflects the degree of spatial difference of the corresponding points in the point cloud registration process; σ is the mean square error of the residual, which reflects the degree of discreteness of the residual; k0, k1 are the set threshold parameters;

[0072] In the process of finding corresponding points by the ICP algorithm, the reweighted iterative least squares constraint can gradually eliminate the registration points with large distance differences, and retain the true invariant points in the two phases of point clouds;

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

[0074]

[0075] Where ρ is the scaling parameter and R is the rotation matrix. By calculating the rotation matrix and scaling coefficient based on the corresponding relationship, 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:

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

[0078]

[0079] The parameter estimation criteria are:

[0080]

[0081] Where Y represents the observation vector; A represents the coefficient matrix, E A represents the error matrix of the coefficient matrix A; X represents the parameter vector to be estimated; vecE A Indicates that the 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 coefficient matrix error is the same as Q A are inverse matrices of each other;

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

[0083] Furthermore, the specific method of establishing the elevation change contour map of the deformation area in S4 is:

[0084] Based on the registered point cloud of the monitoring area, the elevation difference of the corresponding point of the other point cloud is projected onto the image with one phase as the benchmark, and the pixel value is normalized to the elevation difference in the range of 0 to 255; the pixel position of the image represents the plane 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-period point cloud registration system in 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 ground-based 3D laser scanner, and cut 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; and obtain a preliminary registered point cloud based on the same-name points corresponding to the same descriptors;

[0087] The key point screening unit is used to compare the features of the same-name 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 top N points with the smallest change in the two phases of point clouds as global invariant points or slight change points, 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 weighted by the IGGIII method. The weight function factor is determined by calculating the Euclidean distance between corresponding points in the registration point cloud. During the iteration process, points with large distance differences are eliminated, and the truly unchanged points are retained. Then, the spatial transformation relationship between the two phases of point clouds is calculated to achieve fine registration.

[0089] The deformation monitoring and analysis unit is used to establish an elevation change contour map of the deformed area based on the registered point cloud, render different distances represented by the contour map, display the deformation of the monitored area, and analyze the deformation trend according to 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-period point cloud registration method in the ground station scenario.

[0091] In general, the above technical solutions conceived by the present invention can achieve the following beneficial effects compared with the prior art:

[0092] 1. The present invention's multi-period point cloud registration method for ground station scenarios divides monitoring and non-monitoring areas. By calculating local attributes and searching for generalized feature points in the non-monitoring areas, the method then calculates high-level feature descriptors for preliminary registration. This effectively addresses the difficulty of directly determining point cloud correspondences in complex real-world scenarios. This preliminary registration lays the foundation for subsequent precise registration, enabling approximate alignment of point clouds from different periods. This overcomes the initial disorder of point cloud data, provides a viable data foundation for subsequent screening of stable points and precise registration, and improves the efficiency and accuracy of point cloud registration.

[0093] 2. The present invention's multi-period point cloud registration method for ground station scenarios compares the features of initially registered point clouds with similar names and constructs a normalized weight function to select globally invariant points or points with minor changes. This addresses the issue of point position shifts caused by factors such as environmental vibration. These selected key points are less susceptible to environmental interference and can serve as a reliable registration benchmark, enabling stable correspondences between multi-period point clouds even in complex environments. This significantly enhances the reliability of point cloud registration, effectively preventing the accumulation of registration errors caused by environmental factors, and providing a strong guarantee for accurate deformation assessment.

[0094] 3. The multi-period point cloud registration method of the present invention, in a ground station scenario, achieves precise registration of two point clouds by employing an iteratively reweighted least squares constrained ICP registration algorithm with IGGIII weights. It then generates deformation monitoring contour maps based on the registered point clouds. Precise registration ensures precise spatial alignment of the point clouds, while the deformation monitoring contour maps visually display deformation conditions, facilitating analysis of deformation trends. This complete process enables complete and accurate monitoring from point cloud registration to deformation analysis, meeting the practical engineering requirements for high-precision and high-reliability deformation monitoring. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0096] Figure 2 is a specific flow chart of an embodiment of the present invention;

[0097] Figure 3 2 is a diagram of system units according to an embodiment of the present invention. DETAILED DESCRIPTION

[0098] In order to make the objectives, technical solutions and advantages of the present invention more clearly understood, the present invention is further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only intended to illustrate the present invention and are not intended to limit the present invention. In addition, the technical features involved in the various embodiments of the present invention described below may 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-period 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 the monitored area, leaving the remaining area as the non-monitored area. Calculate local properties of the point cloud in the non-monitored area and search for generalized feature points. Calculate high-level feature descriptors for the searched generalized feature points. Preliminary registration of point clouds is performed based on the corresponding points with the same descriptors.

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

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

[0104] S4. Based on the registered point cloud, establish an elevation change contour map of the deformed area, render 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 explains the above steps:

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

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

[0108] Collect point cloud data at different times using a ground-based 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 ground-based 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 by a single station or the data after splicing the original point cloud data collected by multiple stations at the same time.

[0110] 1.2 Search for generalized feature points

[0111] In the non-monitoring area, calculate the local attributes of the point cloud and search for generalized feature points. 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.

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

[0113] For any point P in the non-monitoring area point cloud, 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 [[ID=~19]]‖ < r}, there are m points in total, and the calculation expression of its local point cloud density M is: [[ID=2~0]]

[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 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, and its statistical principle is:

[0134] 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;

[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 angle is:

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

[0143] These relationships are counted into a simplified histogram SPFH, and then the SPFH of each neighborhood point is weighted and accumulated to obtain the FPFH feature:

[0144]

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

[0146] The SHOT feature descriptor (Signature of Histograms of Orientations, SHOT) is a local feature descriptor for 3D point cloud data that can capture 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 the point. Its statistical principle is:

[0147] Construct a local reference frame for each key point so that the feature is invariant to rotation. Determine the x, y, and z axes of the local reference frame by calculating the covariance matrix around the key point and performing eigenvalue decomposition.

[0148] Divide the neighborhood into multiple solid angle cubes (bins), count the normal vector direction and color information in each cube to construct a descriptor histogram;

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

[0150] Due to the above-mentioned local reference frame constructed by the SHOT feature, the spherical support area is divided into 32 small areas, 2 in the radial direction, 2 in the latitude direction, and 8 in the longitude direction. A local histogram is constructed for each small area, and 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, so the SHOT descriptor is 32×11=352 dimensions.

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

[0152] (2) Screening global invariant points or minor change points

[0153] Compare the features of the points with the same name in the preliminary registration point clouds, construct the normalized weight function of the overall feature descriptor, determine the change of the point position according to the value of the normalized weight function, sort the weight values, and screen out the top N points with the smallest change in the two phases of point clouds. These points are global invariant points or points with slight change, and are used as key points for point cloud registration.

[0154] For any point (x j ,y j ,z j ), the local description features in the k-th measurement point cloud are Where t is the number of features, the following feature vector is constructed Indicates the comprehensive situation of the local description features of the point with serial number j in the kth measurement point cloud:

[0155]

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

[0157]

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

[0159]

[0160] Where, Represents the normalized weight function value of the overall feature descriptor obtained after weighting and normalizing the residual vector of the residual point; a u is the weight when weighting the u-th item in the residual vector, and t represents the number of features in the residual vector.

[0161] (3) Two-phase point cloud registration

[0162] For the two-phase point cloud Ps and P t After the operations S1-S2, the corresponding global invariant points or slight change points are obtained, and there are still slight differences between the corresponding feature points. In order to improve the robustness of the registration, the iterative reweighted least squares method based on the IGGIII weighting is used to constrain the ICP registration algorithm. The regression model of the weighted total least squares is:

[0163]

[0164] The parameter estimation criteria are:

[0165]

[0166] Where Y represents the observation vector; A represents the coefficient matrix, E A represents the error matrix of the coefficient matrix A; X represents the parameter vector to be estimated; vecE A Indicates that the 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 coefficient matrix error is the same as Q A are inverse matrices of each other;

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

[0168] For the two-phase point cloud P s and P t There are t point cloud registration key points in P, where P t is the target point cloud, P s It is the source point cloud, and the source point cloud needs to be registered to the target point cloud; then the point cloud P t As a reference, calculate the point cloud P s With the reference point cloud P t The Euclidean distance corresponding to the registration key points is:

[0169]

[0170] Where, (X s ,Y s ,Z s ) is the source point cloud P s The key point coordinates in (X t ,Y t ,Z t ) is the target point cloud P t The corresponding key point coordinates in ;

[0171] After obtaining the distance of j key point pairs, the distance is normalized to obtain the normalized distance residual value V j , the calculation method is as follows:

[0172]

[0173] Where, min d E and max d E are the maximum and minimum values ​​of the distance between the registration key points respectively; ε is a small amount that prevents the normalized distance residual from being zero;

[0174] Using the above normalized distance as the basis for determining the weight of the IGGⅢ method, the weight function factor calculation formula is:

[0175]

[0176] Where w j is the weight function factor, which is used to determine the weight of the corresponding data point in the iterative reweighted least squares method; V j is the normalized distance residual value between the corresponding points of the above-mentioned registration point cloud, which reflects the degree of spatial difference of the corresponding points in the point cloud registration process; σ is the mean square error of the residual, which reflects the degree of discreteness of the residual; k0, k1 are the set threshold parameters;

[0177] In the process of finding corresponding points by the ICP algorithm, the reweighted iterative least squares constraint can gradually eliminate the registration points with large distance differences, and retain the true invariant points in the two phases of point clouds;

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

[0179]

[0180] Where ρ is the scaling parameter and R is the rotation matrix. By calculating the rotation matrix and scaling coefficient based on the corresponding relationship, 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 the deformed area: Based on the registered monitoring area point cloud, the elevation differences of corresponding points in one phase of the point cloud are projected onto an image, with the pixel values ​​normalized to a range of 0 to 255. The pixel position in the image represents the plane position of the point cloud, and the pixel value on the image represents the deformation of the point cloud at the corresponding location.

[0183] The generated contour map can be used to determine deformation in the monitored area. Normalized pixel values ​​are converted to color gradients using a specified color scale for intuitive rendering. Areas with dense elevation changes are analyzed. These areas correspond to dense contour lines, indicating significant deformation at these locations in the two point clouds and requiring special attention.

[0184] Example 2

[0185] Please refer to Figure 3 This embodiment 2 provides a multi-period point cloud registration system in 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 ground-based 3D laser scanner, and cut 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; and obtain a preliminary registered point cloud based on the same-name points corresponding to the same descriptors;

[0187] The key point screening unit is used to compare the features of the same-name 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 top N points with the smallest change in the two phases of point clouds as global invariant points or slight change points, 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 weighted by the IGGIII method. The weight function factor is determined by calculating the Euclidean distance between corresponding points in the registration point cloud. During the iteration process, points with large distance differences are eliminated, and the truly unchanged points are retained. Then, the spatial transformation relationship between the two phases of point clouds is calculated to achieve fine registration.

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

[0190] Example 3

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

[0192] The computer-readable storage medium may include: a USB flash drive, a mobile hard disk, a read-only memory (ROM), a random access memory (RAM), a magnetic disk or an optical disk, and other media that can store program codes.

[0193] For an introduction to the computer-readable storage medium provided in this application, please refer to the above method embodiment, and this application will not go into details here.

[0194] It will be easily understood by those skilled in the art that the above description is only 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 in the scope of protection of the present invention.

Claims

1. A multi-period 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 the monitored area, leaving the remaining area as the non-monitored area. Calculate local properties of the point cloud in the non-monitored area and search for generalized feature points. Calculate high-level feature descriptors for the searched generalized feature points. Preliminary registration of point clouds is performed based on the corresponding points with the same descriptors. S2. Compare the features of the same-name points in the preliminary registration point clouds, construct a normalized weight function for the overall feature descriptor, determine the point position change based on the function value, sort the weight values, and select the top N points with the smallest change between the two point clouds as the global invariant points or the points with slight change, which are used as the key points for point cloud registration; S3. Based on the selected key points, the ICP registration method is constrained using the iterative reweighted least squares method weighted by the IGGIII method. The weighting function is determined by calculating the Euclidean distance between corresponding points in the registered point clouds. During the iteration process, points with large distance differences are eliminated, and the truly invariant points are retained. The spatial transformation relationship between the two point clouds is then calculated to achieve precise registration. S4. Based on the registered point cloud, establish an elevation change contour map of the deformed area, render 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-period point cloud registration method in a ground station scenario according to claim 1 is characterized in that: In S1, the local attributes of the point cloud are calculated in the non-monitoring area, and the specific method for searching for generalized feature points is: The local attributes of the point cloud include information including 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; The local normal vector V is calculated as follows: 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: Translate all points to the centroid and obtain the translated point set N ′ (P) = {P i ′ |P i ∈N(P)}, where The covariance matrix of the point set after translation is: Perform eigenvalue decomposition on the covariance matrix: Cv i =λ i v i Wherein, the dimension of the covariance matrix C is 3×3, λ1, λ2, λ3 are the eigenvalues ​​of the covariance matrix, and λ1≤λ2≤λ3; v1, v2, v3 are the corresponding eigenvectors; Then the local normal vector V = v1; The calculation method of the principal curvature K of the local point cloud is: The calculation method of the local point cloud density M and the local point cloud intensity I is: For any point P in the point cloud of the non-monitoring area and its neighborhood point set N(P) = {P i |‖P - P i ‖ < r}, there are m points in total, and the calculation expression of its local point cloud density M is: Where π is the ratio of circumference to diameter; The calculation expression of its local point cloud intensity I is: According to the above calculation results, corner points, inflection points, wall points, high reflection intensity points, and high density points in the non-monitoring area are searched as the generalized feature points.

3. The multi-period point cloud registration method in a ground station scenario according to claim 1 is characterized in that: The specific method for calculating the high-level feature descriptor of the searched generalized feature points in S1 is as follows: Calculating the FPFH features and SHOT feature descriptors of 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 neighboring points. Its statistical principle is: For any point P and its neighborhood point Q in the generalized feature region point cloud, their respective normal vectors are V P and V Q , then there are three angles and one distance relationship between the two points; The vector between two points is Normalize d: Normal vector V P The angle between d and d is: In V Q The projection angle on is: V P and V Q The angle is: θ=arccos(V P ·V Q , These relationships are counted into a simplified histogram SPFH, and then the SPFH of each neighborhood point is weighted and accumulated to obtain the FPFH feature: Where N k is the number of SPFHs in the neighborhood; W PQ It is the weight between point P and the neighboring point Q, reflecting the contribution of the SPFH feature of the neighboring point Q in calculating the FPFH feature of point P.

4. The multi-period point cloud registration method in a ground station scenario according to claim 3 is characterized in that: The SHOT feature descriptor is specifically: The SHOT feature descriptor is a local feature descriptor for 3D point cloud data that can capture 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 the point. Its statistical principle is: Construct a local reference frame for each key point so that the feature is invariant to rotation; determine the x, y, and z axes of the local reference frame by calculating the covariance matrix around the key point and performing eigenvalue decomposition; Divide the neighborhood into multiple solid angle cubes, count the normal vector direction and color information in each cube to construct a descriptor histogram; In each cube, the normal vector direction and color information are statistically analyzed to construct a descriptor histogram; by normalizing the histogram, the robustness to point cloud density and noise is enhanced.

5. The multi-period point cloud registration method in a ground station scenario according to claim 1 is characterized in that: The specific method for constructing the normalized weight function of the overall feature descriptor in S2 is: For any point (x j ,y j ,z j ), the local description features in the k-th measurement point cloud are Where t is the number of features, then construct the following feature vector V j k , which represents the comprehensive situation of the local description features of the point with serial number j in the kth measurement point cloud: Then, in the two measurements, the residual vectors of period k and period k+1 are expressed as: △V j =V j k -V j k+1 Each item in the current point residual vector is weighted and normalized to obtain the normalized weight function of the overall feature descriptor as follows: Where, Represents the normalized weight function value of the overall feature descriptor obtained after weighting and normalizing the residual vector of the residual point; a u is the weight when weighting the u-th item in the residual vector, and t represents the number of features in the residual vector.

6. The multi-period point cloud registration method in a ground station scenario according to claim 1 is characterized in that: The iterative reweighted least squares method constrained ICP registration method of the IGGIII legal weights in S3 is specifically: For the two-phase point cloud P s and P t There are t point cloud registration key points in P, where P t is the target point cloud, P s It is the source point cloud, and the source point cloud needs to be registered to the target point cloud; then the point cloud P t As a reference, calculate the point cloud P s With the reference point cloud P t The Euclidean distance corresponding to the registration key points is: Where, (X s ,Y s ,Z s ) is the source point cloud P s The key point coordinates in (X t ,Y t ,Z t ) is the target point cloud P t The corresponding key point coordinates in ; After obtaining the distance of j key point pairs, the distance is normalized to obtain the normalized distance residual value V j , the calculation method is as follows: In the formula, mind E and maxd E are the maximum and minimum values ​​of the distance between the registration key points respectively; ε is a small amount that prevents the normalized distance residual from being zero; Use the above normalized distance residual value V j As the basis for determining weights in the IGG III method, the weight function factor calculation formula is: Where w j is the weight function factor, which is used to determine the weight of the corresponding data point in the iterative reweighted least squares method; V j is the normalized distance residual value between the corresponding points of the above-mentioned registration point cloud, which reflects the degree of spatial difference of the corresponding points in the point cloud registration process; σ is the mean square error of the residual, which reflects the degree of discreteness of the residual; k0, k1 are the set threshold parameters; In the process of finding corresponding points by the ICP algorithm, the reweighted iterative least squares constraint can gradually eliminate the registration points with large distance differences, and retain the true invariant points in the two phases of point clouds; Finally, calculate the spatial transformation relationship between the two point clouds: Where ρ is the scaling parameter and R is the rotation matrix. By calculating the rotation matrix and scaling coefficient based on the corresponding relationship, the precise registration of the two phases of point cloud data can be achieved.

7. The multi-period 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: The weighted total least squares regression model is: The parameter estimation criteria are: Where Y represents the observation vector; A represents the coefficient matrix, E A represents the error matrix of the coefficient matrix A; X represents the parameter vector to be estimated; vecE A Indicates that the 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 coefficient matrix error is the same as Q A are inverse matrices of each other; By minimizing To estimate the parameter X, the observation error and the coefficient matrix error are minimized in a weighted sense.

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

9. A multi-period point cloud registration system in a ground station scenario, characterized in that: include: The data preprocessing and preliminary registration unit is used to collect point cloud data at different times using a ground-based 3D laser scanner, and cut 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; and obtain a preliminary registered point cloud based on the same-name points corresponding to the same descriptors; The key point screening unit is used to compare the features of the same-name 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 top N points with the smallest change in the two phases of point clouds as global invariant points or slight change points, which serve as key points 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 weighted by the IGGIII method. The weight function factor is determined by calculating the Euclidean distance between corresponding points in the registration point cloud. During the iteration process, points with large distance differences are eliminated, and the truly unchanged points are retained. Then, the spatial transformation relationship between the two phases of point clouds is calculated to achieve fine registration. The deformation monitoring and analysis unit is used to establish an elevation change contour map of the deformed area based on the registered point cloud, render different distances represented by the contour map, display the deformation of the monitored area, and analyze the deformation trend according to 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 to implement the multi-period point cloud registration method in a ground station scenario as described in any one of claims 1-8.

Images