Desert area mining subsidence three-dimensional deformation data extraction method and device

Abstract

The invention discloses a desert area mining subsidence three-dimensional deformation data extraction method and a desert area mining subsidence three-dimensional deformation data extraction device, and relates to the technical field of earth observation. Normal vector included angle calculation and tangent plane projection analysis are performed on each point in point cloud data to form a plurality of pairs of point cloud contour sets, and the whole point cloud contour sets are based on local geometric characteristics of point clouds; feature points reflecting the contour shape of the point cloud are accurately extracted; then, a plurality of pairs of point cloud contour sets are used as registration objects, K-4PCS and VGICP registration algorithms are adopted to carry out registration on the three-stage point cloud data so as to accurately align the point cloud data, point clouds and a transformation matrix after registration are obtained, and earth surface three-dimensional deformation quantity is extracted from transformation parameters in the transformation matrix; in the process, sinking extraction based on a point DEM and plane displacement extraction based on rigid body matching do not need to be carried out, and desert area plane displacement extraction can be accurately carried out, so that earth surface three-dimensional deformation data caused by desert area mining can be accurately obtained.

Description

Technical Field

[0001] The present invention relates to the technical field of earth observation, in particular to a method, device, equipment and medium for extracting three-dimensional deformation data of mining subsidence in desert areas. Background Technique

[0002] As an important energy source and raw material for the development of China's economic industry, coal resources are the key elements to ensure the vitality of the national economic development, the "ballast stone" for the stable supply of national energy security, and an important starting point for implementing the "dual carbon" goal; while the large-scale mining of coal resources meets the energy demand, it also causes geological environment problems such as surface subsidence, crop yield reduction, soil erosion, and vegetation degradation, which have caused great harm and troubles to the production and life of the people in the mining area. Therefore, the secondary disaster problems caused by coal mining subsidence have become an important factor affecting the stability and economic development of the mining area, and 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 at home and abroad have been innovating in the monitoring means of surface subsidence characteristics, and fruitful results have been achieved in the research of subsidence mechanism. Traditional subsidence monitoring in mining areas mainly relies on technical means such as level instruments, total station instruments, and Real - time kinematic (RTK). The above - mentioned methods require a large amount of manpower and material resources and are too difficult to implement, so they have gradually been replaced by emerging measurement techniques. Some researchers have proposed the SGI - SF method for extracting three - dimensional displacement results from InSAR datasets. The results show that the average root mean square error (RMSE) of three - dimensional deformation in the vertical and horizontal directions of this method is about 9.28 mm and 13.10 mm respectively, which is much smaller than the displacement caused by mining. Some researchers have used SBAS - InSAR to invert the advance influence angle of mining subsidence, providing important reference significance for the dynamic protection of surface buildings and other facilities in western mining areas. Some researchers have proposed a monitoring method combining InSAR and the surface displacement vector attenuation angle model to obtain more detailed and accurate deformation information of the entire basin. Some researchers have combined D - InSAR and UAV technology to monitor the deformation of mining areas, obtaining surface movement and deformation, and the error of the fused data is within 10 cm. The research conducted by the above - described researchers all focuses on InSAR monitoring. InSAR monitoring has significant advantages in obtaining small - magnitude deformations, but it often has many limitations due to the decorrelation phenomenon when facing large deformations of the m - level caused by coal mining. Therefore, airborne photogrammetry and lidar technology have also been introduced into coal - mining subsidence monitoring. Some researchers have obtained ground point cloud data through UAV - LiDAR and used the Kriging interpolation method to differentially process the DEM to obtain ground subsidence data, achieving a ground DEM accuracy of 15 mm and a root mean square error of the ground subsidence basin of 39 mm. However, UAV - LiDAR technology generates a large amount of point cloud data, and the processing and analysis of these data require high - performance computing resources. Moreover, the accuracy of the Kriging interpolation method highly depends on the distribution and quantity of data points. If the data points are unevenly distributed or insufficient in quantity, it will lead to a decrease in the accuracy of the interpolation result. At the same time, the Kriging interpolation method needs to calculate the semi - variance function model and perform interpolation calculations based on this function, but the selection of the semi - variance function model often requires subjective judgment according to the characteristics and distribution of the data, which increases the subjectivity and uncertainty of model selection.

[0004] In view of the problems existing in UAV-LiDAR technology and Kriging interpolation method, various point cloud physical extraction algorithms have emerged one after another, which promotes the application of laser scanning data in mining monitoring. Some researchers have proposed a clustering segmentation algorithm that optimizes the concavity and convexity of supervoxels for the segmentation of ground features with complex terrain characteristics in mining areas, improving the accuracy of ground feature segmentation; some researchers have also proposed a rigid body registration algorithm based on connected domain segmentation to synchronously obtain the three-dimensional movement of the ground surface, solving the problems caused by insufficient rigid bodies. Although such methods have successfully solved the problems of subsidence extraction based on point DEM and planar displacement extraction based on rigid body matching, as the focus of mining activities has shifted westward, a large number of mining activities have entered sparsely populated desert areas. There are few artificial buildings and scarce or even missing rigid bodies in the area. Therefore, the current methods are difficult to extract planar displacement in desert areas to obtain three-dimensional surface deformation data caused by mining in desert areas. Summary of the Invention

[0005] An embodiment of the present invention provides a method and device for extracting three-dimensional deformation data of mining subsidence in desert areas, which can solve the problem that in the prior art, the current method is difficult to extract planar displacement in desert areas to obtain three-dimensional surface deformation data caused by mining in desert areas.

[0006] An embodiment of the present invention provides a method for extracting three-dimensional deformation data of mining subsidence in desert areas, including the following steps: Obtain three-phase surface point cloud data corresponding to different shooting time points in the desert area, and extract the overlapping areas in the three-phase surface point cloud data; wherein, the surface point cloud data corresponding to the overlapping area is composed of ground points and non-ground points; Obtain the normal vector angle between each point in the point cloud data of ground points in each phase of surface point cloud data in the overlapping area and the corresponding points in the other two phases, and the angle between each point and the tangent plane projection of its own point normal vector; Set the points with normal vector angles higher than the preset value as contour feature points; and among the three points formed by the points with normal vector angles higher than the preset value and the corresponding points in the other two phases, set the point with the largest angle between the tangent plane projection of its own point normal vector as the contour feature point; wherein, the points with normal vector angles higher than the preset value and the points with the largest angle between the tangent plane projection of their own point normal vectors are used as a pair of point cloud contours, and multiple pairs of point cloud contour sets are obtained; Taking one-phase point cloud data as the original point cloud, for each point cloud data in the original point cloud, register the point cloud data in the other two phases that belongs to the same pair of point cloud contours as each point cloud data in the original point cloud, and obtain three-dimensional deformation data representing the displacement of the point cloud during the registration process.

[0007] Preferably, the division of the ground points and non-ground points includes: After extracting the ground point cloud data in the overlapping area of the three-phase surface point cloud data, the progressive encrypted triangular mesh filtering algorithm IPTD is used to divide the ground point cloud data in the overlapping area into ground points and non-ground points.

[0008] Preferably, after dividing the ground point cloud data corresponding to the overlapping area into ground points and non-ground points, the ground point cloud data of the overlapping area is constructed into a KD-tree tree structure, including: For the point cloud data of the ground points corresponding to the three-phase overlapping area, a cubic bounding box containing all the point clouds is established according to the global coordinate system of the point cloud. For the cube containing more than 1 point, a splitting plane is constructed; The divided subspaces and the points on the splitting plane form branches and connection points to construct the KD-tree data structure.

[0009] Preferably, obtaining the normal vector angle between each point in the point cloud data of the ground points in each phase of the overlapping area and the corresponding points in the other two phases, and the angle between each point and the tangent plane projection of its own point normal vector, includes: In the kd-tree tree structure, set the search radius r, and record the neighborhood points within the search radius r as the set , that is ; Set the surface equation , , take The corresponding set , calculate The distance to the surface , solve The eigenvector corresponding to the minimum is the normal vector n of the corresponding point; The adjacent point normal vector angle is expressed as ; According to And its normal vector n, make the tangent plane of the corresponding point , and project the normal vector n onto the tangent plane Form the tangent plane projection of the normal vector. In the kd-tree tree structure, set the search radius r, and record the neighborhood points within the search radius r as the set , that is ; Set the surface equation , , take The corresponding set , calculate The distance to the surface , solve The eigenvector corresponding to the minimum is the normal vector n of the corresponding point; The adjacent point normal vector angle is expressed as ; according to Make a tangent plane with the corresponding point and its normal vector n , and project the normal vector n onto the tangent plane Form the tangent plane projection of the normal vector.

[0010] Preferably, forming multiple pairs of point cloud contour sets includes: Set the points where the angle between the normal vectors of adjacent points is greater than 35° as contour feature points; according to Make a tangent plane to the point with its normal vector n , will gather Projecting the points inside onto the tangent plane , recorded as ,exist Pick one ,by is the u axis, n is axis, is the v-axis, Construct a local coordinate system for the coordinate center, expressed as , respectively calculate the set Other points Arrive Vector The clockwise angle with the coordinate axis u , take the difference between two adjacent angles to get the angle set ,in , will gather Arrange the elements in descending order and find the largest angle , the largest angle Set as contour feature point; Multiple pairs of point cloud contour sets are formed based on two types of contour feature points.

[0011] Preferably, the acquisition of the three-dimensional deformation data includes: Set the target point cloud as the point set to be registered P, and the point cloud to be registered as the reference point set Q; For the point set P to be registered and the reference point set Q, a voxel-based grid filter is used for sampling, and then the key point set is extracted by the Gaussian difference DoG key point detector. 、 ; expressed as: ; in: Indicates DoG response, Indicates the blur level, x, y, z indicate the three-dimensional coordinates; In the key point set A set of four coplanar points is randomly selected from , and the corresponding set is searched using affine invariance to obtain an initial transformation matrix; Based on the initial transformation matrix, the similarity between point clouds is calculated through voxelization and generalized distance metric strategy of the VGICP algorithm, 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: ; The final transformation matrix includes a rotation matrix of R×3 and a translation vector of 3×1 in the point cloud coordinate system; among them, the rotation matrix R represents the direction change of the point cloud in three-dimensional space, and the translation vector represents the moving distances of the point cloud in the X, Y, and Z directions; According to the moving distances of the point cloud in the X, Y, and Z directions, three-dimensional deformation data of mining subsidence in desert areas is formed.

[0012] An embodiment of the present invention also provides a device for extracting three-dimensional deformation data of mining subsidence in desert areas, including: A data acquisition module for acquiring three-phase surface point cloud data corresponding to different shooting time points in the desert area and extracting the overlapping area in the three-phase surface point cloud data; among them, the surface point cloud data corresponding to the overlapping area is composed of ground points and non-ground points; A contour feature point determination module for obtaining the normal vector angle between each point in the point cloud data of the ground points in each phase of the surface point cloud data in the overlapping area and the corresponding points in the other two phases, and the angle between each point and the tangent plane projection of its own 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 the preset value and the corresponding points in the other two phases, the point with the largest angle between the tangent plane projection of its own point normal vector is set as a contour feature point; among them, the points with a normal vector angle higher than the preset value and the points with the largest angle between the tangent plane projection of their own point normal vector are used as a pair of point cloud contours, and multiple sets of point cloud contours are obtained; A registration module for using the point cloud data of one phase as the original point cloud, and for each point cloud data in the original point cloud, registering 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 to obtain three-dimensional deformation data representing the displacement of the point cloud during the registration process.

[0013] An embodiment of the present invention also provides an electronic device, including a memory and a processor; The memory is used to store a computer program; 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 a desert area as described above.

[0014] An embodiment of the present invention also provides a computer-readable storage medium for storing a computer program. When the computer program is executed by a processor, it implements the steps of the method for extracting three-dimensional deformation data of mining subsidence in a desert area as described above.

[0015] An embodiment of the present invention provides a method and device for extracting three-dimensional deformation data of mining subsidence in a desert area. Compared with the prior art, the beneficial effects are as follows: In the present invention, by collecting the surface point cloud data of three different time points in the research area, extracting the point cloud data of the ground points in the overlapping area of the three phases, respectively constructing the KD-tree tree structure for the point cloud data of the ground points in the three phases, calculating the normal vector angle and performing tangent plane projection analysis on each point in the point cloud data to form multiple pairs of point cloud contour sets; then, using the K-4PCS and VGICP registration algorithms to register the point cloud data of the three phases with the multiple pairs of point cloud contour sets as the registration objects, obtaining the registered point cloud and the transformation matrix, and extracting the three-dimensional deformation amount of the ground surface from the transformation parameters in the transformation matrix; in this process, when calculating the normal vector angle and performing tangent plane projection analysis on each point in the point cloud data, it is based on the local geometric characteristics of the point cloud to accurately extract the feature points reflecting the shape of the point cloud contour, and registering based on the 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 to perform subsidence extraction based on point DEM and plane displacement extraction based on rigid body matching, and can accurately extract the plane displacement in the desert area to accurately obtain the three-dimensional deformation data of the ground surface caused by mining in the desert area. Description of the Drawings

[0016] Figure 1 It is a schematic diagram of a method for extracting three-dimensional deformation data of mining subsidence in a desert area provided by an embodiment of the present invention; Figure 2 It is a schematic diagram of the UAV flight path during the acquisition process of a method for extracting three-dimensional deformation data of mining subsidence in a desert area provided by an embodiment of the present invention; Figure 3 It is a schematic diagram of the original point cloud data collected on-site of a method for extracting three-dimensional deformation data of mining subsidence in a desert area provided by an embodiment of the present invention; Figure 4Schematic diagram of ground points and non-ground points obtained by classifying point cloud data of each period using an incremental encryption triangular mesh filtering algorithm for a method for extracting three-dimensional deformation data of mining subsidence in desert areas provided by an embodiment of the present invention; (a) therein represents the filtering result; (b) therein represents the ground points; (c) therein represents the non-ground points; Figure 5 Schematic diagram of the registration process for a method for extracting three-dimensional deformation data of mining subsidence in desert areas provided by an embodiment of the present invention; Figure 6 Partial point cloud schematic diagram for a method for extracting three-dimensional deformation data of mining subsidence in desert areas provided by an embodiment of the present invention; Figure 7 KD-tree point cloud schematic diagram for a method for extracting three-dimensional deformation data of mining subsidence in desert areas provided by an embodiment of the present invention; Figure 8 Overall point cloud schematic diagram for a method for extracting three-dimensional deformation data of mining subsidence in desert areas provided by an embodiment of the present invention; Figure 9 Point cloud segmentation schematic diagram for a method for extracting three-dimensional deformation data of mining subsidence in desert areas provided by an embodiment of the present invention; Figure 10 Schematic diagram of contour feature points for a method for extracting three-dimensional deformation data of mining subsidence in desert areas provided by an embodiment of the present invention; Figure 11 Point cloud contour set schematic diagram for a method for extracting three-dimensional deformation data of mining subsidence in desert areas provided by an embodiment of the present invention; Figure 12 Schematic diagram of three-dimensional surface deformation for a method for extracting three-dimensional deformation data of mining subsidence in desert areas provided by an embodiment of the present invention; (a1) therein represents early ground subsidence; (a2) therein represents early working face dip displacement; (a3) therein represents early working face strike displacement; (b1) therein represents middle ground subsidence; (b2) therein represents middle working face dip displacement; (b3) therein represents middle working face strike displacement; (c1) therein represents late ground subsidence; (c2) therein represents late working face dip displacement; (c3) therein represents late working face strike displacement; Figure 13 Isoline schematic diagram of the probability integral method and the research method of the present invention for a method for extracting three-dimensional deformation data of mining subsidence in desert areas provided by an embodiment of the present invention; (a) therein represents the displacement in the Y direction; (b) therein represents the displacement in the X direction; (c) therein represents the displacement in the Z direction; Figure 14 Schematic diagram of extracting numerical values by taking one point every 5 m for the displacement in the Y direction for a method for extracting three-dimensional deformation data of mining subsidence in desert areas provided by an embodiment of the present invention; Figure 15 Schematic diagram of extracting numerical values by taking a point every 5 m in the X - direction displacement for a method of extracting three - dimensional deformation data of mining subsidence in desert areas provided by an embodiment of the present invention. Detailed implementation manners

[0017] In order to make the above - mentioned objects, features and advantages of the present invention more obvious and understandable, the following detailed description of the specific implementation manners of the present invention will be given with reference to the accompanying drawings. Many specific details are set forth in the following description in order to fully understand the present invention. However, the present invention can be implemented in many other ways different from those described herein, and those skilled in the art can make similar improvements without departing from the connotation of the present invention. Therefore, the present invention is not limited by the specific embodiments disclosed below.

[0018] Refer to Figure 1 , an embodiment of the present invention provides a method for extracting three - dimensional deformation data of mining subsidence in desert areas. The observation of moving deformation in the mining subsidence area has always been a hot topic of research. However, conventional observations can only obtain non - synchronous information in the horizontal and subsidence directions. It is particularly important to monitor the real - time three - dimensional dynamic change process of the subsidence area. Based on the new observation means of UAV - LiDAR, the present invention proposes a method for extracting three - dimensional deformation of mining subsidence by using the point - cloud contour registration method; the main steps include: (1) pre - processing the LiDAR point cloud by block to obtain the ground point cloud in the same area; (2) obtaining the point - cloud contour by using an improved algorithm for extracting contour feature points; (3) taking the point - cloud contour as the registration object, registering and extracting three - dimensional deformation from the transformation matrix; the present invention takes the 150202 working face of Zaoquan Coal Mine as the research area for application, and obtains the results: the maximum value of ΔY deformation is 1.22 and - 1.36 m, the maximum value of movement in the ΔX direction is 1.24 and - 0.87 m, and the maximum value of ΔZ deformation is - 5.10 m; this area is semi - desert and a nature reserve, and traditional monitoring means cannot be carried out. The present invention uses the probability integral method to dynamically predict and analyze the three - dimensional deformation of the mining working face, and the comparative analysis preliminarily proves the rationality and feasibility of the research method proposed by the present invention.

[0019] Specifically: I. Research area.

[0020] 1. Research area.

[0021] The research area is located in the 150202 working face of Zaoquan Coal Mine, which is located in the Mu Us Desert. The working face runs from north to south. The landform belongs to the low mountain and hilly area, and is mostly covered by sand dunes. The surface vegetation is a small amount of shrubs and herbaceous plants. The strike length of this working face is 1552 meters, the dip width is 212 meters, the average mining depth is 170m, the average dip angle is 29°, and the average coal seam thickness is 8.25m. The longwall mechanized mining method is adopted, and the roof is managed by the caving method. The mining started on October 14, 2023 and is designed to stop mining on November 16, 2024.

[0022] 2. Monitoring method.

[0023] In this invention, the monitoring uses a DJI M600 drone equipped with a Zhihui ARS-450i laser measurement system for data collection; the designed heading overlap rate and side overlap rate are both 30%, the flight altitude of the drone is 75m, and the point cloud density is 9.52 points / m 2 , and the monitoring range extends 200m outside the working face to cover the entire subsidence area; a total of three phases of data collection are carried out before and during the mining of the 150202 working face, and the data collection times are May 6, 2023, April 25, 2024, and August 1, 2024 respectively. The flight line planning and on-site operation conditions are as Figure 2 shown, and some original point clouds are as Figure 3 shown.

[0024] 3. Data preprocessing.

[0025] After the original data is obtained, first register and splice the multi-pass point cloud data obtained in each phase to obtain complete monitoring data, and then place the point cloud and the underground working face in the same coordinate system through coordinate system conversion; for the analysis of this invention, stable ground points that do not produce random fluctuations are required. Therefore, the StatisticalOutlier Removal filter [1, 2] is used for each phase of point cloud data to remove noise points and outliers, and the Improved Progressive TIN Densification (IPTD) filtering algorithm is used to classify the ground points and non-ground points as Figure 4 shown.

[0026] II. Three-dimensional deformation extraction method based on point cloud contour registration.

[0027] 1. Basic idea principle.

[0028] The filtered ground points are divided into blocks, and then the point cloud contour feature points in the surrounding area are used as the original point cloud and the target point cloud for registration to obtain the transformation matrix. The deformation amount extracted from the transformation matrix is the three-dimensional movement deformation amount of the ground object points. The method flow is as Figure 5 shown. The specific steps include:

[0029] ① Crop the point cloud data to ensure sufficient overlap of the point cloud. Use a point cloud segmentation algorithm to divide the two-phase ground point cloud into blocks to obtain multiple pairs of point sets.

[0030] ② Use the method for extracting point cloud contour feature points to obtain multiple pairs of point cloud contour sets.

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

[0032] ④ Extract the three-dimensional surface deformation quantity using the transformation parameters in the transformation matrix.

[0033] 2. Point cloud block division and contour extraction based on KD-tree.

[0034] When dividing the point cloud and extracting the contour, first construct a point cloud KD-tree. Range search and nearest neighbor search of the KD-tree data structure are required for subsequent registration to improve the operation efficiency. The steps for establishing the KD-tree point cloud are as follows: 1) Establish a cubic bounding box containing all point clouds according to the global coordinate system of the point cloud; 2) For the cube containing more than 1 point, construct a splitting plane; 3) The split subspace and the points on the splitting plane form branches and connection points; 4) Split the subspace. If the number of internal points exceeds 1, execute step 2 to continue splitting; as Figure 6 shown in part of the point cloud, as Figure 7 shown in the established KD-tree point cloud. No deletion of the original data is performed during the establishment of the KD-tree point cloud, and the point cloud accuracy and number remain intact.

[0035] The steps for dividing and extracting the contour feature points of the above-described KD-tree point cloud are as follows: 1) Since the acquisition ranges of the three-phase data are inconsistent, first obtain the overlapping area of the point clouds at the two acquisition time points; 2) Set the grid size and perform splitting processing on the overlapping area as Figure 8 and Figure 9 shown; 3) Match the point pairs in the overlapping area 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.

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

[0037] The specific algorithm includes: 1) Input the Point1 and Point2 point clouds and construct a kd-tree; 2) Calculate the common part using the maximum and minimum values of x and y in the point cloud ; 3) Set the grid side length grid and create a grid index ; 4) Obtain the point cloud pairs in the same area using the coordinate range of the index , for in set the search radius r, and record the neighborhood points within the search radius r as a set , that is ; 5) Set the surface equation , , take corresponding set , calculate the distance to the surface , solve the eigenvector corresponding to the minimum is the normal vector n of this point; 6) Principle for determining boundary points: ① The included angle between the normal vectors of adjacent points , obtain the points that meet the requirements ; ② According to and its normal vector n, make the tangent plane of this point , project the points in the set onto the tangent plane , denoted as , in select a point , with as the u-axis, n as axis, as the v-axis, and as the coordinate center to construct a local coordinate system, denoted as , calculate the vectors of other points in the set to the point respectively and the clockwise included angle with the u-axis of the coordinate axis, make the difference between adjacent included angles pairwise to obtain an included angle set , where , sort the elements in the set in descending order, find the largest included angle , when is greater than the threshold, determine this point as an edge point.

[0038] The extraction results of the contour feature points are as Figure 10 shown, where the red points are the contour feature points and the green points are the original point cloud; A certain point cloud contour set pair obtained is as Figure 11 shown.

[0039] 3. 3D Deformation Extraction Method of K-4PCS Combined with VGICP Registration Algorithm

[0040] K-4PCS is a rough registration algorithm based on the four-point consistency principle of key points. 3D-DoG key points are extracted from the point cloud as a point set, and the rough registration is performed using radiometric invariance to obtain the initial transformation matrix. The voxelized generalized iterative closest point (VGICP) algorithm is an algorithm for accurately registering slightly moving point clouds. However, if the object moves significantly, the registration effect will be poor. Therefore, the K-4PCS registration algorithm adopted in the present invention is used for rough registration to obtain the initial transformation matrix, and under the initial transformation matrix, the VGICP algorithm is used to obtain the registered point cloud and the superimposed transformation matrix.

[0041] The specific algorithm includes:[[]] 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 grid-based filter, and then the key point set is extracted through a Difference of Gaussian (DoG) key point detector. 、 ; It is expressed as: ; Among them: represents the DoG response, represents the blur level, and x, y, z represent three-dimensional coordinates.

[0042] 2) Randomly extract a coplanar point set from the key point set , and use affine invariance to search for the corresponding set to obtain the initial transformation matrix and the target point cloud after transformation .

[0043] 3) The original point cloud , the target point cloud , the transformation matrix , each sampling point comes from a Gaussian distribution: , , The distance between it and its adjacent points: , The distribution of ; ; ; Estimate through the maximum likelihood estimation method, then there is: ; ; ; Wherein: is the number of adjacent points.

[0044] The transformation matrix can be calculated from the above formula, and the initial transformation matrix can be superimposed with it to obtain the deformation amount.

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

[0046] .

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

[0048] III. Research area data processing and analysis.

[0049] Extract the contour feature points of the three-phase data according to the above method. After multiple attempts of grid division, the grid size of 30m×30m has the best effect. The preliminary reason is related to the contour feature points being at the boundary of the point cloud block and the size of the three-dimensional deformation amount; obtain the three-dimensional deformation amounts of each contour position along the working face strike, working face dip, and subsidence direction according to the method of obtaining the transformation matrix by contour feature matching; use Kriging interpolation [01, 1] to interpolate the deformation amounts of discrete points. The discrete point weight function is distributed according to the reciprocal of the first power of the distance from the interpolation point. The farther the distance, the smaller the weight; this interpolation method is applicable to dense point clouds and can reduce the influence of noise on interpolation at the same time. The value of the interpolation point is equal to the weighted average of the nearby discrete points; draw the ground movement and deformation cloud map using the interpolated three-dimensional deformation amounts, and the results are as Figure 12 shown.

[0050] From Figure 12 it can be seen that: (1) The extraction results of the ground subsidence in each phase all show a good correspondence with the working face mining situation. The entire subsidence basin shows the characteristic of shifting towards the down-hill direction; the maximum subsidence value from May 2023 to April 2024 is -4.20m, the maximum subsidence value from May 2023 to August 2024 is -5.06m, and the maximum subsidence value from April 2023 to August 2024 is -5.10m.

[0051] (2) In the dip direction, the planar displacements along the working face in each period also show a good correspondence with the working face. The points in the downhill direction move towards the uphill direction, and the points in the uphill direction move towards the downhill direction. The symmetric position of the movement is located towards the downhill direction of the center of the working face, showing good consistency with the subsidence basin. Among them, the maximum dip displacement from May 2023 to April 2024 is -1.15 m, the maximum dip displacement from May 2023 to August 2024 is -1.48 m, and the maximum dip displacement from April 2023 to August 2024 is -1.36 m. The ratios of the maximum displacement to the maximum subsidence are approximately 0.35, 0.29, and 0.27 respectively, which conform to the general law of mining subsidence.

[0052] (3) In the strike direction, the planar displacement at the advancing position of the working face shows a direction pointing towards the center of the working face, which is consistent with the law of mining subsidence. However, there is no obvious planar displacement at the open-off cut position. Field investigation found that due to the overly shallow mining depth at the open-off cut position, the minimum is only 27 m. During the mining process of the working face, collapse pits appeared at this position, and the collapse pits were subsequently backfilled manually, resulting in the artificial disruption of the surface subsidence law in this area. Therefore, when analyzing the strike data, only the dynamic monitoring results between April 2024 and August 2024 are considered, and the data at the open-off cut position of the working face is no longer considered. Among the surface strike planar displacements from April 2024 to August 2024, the maximum value is 1.24 m, and the ratio of the maximum displacement to the maximum subsidence is approximately 0.24, which also conforms to the general law of mining subsidence.

[0053] To further analyze the reliability of the above-extracted data results, taking the data from April 2024 to August 2024 as an example, a theoretical calculation method is used for comparative analysis.

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

[0055] The random medium theory has developed into the probability integral method. This method believes that the laws of strata and surface movement caused by mining are macroscopically similar to the granular medium model as a random medium. According to the principle of the probability integral method, the settlement value of any point on the ground caused by mining can be expressed as: Among them: .

[0056] Among them: , respectively represent the subsidence distribution coefficients of the projection points of the point to be determined on the strike and dip main sections; , respectively represent the calculated lengths of the working face in the strike and dip directions; represents the mining unit of the probability integral; represents the coal seam thickness; denotes the subsidence coefficient; denotes the dip angle of the coal seam; denotes the main influence radius, , where denotes the average mining depth, denotes the main influence tangent angle.

[0057] The dynamic prediction of subsidence conducts an overall prediction and display of mining subsidence through parameters such as the working face coordinates, mining depth, and mining thickness. The formula is expressed as: Where: denotes the area of the region; denotes the maximum subsidence value under geological and mining conditions; denotes the subsidence at the mining time t.

[0058] Because the area is desertified and the superposition effect of the three-dimensional deformation obtained by the probability integral method and the research method of this study is as Figure 13 shown, and numerical values are extracted at a point every 5m on the horizontal displacement strike line and dip line as Figure 14 and Figure 15 shown.

[0059] The maximum displacement values in the Y direction of the dynamic prediction results of the probability integral method are 1.02m and -0.92m, the maximum displacement values in the X direction are 1.1m and -0.90m, and the maximum subsidence value is -5.10m. From the isogram and data comparison analysis, it can be obtained that the measured results are basically consistent with the theoretical calculation results, initially proving the rationality and feasibility of the research method proposed by the present invention; on the other hand, the main reasons for the deviation between the two include: (1) The terrain of this 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 impacts: this area is in a desertified nature reserve and is regularly treated to prevent water loss; the open-off cut has a large area of collapse and is filled.

[0060] (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.

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

[0062] (4) In addition to natural factors such as weather, the errors in the originally collected point cloud data are also affected by certain errors in the positioning and attitude determination of the unmanned aerial vehicle during its movement.

[0063] For special landforms such as deserts, mountains, and densely vegetated areas, as well as sensitive areas such as desertification areas and nature reserves, the present invention proposes a method for monitoring three-dimensional deformation of a mining area using point cloud contour registration to obtain the three-dimensional deformation of the surface in the X, Y, and Z directions, providing more diverse data for the analysis of surface subsidence laws and better guiding related activities such as mining in the mining area or the restoration of the ecological environment, etc.

[0064] Combined with the method of the present invention, the three-dimensional deformation data of the surface caused by mining in the research area is obtained. The maximum value of the ΔY deformation is 1.22 and -1.36 m, the maximum value of the movement in the ΔX direction is 1.24 and -0.87 m, and the maximum value of the ΔZ deformation is -5.10 m. It is analyzed that artificial landfill and relatively shallow coal seams result in no obvious movement and deformation law near the open-off cut from May 2023 to April 2024 and from May 2023 to August 2024.

[0065] The traditional monitoring methods cannot be carried out in the semi-desert research area mentioned above, which is also a nature reserve. The three-dimensional deformation caused by the mining activities from April 2024 to August 2024 is dynamically predicted using the probability integral method, and the comparison with the results obtained by the research method of the present invention proves its rationality and feasibility.

[0066] The above-described embodiments merely represent several implementation manners of the present invention, and their descriptions are relatively specific and detailed, but should not be construed as limiting the scope of the invention patent. It should be noted that for those of ordinary skill in the art, without departing from the concept of the present invention, several modifications and improvements can still be made, and these all belong to the protection scope of the present invention. Therefore, the protection scope of the present invention patent shall be subject to the appended claims.

Claims

1. A method for extracting three-dimensional deformation data of mining subsidence in desert areas, characterized in that, Including the following steps: Obtain three-phase surface point cloud data corresponding to different shooting time points in the desert area, and extract the overlapping areas in the three-phase surface point cloud data; wherein, the surface point cloud data corresponding to the overlapping areas consists of ground points and non-ground points; Obtain the normal vector angles between each point in the point cloud data of the ground points in each phase of the surface point cloud data in the overlapping area and the corresponding points in the other two phases, as well as the angles between each point and the tangent plane projection of its respective point normal vector; Set the points with normal vector angles higher than the preset value as contour feature points; and among the three points formed by the points with normal vector angles higher than the preset value and the corresponding points in the other two phases, set the point with the largest angle between the tangent plane projection of its respective point normal vector as the contour feature point; wherein, the points with normal vector angles higher than the preset value and the points with the largest angles between the tangent plane projections of their respective point normal vectors are used as a pair of point cloud contours, and multiple pairs of point cloud contour sets are obtained; Taking one-phase point cloud data as the original point cloud, for each point cloud data in the original point cloud, register the point cloud data in the other two-phase point cloud data that belongs to the same pair of point cloud contours as each point cloud data in the original point cloud, and obtain three-dimensional deformation data representing the displacement of the point cloud during the registration process.

2. The three-dimensional deformation data extraction method for mining subsidence in desert areas according to claim 1, wherein, The division of the ground points and non-ground points includes: After extracting the surface point cloud data in the overlapping area in the three-phase surface point cloud data, use the progressive densification triangular mesh filtering algorithm IPTD to divide the surface point cloud data in the overlapping area into ground points and non-ground points.

3. A method for extracting three-dimensional deformation data of mining subsidence in desert areas according to claim 2, characterized in that, After dividing the surface point cloud data corresponding to the overlapping area into ground points and non-ground points, construct a KD-tree tree structure for the surface point cloud data of the ground points corresponding to the three-phase overlapping area, including: For the point cloud data of the ground points corresponding to the three-phase overlapping area, establish a cubic bounding box containing all the point clouds according to the global coordinate system of the point cloud, and for the cubes containing more than 1 point, construct a splitting plane; The divided subspaces and the points on the splitting plane form branches and connection points to construct a KD-tree data structure.

4. A method for extracting three-dimensional deformation data of mining subsidence in desert areas according to claim 3, characterized in that, The obtaining of the normal vector angles between each point in the point cloud data of the ground points in each phase of the surface point cloud data in the overlapping area and the corresponding points in the other two phases, as well as the angles between each point and the tangent plane projection of its respective point normal vector, includes: In the kd-tree structure, set the search radius r, and denote the neighboring points within the search radius r as a set , that is ; Set the surface equation , , take The corresponding set , calculate The distance to the surface , solve The eigenvector corresponding to the minimum is the normal vector n of the corresponding point; The included angle of the normal vectors of adjacent points is expressed as ; According to and its normal vector n, a tangent plane of corresponding points is made , and the normal vector n is projected onto the tangent plane to form the tangent plane projection of the normal vector.

5. A method for extracting three-dimensional deformation data of mining subsidence in desert areas according to claim 4, characterized in that, The forming of multiple pairs of point cloud contour sets includes: Set the points with the adjacent point normal vector angles greater than 35° as contour feature points; According to and its normal vector n, construct the tangent plane at this point , project the points in the set onto the tangent plane , denoted as , select a point in , with as the u-axis, n as the axis, as the v-axis, and as the coordinate center to construct a local coordinate system, denoted as , calculate the vector from other points in the set to the point and the clockwise angle between the vector and the u-axis, and calculate the differences between adjacent angles pairwise to obtain the angle set , where , sort the elements in the set in descending order and find the largest angle , set the largest angle as the contour feature point; Form multiple pairs of point cloud contour sets according to two types of contour feature points.

6. The three-dimensional deformation data extraction method for mining subsidence in desert areas according to claim 1, characterized in that The obtaining of the three-dimensional deformation data includes: Set the target point cloud as the point set P to be registered, and set the point cloud to be registered as the reference point set Q; For the point set P to be registered and the reference point set Q, sampling is performed using a filter based on a voxelized grid, and then a key point set is extracted through a Difference of Gaussians (DoG) key point detector. , ; It is expressed as: ； Wherein: represents the DoG response, represents the blur level, and x, y, and z represent three-dimensional coordinates; Among the key point sets randomly extract a coplanar point set of 4 points , and use affine invariance to search for the corresponding set to obtain the initial transformation matrix; Based on the initial transformation matrix, calculate the similarity between the point clouds through the voxelization and generalized distance metric strategy of the VGICP algorithm, gradually align the point clouds, obtain the transformation matrix during the alignment process, and superimpose the transformation matrix during the alignment process with the initial transformation matrix to obtain the final transformation matrix, expressed as: ； The final transformation matrix includes the rotation matrix of R×3 in the point cloud coordinate system and the translation vector of 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 distances of the point cloud in the X, Y, and Z directions; Three-dimensional deformation data of mining subsidence in desert areas is formed according to the moving distances of the point cloud in the X, Y, and Z directions.

7. A three-dimensional deformation data extraction device for mining subsidence in desert areas, characterized in that, It includes: A data acquisition module, which is used to obtain three-phase surface point cloud data corresponding to different shooting time points in the desert area, and extract the overlapping areas in the three-phase surface point cloud data; among them, the surface point cloud data corresponding to the overlapping area is composed of ground points and non-ground points; A contour feature point determination module, which is used to obtain the normal vector angle between each point in the point cloud data of the ground points in each phase of the surface point cloud data in the overlapping area and the corresponding points in the other two phases, and the angle between each point and the tangent plane projection of its own 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 the preset value and the corresponding points in the other two phases, the point with the largest angle between the tangent plane projection of its own point normal vector is set as a contour feature point; among them, the points with a normal vector angle higher than the preset value and the points with the largest angle between the tangent plane projection of their own point normal vector are used as a pair of point cloud contours, and multiple pairs of point cloud contour sets are obtained; A registration module, which uses the point cloud data of one phase 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 belongs to the same pair of point cloud contours as each point cloud data in the original point cloud is registered to obtain three-dimensional deformation data representing the displacement of the point cloud during the registration process.

8. An electronic device, characterized in that, It includes: A memory and a processor; The memory is used to store computer programs; When the processor is used to execute the computer program stored in the memory, it realizes the steps of a method for extracting three-dimensional deformation data of mining subsidence in desert areas as described in any one of claims 1 to 6.

9. A computer-readable storage medium, characterized in that, For storing a computer program, when the computer program is executed by a processor, it realizes the steps of a method for extracting three-dimensional deformation data of mining subsidence in desert areas as described in any one of claims 1 to 6.

Images