A method for synchronous extraction of 3D surface movement vectors based on point clouds

Abstract

This invention relates to a method for synchronously extracting three-dimensional surface movement vectors based on point clouds, belonging to the field of surface deformation monitoring technology. It solves the problems of complex calculations, large errors, and inability to reflect the synchronicity of deformation in existing mining area deformation monitoring technologies. The method includes the following steps: acquiring a first-phase cloud and a second-phase point cloud of the mining area; extracting ground feature points from the first-phase point cloud and the second-phase point cloud respectively to obtain first-point cloud data and second-point cloud data; extracting feature point clouds from the first-phase point cloud data and the second-phase point cloud data respectively based on an improved connected component segmentation algorithm; registering the feature point clouds from the first-phase point cloud data and the second-phase point cloud data; and obtaining the three-dimensional movement vector of the mining area surface based on the registered feature point clouds. This achieves accurate and efficient synchronous extraction of three-dimensional surface movement.

Description

Technical Field

[0001] This invention relates to the field of surface deformation monitoring technology, and in particular to a method for synchronous extraction of three-dimensional surface movement vectors based on point clouds. Background Technology

[0002] Unmanned Aerial Vehicle (UAV) LiDAR technology boasts advantages such as high precision, high resolution, and accurate scanning of surface morphology, and is increasingly being applied to monitoring surface subsidence in mining areas. UAV-LiDAR monitoring processes high-precision 3D coordinate point clouds collected by laser scanners in the same area to obtain a digital elevation model (DEM) of the surface, and multi-stage DEM differential analysis yields vertical deformation. However, it lacks the crucial indicator of horizontal surface deformation. Uneven subsidence and horizontal movement are sources of surface curvature and strain, and these changes can damage surface objects.

[0003] To determine the corresponding three-dimensional displacement, multi-source data is typically fused. However, the acquisition time and systematic errors of these multi-source data are inconsistent, making it difficult to obtain spatiotemporally synchronized three-dimensional surface subsidence and deformation. For example, the main method for monitoring mining areas in the past was traditional geodetic techniques, including total station surveying and Global Navigation Satellite System (GNSS) surveying. These techniques can monitor high-precision three-dimensional deformation, but they are time-consuming, labor-intensive, and can only acquire deformation data from low-density points. Therefore, traditional geodetic techniques can accurately but only provide limited information on surface three-dimensional deformation. In recent years, with the development of high spatial resolution remote sensing technologies, such as synthetic aperture radar interferometry (DInSAR), ground-based three-dimensional laser scanning, airborne laser scanning, and UAV photogrammetry, it is possible to monitor the surface deformation of the entire mining area. However, most of these methods rely on DEM differential to obtain one-dimensional deformation, which cannot determine the surface three-dimensional deformation vector. These conventional monitoring processes extract subsidence and horizontal movement separately. This method introduces errors and fails to reflect the synchronicity of deformation. Furthermore, the models constructed using these methods are highly complex or require the fusion of multiple data sources, which is not conducive to engineering practice. Summary of the Invention

[0004] In view of the above analysis, the embodiments of the present invention aim to provide a method and system for synchronous extraction of three-dimensional moving vectors on the ground surface based on point cloud, so as to solve the problems of complex calculation, large error and inability to reflect the synchronicity of deformation monitoring in existing mining areas.

[0005] On one hand, embodiments of the present invention provide a method for synchronous extraction of three-dimensional movement vectors on the ground based on point clouds, including the following steps:

[0006] The first phase cloud and the second phase point cloud of the mining area are obtained, and the ground feature points in the first phase point cloud and the second phase point cloud are extracted respectively to obtain the first point cloud data and the second point cloud data.

[0007] Based on the improved connected component segmentation algorithm, feature feature point clouds are extracted from the first point cloud data and the second point cloud data respectively.

[0008] The feature point clouds in the first point cloud data and the feature point clouds in the second point cloud data are registered; based on the registered feature point clouds, the three-dimensional movement vector of the surface of the mining area is obtained.

[0009] Based on further improvements to the above method, the process of extracting feature point clouds from the first point cloud data and the process of extracting feature point clouds from the second point cloud data are the same. The extraction of feature point clouds from the first point cloud data based on the improved connected component segmentation algorithm includes:

[0010] The first point cloud data is divided into an octree;

[0011] Starting from the topmost leaf node in spatial location, each leaf node is marked sequentially;

[0012] Iterate through each marked leaf node. If there is a marked leaf node in the second neighborhood of the current leaf node, then re-mark the current leaf node according to the marked leaf node.

[0013] The point clouds within leaf nodes with the same label constitute a feature point cloud.

[0014] Based on a further improvement to the above method, each leaf node is marked sequentially, starting from the topmost leaf node in terms of spatial location, including:

[0015] Set the current flag value to 1;

[0016] Starting from the topmost leaf node in spatial location, traverse each leaf node sequentially;

[0017] For the current leaf node, if it contains points from the first point cloud data, then determine whether there is a label among the leaf nodes in the first neighborhood of the current leaf node. If there is, take the smallest label as the label of the current leaf node; if there is no label, mark the current leaf node with the current label value and increment the current label value by 1. The leaf nodes in the first neighborhood of the current leaf node are the leaf nodes connected to the front, back, left, right or top surfaces of the current leaf node, as well as the leaf nodes connected to the edges or vertices of the top surface of the current leaf node.

[0018] Based on a further improvement to the above method, if there is a marked leaf node in the second neighborhood of the current leaf node, then the current leaf node is re-marked according to the marked leaf node, including:

[0019] If there are labels among the leaf nodes in the second neighborhood of the current leaf node, then the smallest label value is taken as the label of the current leaf node; where the leaf nodes in the second neighborhood of the current leaf node are leaf nodes connected to the face, edge or point of the current leaf node.

[0020] Based on further improvements to the above method, a three-dimensional movement vector of the mining area surface is obtained from the registered feature point cloud, including:

[0021] For the region containing the i-th registered feature point cloud pair, the 3D deformation vector [ΔX] is calculated according to the following formula. i ΔY i ΔZ i ] T :

[0022]

[0023]

[0024] Among them, Let be the coordinate matrix of the reference point cloud for centering the registered feature point cloud. This represents the matrix after the coordinate matrix of the reference point cloud has been transformed by the transformation matrix. Let represent the rotation matrix in the transformation matrix of the i-th registered feature point cloud pair. Let represent the translation matrix in the transformation matrix of the i-th registered feature point cloud pair.

[0025] Based on further improvements to the above method, the feature point clouds in the first point cloud data and the feature point clouds in the second point cloud data are registered, including:

[0026] Each feature point cloud in the second point cloud data is taken as the current feature point cloud; based on the center position of the current feature point cloud, the corresponding feature point cloud in the first point cloud data is searched to form a feature point cloud pair;

[0027] For each pair of feature point cloud pairs, the initial transformation matrix is ​​calculated based on the fast point feature histogram;

[0028] Based on the initial transformation matrix, the feature point cloud pairs are accurately registered to obtain the registration result for each feature point cloud pair.

[0029] Based on the further improvement of the above method, for each pair of feature point cloud pairs, the initial transformation matrix is ​​calculated based on the fast point feature histogram, including:

[0030] S321. Extract n feature points from the point cloud to be registered of the feature feature point cloud pair;

[0031] S322. Based on the fast point feature histogram of n feature points, find the corresponding point of each feature point in the reference point cloud of the feature feature point cloud pair;

[0032] S323. Calculate the registration transformation matrix and the corresponding distance error based on the reference point and the corresponding point.

[0033] S324. Repeat steps S321 to S323 until the preset number of iterations is reached, and select the registration transformation matrix corresponding to the distance error and the minimum value during the iteration process as the initial transformation matrix.

[0034] Based on a further improvement of the above method, and based on the initial transformation matrix, the feature point cloud pairs are precisely registered to obtain the registration result for each feature point cloud pair, including:

[0035] S331. Transform the point cloud to be registered for the current feature point cloud pair according to the initial transformation matrix;

[0036] S332, According to the formula Calculate the corresponding point in the reference point cloud of the current feature point cloud pair for each point in the transformed point cloud to be registered, forming a registration point pair; where P i (R,t) represents the point in the point cloud to be registered after transformation according to the initial transformation matrix (R,t), where N Q Q represents the number of points in the reference point cloud. cl (i,R,t) represents the corresponding point in the reference point cloud, Q j Let j represent the j-th point of the reference point cloud, and ||·| represent the 1-norm of the matrix;

[0037] S333, Based on the N registration points with the smallest centering distance... P0 The transformation matrix is ​​calculated for each registration point pair and used as the initial transformation matrix;

[0038] S334. If the iteration termination condition is met, the current initial transformation matrix is ​​the final transformation matrix corresponding to the current feature point cloud pair, and the registration result of the current feature point cloud pair is obtained; otherwise, return to step S331.

[0039] Based on further improvements to the above method, ground feature points are extracted from the first phase cloud and the second phase point cloud to obtain first point cloud data and second point cloud data, including:

[0040] S11. Based on the preset size and the point cloud projected onto the grid, perform adaptive grid division on the mining area plane; for each grid, extract the lowest point mapped to that grid as the initial seed point corresponding to that grid; mark all initial seed points as ground points, and mark the remaining points as points to be judged;

[0041] S12. Construct a triangular network based on ground points;

[0042] S13. For each point to be judged in the point cloud data, determine whether the point to be judged is a ground point based on the distance from the point to be judged to the nearest triangle and the angle between the line connecting the point to be judged to the vertex of the nearest triangle and the nearest triangle. If it is, mark the point to be judged as a ground point and return to step S12. Otherwise, return to step S13 to traverse the next point to be judged.

[0043] S14. If no point to be judged in step S13 is marked as a ground point, then extract all non-ground points in the point cloud data as ground feature points.

[0044] Based on a further improvement to the above method, adaptive meshing of the mining area plane is performed according to a preset size and the point cloud projected onto the mesh, including:

[0045] S111. Divide the mining area plane into an initial grid according to the preset size; project the point cloud onto the mining area plane;

[0046] S112. For a grid with more than one projection point, use the RANSAC algorithm to calculate the fitting plane of the grid, calculate the angle between the normal vector of the fitting plane and the horizontal plane. If the angle is greater than the second threshold, divide the current grid into multiple grids. Repeat step S112 until there are no grids that need to be divided.

[0047] This invention extracts feature features through an improved connected component segmentation algorithm, treating feature features in the deformation zone as rigid bodies. By registering feature features before and after deformation in the mining area, a three-dimensional vector of the mining surface is obtained based on the registration result. This invention realizes a method for synchronously extracting three-dimensional movement vectors of the mining surface. This method can synchronously extract the three-dimensional spatial movement vectors of the mining subsidence surface, thereby reducing errors. Furthermore, it does not require the deployment of complex detection equipment or other data sources, and is computationally simple and easy to implement.

[0048] 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

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

[0050] Figure 1This is a flowchart of a method for synchronous extraction of three-dimensional movement vectors on the ground based on point clouds, according to an embodiment of the present invention.

[0051] Figure 2 This is a schematic diagram showing the distance and included angle between the point to be judged and the nearest triangle in an embodiment of the present invention;

[0052] Figure 3 This is a schematic diagram of ground points and non-ground elements in an embodiment of the present invention;

[0053] Figure 4 This is a schematic diagram of the feature point cloud extracted in an embodiment of the present invention;

[0054] Figure 5 This is a schematic diagram of the three-dimensional deformation results of the mining area in an embodiment of the present invention;

[0055] Figure 6 This is a schematic diagram of a two-dimensional horizontal displacement vector in an embodiment of the present invention;

[0056] Figure 7 This is a schematic diagram of three-dimensional spatial displacement vectors in an embodiment of the present invention;

[0057] Figure 8 This is a schematic diagram of the connected domain segmentation result of the high tower in an embodiment of the present invention. Detailed Implementation

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

[0059] A specific embodiment of the present invention discloses a method for synchronous extraction of three-dimensional movement vectors on the ground surface based on point clouds, such as... Figure 1 As shown, it includes the following steps:

[0060] S1. Obtain the first phase point cloud and the second phase point cloud of the mining area, and extract the ground feature points from the first phase point cloud and the second phase point cloud respectively to obtain the first point cloud data and the second point cloud data.

[0061] S2. Based on the improved connected component segmentation algorithm, feature feature point clouds are extracted from the first point cloud data and the second point cloud data respectively.

[0062] S3. Register the feature point clouds in the first point cloud data and the feature point clouds in the second point cloud data; obtain the three-dimensional movement vector of the surface of the mining area based on the registered feature point clouds.

[0063] During implementation, AV-LiDAR is used to monitor surface deformation in the mining area. Point clouds of rigid objects on the surface, such as high-voltage power line towers, utility poles, and surface buildings and structures, are collected. The movement of these objects reflects surface deformation, becoming an important indicator for evaluating surface deformation monitoring. The first phase point cloud is the point cloud collected before deformation, and the second phase point cloud is the point cloud collected after deformation. Feature features are extracted using an improved connected component segmentation algorithm, treating them as rigid bodies in the deformation area. By registering the feature features before and after deformation, a three-dimensional vector of the mining area surface is obtained based on the registration result. This enables a method for synchronously extracting the three-dimensional spatial movement vector of the mining surface. This method can simultaneously extract the three-dimensional spatial movement vector of the mining subsidence surface, thereby reducing errors. Furthermore, it does not require complex detection equipment or other data sources, is computationally simple, and easy to implement.

[0064] After acquiring the first and second phase point clouds, it is necessary to extract ground feature points to facilitate subsequent registration calculations. In practice, this application identifies ground points and non-ground points, thereby removing ground points to extract ground feature points. Specifically, step S1 extracts ground feature points from the first and second phase point clouds to obtain first and second point cloud data, including:

[0065] S11. Based on the preset size and the point cloud projected onto the grid, perform adaptive grid division on the mining area plane; for each grid, extract the lowest point mapped to that grid as the initial seed point corresponding to that grid; mark all initial seed points as ground points, and mark the remaining points as points to be judged;

[0066] S12. Construct a triangular network based on ground points;

[0067] S13. For each point to be judged in the point cloud data, determine whether the point to be judged is a ground point based on the distance from the point to be judged to the nearest triangle and the angle between the line connecting the point to be judged to the vertex of the nearest triangle and the nearest triangle. If it is, mark the point to be judged as a ground point and return to step S12. Otherwise, return to step S13 to traverse the next point to be judged.

[0068] S14. If no point to be judged in step S13 is marked as a ground point, then extract all non-ground points in the point cloud data as ground feature points.

[0069] Step S11 involves adaptively dividing the mining area plane into grids based on preset dimensions and the point cloud projected onto the grid, including:

[0070] S111. Divide the mining area plane into an initial grid according to the preset size; project the point cloud onto the mining area plane;

[0071] S112. For a grid with more than one projection point, use the RANSAC algorithm to calculate the fitting plane of the grid, calculate the angle between the normal vector of the fitting plane and the horizontal plane. If the angle is greater than the second threshold, divide the current grid into multiple grids. Repeat step S112 until there are no grids that need to be divided.

[0072] During implementation, the first and second thresholds can be set according to the requirements of calculation accuracy and speed.

[0073] During implementation, the corresponding ground feature points are extracted from the first phase point cloud and the second phase point cloud according to the process of steps S11-S14.

[0074] The following explanation uses the second phase point cloud as an example. Initially, the preset size is represented as D. First, the mining area ground (XY plane) is divided according to size D. Then, the second phase point cloud is projected onto the XY plane. Due to varying terrain slopes and different building and tree coverage, some grids have more projected points than others. To ensure computational efficiency while accurately extracting ground features, if the number of projected points in a grid exceeds a first threshold, the RANSAC algorithm is first used to fit the corresponding fitting plane for that grid. It should be noted that this plane fitting is performed based on the original points projected into the grid. After obtaining the fitting plane, the normal vector of the fitting plane is calculated. If the angle between the normal vector of the fitting plane and the horizontal plane is greater than a second threshold, it indicates that the terrain slope within the current grid varies significantly, and the grid should be subdivided for accurate identification. In practice, the grid is divided into four equal grids, and the same method is used to determine whether each subdivided grid needs further subdivision. If all grids either have a number of projection points less than or equal to the first threshold, or the included angle calculated according to step S112 is less than the second threshold, meaning there are no grids that need to be divided, then grid division is stopped.

[0075] Adaptive mesh generation makes the extraction of initial seed points more accurate, thus providing a basis for subsequent identification of ground and non-ground points and improving the accuracy of the identification.

[0076] After the mesh is generated, the lowest point mapped to each mesh is extracted and used as the initial seed point for that mesh. The initial seed point is marked as a ground point, and the remaining points are used as points to be evaluated. A triangular mesh is constructed based on the ground points. In implementation, a TIN irregular triangular mesh can be constructed.

[0077] For each point p to be judged in the second point cloud, see [link to relevant documentation]. Figure 2Find the nearest triangle (V1, V2, V3 are the vertices of the triangle), calculate the distance d from the point to be judged to the triangle, and the three angles α, β, γ between the line connecting the point to be judged to the three vertices of the triangle and the plane containing the triangle. If d is less than a preset distance threshold, and α, β, γ are all less than preset angle thresholds, then mark point p as a ground point, and return to step S12 to update the triangle network. If p is not a ground point, then traverse the next point to be judged.

[0078] If no point to be judged is marked as a ground point in step S13, the iteration ends, and all non-ground points in the point cloud data are extracted as ground feature points.

[0079] The first point cloud data is obtained by extracting ground feature points from the first phase, and the second point cloud data is obtained by extracting ground feature points from the second phase. Feature point clouds can be obtained through connected component analysis. Since UAV-LiDAR acquires point clouds of ground features from the air, the number of points reflected from the upper surface of the object is greater, while the number of point clouds in the middle of the object is relatively smaller. Conventional connected component segmentation can easily divide it into multiple parts. Taking a high-voltage tower as an example, conventional connected component segmentation can easily divide the high-voltage tower into two point cloud sets, such as... Figure 8 As shown in (a) above. To avoid this situation, it is necessary to perform connected component analysis on the point cloud of the upper surface of the object, treating it as a connected whole from top to bottom. Based on this, this invention modifies the conventional connected component algorithm, setting the main rule to analyze connectivity from top to bottom, so that it can extract an object completely from the point cloud, as shown in (a). Figure 8 As shown in (b). In practice, the process of extracting feature point clouds from the first point cloud data and extracting feature point clouds from the second point cloud data based on the improved connected component segmentation algorithm is the same.

[0080] Specifically, feature point clouds of ground objects are extracted from the first point cloud data based on an improved connected component segmentation algorithm, including:

[0081] S21. Perform octree partitioning on the first point cloud data;

[0082] S22. Starting from the topmost leaf node in spatial location, mark each leaf node in turn;

[0083] S23. Traverse each marked leaf node. If there is a marked leaf node in the second neighborhood of the current leaf node, then re-mark the current leaf node according to the marked leaf node.

[0084] S24. Point clouds within leaf nodes with the same label constitute a feature point cloud.

[0085] Connectivity analysis cannot be directly applied to unordered point clouds; an octree is needed to convert the unordered point cloud into an ordered point cloud (octree point cloud). Each segmented cube is a leaf node of the octree.

[0086] To analyze connectivity from top to bottom, the nodes are labeled sequentially starting from the topmost leaf node in spatial location. Step S22 specifically includes:

[0087] Set the current flag value to 1;

[0088] Starting from the topmost leaf node in spatial location, traverse each leaf node sequentially;

[0089] For the current leaf node, if it contains points from the first point cloud data, then determine whether there is a label among the leaf nodes in the first neighborhood of the current leaf node. If there is, take the smallest label as the label of the current leaf node; if there is no label, mark the current leaf node with the current label value and increment the current label value by 1. The leaf nodes in the first neighborhood of the current leaf node are the leaf nodes connected to the front, back, left, right or top surfaces of the current leaf node, as well as the leaf nodes connected to the edges or vertices of the top surface of the current leaf node.

[0090] Initially, the current label value is set to 1. When traversing to the first leaf node containing a point (the cube corresponding to the leaf node contains a point), the leaf node containing the first point cloud is labeled with 1, and the current label value is set to 2. When traversing to the second leaf node containing a point, it is checked whether there is a labeled leaf node in its first neighborhood. If so, the smallest label is taken as the label of the current leaf node. Since the smallest label in its first neighborhood is 1, the leaf node containing the second point is labeled with 1. If there is no labeled leaf node in the first neighborhood of the leaf node containing the second point, the leaf node containing the second point is labeled with 2, and the current label value is set to 3. This process continues until all leaf nodes have been traversed from top to bottom, thus achieving preliminary connected component analysis.

[0091] After the first traversal, there's a problem of labeling leaf nodes within the same connected component with different tags. Therefore, it's necessary to merge adjacent regions with different tags into the same connected component. Thus, for each labeled leaf node, if there's a labeled leaf node in its second neighborhood, the current leaf node is re-labeled based on that labeled leaf node. Specifically, this includes:

[0092] If there are labels among the leaf nodes in the second neighborhood of the current leaf node, then the smallest label value is taken as the label of the current leaf node; where the leaf nodes in the second neighborhood of the current leaf node are leaf nodes connected to the face, edge or point of the current leaf node.

[0093] By performing two traversals, adjacent regions with different labels are merged into a single label. The point cloud within the leaf node with the same label is a feature point cloud, thus achieving accurate and complete feature point cloud extraction.

[0094] After extracting the feature point clouds from the first point cloud data, the same process is followed to extract the feature point clouds from the second point cloud data. Each feature point cloud represents a rigid body. By pairing the rigid bodies from the second phase with those from the first phase, the surface deformation at the corresponding locations is calculated.

[0095] Specifically, in step S3, the feature point clouds in the first point cloud data and the feature point clouds in the second point cloud data are registered, including:

[0096] S31. Sequentially take each feature point cloud in the second point cloud data as the current feature point cloud; based on the center position of the current feature point cloud, find the corresponding feature point cloud in the first point cloud data to form a feature point cloud pair;

[0097] S32. For each pair of feature point cloud pairs, calculate the initial transformation matrix based on the fast point feature histogram.

[0098] S33. Based on the initial transformation matrix, perform precise registration on the feature point cloud pairs to obtain the registration result for each feature point cloud pair.

[0099] During implementation, the feature point clouds in the first point cloud serve as the reference point cloud, while the feature point clouds in the second point cloud serve as the point clouds to be registered. For each feature point cloud in the second point cloud data, its center point is first calculated. Then, the feature point cloud closest to the center point is found in the first point cloud data, forming a feature point cloud pair (P, Q), where P represents the point cloud to be registered and Q represents the reference point cloud, resulting in multiple feature point cloud pairs.

[0100] For each pair of feature point cloud pairs, the initial transformation matrix is ​​calculated based on the fast point feature histogram, specifically including:

[0101] S321. Extract n feature points from the point cloud to be registered of the feature feature point cloud pair;

[0102] S322. Based on the fast point feature histogram of n feature points, find the corresponding point of each feature point in the reference point cloud of the feature feature point cloud pair;

[0103] S323. Calculate the registration transformation matrix and the corresponding distance error based on the reference point and the corresponding point.

[0104] S324. Repeat steps S321 to S323 until the preset number of iterations is reached, and select the registration transformation matrix corresponding to the distance error and the minimum value during the iteration process as the initial transformation matrix.

[0105] During implementation, firstly, a fast point feature histogram is calculated for each point in the point cloud to be registered, and then a fast point feature histogram is calculated for each point in the point cloud to be registered. Next, n feature points are extracted from the point cloud to be registered to solve for the transformation matrix. To make the registration more accurate, the distance between the n feature points extracted from the point cloud to be registered should be greater than a preset distance threshold.

[0106] Based on a fast point feature histogram of n feature points, the corresponding point for each feature point is found in the reference point cloud of the feature feature point cloud pair. This includes: finding the k points in the reference point cloud that are most similar to the fast point feature histogram of the current feature point, and randomly selecting one of these points as the corresponding point of the current feature point. In practice, the similarity can be measured using cosine similarity.

[0107] After obtaining n pairs of feature points and their corresponding points, the registration transformation matrix can be solved using singular value decomposition. The distance error and M can then be calculated using the following formula:

[0108]

[0109]

[0110] Among them, l i Let m represent the distance error between the i-th feature point and its corresponding point calculated using the registration transformation matrix, m represent the preset threshold, and ||·|| represent the 1-norm of the matrix.

[0111] Repeat steps S321 to S323 until the preset number of iterations is reached. Select the registration transformation matrix corresponding to the minimum distance error from the set number of iterations. This matrix is ​​the initial transformation matrix (R,t), where R represents the rotation matrix and t represents the translation matrix.

[0112] The initial transformation matrix is ​​not precise and requires further precise registration to obtain the final registration result. Specifically, based on the initial transformation matrix, precise registration is performed on the feature feature point cloud pairs to obtain the registration result for each pair, including:

[0113] S331. Transform the point cloud to be registered for the current feature point cloud pair according to the initial transformation matrix;

[0114] S332, According to the formula Calculate the corresponding point in the reference point cloud of the current feature point cloud pair for each point in the transformed point cloud to be registered, forming a registration point pair; where P i(R,t) represents the point in the point cloud to be registered after transformation according to the initial transformation matrix (R,t), where N Q Indicates the number of points in the reference point cloud; Q cl (i,R,t) represents the corresponding point in the reference point cloud; Q j Let j represent the j-th point of the reference point cloud; ||·| represents the 1-norm of the matrix;

[0115] S333, Based on the N registration points with the smallest centering distance... P0 The transformation matrix is ​​calculated for each registration point pair and used as the initial transformation matrix;

[0116] S334. If the iteration termination condition is met, the current initial transformation matrix is ​​the final transformation matrix corresponding to the current feature point cloud pair, and the registration result of the current feature point cloud pair is obtained; otherwise, return to step S331.

[0117] After obtaining the initial transformation matrix (R,t), the point cloud P to be registered is first transformed according to the initial transformation matrix. The point representation P of the i-th point in the point cloud P after transformation is obtained. i (R,t). For each transformed point, find the nearest point Q in the reference point cloud according to the formula in step S332. j P i (R,t) and Q j By constructing registration point pairs, multiple sets of registration point pairs are obtained.

[0118] The correspondences between these registration point pairs are not necessarily all correct; incorrect correspondences will affect the final registration result. Therefore, the N points with the smallest distance among the registration point pairs are used to determine the correct registration method. P0 The transformation matrix is ​​calculated for each registration point pair and used as the initial transformation matrix.

[0119] During implementation, N P0 Through N P0 =ξN P Calculate, where ξ represents the overlap, and N P This indicates the number of points in the point cloud to be registered. The overlap can be preset according to the required calculation accuracy.

[0120] During implementation, the distance d between registration point pairs is calculated using the following formula. i (R,t)=||Q cl (i,R,t)-P i (R,t)||. Extract the N values ​​with the smallest distance. P0 The transformation matrix (R, t) of each registration point can be solved using the least squares method and used as the current initial transformation matrix. The N points with the smallest distance values ​​are then selected. P0 The sum of the distances between the registration points is denoted as St.

[0121] During implementation, when the maximum number of iterations is reached, or the current distance value is the smallest N, the implementation will proceed. P0 The distances between the registration points and the N points that have the smallest distance value compared to the previous iteration. P0 If the difference between the sum of the distances of the registration points meets the preset accuracy requirement, then the iteration termination condition is met, and the iteration terminates. The current initial transformation matrix is ​​the final transformation matrix corresponding to the current feature point cloud pair.

[0122] Calculate the final transformation matrix for each pair of feature point clouds. The i-th pair of feature point clouds is represented as (P i Q i The corresponding transformation matrix is ​​represented as (R). i ,t i ).

[0123] Due to mining activity, the spatial positions of the rigid body point clouds in the second phase are obtained by moving (translating and rotating) the same rigid body point clouds from the first phase. The transformation matrix corresponding to the feature point clouds can describe the rigid body movement process of the point clouds, and the transformation matrix includes the point cloud coordinate system. rotation matrix and translation vector

[0124] After obtaining the transformation matrix of the i-th group of feature point clouds, the three-dimensional deformation vector of the region where the i-th group of feature point clouds is located is calculated according to the following formula:

[0125]

[0126] Among them, For the reference point cloud coordinate matrix, The coordinates of the registered point cloud, i.e., the coordinate matrix of the reference point cloud after transformation by the transformation matrix, are equal to... This is the 3D deformation vector of the region where the obtained feature point cloud is located. Wherein, Let be the matrix consisting of the X-axis coordinates of all points in the reference point cloud of the i-th feature point cloud pair. Let be the matrix consisting of the Y-axis coordinates of all points in the reference point cloud of the i-th feature point cloud pair. It is a matrix consisting of the Z-axis coordinates of all points in the reference point cloud of the i-th feature point cloud pair. Let represent the rotation matrix in the transformation matrix of the i-th registered feature point cloud pair. Let represent the translation matrix in the transformation matrix of the i-th registered feature point cloud pair.

[0127] For areas without feature point clouds, the corresponding deformation can be calculated using inverse distance weighted difference, thereby obtaining the three-dimensional surface movement vector of the entire mining area.

[0128] The effectiveness of this method is illustrated using the 3S201 working face at the Wangjiata Coal Mine in the Ordos Coalfield as an example. This working face has a strike length of 1050m, a dip length of 250m, an average mining depth of 300m, an average mining thickness of 3.25m, and an average coal seam dip angle of 2°. Mining began on October 7, 2020, and ceased on January 20, 2021. The area is hilly with significant topographic relief. Figure 7 As shown. UAV-LiDAR was used for scanning and monitoring. The first phase of monitoring was carried out on August 24, 2020, and the second phase of monitoring was carried out on March 16, 2021. The flight altitude was 60m, and the observation range was determined to be 200m outside the working face using the boundary angle. Figure 3 In the diagram, (a) represents the obtained ground points and non-ground points, (b) represents ground points, (c) represents non-ground points, which is the rigid body point cloud, and (d) represents local details.

[0129] The improved connected component algorithm described above is used to segment suitable feature point clouds from non-terrestrial point clouds. The segmented point cloud is shown below. Figure 4 As shown. Figure 4 (a) shows the overall segmentation result, (b) is a top view of the segmented houses, and (c) shows the segmented high-voltage towers and houses, etc. A total of 1233 rigid bodies were segmented, with a total point cloud count of 10,211,942. In the study area, there was at least one point per square meter in 75.8% of the area, and in 95.2% of the area, the point cloud density was greater than 8 points / m². 2 .

[0130] The final three-dimensional deformation result is as follows Figure 5 As shown in the diagram, the deformation component Δx indicates that the surface mainly moves in the positive X-axis direction, which is also the direction of tunneling in the mining area. Near the stop-mining line, the movement direction is reversed. The maximum X-axis movement occurs at the lower left and upper right corners of the mining area, at 0.53m and -0.65m respectively. The deformation component Δy roughly approaches the strike line, with the maximum deformation values ​​occurring near the mining boundary at 0.65m and -0.75m respectively. The maximum Z-axis movement is -2.723m, with the maximum value located in the center of the basin.

[0131] Let the positive directions of the X-axis and Y-axis be the positive directions of Δx and Δy, respectively, and let the downward direction of the z-axis be the positive direction of Δz. The calculated horizontal displacement value is the composite two-dimensional horizontal displacement vector of the vector magnitudes. Figure 6 ,by The calculated value is the composite three-dimensional spatial displacement vector of vector magnitudes. Figure 7 The white arrow indicates the direction of horizontal displacement; from Figure 6 It can be seen that the horizontal movement vector points towards the center of the basin, with a maximum horizontal movement vector of 823 mm, located at the mining boundary. The horizontal displacement is relatively small in the central basin area of ​​the working face. Figure 7 It can be seen that the three-dimensional spatial movement vector of the surface points to the center of the goaf, with the largest movement at the center of the basin and the smallest movement at the boundary. The overall surface movement in the mining area conforms to the mining subsidence pattern, indicating the correctness of extracting the surface spatial movement vector.

[0132] The three-dimensional movement of multiple points was measured using GPS, and this data was used as a benchmark to evaluate the accuracy of the three-dimensional movement deformation extracted in the experiment. The data and error calculation results are shown in Table 1. The mean square errors in the X, Y, and Z directions are 54 mm, 83 mm, and 36 mm, respectively, indicating that the method proposed in this invention is feasible and the obtained three-dimensional movement vector data is reliable. Since the vertical accuracy of UAV-LiDAR is higher than its horizontal accuracy, the accuracy of the surface deformation Δz in the experiment is also significantly higher than its horizontal accuracy.

[0133] Table 1 Measured Deformation and Registration Deformation

[0134]

[0135]

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

[0137] 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 synchronously extracting ground surface three-dimensional movement vectors based on point clouds, characterized in that, The method comprises the following steps: obtaining first period point cloud and second period point cloud of a mining area, and extracting feature ground points in the first period point cloud and the second period point cloud to obtain first point cloud data and second point cloud data; extracting feature ground points in the first point cloud data and the second point cloud data based on an improved connected domain segmentation algorithm; registering the feature ground points in the first point cloud data and the feature ground points in the second point cloud data, and obtaining three-dimensional movement vectors of the ground surface of the mining area based on the registered feature ground points; the process of extracting the feature ground points in the first point cloud data and the feature ground points in the second point cloud data is the same, and the feature ground points in the first point cloud data are extracted based on the improved connected domain segmentation algorithm, which comprises: performing octree segmentation on the first point cloud data; starting from the uppermost leaf node in the spatial position and sequentially marking each leaf; traversing each marked leaf node, and if there is a marked leaf node in the second neighborhood of the current leaf node, re-marking the current leaf node according to the marked leaf node; the point cloud in the leaf nodes with the same mark constitutes a feature ground point cloud; starting from the uppermost leaf node in the spatial position and sequentially marking each leaf, which comprises: setting the current mark value to 1; traversing each leaf node starting from the uppermost leaf node in the spatial position; for the current leaf node, if the points in the first point cloud data are contained therein, judging whether there is a mark in the leaf nodes in the first neighborhood of the current leaf node, if there is, taking the smallest mark as the mark of the current leaf node, if not, marking the current leaf node as the current mark value, and increasing the current mark value by 1; wherein the leaf nodes in the first neighborhood of the current leaf node are the leaf nodes connected with the front surface, rear surface, left surface, right surface or upper surface of the current leaf node, and the leaf nodes connected with the edge or vertex of the upper surface of the current leaf node; if there is a marked leaf node in the second neighborhood of the current leaf node, re-marking the current leaf node according to the marked leaf node, which comprises: if there is a mark in the leaf nodes in the second neighborhood of the current leaf node, taking the smallest mark value as the mark of the current leaf node; wherein the leaf nodes in the second neighborhood of the current leaf node are the leaf nodes connected with the surface, edge or point of the current leaf node.

2. The method for extracting ground surface three-dimensional movement vector synchronously based on point cloud according to claim 1, characterized in that, obtaining three-dimensional movement vectors of the ground surface of the mining area based on the registered feature ground points, which comprises: For the region where the i-th set of registered feature object point clouds is located, the three-dimensional deformation vector is calculated according to the following formula : ； ； wherein, wherein is a coordinate matrix of the reference point cloud in the i-th set of registered pairs of feature ground point clouds, denotes a matrix of the coordinate matrix of the reference point cloud transformed by the transformation matrix, denotes a rotation matrix in the transformation matrix of the i-th set of registered pairs of feature ground point clouds, denotes a translation matrix in the transformation matrix of the i-th set of registered pairs of feature ground point clouds. 3.The method of claim 1, wherein, registering the feature ground points in the first point cloud data and the feature ground points in the second point cloud data, which comprises: sequentially taking each feature ground point cloud in the second point cloud data as a current feature ground point cloud; finding the corresponding feature ground point cloud in the first point cloud data according to the center position of the current feature ground point cloud to form a feature ground point cloud pair; for each feature ground point cloud pair, calculating an initial transformation matrix based on a fast point feature histogram; based on the initial transformation matrix, accurately registering the feature ground point cloud pair to obtain the registration result of each feature ground point cloud pair.

4. The method of claim 3, wherein, For each pair of feature object point cloud, an initial transformation matrix is calculated based on a fast point feature histogram, comprising: S321, extracting n feature points from the to-be-registered point cloud of the feature object point cloud pair; S322, finding the corresponding point of each feature point from the reference point cloud of the feature object point cloud pair based on the fast point feature histogram of the n feature points; S323, calculating a registration transformation matrix and a corresponding distance error sum according to the reference point and the corresponding point; S324, repeating steps S321 to S323 until a preset iteration number is reached, and selecting the registration transformation matrix corresponding to the minimum distance error sum in the iteration process as the initial transformation matrix.

5. The method of claim 3, wherein, Based on the initial transformation matrix, the feature object point cloud pair is accurately registered to obtain the registration result of each pair of feature object point cloud, comprising: S331, transforming the to-be-registered point cloud of the current feature object point cloud pair according to the initial transformation matrix; S332、According to the formula Each point in the transformed to-be-registered point cloud is calculated in the corresponding point in the reference point cloud of the current feature ground point cloud pair, forming a registration point pair; wherein, represents the i-th point in the to-be-registered point cloud after transformation according to the initial transformation matrix (R, t), represents the number of points of the reference point cloud, represents the corresponding point in the reference point cloud, represents the j-th point of the reference point cloud, represents the 1-norm of the matrix, R represents the rotation matrix, and t represents the translation matrix; S333、according to the registration point pair whose distance is the smallest calculating a transformation matrix as the initial transformation matrix for the registration point pair S334, if the iteration termination condition is met, the current initial transformation matrix is the final transformation matrix corresponding to the current feature object point cloud pair, and the registration result of the current feature object point cloud pair is obtained; otherwise, return to step S331.

6. The method of claim 1, wherein, The feature object points in the first period point cloud and the second period point cloud are extracted to obtain first point cloud data and second point cloud data, comprising: S11, performing adaptive grid division on the mine plane according to a preset size and the point cloud projected on the grid; for each grid, extracting the lowest point mapped to the grid as the initial seed point corresponding to the grid; marking all initial seed points as ground points, and marking the remaining points as to-be-judged points; S12, constructing a triangular mesh based on the ground points; S13, for each to-be-judged point in the point cloud data, determining whether the to-be-judged point is a ground point according to the distance of the to-be-judged point to the nearest triangle and the included angle of the line connecting the to-be-judged point and the vertex of the nearest triangle with the nearest triangle, if yes, marking the to-be-judged point as a ground point, returning to step S12, otherwise returning to step S13 to traverse the next to-be-judged point; S14, if no to-be-judged point is marked as a ground point in step S13, extracting all non-ground points in the point cloud data as feature object points.

7. The method of claim 6, wherein, Adaptive grid division is performed on the mine plane according to a preset size and the point cloud projected on the grid, comprising: S111, performing initial grid division on the mine plane according to a preset size; projecting the point cloud on the mine plane; S112, for a grid with a number of projected points greater than a first threshold, using a ransac algorithm to calculate a fitting plane of the grid, calculating the included angle of the normal vector of the fitting plane with the horizontal plane, and if the included angle is greater than a second threshold, the current grid is divided into multiple grids, and step S112 is repeatedly executed until there is no grid that needs to be divided.

Images