An information extraction method for mining-induced surface subsidence based on unmanned aerial vehicle remote sensing images

Abstract

A kind of mining surface subsidence information extraction method based on unmanned aerial vehicle remote sensing image, comprising: determining monitoring area;Image map of monitoring area is obtained by unmanned aerial vehicle, and difference value DEM is generated according to the image map of monitoring area in different periods;The initial coordinate information of each sampling point in difference value DEM is obtained;According to the inclination of each sampling point, the subsidence value of each sampling point is corrected for the first time to eliminate accidental error;System error value is calculated by the data of the corresponding area in difference value DEM not affected by mining, and the subsidence value of each sampling point is corrected for the second time according to system error value to eliminate system error;Finally, the corrected subsidence value is generated to obtain the corrected difference value DEM, and the mining surface subsidence information is obtained.The design is not only low in cost, real-time, easy to obtain the change process of entire mining area, and the monitoring precision of subsidence information is very high.

Description

Technical Field

[0001] This invention relates to a method for extracting settlement information, and more particularly to a method for extracting settlement information of mined-affected surfaces based on UAV remote sensing imagery, specifically applicable to obtaining settlement information of mined-affected surfaces. Background Technology

[0002] Changes in surface plane and elevation during mining operations are crucial data for revealing the patterns of surface movement and deformation, constructing methods for predicting surface subsidence, and designing for surface subsidence disaster prevention and control. Traditional methods for monitoring surface elevation changes typically involve deploying ground-based mobile observation stations and employing techniques such as levels, total stations, and GNSS to measure changes in the coordinates of monitoring points during mining operations, thereby calculating changes in surface elevation coordinates. However, traditional monitoring methods suffer from drawbacks such as high labor intensity, long monitoring times, and insufficient real-time performance. Furthermore, they only monitor changes in the coordinates of individual monitoring points, failing to capture the entire process of change within the moving basin. Additionally, monitoring points deployed using traditional methods are easily damaged, making protection difficult.

[0003] In recent years, with the rapid development of modern surveying and mapping technology, new technologies such as UAV photogrammetry, 3D laser scanning, and differential interferometric radar have been increasingly applied to surface subsidence monitoring. These technological advancements have significantly improved the acquisition of deformation data in subsidence basins. UAV photogrammetry, with its advantages of low cost and high real-time performance, has gradually become a primary technical means for monitoring surface movement and deformation. UAV photogrammetry for monitoring surface subsidence mainly utilizes images to construct a digital elevation model (DEM) of the monitoring area, and then calculates the difference between DEMs from different periods to obtain the surface subsidence. However, the accuracy of UAV measurements is affected by many factors, such as the location, number, and quality of control points, flight parameters, and camera parameters. Currently, the horizontal accuracy of UAV surveying is approximately 2-3 cm, and the vertical accuracy is approximately 5-10 cm. Improving monitoring accuracy, especially vertical accuracy, is a major problem that UAV surveying needs to solve. Since monitoring the boundaries of subsidence basins generally requires a vertical accuracy better than 1 cm, the current accuracy of UAV surveying cannot meet the requirements for monitoring surface subsidence caused by mining. Furthermore, improving measurement accuracy remains a driving force for the development of surveying and mapping work.

[0004] Chinese Patent Publication No. CN106969751A, published on July 21, 2017, discloses a method for monitoring and calculating surface subsidence in coal mining based on UAV remote sensing. This method establishes control points based on the study area's scope and geographical features, generates point cloud data from UAV-captured images, and derives the three-dimensional coordinates of the study area's surface. Then, based on ground control survey results, a 4-parameter surface fitting correction model is used to correct the three-dimensional coordinates obtained from the UAV aerial survey, comparing them with the elevation values ​​from the pre-mining topographic map to calculate the post-mining surface subsidence. Chinese Patent Publication No. CN114279398A, published on April 5, 2022, demonstrates a method for monitoring surface subsidence in metal mining based on UAV aerial survey technology. This method generates a digital surface model based on UAV image information, then extracts surface subsidence monitoring points with fixed boundaries, and draws contour maps of the surface subsidence monitoring points. Finally, a 3D model of surface subsidence for this monitoring period was established based on contour maps. EPS software was used to monitor surface subsidence cracks, collecting the three-dimensional coordinates of the cracks and plotting them. Chinese Patent Publication No. CN114612806A, published on June 10, 2022, discloses a method for improving the accuracy of consumer-grade UAV DEM products. This method first constructs a gradient-cloth filtering model to filter the obtained UAV 3D point cloud to obtain a ground point cloud. A ground seed DEM is then constructed using this ground point cloud. Next, an elevation anomaly surface model is constructed based on the elevation difference between GNSS RTK measurements and UAV measurements. Finally, the elevation anomaly surface model is used to compensate and correct the ground seed DEM, thereby improving DEM accuracy. The above-mentioned existing technologies all monitor surface subsidence through UAV photography, but the following problems still exist:

[0005] 1. It is difficult to simultaneously capture the detailed features of the terrain and achieve a smooth effect. For example, in the patent with publication number CN106969751A, surface fitting is used to correct the coordinates of three-dimensional points monitored by UAVs to improve measurement accuracy. However, the actual ground surface is not a smooth curved surface, and filtering will inevitably eliminate the detailed features of the ground surface.

[0006] 2. Significant measurement errors still exist when using UAV photogrammetry to extract land subsidence. Current technologies for improving the accuracy of UAV photogrammetry mainly focus on data acquisition, such as selecting more suitable image control point layouts, but do not address the inherent errors of the method itself. Therefore, the current measurement accuracy for actual land subsidence is approximately 5-10 cm, which is insufficient for calculating land tilt values ​​and subsidence boundary ranges. Summary of the Invention

[0007] The purpose of this invention is to overcome the problems of traditional surface subsidence monitoring methods in the prior art, such as high labor intensity, long monitoring time, insufficient real-time performance, and difficulty in obtaining the change process of the entire mining area. On the other hand, monitoring surface subsidence using UAV photogrammetry technology cannot achieve the required monitoring accuracy. This invention provides a low-cost, real-time, and high-accuracy method for extracting surface subsidence information based on UAV remote sensing images, which can obtain the change process of the entire mining area.

[0008] To achieve the above objectives, the technical solution of the present invention is:

[0009] A method for extracting surface subsidence information caused by mining based on UAV remote sensing imagery, wherein the extraction method specifically includes:

[0010] Step 1: Collect images of the monitoring area at different times using drones, and generate a difference DEM based on the images of the monitoring area at different times;

[0011] Step 2: Set the sampling resolution of the difference DEM, and obtain the planar coordinates and initial subsidence value of the sampling points in the difference DEM at this sampling resolution;

[0012] Step 3: Calculate the inclination value of the sampling point based on its planar coordinates and initial subsidence value. Delineate the corresponding sampling area around the sampling point based on its inclination value. Use the average of the initial subsidence values ​​of all sampling points within the sampling area as the first corrected subsidence value of the corresponding sampling point.

[0013] Step 4: Subtract the system error value from the first corrected subsidence value of the sampling point to obtain the second corrected subsidence value of the sampling point;

[0014] Step 5: Generate a corrected difference DEM based on the plane coordinates of the sampling points and the second corrected subsidence value to obtain the mining-induced surface subsidence information.

[0015] In step three, calculating the tilt value of any sampling point N based on the planar coordinates of the sampling point and the initial subsidence value specifically includes:

[0016] Obtain the planar coordinates and initial subsidence values ​​of all sampling points within a region of radius L centered at sampling point N in the difference DEM. Use linear regression analysis to fit the relationship between the obtained planar coordinates and the obtained initial subsidence values ​​of the sampling points to obtain a linear regression equation. Use the slope of the linear regression equation as the inclination value of sampling point N.

[0017] The number of sampling points in the region with a radius of L centered on sampling point N is not less than 10.

[0018] The value of L is 0.5m.

[0019] In step three, calculating the tilt value of any sampling point N based on the planar coordinates of the sampling point and the initial subsidence value specifically includes:

[0020] A1. Calculate the tilt value of sampling point N in the east-west direction and the tilt value of sampling point N in the north-south direction;

[0021] A2. Calculate the vector sum of the tilt value of sampling point N in the east-west direction and the tilt value of sampling point N in the north-south direction to obtain the tilt value of sampling point N.

[0022] In step A1, calculating the tilt value of sampling point N in the east-west direction and the tilt value of sampling point N in the north-south direction specifically includes:

[0023] Sampling points N1 and N2 are two sampling points adjacent to sampling point N in the east-west direction. The difference between the initial sinking values ​​of sampling point N1 and N2 is calculated. This difference is then divided by the distance between sampling points N1 and N2 to obtain the east-west tilt value i of sampling point N. x ;

[0024] Sampling points N3 and N4 are two sampling points adjacent to sampling point N in the north-south direction. The difference between the initial subsidence values ​​of sampling point N3 and N4 is calculated. This difference is then divided by the distance between sampling points N3 and N4 to obtain the north-south tilt value i of sampling point N. y .

[0025] In step A2, calculating the vector sum of the tilt value of sampling point N in the east-west direction and the tilt value of sampling point N in the north-south direction to obtain the tilt value of sampling point N specifically includes:

[0026] The tilt value i of sampling point N is calculated according to formula (1):

[0027]

[0028] In step three, defining the corresponding sampling area around any sampling point N based on its tilt value specifically includes:

[0029] B1. Obtain the tilt value i of sampling point N, and calculate the side length l based on the tilt value i of sampling point N:

[0030]

[0031] In the formula, Δd represents the elevation measurement accuracy of the UAV, and i represents the tilt value of sampling point N;

[0032] B2. In the difference DEM, a square region with a side length of l is defined with the sampling point N as the center. The square region with a side length of l is the sampling region corresponding to the sampling point N.

[0033] The monitoring area includes areas affected by mining and areas not affected by mining. The system error value is calculated based on the data of the areas not affected by mining in the difference DEM.

[0034] Step one, which involves collecting images of the monitoring area at different times using a drone and generating a difference DEM based on these images, specifically includes:

[0035] C1. At time T1, image maps of the monitoring area are acquired using drones, and a first-phase original DEM is generated based on the image maps of the monitoring area acquired at time T1. At time T2, image maps of the monitoring area are acquired using drones, and a second-phase original DEM is generated based on the image maps of the monitoring area acquired at time T2.

[0036] C2. Calculate the difference between the original DEM of Phase I and the original DEM of Phase II, and generate the difference DEM.

[0037] In step C1, acquiring imagery of the monitoring area via a drone at time T1 and generating a primary original DEM based on the acquired imagery of the monitoring area at time T1 specifically includes:

[0038] C11. Set up image control points in the monitoring area by arranging image control boards, and measure the plane coordinates and elevation coordinates of the image control points;

[0039] C12. Use drones to photograph the monitored area;

[0040] C13. Based on the captured images, the planar coordinates and elevation coordinates of the control points, perform aerial triangulation to generate point cloud data;

[0041] C14. Remove non-ground points from the point cloud data and rasterize the point cloud data after removing non-ground points.

[0042] In step C1, the method of acquiring image maps of the monitoring area using a drone at time T2 and generating a second-phase original DEM based on the image maps of the monitoring area acquired at time T2 is the same as the method of acquiring image maps of the monitoring area using a drone at time T1 and generating a first-phase original DEM based on the image maps of the monitoring area acquired at time T1.

[0043] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0044] 1. In the present invention, a method for extracting mining-induced surface subsidence information based on UAV remote sensing imagery, the tilt value of each sampling point is first calculated based on the planar coordinates and initial subsidence value of each sampling point. Then, a first corrected subsidence value is calculated based on the tilt value of each sampling point to eliminate the influence of random errors on the subsidence value of each sampling point in the differential DEM. Subsequently, a systematic error value is calculated using data from areas in the differential DEM unaffected by mining. A second corrected subsidence value is then calculated based on the systematic error value and the first corrected subsidence value of each sampling point to eliminate the influence of systematic errors on the subsidence value of each sampling point in the differential DEM. Finally, the planar coordinates and second corrected subsidence value of the sampling points are imported into Global Mapper software to generate a corrected differential DEM, thus obtaining mining-induced surface subsidence information. The corrected differential DEM generated after eliminating the influence of random and systematic errors has higher accuracy, resulting in more accurate mining-induced surface subsidence information and significantly improving the accuracy of subsidence monitoring. Therefore, this design eliminates the influence of random and systematic errors on the measured surface subsidence information through filtering, which can improve the accuracy of subsidence monitoring to the millimeter level, significantly improving the accuracy of subsidence monitoring compared with traditional monitoring methods.

[0045] 2. In the method for extracting surface subsidence information based on UAV remote sensing imagery of the present invention, when calculating the first corrected subsidence value of a sampling point, the initial subsidence values ​​of all sampling points within a certain area centered on that sampling point are averaged, and the resulting average value is assigned as the first corrected subsidence value of that sampling point. Utilizing the property of accidental error compensation, averaging the initial subsidence values ​​of multiple sampling points that are relatively close together can cancel out accidental errors, thereby enabling the first corrected subsidence value to more accurately reflect the surface subsidence situation of the sampling point and improving the monitoring accuracy of the subsidence value. Therefore, this design utilizes the property of accidental error compensation, combined with the measurement accuracy of the UAV, to average the initial subsidence values ​​of multiple sampling points that are relatively close together, and performs spatial domain noise reduction on the subsidence value data of sampling points in two DEM phases, eliminating the influence of accidental errors and improving the monitoring accuracy.

[0046] 3. In the method for extracting surface subsidence information based on UAV remote sensing imagery of the present invention, when calculating the first corrected subsidence value of a sampling point, the initial subsidence values ​​of all sampling points within a certain area centered on that sampling point are averaged, and the resulting average value is used as the first corrected subsidence value of the corresponding sampling point. The size of the area used to calculate the first corrected subsidence value is jointly determined by the tilt value of the corresponding sampling point and the UAV's shooting accuracy. By utilizing the matching characteristics between local surface tilt values ​​and UAV measurement accuracy, measurement points with equal elevation and the same accuracy are determined to ensure the accuracy of monitoring. Therefore, in this design, the size of the area used to calculate the first corrected subsidence value is jointly determined by the tilt value of the corresponding sampling point and the UAV's shooting accuracy to ensure the accuracy of monitoring. Attached Figure Description

[0047] Figure 1 This is a flowchart of a method for extracting information on ground subsidence caused by land surveying.

[0048] Figure 2 This is a schematic diagram of the original DEM of Phase I and Phase II.

[0049] Figure 3 This is a schematic diagram of the subsidence data of the area in the differential DEM that is not affected by mining.

[0050] Figure 4 This is a schematic diagram of the difference DEM.

[0051] Figure 5 This is a schematic diagram of the corrected difference DEM.

[0052] Figure 6 This is a schematic diagram of the path profile of the difference DEM.

[0053] Figure 7 This is a schematic diagram of the corrected difference DEM path profile.

[0054] Figure 8 It is a contour map of the difference DEM.

[0055] Figure 9 This is the corrected difference DEM contour map. Detailed Implementation

[0056] The present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0057] See Figures 1 to 9 A method for extracting surface subsidence information based on UAV remote sensing imagery, the extraction method comprising the following steps:

[0058] Step 1: Collect images of the monitoring area at different times using drones, and generate a difference DEM based on the images of the monitoring area at different times. Specifically, this includes:

[0059] A monitoring area is selected, encompassing both areas affected by mining activities and those unaffected. According to the "Digital Aerial Photogrammetry" mapping specifications, image control points (ADCs) are set up above the mined-out area of ​​the UAV within the monitoring area using image control boards. The spacing between the ADCs can be set to 50m. The planar and elevation coordinates of the ADCs are measured. The planar coordinates of the ADCs can be obtained using RTK (Real-Time Kinematic) technology, and the elevation coordinates can be obtained using third-order leveling.

[0060] Select the appropriate drone model and flight parameters based on the shooting requirements. Use the drone to film the monitored area. In actual shooting, a Phantom 4 drone can be selected. Set the drone's flight altitude, forward overlap rate, and lateral overlap rate, and control the drone to fly along a pre-set route. Control the drone to fly along the pre-set route at times T1 and T2 to film the monitored area. The drone's flight altitude can be set to 50m, the forward overlap rate to 70%, and the lateral overlap rate to 80%. Based on experimental experience, the drone's elevation measurement accuracy is 3cm under the above measurement equipment and methods.

[0061] The captured images are then imported into photogrammetry software, such as Pix4dmapper. In Pix4dmapper, the CGCS2000 coordinate system is selected, and the planar and elevation coordinates of each control point are imported into the Pix4dmapper software. After point piercing and aerial triangulation calculations, point cloud data is generated.

[0062] The acquired point cloud data is preprocessed to remove various noise points and non-ground points such as vegetation. The point cloud data after removing non-ground points is then imported into map drawing software for rasterization processing, and finally a high-precision original DEM is generated. The map drawing software mentioned above can be Global Mapper.

[0063] At time T1, images are taken of the monitored area according to the above steps, and the acquired images are processed to obtain a primary raw DEM. At time T2, images are taken of the monitored area according to the above steps, and the acquired images are processed to obtain a secondary raw DEM. The primary and secondary raw DEMs are as follows: Figure 2 As shown.

[0064] After obtaining the original DEM of Phase I and Phase II, use map drawing software to generate the difference DEM. The specific steps are as follows: import the original DEM of Phase I and Phase II into the map drawing software Global mapper, select the "Combined Comparison Topographic Layer" tool to calculate the difference between the original DEM of Phase I and Phase II, and generate the difference DEM.

[0065] Step 2: Set the sampling resolution of the difference DEM, and obtain the planar coordinates and initial subsidence value of the sampling points in the difference DEM at this sampling resolution. The specific steps are as follows:

[0066] After generating the difference DEM values ​​using the map-drawing software Global Mapper, the sampling resolution of the difference DEM was set to 10 cm. Since the sampling resolution of the difference DEM is 10 cm, the interval between two adjacent sampling points within the study area is 10 cm. The initial subsidence values ​​of all sampling points within the study area were extracted, obtaining the planar coordinates and initial subsidence values ​​of all sampling points in the difference DEM. The initial subsidence value of each sampling point is its elevation value in the difference DEM.

[0067] Since only a portion of the monitoring area needs to be studied, the study area can be delineated in the difference DEM. In step two, only the planar coordinates and initial subsidence values ​​of the sampling points within the study area are extracted and processed. Similarly, if necessary, the planar coordinates and initial subsidence values ​​of the sampling points throughout the entire monitoring area can also be extracted and processed.

[0068] Step 3: Calculate the inclination value of the sampling point based on its planar coordinates and initial subsidence value. Delineate a corresponding sampling area around the sampling point based on its inclination value. Use the average of the initial subsidence values ​​of all sampling points within the sampling area as the first corrected subsidence value for the corresponding sampling point. Specifically, the steps for calculating the inclination value of any sampling point N based on its planar coordinates and initial subsidence value are as follows:

[0069] The planar coordinates and initial subsidence values ​​of all sampling points within a radius L centered at sampling point N in a planar coordinate system are obtained. Linear regression analysis is then used to fit the relationship between the obtained planar coordinates and subsidence values ​​of the sampling points, resulting in a linear regression equation. The slope of this linear regression equation is the inclination value of sampling point N. During the fitting process, the planar coordinates of the sampling points (in the difference DEM, the planar coordinates of the sampling points can be geographic coordinates or coordinates in other coordinate systems) are used as the independent variable, and the initial subsidence value of the sampling points is used as the dependent variable.

[0070] The radius L of the region can be 0.5m, and the number of sampling points in the region with radius L centered on sampling point N is not less than 10.

[0071] When calculating the inclination value of a sampling point based on its planar coordinates and initial subsidence value, the following method can also be used to calculate the inclination value of any sampling point N:

[0072] Export the planar coordinates and initial subsidence values ​​of the sampling points in the difference DEM along the east-west and north-south directions respectively, to obtain initial subsidence value files sorted along two different directions.

[0073] Based on the exported initial sinking value files sorted along two different directions, relevant code is written to iterate through each sampling point in the initial sinking value files sorted along the east-west direction. The difference in initial sinking values ​​between sampling points N1 and N2, which are adjacent to sampling point N in the east-west direction, is calculated. The difference between the initial sinking values ​​of sampling point N1 and N2, divided by the horizontal distance between sampling points N1 and N2, yields the tilt value i of sampling point N in the east-west direction. x By writing relevant code, the sampling points in the initial sinking value file sorted along the north-south direction are traversed step by step, and the difference between the initial sinking values ​​of two sampling points N3 and N4 that are adjacent to sampling point N in the north-south direction is calculated. Then, the difference between the initial sinking values ​​of sampling point N3 and N4 is divided by the horizontal distance between sampling point N3 and sampling point N4 to obtain the tilt value i of sampling point N in the north-south direction. y .

[0074] The tilt value i of sampling point N in the east-west direction is calculated according to formula (1). x The tilt value i of sampling point N in the north-south direction y The vector sum is used to obtain the tilt value i of sampling point N:

[0075]

[0076] When calculating the tilt value i of sampling point N, the global maximum tilt value can be used to replace the tilt value i of sampling point N.

[0077] After calculating the tilt values ​​of all sampling points within the study area, the tilt values ​​of all sampling points are compiled into a single file and saved.

[0078] In step three above, defining the corresponding sampling area around the sampling point N based on its tilt value includes the following steps:

[0079] Obtain the tilt value of sampling point N, and calculate the side length l based on the tilt value of sampling point N:

[0080]

[0081] In formula (2), Δd is the accuracy of the image map of the monitoring area at time T1 and the image map of the monitoring area at time T2, and i is the tilt value of sampling point N.

[0082] Since the elevation measurement accuracy of the UAV is 3cm under the selected measurement equipment and measurement method, if the tilt value of the sampling point N is 30mm / m, the side length l can be calculated to be 1m according to formula (2).

[0083] In the planar coordinate system of the difference DEM, a square region with a side length of l is defined centered on sampling point N. Since the side length l is 1m, a 1m*1m region is defined in the difference DEM centered on sampling point N; this region is the sampling area corresponding to sampling point N. Because the sampling resolution of the difference DEM is 10cm, there are 121 sampling points within this 1m*1m region. The first corrected settlement value of sampling point N is obtained by calculating the average of the initial settlement values ​​of these 121 sampling points within the 1m*1m region.

[0084] When calculating the first corrected subsidence value at sampling point N, the delineated area used to calculate the first corrected subsidence value can be a square area with a side length of l, or an area of ​​other shapes centered on sampling point N. The size of the delineated area is determined by the tilt value of sampling point N and the elevation measurement accuracy of the UAV.

[0085] The first corrected subsidence value of all sampling points is calculated using the method described above, and all the calculated first corrected subsidence values ​​are saved.

[0086] Step 4: Subtract the system error value from the first corrected subsidence value of the sampling point to obtain the second corrected subsidence value of the sampling point. The system error value is calculated based on the data of the corresponding area in the difference DEM that is not affected by the mining.

[0087] Areas unaffected by mining should not be affected by mining operations, and no subsidence should occur in these areas. For example... Figure 3 As shown, data from the differential DEM corresponding to the area unaffected by mining activity reveals that the elevation value within this unaffected area is 10 mm. This indicates that due to systematic error, a 10 mm settlement is also observed within the unaffected area of ​​the differential DEM, thus confirming the systematic error value as 10 mm. Therefore, subtracting the systematic error value of 10 mm from the first corrected settlement value of each sampling point yields the second corrected settlement value. After correcting the settlement values ​​of the sampling points using the systematic error value, the second corrected settlement value for most sampling points within the unaffected area of ​​the differential DEM is zero or less than a preset threshold, which can be set to 5 mm.

[0088] Step 5: Generate a corrected difference DEM based on the plane coordinates of the sampling points and the second corrected subsidence value to obtain the mining-induced surface subsidence information.

[0089] Generating the corrected difference DEM based on the planar coordinates and second corrected subsidence values ​​of each sampling point specifically involves: importing the planar coordinates and second corrected subsidence values ​​of all sampling points into Global Mapper; using the second corrected subsidence values ​​of the sampling points as elevation values; constructing an elevation triangulation network; and calculating and regenerating the corrected difference DEM based on the elevation triangulation network. Specifically, importing the planar coordinates and second corrected subsidence values ​​of all sampling points into Global Mapper software, and using the "Analyze"—"Create Elevation Grid from 3D Vector Data" tool in Global Mapper, the spatially filtered difference DEM, i.e., the corrected difference DEM, can be generated. The difference DEM and the corrected difference DEM are respectively as follows... Figure 4 , Figure 5 As shown in the figure.

[0090] Path analysis was performed on both the original and corrected difference DEMs: In Global Mapper software, the "Path" profile tool was selected, and a path was chosen. The sink values ​​of this path in both the original and corrected difference DEMs were extracted. The sink values ​​of the selected path in the original and corrected difference DEMs are shown below. Figure 6 , Figure 7 As shown, analysis Figure 6 , Figure 7 As can be seen from the two profile curves, the corrected difference DEM not only retains the overall downward trend of the original data and significantly reduces the impact of measurement errors, but also enhances the continuity and correlation between adjacent data. Therefore, averaging the downward values ​​within a certain range based on the magnitude of the tilt value can effectively reduce random measurement errors and improve monitoring accuracy.

[0091] Contour analysis was performed on the surface subsidence information caused by mining: Initial subsidence data and second corrected subsidence data were imported into Origin software, and differential DEM contour maps and corrected differential DEM contour maps were generated in Origin software at 0.1m intervals based on the initial and second corrected subsidence data, respectively. The differential DEM contour maps and corrected differential DEM contour maps are shown below. Figure 8 , Figure 9 As shown, analysis Figure 8 , Figure 9 As can be seen from the contour lines, the contour map generated based on the second corrected subsidence value data eliminates most of the abrupt deformation areas, enhances the coherence between contour lines, and makes the monitoring data more reliable. This is of great significance for determining the extent and degree of deformation of the subsidence basin.

[0092] The principle of this invention is explained as follows:

[0093] The term "mining" refers to mining and excavation activities such as coal mining, oil mining, and salt mining. After mining activities such as coal mining, the surface points will subside and move horizontally. Uneven subsidence and horizontal movement will cause the surface to tilt. During this process, it is necessary to observe and monitor the changes in the surface.

[0094] DEM is an abbreviation for Digital Elevation Model.

[0095] When using drones to monitor ground subsidence caused by mining, under the same measurement mode, the measurement error of the drone monitoring point coordinates includes both systematic and random errors. Systematic errors are relatively fixed, mainly related to the drone's equipment and flight parameters, and can be considered stable within a given data acquisition cycle. Random errors, on the other hand, follow a normal distribution. Both systematic and random errors affect the accuracy of the monitoring results; therefore, it is necessary to minimize their impact during data processing.

[0096] Random errors conform to a normal distribution and possess four basic characteristics: density, symmetry, compensation, and boundedness. Because random errors are compensatory—meaning that for a set of observations with the same precision, random errors can be positive or negative, and the probabilities of positive and negative values ​​are roughly equal—they can cancel each other out when a large amount of observation data is averaged. This ultimately eliminates the impact of random errors.

[0097] When using drones to monitor surface subsidence caused by mining, the difference between two monitoring digital elevation models (DEMs) is typically used to obtain a difference DEM, which is then used to assess surface subsidence. The elevation difference between any two sampling points in the difference DEM contains both measurement error and surface subsidence. When the tilt values ​​of adjacent sampling points are small, the subsidence difference between them is also small. When the difference in subsidence values ​​between adjacent sampling points is less than the measurement error of the drone, it is assumed that the drone measurement cannot accurately measure the subsidence difference between these two sampling points, and the subsidence of these two points can be approximated as the same. Therefore, when the tilt values ​​of a sampling point and all sampling points within a certain range are less than a certain threshold, it can be considered that the subsidence of this sampling point is the same as that of multiple sampling points within a certain range. In this case, the elevation data of these sampling points only contain error and the same surface subsidence. Therefore, averaging the elevation values ​​of these sampling points as the elevation of that point eliminates the influence of random errors on the elevation of that point. In this design, when calculating the first corrected sinking value of a certain sampling point N, an area is delineated around the sampling point N based on the tilt value of the sampling point N and the shooting accuracy of the UAV. The average of the initial sinking values ​​of all sampling points within the delineated area is used as the first corrected sinking value of the sampling point N, which can effectively eliminate the influence of random errors on the sinking value of the sampling point.

[0098] Systematic errors, primarily related to the UAV's equipment and flight parameters, can be considered stable within a single data acquisition and can be calculated using settlement results from areas unaffected by mining. In areas unaffected by settlement, the theoretical difference between two measurements should be zero. However, the actual difference is usually not zero. Therefore, the difference in settlement values ​​monitored in unaffected areas is the systematic error. Subtracting the systematic error value from all measurements within the entire monitoring area eliminates its influence. In this design, calculating the systematic error value using data from unaffected areas and correcting the settlement data at each sampling point based on this value effectively eliminates the impact of systematic errors.

[0099] Example 1:

[0100] A method for extracting surface subsidence information caused by mining based on UAV remote sensing imagery, wherein the extraction method specifically includes:

[0101] Step 1: Collect images of the monitoring area at different times using drones, and generate a difference DEM based on the images of the monitoring area at different times;

[0102] Step 2: Set the sampling resolution of the difference DEM, and obtain the planar coordinates and initial subsidence value of the sampling points in the difference DEM at this sampling resolution;

[0103] Step 3: Calculate the inclination value of the sampling point based on its planar coordinates and initial subsidence value. Delineate the corresponding sampling area around the sampling point based on its inclination value. Use the average of the initial subsidence values ​​of all sampling points within the sampling area as the first corrected subsidence value of the corresponding sampling point.

[0104] Step 4: Subtract the system error value from the first corrected subsidence value of the sampling point to obtain the second corrected subsidence value of the sampling point;

[0105] Step 5: Generate a corrected difference DEM based on the plane coordinates of the sampling points and the second corrected subsidence value to obtain the mining-induced surface subsidence information.

[0106] The monitoring area includes areas affected by mining and areas not affected by mining. The system error value is calculated based on the data of the areas not affected by mining in the difference DEM.

[0107] Step one, which involves collecting images of the monitoring area at different times using a drone and generating a difference DEM based on these images, specifically includes:

[0108] C1. At time T1, image maps of the monitoring area are acquired using drones, and a first-phase original DEM is generated based on the image maps of the monitoring area acquired at time T1. At time T2, image maps of the monitoring area are acquired using drones, and a second-phase original DEM is generated based on the image maps of the monitoring area acquired at time T2.

[0109] C2. Calculate the difference between the original DEM of Phase I and the original DEM of Phase II, and generate the difference DEM.

[0110] In step C1, acquiring imagery of the monitoring area via a drone at time T1 and generating a primary original DEM based on the acquired imagery of the monitoring area at time T1 specifically includes:

[0111] C11. Set up image control points in the monitoring area by arranging image control boards, and measure the plane coordinates and elevation coordinates of the image control points;

[0112] C12. Use drones to photograph the monitored area;

[0113] C13. Based on the captured images, the planar coordinates and elevation coordinates of the control points, perform aerial triangulation to generate point cloud data;

[0114] C14. Remove non-ground points from the point cloud data and rasterize the point cloud data after removing non-ground points.

[0115] In step C1, the method of acquiring image maps of the monitoring area using a drone at time T2 and generating a second-phase original DEM based on the image maps of the monitoring area acquired at time T2 is the same as the method of acquiring image maps of the monitoring area using a drone at time T1 and generating a first-phase original DEM based on the image maps of the monitoring area acquired at time T1.

[0116] Example 2:

[0117] Example 2 is basically the same as Example 1, except that:

[0118] In step three, defining the corresponding sampling area around any sampling point N based on its tilt value specifically includes:

[0119] B1. Obtain the tilt value i of sampling point N, and calculate the side length l based on the tilt value i of sampling point N:

[0120]

[0121] In the formula, Δd represents the elevation measurement accuracy of the UAV, and i represents the tilt value of sampling point N;

[0122] B2. In the difference DEM, a square region with a side length of l is defined with the sampling point N as the center. The square region with a side length of l is the sampling region corresponding to the sampling point N.

[0123] Example 3:

[0124] Example 3 is basically the same as Example 2, except that:

[0125] In step three, calculating the tilt value of any sampling point N based on the planar coordinates of the sampling point and the initial subsidence value specifically includes:

[0126] Obtain the planar coordinates and initial subsidence values ​​of all sampling points within a region of radius L centered at sampling point N in the difference DEM. Use linear regression analysis to fit the relationship between the obtained planar coordinates and the obtained initial subsidence values ​​of the sampling points to obtain a linear regression equation. Use the slope of the linear regression equation as the inclination value of sampling point N.

[0127] The number of sampling points in the area with a radius of L centered on sampling point N is not less than 10, and the value of L is 0.5m.

[0128] Example 4:

[0129] Example 4 is basically the same as Example 2, except that:

[0130] In step three, calculating the tilt value of any sampling point N based on the planar coordinates of the sampling point and the initial subsidence value specifically includes:

[0131] A1. Calculate the tilt value of sampling point N in the east-west direction and the tilt value of sampling point N in the north-south direction;

[0132] A2. Calculate the vector sum of the tilt value of sampling point N in the east-west direction and the tilt value of sampling point N in the north-south direction to obtain the tilt value of sampling point N.

[0133] In step A1, calculating the tilt value of sampling point N in the east-west direction and the tilt value of sampling point N in the north-south direction specifically includes:

[0134] Sampling points N1 and N2 are two sampling points adjacent to sampling point N in the east-west direction. The difference between the initial sinking values ​​of sampling point N1 and N2 is calculated. This difference is then divided by the distance between sampling points N1 and N2 to obtain the east-west tilt value i of sampling point N. x ;

[0135] Sampling points N3 and N4 are two sampling points adjacent to sampling point N in the north-south direction. The difference between the initial subsidence values ​​of sampling point N3 and N4 is calculated. This difference is then divided by the distance between sampling points N3 and N4 to obtain the north-south tilt value i of sampling point N. y .

[0136] In step A2, calculating the vector sum of the tilt value of sampling point N in the east-west direction and the tilt value of sampling point N in the north-south direction to obtain the tilt value of sampling point N specifically includes:

[0137] The tilt value i of sampling point N is calculated according to formula (1):

[0138]

[0139] The above description is only a preferred embodiment of the present invention. The scope of protection of the present invention is not limited to the above embodiments. Any equivalent modifications or changes made by those skilled in the art based on the content disclosed in the present invention should be included within the scope of protection set forth in the claims.

Claims

1. A method for extracting surface subsidence information based on UAV remote sensing imagery, characterized in that: The method for extracting surface subsidence information caused by mining specifically includes: Step 1: Collect images of the monitoring area at different times using drones, and generate a difference DEM based on the images of the monitoring area at different times; Step 2: Set the sampling resolution of the difference DEM, and obtain the planar coordinates and initial subsidence value of the sampling points in the difference DEM at this sampling resolution; Step 3: Calculate the inclination value of the sampling point based on its planar coordinates and initial subsidence value. Delineate the corresponding sampling area around the sampling point based on its inclination value. Use the average of the initial subsidence values ​​of all sampling points within the sampling area as the first corrected subsidence value of the corresponding sampling point. Step 4: Subtract the system error value from the first corrected subsidence value of the sampling point to obtain the second corrected subsidence value of the sampling point; Step 5: Generate a corrected difference DEM based on the plane coordinates of the sampling points and the second corrected subsidence value to obtain the mining-induced surface subsidence information.

2. The method for extracting surface subsidence information based on UAV remote sensing imagery according to claim 1, characterized in that: In step three, calculating the tilt value of any sampling point N based on the planar coordinates of the sampling point and the initial subsidence value specifically includes: Obtain the planar coordinates and initial subsidence values ​​of all sampling points within a region of radius L centered at sampling point N in the difference DEM. Use linear regression analysis to fit the relationship between the obtained planar coordinates and the obtained initial subsidence values ​​of the sampling points to obtain a linear regression equation. Use the slope of the linear regression equation as the inclination value of sampling point N. The number of sampling points in the region with a radius of L centered on sampling point N is not less than 10.

3. The method for extracting surface subsidence information based on UAV remote sensing imagery according to claim 2, characterized in that: The value of L is 0.5m.

4. The method for extracting surface subsidence information based on UAV remote sensing imagery according to claim 1, characterized in that: In step three, calculating the tilt value of any sampling point N based on the planar coordinates of the sampling point and the initial subsidence value specifically includes: A1. Calculate the tilt value of sampling point N in the east-west direction and the tilt value of sampling point N in the north-south direction; A2. Calculate the vector sum of the tilt value of sampling point N in the east-west direction and the tilt value of sampling point N in the north-south direction to obtain the tilt value of sampling point N.

5. The method for extracting surface subsidence information based on UAV remote sensing imagery according to claim 4, characterized in that: In step A1, calculating the tilt value of sampling point N in the east-west direction and the tilt value of sampling point N in the north-south direction specifically includes: Sampling points N1 and N2 are two sampling points adjacent to sampling point N in the east-west direction. The difference between the initial sinking values ​​of sampling point N1 and N2 is calculated. This difference is then divided by the distance between sampling points N1 and N2 to obtain the east-west tilt value i of sampling point N. x ; Sampling points N3 and N4 are two sampling points adjacent to sampling point N in the north-south direction. The difference between the initial subsidence values ​​of sampling point N3 and N4 is calculated. This difference is then divided by the distance between sampling points N3 and N4 to obtain the north-south tilt value i of sampling point N. y .

6. The method for extracting surface subsidence information based on UAV remote sensing imagery according to claim 5, characterized in that: In step A2, calculating the vector sum of the tilt value of sampling point N in the east-west direction and the tilt value of sampling point N in the north-south direction to obtain the tilt value of sampling point N specifically includes: The tilt value i of sampling point N is calculated according to formula (1):

7. A method for extracting surface subsidence information based on UAV remote sensing imagery according to any one of claims 1-6, characterized in that: In step three, defining the corresponding sampling area around any sampling point N based on its tilt value specifically includes: B1. Obtain the tilt value i of sampling point N, and calculate the side length l based on the tilt value i of sampling point N: In the formula, Δd represents the elevation measurement accuracy of the UAV, and i represents the tilt value of sampling point N; B2. In the difference DEM, a square region with a side length of l is defined with the sampling point N as the center. The square region with a side length of l is the sampling region corresponding to the sampling point N.

8. The method for extracting surface subsidence information based on UAV remote sensing imagery according to claim 1, characterized in that: The monitoring area includes areas affected by mining and areas not affected by mining. The system error value is calculated based on the data of the areas not affected by mining in the difference DEM.

9. A method for extracting surface subsidence information based on UAV remote sensing imagery according to claim 8, characterized in that: Step one, which involves collecting images of the monitoring area at different times using a drone and generating a difference DEM based on these images, specifically includes: C1. At time T1, image maps of the monitoring area are acquired using drones, and a first-phase original DEM is generated based on the image maps of the monitoring area acquired at time T1. At time T2, image maps of the monitoring area are acquired using drones, and a second-phase original DEM is generated based on the image maps of the monitoring area acquired at time T2. C2. Calculate the difference between the original DEM of Phase I and the original DEM of Phase II, and generate the difference DEM.

10. A method for extracting surface subsidence information based on UAV remote sensing imagery according to claim 9, characterized in that: In step C1, acquiring imagery of the monitoring area via a drone at time T1 and generating a primary original DEM based on the acquired imagery of the monitoring area at time T1 specifically includes: C11. Set up image control points in the monitoring area by arranging image control boards, and measure the plane coordinates and elevation coordinates of the image control points; C12. Use drones to photograph the monitored area; C13. Based on the captured images, the planar coordinates and elevation coordinates of the control points, perform aerial triangulation to generate point cloud data; C14. Remove non-ground points from the point cloud data and rasterize the point cloud data after removing non-ground points. In step C1, the method of acquiring image maps of the monitoring area using a drone at time T2 and generating a second-phase original DEM based on the image maps of the monitoring area acquired at time T2 is the same as the method of acquiring image maps of the monitoring area using a drone at time T1 and generating a first-phase original DEM based on the image maps of the monitoring area acquired at time T1.

Images