A mine surface three-dimensional movement extraction method and system based on improved CPD non-rigid point cloud registration

Abstract

The present application relates to a kind of based on improved CPD nonrigid point cloud registration mining area surface three-dimensional movement extraction method and system, belong to surface detection technical field, solve the problem of poor adaptability of nonrigid registration in complex sparse area, low computing efficiency in prior art.The specific steps include: obtaining the nonrigid point cloud in the two-period unmanned aerial vehicle laser point cloud data of the mining area to be measured, respectively on the nonrigid point cloud of two periods multi-layer feature point extraction, obtain the first candidate point cloud set and the second candidate point cloud set;According to the set registration window and the moving step on the first candidate point cloud set and the second candidate point cloud set, obtain multiple local unit point cloud pairs;Using direction constrained CPD algorithm, sequentially nonrigid registration is carried out to each local unit point cloud pair, obtain each local unit deformation field;Superimpose each local unit deformation field, obtain global deformation field, significantly improve the efficiency and precision of mining area surface three-dimensional movement extraction.

Description

Technical Field

[0001] This invention relates to the field of surface exploration technology, and in particular to a method and system for extracting three-dimensional movement of the surface in mining areas based on improved CPD non-rigid point cloud registration. Background Technology

[0002] Surface deformation in mining areas is an unavoidable geological phenomenon during mining operations. It not only affects the stability of the geological environment of the mining area but may also cause a series of ecological and safety problems such as damage to surface buildings and loss of surface water resources. With the continuous expansion of coal mining scale in my country, monitoring surface subsidence and three-dimensional movement in mining areas has become an important issue for mine safety and ecological protection.

[0003] With the rapid development of computer technology, drones are being applied in more scenarios. UAV-LiDAR technology, with its high precision, high resolution, and ability to quickly acquire three-dimensional information about the Earth's surface, provides massive amounts of point cloud data for monitoring surface deformation in mining areas. However, how to efficiently and accurately extract the three-dimensional deformation field from these massive point clouds, especially for the identification and quantification of non-rigid deformations such as local deformations like surface tension and compression, remains a challenge in current research. In non-rigid point cloud registration, early studies were mostly based on feature matching and physical models, such as the energy function method based on thin plate splines (TPS) and the improved Iterative Closest Point (ICP) algorithm based on local affine transformations. These methods have improved the adaptability of non-rigid registration to some extent, but still suffer from problems such as sensitivity to initial position and low computational efficiency.

[0004] With the development of deep learning, unsupervised registration methods based on neural networks have gradually become a research hotspot. The Consistent Point Shift (CPD) algorithm, a non-rigid registration method based on Gaussian mixture models, has been widely used for matching medical images with 3D models. However, the CPD algorithm has low computational efficiency when processing large-scale point clouds, limiting its application in practical engineering. It is evident that existing non-rigid registration methods still have shortcomings in monitoring surface deformation in mining areas, such as strong dependence on point cloud density and features, poor adaptability in complex and sparse regions, and low computational efficiency. Summary of the Invention

[0005] Based on the above analysis, the embodiments of the present invention aim to provide a method and system for three-dimensional motion extraction of the surface of mining areas based on improved CPD non-rigid point cloud registration, so as to solve the problems of poor adaptability and low computational efficiency of non-rigid registration in complex and sparse regions in the prior art.

[0006] The objective of this invention is mainly achieved through the following technical solutions:

[0007] On one hand, embodiments of the present invention provide a method for three-dimensional movement extraction of the surface of a mining area based on improved CPD non-rigid point cloud registration, comprising the following steps: Non-rigid point clouds were obtained from the UAV laser point cloud data of the mining area before and after two periods. Multi-layer feature point extraction was performed on the non-rigid point clouds of the two periods to obtain the first candidate point cloud set and the second candidate point cloud set. According to the set registration window and moving step size, the system moves on the first phase candidate point cloud set and the second phase candidate point cloud set to obtain multiple local unit point cloud pairs; using the orientation-constrained CPD algorithm, non-rigid registration is performed on each local unit point cloud pair in sequence to obtain the deformation field of each local unit. The global deformation field is obtained by superimposing the deformation fields of each local unit.

[0008] Furthermore, the deformation field of any local element is obtained, including: The local unit point cloud pairs are initially registered to obtain the first translation vector corresponding to the original coordinate system of the point cloud pairs. Using the vertical direction of the original coordinate system as the main deformation direction, and using the rotation matrix constructed by the first translation vector, the two phase point clouds within the local unit are projected onto the vertical direction respectively, to obtain the first phase candidate point cloud and the second phase candidate point cloud after rotation. Using the first phase candidate point cloud after rotation as the source point set and the second phase candidate point cloud after rotation as the target point set, CPD non-rigid point cloud registration is performed in the vertical direction to obtain the second translation vector. Based on the first translation vector, the second translation vector, and the rotation matrix, the deformation variables are inverted to obtain the local element deformation field.

[0009] Furthermore, the first phase candidate point cloud or the second phase candidate point cloud is obtained, including: Based on elevation and reflection intensity fluctuations in non-rigid point clouds in the early or later stages, ground point cloud sets are extracted in combination with preset thresholds, and the remaining point cloud after extraction is used as the initial screening point cloud set. Key points and radii are determined on the initial screening point cloud to obtain several regions of interest point clouds. Calculate the curvature ratio feature and normal vector ratio feature of each point in the point cloud of each region of interest, extract points whose curvature ratio feature and normal vector ratio feature are greater than the preset curvature threshold and normal vector threshold respectively, and combine them with the ground point cloud to obtain the first phase candidate point cloud or the second phase candidate point cloud.

[0010] Furthermore, the normal vector ratio feature of the point cloud of any region of interest is obtained based on the following steps: Calculate the three-dimensional geometric centroid of the point cloud of the region of interest; Based on the three-dimensional geometric centroid, construct the covariance matrix of each point data in the point cloud of the region of interest, and extract the minimum value in the covariance matrix as the first feature value of the corresponding point. The average of the sum of the angles between the first eigenvalue of each point and the first eigenvalues ​​of other points is obtained to get the rate of change of the local normal vector of the corresponding point; The rate of change of the normal vector within the region of interest is obtained by averaging the rate of change of all local normal vectors in the point cloud of the region of interest. The normal vector ratio feature of the corresponding point is obtained by comparing the ratio of the local normal vector change rate at each point to the normal vector change rate within the region of interest.

[0011] Furthermore, the curvature ratio features of each point in the point cloud of any region of interest are obtained based on the following steps: Based on each point in the point cloud of the region of interest, the surface is fitted using the least squares method to obtain the local curvature change rate corresponding to each point; The rate of curvature change within the region of interest is obtained by averaging the rates of curvature change of all local points in the point cloud of the region of interest. The curvature ratio feature of the corresponding point is obtained based on the ratio of the local curvature change rate at each point to the curvature change rate within the region of interest.

[0012] Furthermore, a Gaussian difference pyramid is constructed using scale-invariant features to detect scale-space extrema and determine the key points.

[0013] Furthermore, using the vertical direction as the rotation axis, and based on the angle between the first translation vector and the vertical direction, the Rodrigues rotation formula is constructed to obtain the rotation matrix.

[0014] Furthermore, singular value decomposition is used to perform preliminary registration on the local unit point cloud pairs.

[0015] Furthermore, the registration window includes a moving window and a buffer window, and the size of the buffer window is the sum of the moving window and the estimated maximum horizontal displacement; the moving window is advanced on the first phase candidate point cloud according to the moving step size, while the buffer window is advanced synchronously on the second phase candidate point cloud, to obtain multiple local unit point cloud pairs.

[0016] On the other hand, embodiments of the present invention provide a three-dimensional moving extraction system for mining area surface based on improved CPD non-rigid point cloud registration, comprising: The data acquisition module is used to acquire non-rigid point clouds from the UAV laser point cloud data of the mining area before and after two periods. The feature extraction module is used to extract multi-layer feature points from the non-rigid point clouds of the two phases respectively, to obtain the first phase candidate point cloud set and the second phase candidate point cloud set; The registration module is used to move on two candidate point cloud sets according to the set registration window and movement step size to obtain multiple local unit point cloud pairs. Using the orientation-constrained CPD algorithm, non-rigid registration is performed on each local unit point cloud pair in sequence to obtain the deformation field of each local unit. The output module is used to superimpose the deformation fields of each local unit to obtain the global three-dimensional deformation field.

[0017] Compared with the prior art, the present invention can achieve at least one of the following beneficial effects: 1. This invention addresses the problem of scarce surface feature points in mining areas and the presence of non-rigid deformations such as stretching and compression. It proposes to improve the accuracy and discriminative power of point cloud feature point extraction by preprocessing two phases of point cloud data through multi-level semantic feature extraction. At the same time, it combines a registration window strategy and directional constraints for efficient CPD registration of non-rigid point clouds, thereby improving the efficiency and accuracy of surface movement extraction in mining areas.

[0018] 2. The orientation-constrained CPD algorithm transforms the point cloud into the vertical direction for registration through rotation optimization, avoiding the complex displacement field regularization calculation in traditional CPD and significantly reducing computational overhead.

[0019] 3. By using a multi-level semantic feature extraction framework, Euclidean distance, key points, and neighborhood semantics are integrated to improve the stability and robustness of feature points, reduce dependence on fixed thresholds, and enhance adaptability in complex terrains.

[0020] In this invention, the above-described technical solutions can be combined with each other to achieve more preferred combinations. Other features and advantages of this invention will be set forth in the following description, and some advantages may become apparent from the description or be learned by practicing the invention. The objects and other advantages of this invention can be realized and obtained from what is particularly pointed out in the description and drawings. Attached Figure Description

[0021] The accompanying drawings are for illustrative purposes only and are not intended to limit the invention. Throughout the drawings, the same reference numerals denote the same parts.

[0022] Figure 1 This is a flowchart of a method for extracting three-dimensional movement of the surface in a mining area based on improved CPD non-rigid point cloud registration, according to an embodiment of the present invention. Figure 2 This is a logic block diagram of the registration method according to an embodiment of the present invention; Figure 3 This is a schematic diagram of the local registration process in an embodiment of the present invention; Figure 4 This is a schematic diagram of movement after directional constraint calibration according to an embodiment of the present invention. Detailed Implementation

[0023] Preferred embodiments of the present invention will now be described in detail with reference to the accompanying drawings, which form part of this application and are used together with the embodiments of the present invention to illustrate the principles of the present invention, but are not intended to limit the scope of the present invention.

[0024] Example 1 A specific embodiment of the present invention discloses a method for extracting three-dimensional movement of the surface in mining areas based on improved CPD non-rigid point cloud registration, such as... Figure 1 As shown, it includes the following steps: Step S1: Obtain non-rigid point clouds from the UAV laser point cloud data of the mining area before and after the test, and extract multi-layer feature points from the non-rigid point clouds of the two periods respectively to obtain the first candidate point cloud set and the second candidate point cloud set. Step S2: Move on the first and second candidate point cloud sets according to the set registration window and movement step size to obtain multiple local unit point cloud pairs; use the orientation-constrained CPD algorithm to perform non-rigid registration on each local unit point cloud pair in sequence to obtain the deformation field of each local unit. Step S3: Superimpose the deformation fields of each local unit to obtain the global deformation field.

[0025] Using the above method, feature extraction is performed on two phases of non-rigid point clouds through multi-level semantics. Based on the two phases of candidate point clouds after feature extraction, CPD non-rigid point cloud registration is performed using the registration window strategy and directional constraints to obtain the local deformation field of the mining area to be tested. After integration, the global deformation field is obtained, which improves the adaptability of non-rigid registration in complex and sparse regions and improves the registration calculation efficiency.

[0026] It's important to note that point cloud data registration can be categorized into rigid registration and non-rigid registration, depending on whether the object itself undergoes deformation. The principle of rigid registration is that two 3D point clouds are represented by the model point cloud X and the target point cloud Y, respectively. Assuming a relationship T exists between the two point clouds such that T: X->Y, rigid registration aims to find the optimal rigid transformation T, making T(X) and Y coincide as much as possible. Due to the deformation characteristics of the object itself, all vertices in rigid registration share a single transformation matrix. In contrast to rigid registration, point clouds or meshes (i.e., point clouds with connections) in non-rigid registration can undergo deformation. Assuming an object consists of n 3D points, non-rigid registration is considered to consist of a weighted average of several rigid baselines. The mathematical principle of non-rigid registration is as follows: the deformation template (source) is represented by X, the registration target is represented by Y, and the deformed surface is represented by X^. If we use X^=Minimize E(fΘ( X ),Y) to represent the deformation of the template, then Θ is the parameter to be optimized. Generally speaking, the larger the number of parameters, the greater the degree of freedom of template deformation (and generally the more accurate the registration), but at the same time, the optimization difficulty and cost are also higher.

[0027] Based on this, and considering the typical terrain characteristics of mining areas—namely, the lack of rigid features such as buildings and roads, and the widespread distribution of non-rigid features such as bare soil, vegetation, and loose deposits, resulting in sparse natural surface features and susceptibility to seasonal changes—this invention proposes a method for synchronously extracting three-dimensional surface movement based on UAV lidar point cloud and non-rigid body registration technology. Figure 2 As shown.

[0028] Specifically, in step S1, the surface of the mining area is mainly a continuous and flexible natural surface, with sparse features and susceptible to seasonal changes. Therefore, to address the difficulties in feature point extraction caused by complex terrain, uneven point cloud density, and noise interference in the mining area, this method constructs a multi-level semantic feature framework, integrating Euclidean distance semantics, keypoint semantics, and neighborhood semantics. This framework extracts features from two phases of non-rigid point clouds, obtaining corresponding candidate point clouds. This effectively captures and enhances weak features such as slight terrain undulations and vegetation edges, improving the feature representativeness in non-rigid terrain environments. The specific process for obtaining the first and second phase candidate point clouds is as follows: S11. Obtain non-rigid point clouds from the UAV laser point cloud data of the mining area before and after the test; At different time points, unmanned aerial vehicle (UAV) lidar was used to scan the same mining area to obtain lidar point cloud data for the corresponding period. For example, the data scanned at the beginning of the monitoring period was used as the first phase of unmanned lidar point cloud data, and the comparative data collected later was used as the second phase of lidar point cloud data. Each phase of point cloud data contains rigid components in stable areas and non-rigid components in deformable areas. The non-rigid point cloud data was extracted for surface movement analysis of the mining area.

[0029] For example, the two original point cloud data are denoised and thinned, and the coordinate system and elevation benchmark are unified to obtain the first and second preprocessed point clouds. The two preprocessed point clouds are rigidly registered. For example, the key point detection (ISS) is used to lock the key point positions of the two preprocessed point clouds. Then, the descriptor (3DSC) is used to generate high-discrimination vectors for these key point positions (i.e., the surrounding space is divided into grids with each ISS point as the center, and the number of points falling into the grid is counted to form a high-dimensional orthogonal vector) to form the nearest neighbor pairs of the two point clouds. Finally, the geometric verification (RANSAC) is used to match all nearest neighbor pairs and remove erroneous vector pairs, so as to quickly and robustly put the two UAV laser point clouds together roughly to obtain two non-rigid point clouds.

[0030] S12. Based on the elevation fluctuation and reflection intensity fluctuation in the non-rigid point cloud in the early or later stage, extract the ground point cloud set in combination with the preset threshold, and use the remaining point cloud after extraction as the initial screening point cloud set. Specifically, this layer calculates the mean elevation, variance elevation, mean intensity, variance intensity, and relative elevation within the neighborhood of the two non-rigid point clouds. The variance elevation reflects elevation fluctuations, and the variance intensity reflects reflection intensity fluctuations. Large elevation fluctuations indicate significant topographic changes or the presence of ground features; large reflection intensity fluctuations suggest the boundary between different ground features. These fluctuations are used to identify abrupt topographic changes and boundaries between ground features and the ground, thus separating non-near-ground points from the non-rigid point cloud.

[0031] For example, for any period of non-rigid point cloud data, based on the Euclidean space attributes of the point cloud, the following Euclidean semantic feature set is constructed: ;in, For point elevation, The mean elevation of the point cloud (i.e., the z-coordinate) For the elevation variance of point clouds, The relative elevation of the point cloud (the height of each laser point relative to the local terrain surface at its location). Let V be the intensity variance of the point cloud. By calculating the corresponding feature values ​​of each point cloud data point, a Euclidean semantic feature set is obtained, and the neighborhood range is defined based on Euclidean distance, such as the point elevation being located within... Points within the specified range are sorted in descending order of their elevation variance and intensity variance. A preset threshold is used to divide the sorted data into two parts. The points in the first half of the sorted data constitute the initial ground point set, while the points in the second half are considered the candidate ground point set (i.e., the initial screening point set), requiring further refinement to obtain the accurate ground point set. For example, if the preset threshold is set to 50%, the first 50% are considered the initial ground points, and the second 50% are considered the initial screening point set.

[0032] S13. Key points and radii are determined on the initial point cloud set to obtain several regions of interest point clouds, thus defining the neighborhood semantic space. Key points are sparse point sets with stability and discriminative power in the point cloud data. Combining them with local feature descriptors allows for the construction of a compact point cloud structure representation, which helps reduce data volume, improve computational efficiency, and suppress noise interference. Setting the radius too small makes it susceptible to noise interference, leading to estimation bias; while a large radius can resist noise, it increases computational load and excessively smooths local features, reducing discriminative power. The radius is determined by balancing feature resolution and computational efficiency.

[0033] For example, a Gaussian difference pyramid can be constructed using the Scale-Invariant Feature Transform (SIFT) operator to detect scale-space extrema and extract stable keypoints. This method exhibits scale and rotation invariance and good robustness to changes in illumination, noise, and viewpoint. For instance, let... For a given key point in the initial screening cluster, calculate the distance between it and its nearest neighbor key points. Obtain the radius R. Draw a spherical neighborhood in space with R as the range to find more feature points with significant geometric correlations.

[0034] S14. Calculate the curvature ratio feature and normal vector ratio feature of each point in the point cloud of each region of interest, extract the points whose curvature ratio feature and normal vector ratio feature are greater than the preset curvature threshold and normal vector threshold respectively, and combine them with the ground point cloud to obtain the first phase candidate point cloud or the second phase candidate point cloud.

[0035] Specifically, features are extracted based on normal vectors and curvature within the region of interest to more effectively identify significant changes in surface structure. On one hand, the normal vector of a point cloud, as an important geometric attribute, can characterize the undulations of the local surface. In flat areas, the direction of the normal vector changes gradually; while in areas with significant surface undulations, the direction of the normal vector changes drastically. Therefore, by quantifying and comparing the degree of change in the normal vector, points with significant geometric features can be identified, such as terrain boundaries, steep slopes, and other abrupt topographic changes. The process of calculating the normal vector feature points for each region of interest is as follows: (a) Calculate the three-dimensional geometric centroid of the point cloud within the region of interest, expressed as: , in, It is the three-dimensional geometric centroid; The spatial coordinates of point i are represented by k; k is the number of point clouds in the region of interest.

[0036] (b) Construct the covariance matrix of each point data in the point cloud of the region of interest based on the three-dimensional geometric centroid, and extract the minimum value in the covariance matrix as the first feature value of the corresponding point; For example, build Find the covariance matrix of the nearest neighbor and solve for the smallest eigenvalue of the matrix. The first eigenvalue of a point is the minimum value in the covariance matrix. The covariance matrix C is expressed as: , Where, x i y i and z i Let P be the three-dimensional spatial coordinates of the i-th point within the region of interest. ox P oy P oz Then the corresponding three-dimensional geometric centroid P o Spatial location coordinates.

[0037] (c) Obtain the average of the sum of the angles between the first eigenvalue of each point within the region of interest and the first eigenvalues ​​of other points, and obtain the rate of change of the local normal vector of the corresponding point, expressed as: , in, Let be the rate of change of the local normal vector at point i. Let be the angle between the first eigenvalues ​​of point i and point j, and K be the number of point clouds within the region of interest.

[0038] (d) Based on the average of the local normal vector change rates of all point clouds in the region of interest, the normal vector change rate within the region of interest is obtained; based on the ratio of the local normal vector change rate at each point to the normal vector change rate within the region of interest, the normal vector ratio feature of the corresponding point is obtained, expressed as: , , in, The ratio characteristic of the normal vector at point i. Let be the rate of change of the normal vector within the region of interest containing the i-th point.

[0039] On the other hand, curvature is a geometric index that quantifies the degree of local bending of a surface. A larger value indicates more significant concavity and convexity changes, commonly found in characteristic areas such as edges and corners. This invention uses a local quadratic surface fitting method to reconstruct the neighborhood geometry, and then estimates the average curvature of each point to obtain the curvature ratio characteristics of each point. The specific implementation process is as follows: (a) Based on each point in the point cloud of the region of interest, fit the surface using the least squares method to obtain the local curvature change rate corresponding to each point; Specifically, by fitting the surface containing each point using the least squares method, a local quadratic surface z is constructed. The curvature of any point at coordinates (x, y) is then the curvature of the local quadratic surface z. The quadratic surface fitting formula is as follows: , We perform a least-squares fit on the above equation. Since the coefficients of z are unknowns, we use mathematical methods, specifically the Lagrange multiplier method, to find the extreme values ​​when the first derivative is zero. The coefficients of the quadratic surface are obtained by solving for the minimum value. Based on the partial derivatives of the surface parametric equations , , , , Obtain the average curvature at this point The calculation formula is as follows: .

[0040] (b) Based on the average of all local curvature change rates of the point cloud of the region of interest, the curvature change rate within the region of interest is obtained; based on the ratio of the local curvature change rate at each point to the curvature change rate within the region of interest, the curvature ratio feature of the corresponding point is obtained, expressed as: , , in, The curvature ratio feature at point i is... Let be the rate of curvature change within the region of interest containing point i.

[0041] Since feature changes within a region of interest are typically gradual, their feature values ​​are often lower than those in local regions. Therefore, each point can be... , The effective feature point cloud is identified by comparing it with the corresponding preset curvature threshold and preset normal vector threshold, such as when , Points with values ​​greater than 1 can be directly considered as valid feature points. All valid feature point sets are merged with the ground point set from step S12 to obtain either the first-phase candidate point set or the second-phase candidate point set.

[0042] By integrating the three semantic features mentioned above and introducing a relative ratio mechanism, efficient and objective feature point extraction is achieved. Compared with traditional methods that rely on fixed thresholds, this strategy not only simplifies the process and reduces dependence on empirical parameters, but also fundamentally improves the automation level and objectivity of feature extraction. It effectively overcomes the limitations of using a single semantic method to miss or misextract features in low curvature regions, significantly improving the accuracy, discriminative power, and robustness of feature point extraction, and providing a stable and reliable feature foundation for subsequent point cloud registration.

[0043] Specifically, in step S2, addressing the problems of computational complexity and limited ability to characterize local deformations in large-scale scenes by global non-rigid registration, this invention proposes a moving registration unit strategy. This strategy divides the overall point cloud into local units for independent registration, thereby transforming complex global deformations into a superposition of a series of local deformation fields. The steps for obtaining the deformation fields of each pair of local units are as follows: S21. Preset the parameters of the moving registration unit, and control the registration process through three key parameters, such as... Figure 3 As shown, the three key parameters are the moving window, the buffer window, and the moving step size. The registration window is composed of the moving window and the buffer window. By moving on the candidate point set in the two phases through the registration window and the moving step size, multiple local unit pairs are obtained.

[0044] In this method, the size of the moving window determines the scale of the registration unit, and the size of the buffer window is set as the sum of the moving window and the estimated maximum horizontal displacement to ensure that the source unit obtains sufficient corresponding points in the target point cloud. The moving step size is usually smaller than the moving window size to ensure that the local registration units after the movement have overlapping areas, thereby avoiding edge errors and verifying the consistency of the results through the overlapping areas. The moving window is advanced on the first candidate point cloud according to the moving step size, while the buffer window is advanced synchronously on the second candidate point cloud, resulting in multiple pairs of local unit point clouds. This partitioning method lays the foundation for subsequent high-precision and high-efficiency non-rigid registration.

[0045] By employing a moving registration unit strategy, global non-rigid deformation is decomposed into superimposed local deformation fields, enabling precise modeling of continuous and non-uniform surface deformation, which is particularly suitable for deformation recovery of flexible media such as land and grassland.

[0046] S22. Using the Coherence Point Drift (CPD) algorithm with directional constraints, non-rigid registration is performed on each local element pair in sequence to obtain the deformation field of each local element pair. It should be noted that the CPD algorithm is a probabilistic registration method based on a Gaussian Mixture Model (GMM). It transforms the registration problem into a probability density estimation problem, where the source point set is used as the movable centroid of the GMM, and the target point set is considered as the observations sampled by the GMM. A non-rigid transformation moves these centroids, ensuring that the GMM defined by the source point set generates the target point set to the greatest extent possible. Simultaneously, a regularization term ensures the smoothness and coherence of the transformation, meaning that the motions of adjacent points should be similar. The specific steps are as follows: First, obtain the input target point set. Source point set N and M represent the number of data points in the corresponding point set, respectively, and D represents the data dimension. A Gaussian mixture model is constructed, assuming that each point in the target point set is generated by a Gaussian distribution centered at the source point, plus a uniform distribution (to handle noise, redundant points, and outliers). A displacement field model is used to represent the non-rigid transformation, where the displacement of each source point is obtained through a linear combination of a set of Gaussian radial basis functions, expressed as: T(Y, W) = Y + G W; where G is a Gaussian kernel matrix (measuring the similarity between points, G_ij = exp(-||y_i - y_j||² / 2β²)); W is the displacement coefficient matrix to be determined; β is the Gaussian kernel width, used to control the local smoothness of the transformation. Then, the likelihood function is maximized using the Expectation-Maximization (EM) algorithm to obtain the optimal model parameters, that is, given the initial model parameters (displacement coefficient matrix W and the variance of the Gaussian distribution components), the probability of observing the target point set is maximized by adjusting the model parameters. Finally, the deformation fields of the two types of data are calculated using the solved optimal model parameters.

[0047] However, the non-rigid registration model of CPD introduces a displacement field regularization term to smooth the displacement field and indirectly constrain the direction of the displacement vector (since adjacent point displacement vectors are similar), which complicates the solution process and incurs huge computational costs. For example, processing a point cloud containing only 10,000 points takes about 3 hours, which severely limits its practical application in large-scale scenarios. Actual surface subsidence and deformation phenomena show that each point in a local spatial range undergoes displacement along a specific direction, and the resulting deformation direction vector has consistent characteristics. Based on this, this invention introduces a direction-constrained CPD algorithm for point cloud calibration. The specific process of applying direction constraints to any local unit point cloud pair using CPD is as follows: S221. Perform preliminary registration on the local unit pairs to obtain the first translation vector corresponding to the original coordinate system of the point cloud pair. Specifically, the Singular Value Decomposition (SVD) algorithm is used to perform global orientation within local cells, and the preliminary three-dimensional global translation vector t1(x) is solved. t1 ,y t1, zt1 This process obtains the first translation vector. This vector represents the dominant trend direction of point cloud movement within the local unit region, essentially reflecting the macroscopic displacement pattern of the surface under the influence of mining: "where the whole is moving." This step quickly captures the overall directional characteristics of deformation, providing a constraint basis for subsequent fine registration.

[0048] S222. Using the vertical direction of the original coordinate system as the main deformation direction, and using the rotation matrix constructed by the first translation vector, the point clouds of the local unit in the previous and next periods are projected to the vertical direction respectively, so as to obtain the first phase candidate point cloud and the second phase candidate point cloud after rotation. For example, to obtain the first translation vector t1(x) t1 ,y t1, z t1 Let t1 be the source vector, and the perpendicular direction of the original coordinate system (i.e., the z-axis) be the target direction. The target unit vector is then represented as k = (0, 0, 1). Construct a rotation matrix R1 such that R1 × t1 = k. The rotation matrix R1 is obtained using Rodriguez's rotation formula. During the solution process, an antisymmetric matrix K is generated based on the unit vector k, and the angle θ between t1 and the k direction is obtained. The rotation matrix R1 is then calculated using the following formula: .

[0049] Based on the obtained rotation matrix R1, this rotation matrix is ​​applied to both phases of the point cloud within the local element pair, causing the two phases of the point cloud to be projected onto the vertical direction, resulting in the rotated first-phase candidate point cloud and the rotated second-phase candidate point cloud. Using the first translation direction as a directional constraint, by transforming the first translation direction to the positive z-axis, during subsequent CPD non-rigid registration, such as... Figure 4 As shown, the originally complex three-dimensional spatial displacement problem is projected into a one-dimensional displacement problem along this specific vertical direction.

[0050] S223. Using the first phase candidate point cloud after rotation as the source point set and the second phase candidate point cloud after rotation as the target point set, perform CPD non-rigid point cloud registration in the vertical direction to obtain the second translation vector.

[0051] When performing non-rigid registration with directional constraints using CPD, after aligning the principal direction t1 with the vertical axis (Z-axis) using the rotation matrix R1, the source and target point clouds are transformed into a new computational coordinate system. In this coordinate system, the complex 3D non-rigid registration problem is cleverly reduced in dimensionality, and the displacement of the point cloud is separated into two relatively independent components: the Z-axis component and the XY-plane component. The Z-axis component represents the primary displacement, carrying the dominant trend of surface movement. Non-rigid registration is performed on this component, adjusting the displacement field model—that is, adjusting the displacement field W in the non-rigid transformation to one dimension—to quickly solve for the displacement z of each point along the Z-axis.t2 The XY plane components represent secondary displacements, perpendicular to the principal directions, reflecting complex deformations such as local lateral stretching, compression, or torsion. Non-rigid registration is performed in the XY plane, corresponding to adjusting the displacement field W to two dimensions, and the displacement quantities (x, y, y) in the plane are solved. t2 ,y t2 To capture details of lateral deformation, the computational complexity of the two-dimensional problem is much lower than that of the three-dimensional problem. Registration is used to obtain the second translation vector t2(x) in three dimensions. t2 ,y t2 ,z t2 ).

[0052] S223. Using the obtained parameters R1, t1, t2, the three-dimensional deformation of the point cloud in the local element in the original coordinate system can be inverted to obtain the deformation field of the local element. Represented as: , Where y is the coordinate of the source and target points, and G is the Gaussian kernel matrix in CPD registration.

[0053] Compared to traditional CPD algorithms, the orientation-constrained CPD non-rigid registration process of this invention compresses the search space, decomposing the 3D search into a "1D + 2D" search, greatly reducing the number of computational parameters and complexity. This simplifies the solution process. Furthermore, a moving registration unit strategy is introduced, dividing the global point cloud into local units and performing orientation-constrained CPD non-rigid registration unit by unit. This strategy not only further reduces the computational complexity of a single registration but also enhances the ability to capture local heterogeneous deformations, enabling the overall method to efficiently process large-scale point clouds while also accurately depicting complex surface deformation details.

[0054] Specifically, in step S3, the deformation fields of each local unit are obtained through step S2, and the deformation fields of each local unit are superimposed to obtain the global deformation field.

[0055] Compared with existing technologies, this embodiment provides a method for extracting 3D movement on the surface of mining areas based on improved CPD non-rigid point cloud registration. First, it constructs a feature point extraction framework based on multi-level semantics by integrating Euclidean distance semantics, key point semantics, and neighborhood semantics. Then, it introduces a movement registration unit strategy to decompose global deformation into local deformation fields for optimization. On this basis, it adds directional constraints to the traditional CPD algorithm to improve the registration efficiency of non-rigid point clouds, thereby achieving high efficiency and high accuracy in extracting 3D movement on the surface of mining areas. This effectively solves the problem of high-precision synchronous extraction of 3D deformation fields in mining environments lacking obvious ground features.

[0056] Example 2 Another specific embodiment of the present invention discloses a three-dimensional moving extraction system for mining area surface based on improved CPD non-rigid point cloud registration, comprising: The data acquisition module is used to acquire non-rigid point clouds from the UAV laser point cloud data of the mining area before and after two periods. The feature extraction module is used to extract multi-layer feature points from the non-rigid point clouds of the two phases respectively, to obtain the first phase candidate point cloud set and the second phase candidate point cloud set; The registration module is used to move on two candidate point cloud sets according to the set registration window and movement step size to obtain multiple local unit point cloud pairs. Using the orientation-constrained CPD algorithm, non-rigid registration is performed on each local unit point cloud pair in sequence to obtain the deformation field of each local unit. The output module is used to superimpose the deformation fields of each local unit to obtain the global three-dimensional deformation field.

[0057] The system can perform surface mobile extraction in the mining area according to any of the methods described in Embodiment 1. Related aspects can be referenced from each other, but are not repeated in this embodiment.

[0058] Compared with existing technologies, this embodiment provides a three-dimensional moving extraction system for mining area surface based on improved CPD non-rigid point cloud registration. Through the coordinated operation of various modules, it achieves high efficiency and high adaptability in the registration of non-rigid point clouds, providing a reliable technical approach for deformation monitoring and disaster early warning of mining subsidence.

[0059] Those skilled in the art will understand that all or part of the processes of the methods described in the above embodiments can be implemented by a computer program instructing related hardware, and the program can be stored in a computer-readable storage medium. The computer-readable storage medium may be a disk, optical disk, read-only memory, or random access memory, etc.

[0060] The above description is only a preferred embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for extracting three-dimensional movement of surface data in mining areas based on improved CPD non-rigid point cloud registration, characterized in that, Includes the following steps: Non-rigid point clouds were obtained from the UAV laser point cloud data of the mining area before and after two periods. Multi-layer feature point extraction was performed on the non-rigid point clouds of the two periods to obtain the first candidate point cloud set and the second candidate point cloud set. According to the set registration window and moving step size, the system moves on the first phase candidate point cloud set and the second phase candidate point cloud set to obtain multiple local unit point cloud pairs; using the orientation-constrained CPD algorithm, non-rigid registration is performed on each local unit point cloud pair in sequence to obtain the deformation field of each local unit. The global deformation field is obtained by superimposing the deformation fields of each local unit.

2. The method according to claim 1, characterized in that, The deformation field of any local element is obtained, including: The local unit point cloud pairs are initially registered to obtain the first translation vector corresponding to the original coordinate system of the point cloud pairs. Using the vertical direction of the original coordinate system as the main deformation direction, and using the rotation matrix constructed by the first translation vector, the two phase point clouds within the local unit are projected onto the vertical direction respectively, to obtain the first phase candidate point cloud and the second phase candidate point cloud after rotation. Using the first phase candidate point cloud after rotation as the source point set and the second phase candidate point cloud after rotation as the target point set, CPD non-rigid point cloud registration is performed in the vertical direction to obtain the second translation vector. Based on the first translation vector, the second translation vector, and the rotation matrix, the deformation variables are inverted to obtain the local element deformation field.

3. The method according to claim 2, characterized in that, Obtain either the first phase candidate point cloud or the second phase candidate point cloud, including: Based on elevation and reflection intensity fluctuations in non-rigid point clouds in the early or later stages, ground point cloud sets are extracted in combination with preset thresholds, and the remaining point cloud after extraction is used as the initial screening point cloud set. Key points and radii are determined on the initial screening point cloud to obtain several regions of interest point clouds. Calculate the curvature ratio feature and normal vector ratio feature of each point in the point cloud of each region of interest, extract points whose curvature ratio feature and normal vector ratio feature are greater than the preset curvature threshold and normal vector threshold respectively, and combine them with the ground point cloud to obtain the first phase candidate point cloud or the second phase candidate point cloud.

4. The method according to claim 3, characterized in that, The normal vector ratio feature of the point cloud of any region of interest is obtained based on the following steps: Calculate the three-dimensional geometric centroid of the point cloud of the region of interest; Based on the three-dimensional geometric centroid, construct the covariance matrix of each point data in the point cloud of the region of interest, and extract the minimum value in the covariance matrix as the first feature value of the corresponding point. The average of the sum of the angles between the first eigenvalue of each point and the first eigenvalues ​​of other points is obtained to get the rate of change of the local normal vector of the corresponding point; The rate of change of the normal vector within the region of interest is obtained by averaging the rate of change of all local normal vectors in the point cloud of the region of interest. The normal vector ratio feature of the corresponding point is obtained by comparing the ratio of the local normal vector change rate at each point to the normal vector change rate within the region of interest.

5. The method according to claim 3, characterized in that, The curvature ratio features of each point in the point cloud of any region of interest are obtained based on the following steps: Based on each point in the point cloud of the region of interest, the surface is fitted using the least squares method to obtain the local curvature change rate corresponding to each point; The rate of curvature change within the region of interest is obtained by averaging the rates of curvature change of all local points in the point cloud of the region of interest. The curvature ratio feature of the corresponding point is obtained based on the ratio of the local curvature change rate at each point to the curvature change rate within the region of interest.

6. The method according to claim 3, characterized in that, The key points are determined by constructing a Gaussian difference pyramid using scale-invariant features to detect scale-space extrema.

7. The method according to claim 2, characterized in that, Using the vertical direction as the rotation axis, and based on the angle between the first translation vector and the vertical direction, the Rodrigues rotation formula is constructed to obtain the rotation matrix.

8. The method according to claim 2, characterized in that, The local unit point cloud pairs are initially registered using singular value decomposition.

9. The method according to any one of claims 1-8, characterized in that, The registration window includes a moving window and a buffer window, and the size of the buffer window is the sum of the moving window and the estimated maximum horizontal displacement. The moving window is advanced on the first candidate point cloud according to the moving step size, while the buffer window is advanced synchronously on the second candidate point cloud, resulting in multiple pairs of local unit point clouds.

10. A three-dimensional moving extraction system for mining area surface based on improved CPD non-rigid point cloud registration, characterized in that, include: The data acquisition module is used to acquire non-rigid point clouds from the UAV laser point cloud data of the mining area before and after two periods. The feature extraction module is used to extract multi-layer feature points from the non-rigid point clouds of the two phases respectively, to obtain the first phase candidate point cloud set and the second phase candidate point cloud set; The registration module is used to move on two candidate point cloud sets according to the set registration window and movement step size to obtain multiple local unit point cloud pairs. Using the orientation-constrained CPD algorithm, non-rigid registration is performed on each local unit point cloud pair in sequence to obtain the deformation field of each local unit. The output module is used to superimpose the deformation fields of each local unit to obtain the global three-dimensional deformation field.

Images