A method and apparatus for extracting three-dimensional deformation data of mining subsidence in desert areas

Abstract

This invention discloses a method and apparatus for extracting three-dimensional deformation data of mining subsidence in desert areas, relating to the field of Earth observation technology. The invention calculates the angle between the normal vectors of each point in the point cloud data and performs tangent plane projection analysis to form multiple pairs of point cloud contour sets. The overall approach is based on the local geometric characteristics of the point cloud, accurately extracting feature points reflecting the shape of the point cloud contours. Then, using these multiple pairs of point cloud contour sets as registration objects, the K-4PCS and VGICP registration algorithms are employed to register the three phases of point cloud data, ensuring precise alignment and obtaining the registered point cloud and transformation matrix. The three-dimensional deformation of the land surface is extracted from the transformation parameters in the transformation matrix. This process eliminates the need for subsidence extraction based on point DEMs and planar displacement extraction based on rigid body matching, enabling accurate extraction of planar displacement in desert areas and precise acquisition of three-dimensional deformation data of the land surface caused by mining in desert regions.

Description

Technical Field

[0001] This invention relates to the field of Earth observation technology, and in particular to a method, apparatus, equipment and medium for extracting three-dimensional deformation data of mining subsidence in desert areas. Background Technology

[0002] Coal resources, as an important energy source and raw material for my country's economic and industrial development, are a key element in ensuring the vitality of national economic development, a "ballast stone" for the stable supply of national energy security, and an important tool for achieving the "dual carbon" target. While large-scale mining of coal resources meets energy demand, it also causes geological and environmental problems such as surface subsidence, crop yield reduction, soil erosion, and vegetation degradation, causing great harm and distress to the production and lives of people in mining areas. Therefore, secondary disasters caused by coal mining subsidence have become an important factor affecting the stability and economic development of mining areas. It is necessary to determine the basic laws of mining subsidence through scientific monitoring and research to provide technical support for disaster prevention and control.

[0003] For a long time, scholars both at home and abroad have been constantly innovating methods for monitoring surface subsidence characteristics and have achieved fruitful results in the study of subsidence mechanisms. Traditional mining area subsidence monitoring mainly relies on leveling instruments, total stations, and real-time dynamic carrier phase differential technology (Real-time... The primary methods used are kinematic and RTK (Real-Time Kinematic) surveying, but these methods require significant manpower and resources and are too difficult to implement, leading to their gradual replacement by emerging measurement technologies. Some researchers have proposed the SGI-SF method to extract three-dimensional displacement results from InSAR datasets. Results show that the average root mean square error (RMSE) of this method in the vertical and horizontal directions is approximately 9.28 mm and 13.10 mm, respectively, far less than the displacement caused by mining. Other researchers have used SBAS-InSAR to invert the advance influence angle of mining subsidence, providing important reference for the dynamic protection of surface buildings and facilities in western mining areas. Some researchers have proposed a monitoring method using InSAR and a surface displacement vector attenuation angle model to obtain more detailed and accurate deformation information for the entire basin. Still others have combined D-InSAR and UAV technology to monitor deformation in mining areas, obtaining surface movement and deformation, with the fused data error within 10 cm. All the research described above focuses on InSAR monitoring. Monitoring has significant advantages in acquiring small-scale deformation data, but it often faces limitations when dealing with large deformations (m-level) caused by coal mining due to correlation issues. Therefore, airborne photogrammetry and lidar technology have been introduced into coal mining subsidence monitoring. Some researchers have acquired ground point cloud data using UAV-LiDAR and applied the Kriging interpolation method to differentially process the DEM to obtain ground subsidence data, achieving a ground DEM accuracy of 15mm and a root mean square error of 39mm for the ground subsidence basin. However, UAV-LiDAR technology generates a large amount of point cloud data, and the processing and analysis of this data requires high-performance computing resources. Furthermore, the accuracy of the Kriging interpolation method is highly dependent on the distribution and quantity of data points. If the data points are unevenly distributed or insufficient in quantity, the accuracy of the interpolation results will be reduced. In addition, the Kriging interpolation method requires the calculation of a semivariance function model and interpolation calculation based on this function. However, the selection of the semivariance function model often requires subjective judgment based on the characteristics and distribution of the data, which increases the subjectivity and uncertainty of model selection.

[0004] To address the problems existing in UAV-LiDAR technology and Kriging interpolation methods, numerous point cloud physical extraction algorithms have emerged, promoting the application of laser scanning data in mining monitoring. Some researchers have proposed clustering segmentation algorithms that optimize the concavity and convexity of supervoxels to improve the accuracy of feature segmentation in complex terrain features of mining areas. Other researchers have proposed rigid body registration algorithms based on connected component segmentation to simultaneously acquire 3D surface movement, solving the problem caused by insufficient rigid bodies. Although these methods have successfully solved the problems of sinking extraction based on point DEMs and planar displacement extraction based on rigid body matching, with the westward shift of mining activities, a large amount of mining activity has entered sparsely populated desert areas. These areas lack artificial structures and have few or no rigid bodies, making it difficult for current methods to extract planar displacement in desert areas to obtain 3D surface deformation data caused by mining. Summary of the Invention

[0005] This invention provides a method and apparatus for extracting three-dimensional deformation data of mining subsidence in desert areas. It can solve the problem that existing methods are difficult to extract planar displacement in desert areas in order to obtain three-dimensional deformation data of the surface caused by mining in desert areas.

[0006] This invention provides a method for extracting three-dimensional deformation data of mining subsidence in desert areas, comprising the following steps:

[0007] Three phases of surface point cloud data corresponding to different shooting time points in a desert area were obtained, and overlapping areas were extracted from the three phases of surface point cloud data; wherein, the surface point cloud data corresponding to the overlapping areas consisted of ground points and non-ground points;

[0008] Obtain the angle between the normal vectors of each point in the point cloud data of the ground points in each period of the overlapping area and the normal vectors of the corresponding points in the other two periods, as well as the angle between each point and the projection of the tangent plane of its own point normal vector.

[0009] Points with a normal vector angle higher than a preset value are set as contour feature points; and among the three points formed by the point with a normal vector angle higher than the preset value and the other two corresponding points, the point with the largest angle with the projection of the tangent plane of its respective point's normal vector is set as a contour feature point; wherein, the point with a normal vector angle higher than the preset value and the point with the largest angle with the projection of the tangent plane of its respective point's normal vector are used as a pair of point cloud contours, and multiple pairs of point cloud contour sets are obtained.

[0010] Using the first phase point cloud data as the original point cloud, for each point cloud data in the original point cloud, the point cloud data in the other two phases that belong to the same pair of point cloud contours as each point cloud data in the original point cloud are registered to obtain three-dimensional deformation data that characterizes the displacement of the point cloud during the registration process.

[0011] Preferably, the division of ground points and non-ground points includes:

[0012] After extracting the surface point cloud data from the overlapping areas in the three phases of surface point cloud data, the progressively encrypted triangular network filtering algorithm (IPTD) is used to divide the surface point cloud data in the overlapping areas into ground points and non-ground points.

[0013] Preferably, after dividing the surface point cloud data corresponding to the overlapping region into ground points and non-ground points, the surface point cloud data of the ground points corresponding to the overlapping region is constructed into a KD-tree structure, including:

[0014] For the point cloud data of ground points corresponding to the overlapping areas of the three phases, a cube bounding box containing all point clouds is established based on the global coordinate system of the point cloud. For cubes containing more than one point, a segmentation plane is constructed.

[0015] The points in the partitioned subspace and the partitioned plane form branches and connection points to construct the KD-tree data structure.

[0016] Preferably, the step of obtaining the angle between the normal vectors of each point in the point cloud data of each period of the overlapping area and the normal vectors of corresponding points in the other two periods, as well as the angle between each point and the projection of the tangent plane of its respective point normal vector, includes:

[0017] In a kd-tree structure, a search radius r is set, and the neighborhood points within the search radius r are denoted as a set. ,Right now ;

[0018] Define the surface equation , ,Pick The corresponding set ,calculate Distance to the surface Solve The eigenvector corresponding to the minimum is the normal vector n of the corresponding point;

[0019] The angle between the normal vectors of adjacent points is expressed as ;

[0020] according to And its normal vector n, make a tangent plane at the corresponding point And project the normal vector n onto the tangent plane. The tangent plane that forms the normal vector is projected onto the kd-tree structure. A search radius r is set, and the neighborhood points within the search radius r are denoted as a set. ,Right now ;

[0021] Define the surface equation , ,Pick The corresponding set ,calculate Distance to the surface Solve The eigenvector corresponding to the minimum is the normal vector n of the corresponding point;

[0022] The angle between the normal vectors of adjacent points is expressed as ;

[0023] according to And its normal vector n, make a tangent plane at the corresponding point And project the normal vector n onto the tangent plane. The projection of the tangent plane onto the normal vector is formed.

[0024] Preferably, forming multiple pairs of point cloud contour sets includes:

[0025] Points with an angle greater than 35° between the normal vectors of adjacent points are set as contour feature points;

[0026] according to Construct a tangent plane at this point along with its normal vector n. , will set Projecting points within the plane onto the tangent plane Above, recorded as ,exist Select a point ,by Let u be the u-axis and n be the n-axis. axis, With the v-axis, Construct a local coordinate system for the coordinate center, denoted as Calculate the sets respectively Other points Time vector The clockwise angle with the coordinate axis u The set of included angles is obtained by subtracting each pair of adjacent included angles. ,in , will set Arrange the elements in descending order and find the largest included angle. , the largest included angle Set as contour feature points;

[0027] Multiple pairs of point cloud contour sets are formed based on the two types of contour feature points.

[0028] Preferably, the acquisition of the three-dimensional deformation data includes:

[0029] Set the target point cloud as the point set to be registered P, and set the point cloud to be registered as the reference point set Q;

[0030] For the registration point set P and the reference point set Q, a voxelized grid-based filter is used for sampling, and then the key point set is extracted by the Difference of Gaussian (DoG) key point detector. , ; indicates as:

[0031] ;

[0032] in: Indicates a DoG response. The fuzziness level is indicated by x, y, and z, which represent the three-dimensional coordinates.

[0033] In the key point set A set of 4 coplanar points randomly selected from the data. Searching for the corresponding set using affine invariance Obtain the initial transformation matrix;

[0034] Based on the initial transformation matrix, the similarity between point clouds is calculated through VGICP algorithm voxelization and generalized distance metric strategy, and the point clouds are gradually aligned to obtain the transformation matrix during the alignment process. The transformation matrix during the alignment process is superimposed with the initial transformation matrix to obtain the final transformation matrix, which is expressed as:

[0035] ;

[0036] The final transformation matrix includes a rotation matrix of R×3 and a translation vector of 3×1 for the point cloud coordinate system; where the rotation matrix R represents the directional change of the point cloud in three-dimensional space, and the translation vector represents the distance the point cloud moves in the X, Y, and Z directions.

[0037] Based on the movement distance of the point cloud in the X, Y, and Z directions, three-dimensional deformation data of mining subsidence in desert areas are generated.

[0038] This invention also provides a device for extracting three-dimensional deformation data of mining subsidence in desert areas, comprising:

[0039] The data acquisition module is used to acquire three phases of surface point cloud data corresponding to different shooting time points in the desert area, and extract the overlapping areas in the three phases of surface point cloud data; wherein, the surface point cloud data corresponding to the overlapping areas consists of ground points and non-ground points;

[0040] The contour feature point determination module is used to obtain the angle between the normal vectors of each point in the point cloud data of the ground points in each period of the overlapping area and the corresponding points in the other two periods, as well as the angle between each point and the projection of the tangent plane of its own point normal vector.

[0041] Points with a normal vector angle higher than a preset value are set as contour feature points; and among the three points formed by the point with a normal vector angle higher than the preset value and the other two corresponding points, the point with the largest angle with the projection of the tangent plane of its respective point's normal vector is set as a contour feature point; wherein, the point with a normal vector angle higher than the preset value and the point with the largest angle with the projection of the tangent plane of its respective point's normal vector are used as a pair of point cloud contours, and multiple pairs of point cloud contour sets are obtained.

[0042] The registration module is used to take the first phase point cloud data as the original point cloud. For each point cloud data in the original point cloud, it registers the point cloud data in the other two phases that belong to the same pair of point cloud contours as each point cloud data in the original point cloud, and obtains three-dimensional deformation data that characterizes the displacement of the point cloud during the registration process.

[0043] This invention also provides an electronic device, including a memory and a processor;

[0044] The memory is used to store computer programs;

[0045] When the processor executes the computer program stored in the memory, it implements the steps of the method for extracting three-dimensional deformation data of mining subsidence in desert areas as described above.

[0046] This invention also provides a computer-readable storage medium for storing a computer program, which, when executed by a processor, implements the steps of a method for extracting three-dimensional deformation data of mining subsidence in desert areas as described above.

[0047] This invention provides a method and apparatus for extracting three-dimensional deformation data of mining subsidence in desert areas. Compared with the prior art, its advantages are as follows:

[0048] This invention collects surface point cloud data from three different time points in a study area, extracts point cloud data of ground points in overlapping areas, constructs KD-tree structures for each of the three point cloud data points, and performs normal vector angle calculation and tangent plane projection analysis on each point in the point cloud data to form multiple pairs of point cloud contour sets. Then, using these multiple pairs of point cloud contour sets as registration objects, the K-4PCS and VGICP registration algorithms are used to register the point cloud data from the three periods, obtaining the registered point clouds and transformation matrices. The three-dimensional surface deformation is then extracted from the transformation parameters in the transformation matrix. This process... In this process, when calculating the angle between the normal vectors and performing tangent plane projection analysis on each point in the point cloud data, the local geometric characteristics of the point cloud are used to accurately extract feature points that reflect the shape of the point cloud contour. Registration is performed based on K-4PCS and VGICP registration algorithms to accurately align the point cloud data and obtain the transformation matrix of the three-dimensional transformation data representing the displacement of the point cloud during the registration process. Without the need for sinking extraction based on point DEM and planar displacement extraction based on rigid body matching, the planar displacement of the desert area can be accurately extracted to accurately obtain the three-dimensional deformation data of the surface caused by mining in the desert area. Attached Figure Description

[0049] Figure 1 A schematic diagram illustrating a method for extracting three-dimensional deformation data of mining subsidence in desert areas, provided by an embodiment of the present invention;

[0050] Figure 2 A schematic diagram of the UAV flight path during the data acquisition process of a method for extracting three-dimensional deformation data of mining subsidence in desert areas, provided in an embodiment of the present invention;

[0051] Figure 3 A schematic diagram of the original point cloud data collected on-site for a method of extracting three-dimensional deformation data of mining subsidence in desert areas, provided in an embodiment of the present invention;

[0052] Figure 4 This is a schematic diagram of ground points and non-ground points obtained by classifying point cloud data from different periods using a progressively encrypted triangular mesh filtering algorithm, which is provided in an embodiment of the present invention for a method of extracting three-dimensional deformation data of mining subsidence in desert areas; where (a) represents the filtering result; (b) represents ground points; and (c) represents non-ground points.

[0053] Figure 5 This is a schematic diagram of the registration process for a method of extracting three-dimensional deformation data of mining subsidence in desert areas, provided by an embodiment of the present invention.

[0054] Figure 6 This is a partial point cloud diagram illustrating a method for extracting three-dimensional deformation data of mining subsidence in desert areas, provided in an embodiment of the present invention.

[0055] Figure 7A KD-tree point cloud diagram illustrating a method for extracting three-dimensional deformation data of mining subsidence in desert areas, provided in an embodiment of the present invention;

[0056] Figure 8 This is a schematic diagram of the overall point cloud of a method for extracting three-dimensional deformation data of mining subsidence in desert areas, provided by an embodiment of the present invention.

[0057] Figure 9 A point cloud block diagram illustrating a method for extracting three-dimensional deformation data of mining subsidence in desert areas, provided in an embodiment of the present invention;

[0058] Figure 10 A schematic diagram of contour feature points for a method for extracting three-dimensional deformation data of mining subsidence in desert areas, provided in an embodiment of the present invention.

[0059] Figure 11 A schematic diagram of the point cloud contour set for a method of extracting three-dimensional deformation data of mining subsidence in desert areas provided by an embodiment of the present invention;

[0060] Figure 12 This is a schematic diagram of three-dimensional surface deformation for a method of extracting three-dimensional deformation data of mining subsidence in desert areas, provided by an embodiment of the present invention; wherein (a1) represents early ground subsidence; (a2) represents early working face dip displacement; (a3) ​​represents early working face strike displacement; (b1) represents mid-term ground subsidence; (b2) represents mid-term working face dip displacement; (b3) represents mid-term working face strike displacement; (c1) represents late-term ground subsidence; (c2) represents late-term working face dip displacement; (c3) represents late-term working face strike displacement.

[0061] Figure 13 The diagram shows the contour lines of the probability integral method and the research method of the present invention for extracting three-dimensional deformation data of mining subsidence in desert areas, provided by an embodiment of the present invention; where (a) represents displacement in the Y direction; (b) represents displacement in the X direction; and (c) represents displacement in the Z direction.

[0062] Figure 14 This is a schematic diagram illustrating the extraction of three-dimensional deformation data of mining subsidence in desert areas, provided by an embodiment of the present invention, where a point is taken every 5m of displacement in the Y direction for numerical extraction.

[0063] Figure 15 This is a schematic diagram illustrating the extraction of three-dimensional deformation data of mining subsidence in desert areas, provided by an embodiment of the present invention, where a point is taken every 5m of displacement in the X direction for numerical extraction. Detailed Implementation

[0064] To make the above-mentioned objects, features, and advantages of the present invention more apparent and understandable, specific embodiments of the present invention will be described in detail below with reference to the accompanying drawings. Many specific details are set forth in the following description to provide a thorough understanding of the present invention. However, the present invention can be practiced in many other ways different from those described herein, and those skilled in the art can make similar modifications without departing from the spirit of the present invention. Therefore, the present invention is not limited to the specific embodiments disclosed below.

[0065] See Figure 1 This invention provides a method for extracting three-dimensional deformation data of mining subsidence in desert areas. Observation of the movement and deformation of mining subsidence areas has always been a hot research topic. However, conventional observation can only obtain asynchronous information in the horizontal and downward directions. Monitoring the real-time three-dimensional dynamic changes of the subsidence area is particularly important. Based on the novel observation method of UAV-LiDAR, this invention proposes a method for extracting three-dimensional deformation of mining subsidence using a point cloud contour registration method. The main steps include: (1) preprocessing the LiDAR point cloud into blocks to obtain ground point clouds of the same area; (2) using improved extraction contour feature points. The algorithm obtains the point cloud contour; (3) The point cloud contour is used as the registration object, and the three-dimensional deformation is extracted from the transformation matrix; This invention is applied to the 150202 working face of Zaoquan Coal Mine as the research area, and the results are as follows: the maximum value of ΔY deformation is 1.22 and -1.36m, the maximum value of ΔX direction movement is 1.24 and -0.87m, and the maximum value of ΔZ deformation is -5.10m; This area is semi-desert and is a nature reserve. Traditional monitoring methods cannot be carried out. This invention uses the probability integral method to dynamically predict and analyze the three-dimensional deformation of the mining working face. The comparative analysis initially proves the rationality and feasibility of the research method proposed in this invention.

[0066] Specifically:

[0067] I. Study Area.

[0068] 1. Study area.

[0069] The study area is located at the 150202 working face of Zaoquan Coal Mine, situated in the Mu Us Desert. The working face runs north-south, and the terrain is characterized by low mountains and hills, mostly covered by sand dunes, with surface vegetation consisting of a small amount of shrubs and herbaceous plants. The working face is 1552 meters long, 212 meters wide, with an average mining depth of 170 meters, an average dip angle of 29°, and an average coal seam thickness of 8.25 meters. Longwall mining is employed, with caving management of the roof. Mining commenced on October 14, 2023, and is scheduled to cease on November 16, 2024.

[0070] 2. Monitoring methods.

[0071] This invention uses a DJI M600 drone equipped with a SmartBeak ARS-450i laser measurement system for data acquisition; the design forward and lateral overlap rates are both 30%, the drone's flight altitude is 75m, and the point cloud density is 9.52 points / m². 2 The monitoring range was extended 200m beyond the working face to cover the entire subsidence area. Three phases of data collection were conducted before and during the mining of the 150202 working face, on May 6, 2023, April 25, 2024, and August 1, 2024, respectively. The route planning and on-site operation details are as follows: Figure 2 As shown, part of the original point cloud is as follows: Figure 3 As shown.

[0072] 3. Data preprocessing.

[0073] After acquiring the raw data, the point cloud data from multiple flights acquired in each period were first registered and stitched together to obtain complete monitoring data. Then, the point cloud and the underground working face were placed in the same coordinate system through coordinate system transformation. This invention requires stable ground points that do not produce random fluctuations for analysis. Therefore, the point cloud data from each period were filtered using StatisticalOutlier Removal [1, 2] to remove noise points and outliers. The Improved Progressive TIN Densification (IPTD) algorithm was used to classify the ground points and non-ground points as follows: Figure 4 As shown.

[0074] II. A 3D Deformation Extraction Method Based on Point Cloud Contour Registration.

[0075] 1. Basic ideas and principles.

[0076] The filtered ground points are divided into blocks, and then the feature points of the surrounding point cloud are used as the original point cloud and the target point cloud for registration to obtain the transformation matrix. The deformation variables extracted from the transformation matrix are the 3D movement deformation variables of the ground features. The method flow is as follows: Figure 5 As shown. The specific steps include:

[0077] ① The point cloud data is cropped to ensure sufficient overlap. A point cloud segmentation algorithm is used to divide the ground point clouds of the two phases into blocks to obtain multiple pairs of point sets.

[0078] ②Use point cloud contour feature point extraction methods to obtain multiple pairs of point cloud contour sets.

[0079] ③ Take the second phase point cloud contour set one by one as the point cloud to be registered, and use the first phase point cloud contour set as the original point cloud to perform registration using the K-4PCS combined with VGICP registration algorithm to obtain the registered point cloud and transformation matrix.

[0080] ④ Extract the three-dimensional shape of the land surface using the transformation parameters in the transformation matrix.

[0081] 2. Point cloud segmentation and contour extraction based on KD-tree.

[0082] In point cloud segmentation and contour extraction, a point cloud KD-tree is first constructed. Subsequent registration requires range search and nearest neighbor search using the KD-tree data structure to improve computational efficiency. The steps for building a KD-tree point cloud are: 1) Construct a cube bounding box containing all points based on the global coordinate system of the point cloud; 2) For cubes containing more than one point, construct segmentation planes; 3) The points on the segmented subspace and segmentation planes form branches and connection points; 4) Segment the subspace. If the number of internal points exceeds one, continue segmenting by executing step 2. Figure 6 The image shown is a partial point cloud, such as... Figure 7 The image shows the KD-tree point cloud that was created. No original data was deleted during the creation of the KD-tree point cloud, and the accuracy and number of points remain intact.

[0083] The above-described KD-tree point cloud is segmented and its contour feature points are extracted. The specific steps include: 1) Due to the inconsistent acquisition ranges of the three data acquisition periods, the overlapping region of the point clouds from the two acquisition times is first obtained; 2) The grid size is set, and the overlapping region is segmented as follows: Figure 8 and Figure 9 As shown; 3) Match overlapping point pairs and set them as the original point cloud and the target point cloud; 4) Calculate the normal vector and direction of each point in the point set.

[0084] The principles for identifying contour feature points include: ① The angle between the normal vectors of adjacent points is greater than 35°; ② Using the point and its normal as a tangent plane, project all points with a radius of 0.5m near the point onto the tangent plane, and take the point with the largest two-dimensional angle.

[0085] The specific algorithms include:

[0086] 1) Input point clouds of Point1 and Point2 and construct a kd-tree;

[0087] 2) Calculate the common part using the maximum and minimum values ​​of x and y in the point cloud. ;

[0088] 3) Set grid side lengths and create a grid index. ;

[0089] 4) Obtain point cloud pairs of the same region using the coordinate range of the index. ,right middle Let the search radius be r, and denote the neighborhood points within the search radius r as a set. ,Right now ;

[0090] 5) Define the surface equation , ,Pick The corresponding set ,calculate Distance to the surface Solve The eigenvector corresponding to the minimum value is the normal vector n of that point;

[0091] 6) Principles for determining boundary points: ① Angle between the normal vectors of adjacent points Get the points that meet the requirements. ;②According to Construct a tangent plane at this point along with its normal vector n. , will set Projecting points within the plane onto the tangent plane Above, recorded as ,exist Select a point ,by Let u be the u-axis and n be the n-axis. axis, With the v-axis, Construct a local coordinate system for the coordinate center, denoted as . Calculate the sets respectively Other points Time vector The clockwise angle with the coordinate axis u The set of included angles is obtained by subtracting each pair of adjacent included angles. ,in , will set Arrange the elements in descending order and find the largest included angle. ,when If the value is greater than the threshold, the point is determined to be an edge point.

[0092] The extraction results of contour feature points are as follows Figure 10 As shown, red dots represent contour feature points, and green dots represent the original point cloud; the obtained point cloud contour set is as follows: Figure 11 As shown.

[0093] 3. A three-dimensional deformation extraction method based on K-4PCS combined with VGICP registration algorithm.

[0094] K-4PCS is a coarse registration algorithm based on the principle of four-point consistency of key points. It extracts 3D-DoG key points from the point cloud as a point set and uses radiation invariance to perform coarse registration to obtain the initial transformation matrix. The Voxelized Generalized Iterative Closest Point (VGICP) algorithm is an algorithm for accurately registering point clouds with small movements, but if the object position moves significantly, the registration effect will be poor. Therefore, this invention uses the K-4PCS registration algorithm to perform coarse registration to obtain the initial transformation matrix. Under the initial transformation matrix, the VGICP algorithm is used to obtain the registered point cloud and the superimposed transformation matrix.

[0095] The specific algorithms include:

[0096] 1) The point set P to be registered (i.e., the target point cloud) and the reference point set Q (i.e., the point cloud to be registered) are sampled using a voxelized mesh-based filter, and then the key point set is extracted using a Difference of Gaussian (DoG) key point detector. , ; indicates as:

[0097] ;

[0098] in: Indicates a DoG response. The fuzziness level is indicated by x, y, and z, which represent the three-dimensional coordinates.

[0099] 2) In the key point set A set of 4 coplanar points randomly selected from the data. Searching for the corresponding set using affine invariance. Obtain the initial transformation matrix and the transformed target point cloud. .

[0100] 3) Original point cloud Target point cloud Transformation matrix Each sampling point comes from a Gaussian distribution: , , Distance between it and its neighboring points: , The distribution is represented as:

[0101] ;

[0102] ;

[0103] ;

[0104] The maximum likelihood estimation method was used to... Estimating, we have:

[0105] ;

[0106] ;

[0107] ;

[0108] in: This represents the number of adjacent points.

[0109] The transformation matrix can be calculated from the above formula, and the deformation variables can be obtained by superimposing the initial transformation matrix on it.

[0110] 9) The transformation matrix includes a rotation matrix of R×3 for the point cloud coordinate system and a translation vector of 3×1.

[0111] .

[0112] The three-dimensional deformation variables ΔX, ΔY, and ΔZ can be obtained from the above transformation matrix.

[0113] III. Data processing and analysis of the study area.

[0114] Based on the above method, contour feature points were extracted from the three phases of data. After multiple attempts at grid division, a grid size of 30m×30m yielded the best results. The preliminary investigation suggests that this is related to the contour feature points being located at the boundaries of point cloud blocks and the magnitude of the 3D deformation. The transformation matrix was obtained by matching contour features to acquire the 3D deformation along the working face direction, dip direction, and subsidence direction for each contour position. Kriging interpolation [0, 1] was used to interpolate the deformation of discrete points. The weight function for discrete points was based on the reciprocal of the first power of the distance to the interpolation point; the greater the distance, the smaller the weight. This interpolation method is suitable for dense point clouds and can reduce the impact of noise on interpolation. The interpolation point value is equal to the weighted average of its nearby discrete points. The interpolated 3D deformation was used to draw a surface movement deformation cloud map, and the results are as follows: Figure 12 As shown.

[0115] from Figure 12 It can be seen from this:

[0116] (1) The ground subsidence extraction results of each period showed a good correspondence with the mining situation of the working face. The entire subsidence basin showed a downward shift in the direction of the mountain. The maximum subsidence value from 202305 to 202404 was -4.20m, the maximum subsidence value from 202305 to 202408 was -5.06m, and the maximum subsidence value from 202304 to 202408 was -5.10m.

[0117] (2) In terms of the dip direction, the planar displacement along the working face in each period also shows a good correspondence with the working face. Points in the downhill direction move towards the uphill direction, and points in the uphill direction move towards the downhill direction. The symmetrical position of the movement is located in the center of the working face, slightly towards the downhill direction, and maintains a good fit with the subsidence basin. Among them, the maximum dip displacement from 202305 to 202404 is -1.15m, the maximum dip displacement from 202305 to 202408 is -1.48m, and the maximum dip displacement from 202304 to 202408 is -1.36m. The ratio of the maximum displacement to the maximum subsidence is about 0.35, 0.29, and 0.27, which is in line with the general law of mining subsidence.

[0118] (3) In terms of strike direction, the planar displacement of the working face advance position points towards the center of the working face, which is consistent with the mining subsidence pattern, but no obvious planar displacement is observed at the opening position. The field survey found that the opening position was too shallow, with a minimum depth of only 27m. During the mining process, a subsidence pit appeared at this position, which was subsequently backfilled artificially, causing the surface subsidence pattern in this area to be artificially disrupted. Therefore, when conducting strike data analysis, only the dynamic monitoring results between 202404 and 202408 are considered, and the data of the opening position of the working face are no longer considered. In the surface strike planar displacement between 202404 and 202408, the maximum value is 1.24m, and the ratio of the maximum displacement to the maximum subsidence is about 0.24, which is also consistent with the general pattern of mining subsidence.

[0119] To further analyze the reliability of the extracted data, a comparative analysis was conducted using theoretical calculations, taking the data from 202404 to 202408 as an example.

[0120] IV. Theoretical analysis of the reliability of data results.

[0121] The theory of random media evolved into the probability integral method, which posits that the movement of rock strata and the earth's surface caused by mining is macroscopically similar to that of a granular medium model, which is considered a random medium. Based on the principle of the probability integral method, the settlement value at any point on the ground caused by mining can be expressed as:

[0122]

[0123] in: .

[0124] in: , denoted as the subsidence distribution coefficients of the projection points of the point to be determined on the strike and dip main cross sections, respectively; , These represent the calculated lengths of the working face, indicating its direction and dip, respectively. A mining unit representing a probability integral; Indicates the thickness of the coal seam; Indicates the subsidence coefficient; Indicates the dip angle of the coal seam; Indicates the radius of main influence. ,in Indicates the average mining depth. This indicates the main influence on the tangent angle.

[0125] Subsidence dynamics prediction uses parameters such as working face coordinates, mining depth, and mining thickness to comprehensively predict and display mining subsidence, expressed by the formula:

[0126]

[0127] in: Indicates the area of ​​the region; This represents the maximum settlement value under geological and mining conditions; This indicates the sinking time t during mining.

[0128] Because the region is desertified, and the superposition effect of the three-dimensional deformation obtained using the probability integral method and the method used in this study is as follows: Figure 13 As shown, values ​​are extracted from a point every 5m along the horizontal displacement trend line and the dip line. Figure 14 and Figure 15 As shown.

[0129] The dynamic prediction results using the probability integral method show that the maximum displacement in the Y direction is 1.02m and -0.92m, the maximum displacement in the X direction is 1.1m and -0.90m, and the maximum subsidence is -5.10m. Analysis of the contour maps and data comparisons shows that the measured results are basically consistent with the theoretical calculations, preliminarily proving the rationality and feasibility of the research method proposed in this invention. On the other hand, the main reasons for the discrepancies between the two methods include:

[0130] (1) The terrain of the area is high in the east and west and low in the middle, which has a certain impact on the measured data. In addition, there are also human factors: the area is located in a desertification nature reserve and is regularly managed to prevent water loss; large-scale collapses occurred in the cut-in eyes and were filled with soil.

[0131] (2) Under different geological conditions and mining methods, the selection of dynamic prediction parameters needs to be adjusted and selected according to the actual situation, which has a certain impact on the results.

[0132] (3) The selection of the segmentation boundary and the search radius of similar feature points in the neighborhood of the contour feature point registration algorithm studied in this invention has a certain impact on the experimental results.

[0133] (4) In addition to natural weather factors, the positioning and attitude determination of the UAV in motion also have certain errors in the original point cloud data, which have a certain impact.

[0134] For special landforms such as deserts, mountains, and densely vegetated areas, as well as sensitive areas such as desertification and nature reserves, this invention proposes a point cloud contour registration method to monitor three-dimensional deformation in mining areas, obtain the X, Y, and Z three-dimensional deformation of the surface, provide more diverse data for the analysis of surface subsidence patterns, and better guide related mining activities or ecological environment restoration.

[0135] Using the method of this invention, the three-dimensional deformation data of the surface caused by mining in the study area were obtained. The maximum values ​​of ΔY deformation were 1.22 and -1.36 m, the maximum values ​​of ΔX displacement were 1.24 and -0.87 m, and the maximum value of ΔZ deformation was -5.10 m. The analysis showed that artificial backfilling and shallow coal seams resulted in no obvious movement and deformation patterns near the cuts of 202305-202404 and 202305-202408.

[0136] The research area mentioned is semi-desert and a nature reserve, where traditional monitoring methods cannot be carried out. The probability integral method was used to dynamically predict the three-dimensional deformation caused by mining activities from 202404 to 202408. The results obtained by comparing with the research method of this invention prove its rationality and feasibility.

[0137] The embodiments described above are merely illustrative of several implementations of the present invention, and while the descriptions are relatively specific and detailed, they should not be construed as limiting the scope of the invention patent. It should be noted that those skilled in the art can make various modifications and improvements without departing from the concept of the present invention, and these all fall within the protection scope of the present invention. Therefore, the protection scope of this invention patent should be determined by the appended claims.

Claims

1. A method for extracting three-dimensional deformation data of mining subsidence in a desert area, characterized in that, The method comprises the following steps: Obtaining three-phase ground point cloud data corresponding to different shooting time points of a desert area, and extracting an overlapping area in the three-phase ground point cloud data; wherein the ground point cloud data corresponding to the overlapping area is composed of ground points and non-ground points; Obtaining the normal vector angle between each point in the point cloud data of the ground points in each phase of the ground point cloud data in the overlapping area and the corresponding points in the other two phases, and the angle of the projection of each point on the tangent plane of the normal vector of the respective point; Setting the points with the normal vector angle higher than a preset value as contour feature points, and setting the point with the largest angle of the projection of each point on the tangent plane of the normal vector of the respective point as a contour feature point among the three points formed by the point with the normal vector angle higher than the preset value and the corresponding points in the other two phases; wherein the point with the normal vector angle higher than the preset value and the point with the largest angle of the projection of each point on the tangent plane of the normal vector of the respective point are a pair of point cloud contours, and a plurality of pairs of point cloud contour sets are obtained; Taking one-phase point cloud data as original point cloud, for each point cloud data in the original point cloud, the point cloud data belonging to the same pair of point cloud contours in the other two phases of point cloud data as the original point cloud is registered, and three-dimensional deformation data representing the displacement of the point cloud in the registration process is obtained.

2. The method according to claim 1, wherein, The division of the ground points and the non-ground points comprises: After the extraction of the ground point cloud data in the overlapping area in the three-phase ground point cloud data, the ground point cloud data in the overlapping area is divided into ground points and non-ground points by using an incremental triangulated irregular network filtering algorithm IPTD.

3. The method according to claim 2, wherein, After the ground point cloud data corresponding to the ground points in the overlapping area is divided into ground points and non-ground points, the ground point cloud data of the ground points corresponding to the overlapping area is constructed into a KD-tree tree structure, comprising: For the point cloud data of the ground points corresponding to the three-phase overlapping area, a cube bounding box containing all point clouds is established according to a point cloud global coordinate system, and a segmentation plane is constructed for the cube containing more than one point; The points on the segmentation sub-space and the segmentation plane form branches and connection points to construct a KD-tree data structure.

4. The method according to claim 3, wherein, The obtaining of the normal vector angle between each point in the point cloud data of the ground points in each phase of the ground point cloud data in the overlapping area and the corresponding points in the other two phases, and the angle of the projection of each point on the tangent plane of the normal vector of the respective point comprises: In the kd-tree tree structure, a search radius r is set, and the neighborhood points within the search radius r are recorded as a set , that is ; Setting the surface equation , , take the corresponding set , calculate the distance to the surface , solve the minimum corresponding eigenvector as the normal vector n of the corresponding point; The angle between the normal vectors of adjacent points is expressed as ; According to and its normal vector n make the tangent plane of the corresponding point , and the normal vector n is projected to the tangent plane The tangent plane projection of the normal vector is formed.

5. The method according to claim 4, wherein, The formation of the plurality of pairs of point cloud contour sets comprises: The points with the normal vector angle between adjacent points greater than 35° are set as contour feature points; according to Construct a tangent plane at this point along with its normal vector n. , will set Projecting points within the plane onto the tangent plane Above, recorded as ,exist Select a point ,by Let u be the u-axis and n be the n-axis. axis, With the v-axis, Construct a local coordinate system for the coordinate center, denoted as Calculate the sets respectively Other points Time vector The clockwise angle with the coordinate axis u The set of included angles is obtained by subtracting each pair of adjacent included angles. ,in , will set Arrange the elements in descending order and find the largest included angle. , the largest included angle Set as contour feature points; The plurality of pairs of point cloud contour sets are formed according to the two types of contour feature points.

6. The method for extracting 3D deformation data of mining subsidence in a desert region according to claim 1, characterized in that, The obtaining of the three-dimensional deformation data comprises: The target point cloud is set as a point set P to be registered, and the point cloud to be registered is set as a point set Q to be registered; For the point set P to be registered and the reference point set Q, a voxelized grid based filter is used for sampling, and then a Difference of Gaussians (DoG) keypoint detector is used to extract the keypoint set , ; is represented as: ； wherein: represents the DoG response, represents the blur level, x, y, z represent three-dimensional coordinates; In the key point set A set of 4 coplanar points randomly selected from the data. Searching for the corresponding set using affine invariance. Obtain the initial transformation matrix; Based on the initial transformation matrix, the similarity between the point clouds is calculated by the VGICP algorithm voxelization and generalized distance measurement strategy, the point clouds are gradually aligned, the transformation matrix in the alignment process is obtained, the transformation matrix in the alignment process and the initial transformation matrix are superimposed, and the final transformation matrix is obtained, which is represented as: ； The final transformation matrix includes a rotation matrix of a point cloud coordinate system R×3 and a translation vector 3×1; wherein the rotation matrix R represents the direction change of the point cloud in the three-dimensional space, and the translation vector represents the moving distance of the point cloud in the X, Y and Z directions; According to the moving distance of the point cloud in the X, Y and Z directions, three-dimensional deformation data of the mining subsidence in the desert area is formed.

7. A device for extracting three-dimensional deformation data of mining subsidence in desert areas, characterized in that, The method comprises the following steps: a data acquisition module is configured to acquire three-stage ground point cloud data corresponding to different shooting time points in the desert area, and extract overlapping areas in the three-stage ground point cloud data; wherein the ground point cloud data corresponding to the overlapping areas is composed of ground points and non-ground points; a contour feature point determination module is configured to acquire the normal vector angle between each point in the point cloud data of the ground points in each stage of ground point cloud data and the corresponding points in the other two stages, and the angle of the projection of each point and the tangent plane of the respective point normal vector; points with a normal vector angle higher than a preset value are set as contour feature points; and among the three points formed by the points with a normal vector angle higher than a preset value and the corresponding points in the other two stages, the point with the largest angle of the projection of each point and the tangent plane of the respective point normal vector is set as a contour feature point; wherein the points with a normal vector angle higher than a preset value and the points with the largest angle of the projection of each point and the tangent plane of the respective point normal vector are a pair of point cloud contours, and a plurality of pairs of point cloud contours are acquired; a registration module is configured to take one stage of point cloud data as original point cloud, for each point cloud data in the original point cloud, register the point cloud data in the other two stages of point cloud data which belongs to a pair of point cloud contours with each point cloud data in the original point cloud, to obtain three-dimensional deformation data representing point cloud displacement in the registration process.

8. An electronic device, comprising: The method comprises the following steps: a memory and a processor; the memory is configured to store a computer program; the processor is configured to execute the computer program stored in the memory, to realize the steps of the method for extracting three-dimensional deformation data of mining subsidence in a desert area according to any one of claims 1-6.

9. A computer-readable storage medium, characterized in that, a computer program is stored, and the computer program is executed by a processor to realize the steps of the method for extracting three-dimensional deformation data of mining subsidence in a desert area according to any one of claims 1-6.

Images