Mining area subsidence dynamic monitoring method, device, equipment, medium and product

Abstract

The invention discloses a mining area subsidence dynamic monitoring method, device and equipment, a medium and a product, and relates to the field of mining area mining subsidence monitoring. The method comprises the following steps: firstly, carrying out data preprocessing on an unmanned aerial vehicle image and laser radar point cloud time sequence data observed in the same period to obtain DOM and DEM data in each period; based on each period of DOM data, extracting east-west and south-north horizontal displacements between any two periods of DOM data by using a normalized cross-correlation matching algorithm; for the newly-added DOM data, dynamically updating historical and newly-added east-west and south-north direction horizontal displacements according to a sequential least square adjustment principle to obtain adjusted east-west and south-north direction horizontal displacements; carrying out plane position correction on DEM data of each period based on the adjusted horizontal displacement in the east-west direction and the south-north direction; and performing difference calculation on the corrected DEM data of each period and the reference DEM data before deformation to obtain time sequence subsidence data of the mining area, thereby realizing dynamic monitoring of mining area subsidence, and improving mining area subsidence monitoring precision and efficiency.

Description

Technical Field

[0001] This application relates to the technical field of mining area subsidence monitoring, and particularly to a method, device, equipment, medium and product for dynamic monitoring of mining area subsidence. Background Art

[0002] Coal mining causes geological disasters such as surface subsidence, damage to buildings (structures), surface cracks and landslides, which seriously threaten the ecological environment of mining areas and the residential safety of residents. Therefore, obtaining high-precision three-dimensional surface deformation data of mining areas is of great significance for ensuring the ecological environment and residential safety of mining areas. Traditional mining area subsidence monitoring methods mainly include leveling, trigonometric leveling and global navigation satellite system measurement. Although these methods have high monitoring accuracy, they have deficiencies such as large workload, high cost and harsh working conditions. In contrast, unmanned aerial vehicle (UAV) photogrammetry and LiDAR (Light Detection and Ranging) technology have been widely used in mining area subsidence monitoring due to their advantages of high efficiency, flexibility and low cost. However, existing UAV photogrammetry and LiDAR technology mainly constructs subsidence basins based on the DEM (Digital Elevation Model) differential method. This method only reflects the ground elevation changes under the same geographical space coordinates and ignores the influence of horizontal displacement on vertical deformation. Summary of the Invention

[0003] Aiming at the problems pointed out in the background art, this application provides a method, device, equipment, medium and product for dynamic monitoring of mining area subsidence, which considers the influence of horizontal displacement on vertical deformation during the calculation of subsidence values, and improves the accuracy and efficiency of mining area subsidence monitoring.

[0004] To achieve the above object, this application provides the following solutions.

[0005] In a first aspect, this application provides a method for dynamic monitoring of mining area subsidence, including:

[0006] Obtaining UAV images and LiDAR point cloud time series data observed in the mining area during the same period and performing data preprocessing to obtain DOM and DEM data for each period; where the same period refers to the same observation moment; the DOM and DEM data for each period refer to the DOM and DEM data at each observation moment;

[0007] Based on the DOM data for each period, using the normalized cross-correlation matching algorithm to extract the horizontal displacements in the east-west and north-south directions between any two periods of DOM data;

[0008] For the newly added DOM data, according to the principle of sequential least squares adjustment, the historical horizontal displacements in the east-west and north-south directions and the newly added horizontal displacements in the east-west and north-south directions are dynamically updated to obtain the adjusted horizontal displacements in the east-west and north-south directions;

[0009] Based on the adjusted horizontal displacements in the east-west and north-south directions, the plane position correction is performed on each period of DEM data to obtain the corrected DEM data for each period;

[0010] The corrected DEM data for each period is differentially calculated with the reference DEM data before deformation to obtain the time-series subsidence data of the mining area;

[0011] Based on the time-series subsidence data of the mining area and the adjusted horizontal displacements in the east-west and north-south directions, the time-series three-dimensional deformation data of the mining area is constructed to realize the dynamic monitoring of the subsidence of the mining area.

[0012] Optionally, the obtaining of the time-series data of UAV images and lidar point clouds for synchronous observations in the mining area and the data preprocessing to obtain the DOM and DEM data for each period specifically include:

[0013] Using Pix4D software to process the UAV image data at each observation moment in the time-series data of UAV images to generate the corresponding DOM data for each period;

[0014] For the lidar point cloud data at each observation moment in the time-series data of lidar point clouds, an improved progressive densification triangular mesh filtering algorithm is used to remove the non-ground points in the lidar point cloud data, extract the ground points, and use the inverse distance weighted interpolation algorithm to interpolate the ground points to generate the corresponding DEM data.

[0015] Optionally, the extracting of the horizontal displacements in the east-west and north-south directions between any two periods of DOM data based on the DOM data for each period by using the normalized cross-correlation matching algorithm specifically includes:

[0016] Convert the selected two periods of DOM data into the corresponding grayscale images;

[0017] Take the grayscale image with an earlier observation moment as the reference image and the grayscale image with a later observation moment as the target image, and construct the corresponding matching windows in the reference image and the target image;

[0018] Calculate the similarity of the image blocks in the two matching windows through the normalized cross-correlation matching algorithm to obtain the normalized cross-correlation coefficient;

[0019] Adopt the winner-takes-all strategy, use the peak point of the normalized cross-correlation coefficient as the matching point, and the corresponding left-right disparity and up-down disparity are the relative displacements in the row and column directions respectively. Then use the sinc interpolation function to interpolate the relative displacements in the row and column directions respectively to obtain the relative displacements in the sub-pixel level row and column directions;

[0020] Combined with the projection coordinate system parameters and the pixel resolution, convert the relative displacements in the sub-pixel level row and column directions into the horizontal displacements in the east-west and north-south directions in the projection coordinate system.

[0021] Optionally, for the newly added DOM data, dynamically update the historical east-west and north-south horizontal displacements and the newly added east-west and north-south horizontal displacements according to the sequential least squares adjustment principle to obtain the adjusted east-west and north-south horizontal displacements, which specifically includes:

[0022] Based on the relationship between the east-west and north-south horizontal displacements between any two periods of DOM data, establish an error equation;

[0023] According to the least squares principle, obtain the corresponding normal equation from the error equation;

[0024] Based on the normal equation, perform robust estimation and obtain the historical east-west and north-south horizontal displacements based on the estimation results;

[0025] When the newly added DOM data is obtained, jointly perform adjustment calculation on the newly added DOM data and the historical DOM data, and establish a sequential adjustment error equation;

[0026] Combined with the sequential adjustment error equation, according to the Bayesian least squares criterion, construct a sequential adjustment normal equation;

[0027] Based on the sequential adjustment normal equation, perform robust estimation and obtain the adjusted east-west and north-south horizontal displacements based on the estimation results.

[0028] Optionally, based on the adjusted east-west and north-south horizontal displacements, perform plane position correction on each period of DEM data to obtain the corrected DEM data for each period, which specifically includes:

[0029] For each period of DEM data, taking the plane positions of each feature point in the pre-deformation DEM data as the reference, combined with the adjusted east-west and north-south horizontal displacements, correct the plane positions of the corresponding feature points in each period of the post-deformation DEM data, so that the plane positions of each corresponding feature point before and after deformation are in the same vertical direction;

[0030] Adopt the inverse distance weighted interpolation algorithm to interpolate each period of the post-deformation DEM data at the plane positions of the corrected feature points to obtain the corrected elevation values of each feature point, and form the corrected DEM data for each period.

[0031] Optionally, the differential calculation of the corrected DEM data of each period and the reference DEM data before deformation is performed to obtain the time-series subsidence data of the mining area, specifically including:

[0032] Directly subtract the corrected DEM data of each period from the reference DEM data before deformation to obtain the DEM differential data of each period, and arrange them according to the time series of each period to form the time-series subsidence data of the mining area.

[0033] In a second aspect, the present application provides a device for dynamically monitoring mining area subsidence, including:

[0034] A data preprocessing module, configured to obtain the time-series data of UAV images and lidar point clouds observed in the same period in the mining area and perform data preprocessing to obtain DOM and DEM data of each period; where the same period refers to the same observation moment; each period refers to each observation moment;

[0035] A horizontal displacement extraction module, configured to extract the horizontal displacements in the east-west and north-south directions between any two periods of DOM data based on the DOM data of each period by using the normalized cross-correlation matching algorithm;

[0036] A horizontal displacement sequential adjustment module, configured to perform dynamic update on the historical east-west and north-south horizontal displacements and the newly added east-west and north-south horizontal displacements according to the principle of sequential least squares adjustment for the newly added DOM data, to obtain the adjusted east-west and north-south horizontal displacements;

[0037] A planar position correction module, configured to perform planar position correction on the DEM data of each period based on the adjusted east-west and north-south horizontal displacements to obtain the corrected DEM data of each period;

[0038] A differential calculation module, configured to perform differential calculation on the corrected DEM data of each period and the reference DEM data before deformation to obtain the time-series subsidence data of the mining area;

[0039] A mining area subsidence dynamic monitoring module, configured to construct the time-series three-dimensional deformation data of the mining area based on the time-series subsidence data of the mining area and the adjusted east-west and north-south horizontal displacements, to realize the dynamic monitoring of the mining area subsidence.

[0040] In a third aspect, the present application provides a computer device, including: a memory, a processor, and a computer program stored on the memory and executable on the processor, where the processor executes the computer program to implement the mining area subsidence dynamic monitoring method.

[0041] In a fourth aspect, the present application provides a computer-readable storage medium, on which a computer program is stored, and when the computer program is executed by a processor, the mining area subsidence dynamic monitoring method is implemented.

[0042] In a fifth aspect, the present application provides a computer program product, including a computer program which, when executed by a processor, implements the dynamic monitoring method for mining area subsidence.

[0043] According to the specific embodiments provided by the present application, the following technical effects are disclosed.

[0044] In the dynamic monitoring method, device, equipment, medium and product for mining area subsidence provided by the present application, on the one hand, the influence of horizontal displacement in the east-west and north-south directions on vertical deformation is considered in the calculation process of time-series subsidence data of the mining area. The plane position of the DEM data is corrected by using the horizontal displacement, and the dynamic monitoring of mining area subsidence is carried out based on the corrected DEM data, which can significantly improve the accuracy of mining area subsidence monitoring; on the other hand, for historical and new DOM data, the high-efficiency and dynamic update of historical and new horizontal displacements in the east-west and north-south directions can be realized, thereby improving the efficiency of mining area subsidence monitoring. BRIEF DESCRIPTION OF THE DRAWINGS

[0045] In order to more clearly illustrate the technical solutions in the embodiments of the present application or the prior art, the following will briefly introduce the drawings required in the embodiments. Obviously, the drawings in the following description are only some embodiments of the present application. For those of ordinary skill in the art, other drawings can be obtained based on these drawings without creative efforts.

[0046] Figure 1 is a schematic flowchart of a dynamic monitoring method for mining area subsidence according to the present application;

[0047] Figure 2 is a schematic diagram of the main process of dynamic monitoring of mining area subsidence by using the method of the present application;

[0048] Figure 3 is a schematic diagram of the process of plane position correction for DEM data. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0049] The following will clearly and completely describe the technical solutions in the embodiments of the present application with reference to the drawings in the embodiments of the present application. Obviously, the described embodiments are only some embodiments of the present application, rather than all embodiments. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present application without creative efforts belong to the scope of protection of the present application.

[0050] The present application proposes a dynamic monitoring method, device, equipment, medium and product for mining area subsidence, aiming to establish a dynamic monitoring system for mining area subsidence that takes into account horizontal displacement and is based on time-series DOM / DEM data of unmanned aerial vehicles, and considers the influence of horizontal displacement on vertical deformation in the calculation process of subsidence values, so as to improve the accuracy and efficiency of mining area subsidence monitoring.

[0051] To make the above objects, features, and advantages of the present application more obvious and understandable, the present application will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0052] In an exemplary embodiment, as Figure 1 shown, the present application provides a method for dynamic monitoring of mining area subsidence, including the following steps 1 to 6.

[0053] Step 1: Obtain the UAV image and LiDAR point cloud time-series data for synchronous observation of the mining area and perform data preprocessing to obtain DOM and DEM data for each period.

[0054] Figure 2 shows the main process of dynamic monitoring of mining area subsidence using the method of the present application. Refer to Figure 2 , the present application uses UAV photogrammetry and LiDAR technology to observe the target mining area and obtain the UAV image and LiDAR point cloud time-series data for synchronous observation. Among them, the UAV image time-series data includes UAV image data at a total of n observation times (also simply referred to as times) T0, T1,..., T (n-1) . Each observation time is also called each period. Correspondingly, the LiDAR point cloud time-series data also includes LiDAR point cloud data at a total of n observation times T0, T1,..., T (n-1) . Where n is usually an integer greater than 3. For example, the UAV image data observed at time T (n-1) and the LiDAR point cloud data observed at time T (n-1) are the synchronous observation data for the (n - 1)-th period. The purpose of the data preprocessing stage is to generate DOM (Digital Orthophoto Map) and DEM data for times T0, T1,..., T (n-1) based on the UAV image and LiDAR point cloud time-series data respectively.

[0055] On the one hand, for the UAV image time-series data, use Pix4D software to process the UAV image data at these n observation times T0, T1,..., T (n-1) in the UAV image time-series data to generate the corresponding DOM data for each period, that is, generate a total of n periods of DOM data T0, T1,..., T (n-1) .

[0056] On the other hand, for the lidar point cloud data at each observation moment in the lidar point cloud time series data, an improved progressive TIN densification filtering algorithm is used to remove non-ground points such as vegetation and buildings, extract ground points, and use the inverse distance weighted interpolation algorithm to interpolate these ground points to generate corresponding DEM data. Among them, the improved progressive TIN densification filtering algorithm refers to Improved progressive TIN densification filtering algorithm for airborne LiDAR data in forested areas (Zhao et al., 2016).

[0057] Further, the generated DOM data and DEM data at times T0, T1, …, T (n-1) are georegistered to ensure that the coordinate systems and spatial resolutions of each period of DOM and DEM are consistent. The coverage ranges of the generated DOM and DEM are usually not exactly the same. Crop the common areas of each period of DOM and DEM to obtain DOM and DEM data for each period with consistent spatial coverage ranges and consistent planar positions of corresponding ground feature points.

[0058] Step 2: Based on the DOM data for each period, use the normalized cross-correlation matching algorithm to extract the horizontal displacements in the east-west and north-south directions between any two periods of DOM data.

[0059] Based on the DOM data at these n observation times T0, T1, …, T (n-1) , use the normalized cross-correlation matching algorithm to extract the horizontal displacements in the east-west and north-south directions between the DOM data observed at any two observation times (such as T 0,1 , T 0,2 , T 0,3 , …, T 0,(n-1) , T 1,2 , T 1,3 , …, T 1,(n-1) , …, T (n-2),(n-1) ). Among them, T (n-2),(n-1) represents the two observation times n - 2 and n - 1. See Figure 2 . The specific steps of step 2 include the following steps 2.1 to 2.5.

[0060] Step 2.1: Convert any two selected periods of DOM data into corresponding grayscale images.

[0061] The conversion formula for each period of DOM data is as follows:[[]]

[0062] I gray = 0.299R + 0.587G + 0.114B (1)

[0063] In the formula, I gray is the converted grayscale image; R, G, and B are the pixel values of the red, green, and blue channels of the DOM, respectively.

[0064] Step 2.2: Use the grayscale image with an earlier observation time as the reference image, and the grayscale image with a later observation time as the target image, and construct corresponding matching windows in the reference image and the target image.

[0065] Select the grayscale image corresponding to the DOM with an earlier acquisition time as the reference image I1, and the grayscale image corresponding to the DOM with a later acquisition time as the target image I2. In the reference image I1, construct a matching window with a size of h×h centered on the position of the pixel to be matched. In the target image I2, construct a matching window with the same size centered on the position of the target pixel.

[0066] Step 2.3: Calculate the similarity of the image blocks within the two matching windows through the normalized cross-correlation matching algorithm to obtain the normalized cross-correlation coefficient.

[0067] Calculate the similarity of the image blocks within these two matching windows through the normalized cross-correlation matching algorithm, and its calculation formula is:

[0068]

[0069] In the formula, NCC(u,v) represents the normalized cross-correlation coefficient, and its value range is [-1,1]. I1(x,y) is the pixel value at the pixel to be matched (x,y) in the reference image I1, and the mean value of all pixels within its matching window is W p is the matching window centered on the pixel to be matched. I2(x+u,y+v) is the pixel value of the central pixel (i.e., the target pixel) of the matching window in the target image I2, and the mean value of its corresponding matching window is (x+u,y+v) are the coordinates of the corresponding target pixel in the target image I2. u and v are the relative displacements in the row and column directions during the traversal process, and subsequently, the u and v corresponding to the peak point of the normalized cross-correlation function are used as the relative displacements u' and v' in the row and column directions of this matching point.

[0070] Step 2.4: Adopt the winner-takes-all strategy, use the peak point of the normalized cross-correlation coefficient as the matching point, and the corresponding left-right disparity and up-down disparity are the relative displacements in the row and column directions, respectively. Then use the sinc interpolation function to interpolate the relative displacements in the row and column directions respectively to obtain the relative displacements at the sub-pixel level in the row and column directions.

[0071] Adopt the winner-takes-all strategy, and use the peak point of the normalized cross-correlation coefficient NCC(u, v) as the matching point. The corresponding left-right disparity and up-down disparity at the matching point are the relative displacements u and v in the row and column directions respectively. Since the actual displacement usually does not fall on the integer pixel position, the peak of the corresponding normalized cross-correlation function often lies at the sub-pixel position. To accurately locate this peak, the sinc interpolation function is used to interpolate the relative displacements u and v in the row and column directions respectively, so as to achieve sub-pixel-level disparity estimation. Denote the relative displacement in the row direction at the sub-pixel level as u', and the relative displacement in the column direction at the sub-pixel level as v'.

[0072] Step 2.5: Combine the projection coordinate system parameters and the pixel resolution to convert the relative displacements in the row and column directions at the sub-pixel level into the horizontal displacements in the east-west and north-south directions in the projection coordinate system.

[0073] Combine the projection coordinate system parameters and the pixel resolution, and convert the relative displacements in the row and column directions at the sub-pixel level calculated in pixel units into the horizontal displacements in the east-west and north-south directions in the projection coordinate system according to the following formula (3):

[0074]

[0075] In the formula, u' and v' are the relative displacements in the row and column directions at the sub-pixel level respectively. The pixel resolution is Δ x *Δ y . D EW is the horizontal displacement in the east-west direction, and D NS is the horizontal displacement in the north-south direction. Subsequently, the horizontal displacements in the east-west and north-south directions between any two periods of DOM data, for example, denote the horizontal displacement in the east-west direction between the (n - 2)th and (n - 1)th periods of DOM data as Denote the horizontal displacement in the north-south direction between the (n - 2)th and (n - 1)th periods of DOM data as

[0076] Step 3: For the newly added DOM data, dynamically update the historical horizontal displacements in the east-west and north-south directions and the newly added horizontal displacements in the east-west and north-south directions according to the principle of sequential least squares adjustment to obtain the adjusted horizontal displacements in the east-west and north-south directions.

[0077] As mentioned above, this application has pre-collected the DOM data at n observation times T0, T1, …, T (n-1) , which are also called historical DOM data. Then for the DOM data at the subsequent nth observation time T n , this application refers to it as the newly added DOM data. Correspondingly,[[]] and are called the historical horizontal displacements in the east-west and north-south directions; and It is called the new horizontal displacement in the east-west and north-south directions. Next, based on the principle of sequential least squares adjustment, the historical horizontal displacements in the east-west and north-south directions and the new horizontal displacements in the east-west and north-south directions are dynamically updated, and the M-estimation method is introduced to reduce the influence of gross errors on the parameter estimation results, thereby improving the accuracy and reliability of parameter calculation. See Figure 2 Step 3 specifically includes the following steps 3.1 to 3.6.

[0078] Step 3.1: Based on the relationship between the horizontal displacements in the east-west and north-south directions between any two periods of DOM data, establish an error equation.

[0079] Based on the obtained T 0,1 , T 0,2 , …, T (n-2),(n-1) The relationship between the horizontal displacement observations in the east-west and north-south directions, the error equation is established as follows:

[0080]

[0081] The extended form of the error equation is as follows:

[0082]

[0083] In the formula, is the matrix of correction values of the historical horizontal displacement parameters in the east-west or north-south direction, and the element represents the correction value of the historical horizontal displacement parameter in the east-west or north-south direction between the 0th and (n - 1)th periods of DOM. is the matrix composed of the difference between the historical horizontal displacement in the east-west or north-south direction and the approximate value of the horizontal displacement in the east-west or north-south direction, and the element represents the difference between the historical horizontal displacement in the east-west or north-south direction between the (n - 2)th and (n - 1)th periods of DOM and its approximate value of the horizontal displacement in the east-west or north-south direction. is the matrix of corrections of the historical horizontal displacement in the east-west or north-south direction, and the element represents the correction of the historical horizontal displacement in the east-west or north-south direction between the (n - 2)th and (n - 1)th periods of DOM.

[0084] is a coefficient matrix of (n(n - 1) / 2)*(n - 1), which is related to the horizontal displacement observations in the east-west or north-south direction extracted from any two periods of DOM. The generation rule of each row of the matrix is as follows. For the horizontal displacement in the east-west or north-south direction extracted from the i, j periods of DOM, where 0 ≤ i < j ≤ (n - 1); when i = 0, the element in the jth column of the corresponding row in the matrix is 1, and the remaining elements are 0; when i > 0, the matrix The element in the \(i\)-th column of the corresponding row is -1, the element in the \(j\)-th column is 1, and the remaining elements are 0. For example, when \(i = 0\) and \(j = 1\), the corresponding correction is matrix The element in the first column of the corresponding one is 1, and the remaining elements are 0. When \(i = 0\) and \(j = 2\), the corresponding correction is matrix The element in the second column of the corresponding row is 1, and the remaining elements are 0. When \(i = 1\) and \(j = 2\), the corresponding correction is matrix The element in the first column of the corresponding row is -1, the element in the second column is 1, and the remaining elements are 0. And so on.

[0085] For more convenient understanding, the meanings of the elements in the entire error equation (5) are illustrated by examples below. For example, assume that the horizontal displacement in the east-west direction between the DOM of the 3rd period and the DOM of the 4th period extracted is Due to the existence of errors, the observed value needs to be adjusted, and the horizontal displacement in the east-west direction after adjustment is Then According to the indirect adjustment principle, the horizontal displacements in the east-west direction calculated at the selected observation times \(T\) 0,1 , \(T\) 0,2 , \(T\) 0,3 , …, \(T\) 0,(n-1) are taken as parameters During adjustment, generally an approximate value \(D'\) is taken for the parameter EW , for example, the approximate value of Then When adjusting the horizontal displacement in the north-south direction, the principle is the same.

[0086] Step 3.2: According to the least squares principle, obtain the corresponding normal equation from the error equation.

[0087] According to the least squares principle, the normal equation can be obtained from formulas (4) and (5) as follows:

[0088]

[0089] where is the cofactor matrix of the historical horizontal displacement in the east-west or north-south direction, is the weight matrix of the observed values of the historical horizontal displacement in the east-west or north-south direction.

[0090] Step 3.3: Perform robust estimation based on the normal equation, and obtain the historical horizontal displacements in the east-west and north-south directions based on the estimation results.

[0091] During the robust estimation process, the M-estimation method is introduced to reduce the influence of gross errors on the estimation of the correction values of the horizontal displacement parameters in the east-west or north-south directions. The weighted iterative method is used to solve the correction values of the historical horizontal displacement parameters in the east-west or north-south directions. The iterative calculation formula is as follows:

[0092]

[0093] In the formula, are the correction value matrix of the historical horizontal displacement parameters in the east-west or north-south directions, the correction number matrix of the historical horizontal displacement in the east-west or north-south directions, the cofactor matrix of the historical horizontal displacement in the east-west or north-south directions, and the equivalent weight matrix of the historical horizontal displacement observations in the east-west or north-south directions in the k-th iterative calculation, respectively. w(k - 1) is the weight factor matrix in the (k - 1)-th iterative calculation.

[0094] The weight factor matrix w(k) in the k-th iterative calculation can be obtained by constructing the Huber weight function:

[0095]

[0096] In the formula, w mm (k) is the weight factor obtained at the (m, m) position in the weight factor matrix w(k) in the k-th iterative calculation; that is to say, w mm (k) is the weight factor for iterative calculation at the pixel (m, m) position, and w(k) is the weight factor matrix composed of the results of iterative calculation for each pixel. is the standardized residual, v m (k) is the corresponding residual, and σ m (k) is the corresponding standard deviation. c is a constant, generally selected as 2.0 - 2.5.

[0097] During the iterative process, when the difference between the correction values of the historical horizontal displacement parameters in the east-west or north-south directions obtained from two consecutive calculations satisfies the preset iterative convergence criterion, the iterative process terminates. Substituting the obtained by iterative calculation into Equation (4), the correction number matrix can be obtained. Since the elements in By adding the horizontal displacement calculated in Step 2 to the correction number calculated in Step 3.3, a more accurate horizontal displacement can be obtained, that is, the adjusted historical horizontal displacement in the east-west or north-south directions where, is the horizontal displacement in the east-west or north-south directions between the DOM of the (n - 2)-th period and the DOM of the (n - 1)-th period, is the corresponding adjusted east-west or north-south horizontal displacement. Since it is for historical DOM data, it is also called the adjusted historical east-west or north-south horizontal displacement.

[0098] Step 3.4: When new DOM data is obtained, jointly perform adjustment calculation on the new DOM data and historical DOM data to establish a sequential adjustment error equation.

[0099] When there is new DOM data, use the method in Step 2 to extract the east-west or north-south horizontal displacement observation values between the new DOM data and the previous DOM data (referred to as historical DOM data), and jointly perform adjustment calculation on the new observation values and the existing observation values. The expression of the established sequential adjustment error equation is as follows:

[0100]

[0101] In the formula, is the correction matrix of the new and historical east-west or north-south horizontal displacements; is the coefficient matrix related to the historical east-west or north-south horizontal displacement; A is the coefficient matrix related to the new east-west or north-south horizontal displacement; is the matrix composed of the difference between the approximate values of the new east-west or north-south horizontal displacement and the historical east-west or north-south horizontal displacement. is the matrix composed of the estimated values of the parameter corrections of the sequential adjustment. Among them, is the estimated value matrix of the parameter corrections of the new horizontal displacement for sequential adjustment based on the new horizontal displacement. is the estimated value matrix of the parameter corrections for updating the historical horizontal displacement when jointly performing sequential adjustment on the adjusted (after k iterations) historical horizontal displacement and the new horizontal displacement after adding the new horizontal displacement data. In the parameter subscripts, EW / NS both refer to the east-west or north-south horizontal displacement, and in some explanations of parameter meanings, they are not specifically distinguished. When performing east-west or north-south horizontal displacement adjustment, substitute the corresponding east-west or north-south direction calculated values.

[0102] Step 3.5: Combine the sequential adjustment error equation and construct a sequential adjustment normal equation according to the Bayesian least squares criterion.

[0103] Combine the sequential adjustment error equation in Equation (9), construct the objective function as shown in Equation (10) below according to the Bayesian least squares criterion, and further construct the corresponding sequential adjustment normal equation (similar to the construction method of the normal equation (6)).

[0104]

[0105] Among them, min means taking the minimum value of the right side of the equal sign. It is the weight matrix for the newly added horizontal displacement in the east-west or north-south direction. It is the correction value matrix of the historical horizontal displacement parameters obtained after k iterations by performing adjustment only using the historical horizontal displacement. And It is the estimated value matrix of the correction number of the parameters for updating the historical horizontal displacement by sequential adjustment after adding the new horizontal displacement.

[0106] Step 3.6: Conduct robust estimation based on the sequential adjustment normal equations, and obtain the adjusted horizontal displacements in the east-west and north-south directions based on the estimation results.

[0107] Combined with the error equations, according to the Bayesian least squares criterion, construct an objective function to solve for the estimated value of the correction value of the parameters for sequential adjustment and the cofactor matrix, so as to realize the dynamic update of the horizontal displacements in the east-west and north-south directions in time series after adjustment. Specifically, by solving formulas (9) and (10), the estimated value matrix of the correction value of the parameters for sequential adjustment can be obtained as:

[0108]

[0109] After obtaining the initial sequential solution, in order to resist the influence of outliers in the newly added horizontal displacement observations in the east-west or north-south direction on parameter estimation, the M-estimation method is further introduced to effectively suppress the influence of outliers on parameter estimation. Solve the correction value matrix of the parameters of the newly added horizontal displacement in the east-west or north-south direction through the weight selection iteration method, and determine the equivalent weight matrix of the newly added horizontal displacement in the east-west or north-south direction through iterative solution.

[0110]

[0111] Among them, the intermediate parameters in the calculation process

[0112]

[0113] When performing joint adjustment calculation for the newly added horizontal displacement, the cofactor matrix of the sequential adjustment used is as follows:

[0114]

[0115] Among them,

[0116]

[0117] So far, by applying formulas (12) to (15), the dynamic update of the horizontal displacements in the east-west and north-south directions in time series after adjustment can be realized. This process is based on the recursive principle of sequential adjustment, and the estimated value matrix Ensure the global consistency and dynamic optimality of the solution results, and substitute the updated parameter correction value estimates into Equation (9) to obtain the correction matrix of the new and historical horizontal displacements in the east-west or north-south directions. Correct the extracted horizontal displacements in the east-west or north-south directions to obtain more accurate horizontal displacements in the east-west or north-south directions. Finally, the adjusted new and historical horizontal displacements in the east-west and north-south directions are collectively referred to as the adjusted horizontal displacements in the east-west and north-south directions.

[0118] When new horizontal displacement data is added next time, (including) the horizontal displacement data updated this time will be used as historical horizontal displacement data. For sequential adjustment in combination with the newly added horizontal displacement next time, with each addition of horizontal displacement data, the estimated matrix of the correction values of the historical horizontal displacement parameters obtained by Equation (12) and the estimated matrix of the correction values of the newly added horizontal displacement parameters will be dynamically updated.

[0119] As the mining face in the mining area continues to advance, the flight cycle of the unmanned aerial vehicle and the amount of observation data gradually increase. These large amounts of observation data provide valuable opportunities for revealing the deformation laws in the mining area, but at the same time pose new challenges for quickly and efficiently obtaining more accurate deformation data. With the increase in the data collection frequency of the unmanned aerial vehicle, DOM data gradually accumulates, and the horizontal displacement data between every two periods of DOM doubles, resulting in a significant increase in the amount of horizontal displacement data. When the traditional overall adjustment calculation method faces a large amount of historical horizontal displacement data and newly added horizontal displacement data, the calculation amount increases sharply, resulting in a significant reduction in the overall adjustment efficiency, thereby affecting the real-time dynamic update of horizontal displacement data. By using the adjustment method proposed in Step 3 of this application, efficient dynamic update of horizontal displacement data in the east-west and north-south directions can be achieved.

[0120] Step 4: Based on the adjusted horizontal displacements in the east-west and north-south directions, perform plane position correction on each period of DEM data to obtain the corrected DEM data for each period.

[0121] This application uses the adjusted horizontal displacement data to correct the plane positions of each period of deformed DEM. Therefore, when performing data preprocessing in Step 1, it is necessary to crop each period of DOM and DEM data to keep their coverage ranges the same. The DEM data at the observation time T0 is called the pre-deformation DEM data, and the subsequent T1, T2,..., T nThe DEM data at these n observation times is called the deformed DEM data. For each period of DEM data in the deformed DEM data, taking the planar positions of each feature point in the pre-deformation DEM data as the reference, and combining the horizontal displacements in the east-west and north-south directions after adjustment, the planar positions of the corresponding feature points in each period of the deformed DEM data are corrected to be consistent with the planar positions of the corresponding feature points in the pre-deformation DEM data. The ultimate goal of this step is to make the planar positions of each corresponding feature point before and after deformation in the same vertical direction.

[0122] Furthermore, using the inverse distance weighted interpolation algorithm, interpolation is performed on each period of the deformed DEM data at the planar positions of the corrected feature points to obtain the corrected elevation values of each feature point, constituting each period of the corrected DEM data.

[0123] The specific correction process is as Figure 3 shown, Figure 3 in which point A and point A' respectively represent the same feature point before and after horizontal displacement, and D EW and D NS respectively represent the horizontal displacements in the east-west and north-south directions from point A to point A'. In order to accurately obtain the displacement amount (elevation change) of point A in the vertical direction, first, point A' is corrected to the position of point A in the deformed DEM using the horizontal displacement information in the east-west and north-south directions after adjustment. Then, taking the planar position of point A as the reference, and combining the horizontal displacements in the east-west and north-south directions after adjustment, through the inverse distance weighted interpolation algorithm, the elevation value stored in the pixel of point A' is assigned to the pixel of point A, so that the elevation values of the pixels where point A was before deformation and point A' is after deformation are in the same vertical direction. The same method is used for other feature points to construct each period of the deformed DEM corresponding to each feature point in the pre-deformation DEM, so that each period of the deformed DEM after correction and each corresponding feature point in the pre-deformation DEM are in the same vertical direction.

[0124] Step 5: Perform differential calculation on each period of the corrected DEM data and the reference DEM data before deformation to obtain the time-series subsidence data of the mining area.

[0125] During the correction process, the reference DEM data before deformation usually refers to the DEM data obtained by observation at time T0, and the deformed DEM data is usually the DEM data obtained by observation at times T1, T2, …, T n times. Taking Figure 3 point A in it as an example, assuming that when observing at time T1, point A has a horizontal displacement and moves to the horizontal position of point A', where D EW and D NSIt is the horizontal displacement in the east-west and north-south directions when point A moves to point A' during the time period from T0 to T1. Using the horizontal displacements in the east-west and north-south directions after adjustment in step 3, the A' point of the DEM at T1 is corrected back to the position corresponding to point A at T0, and through the inverse distance weighted algorithm, the elevation at the position of point A' at T1 is assigned to the position corresponding to point A in the DEM at T1 after correction, so as to perform differential calculation through two-phase DEMs subsequently.

[0126] Specifically, through the correction process in step 4, it has been ensured that each homologous ground feature point before and after deformation is in the same vertical direction. At this time, the corrected DEM data of each period is directly subtracted from the reference DEM data before deformation, and the DEM differential data of each period can be obtained. Then, the DEM differential data of each period is arranged according to the time sequence of each period, and the time-sequence subsidence data of the mining area can be formed.

[0127] Step 6: Construct the time-sequence three-dimensional deformation data of the mining area based on the time-sequence subsidence data of the mining area and the horizontal displacements in the east-west and north-south directions after adjustment, so as to realize the dynamic monitoring of the subsidence of the mining area.

[0128] Taking the reference DEM before deformation as a reference, the DEMs of each period after correction are differentiated from the reference DEM to obtain the time-sequence subsidence data of the mining area at each observation moment after deformation relative to before deformation. The horizontal displacements in the east-west and north-south directions after adjustment are arranged according to the time sequence, and the time-sequence horizontal displacement data after adjustment can be obtained. The time-sequence subsidence data of the mining area combined with the time-sequence horizontal displacement data after adjustment can obtain the time-sequence three-dimensional deformation data of the mining area. Based on the time-sequence three-dimensional deformation data of the mining area, a time-sequence three-dimensional deformation model of the subsidence basin in the mining area is constructed, thereby realizing the dynamic monitoring of the subsidence of the mining area.

[0129] The method of the present application solves the problem that the influence of horizontal displacement on the calculation of subsidence value is not considered in the existing subsidence basin construction technology, and avoids the problem of error caused by horizontal displacement in the calculation of subsidence value in areas with large terrain undulations. And it realizes the efficient dynamic update of the time-sequence observation data of horizontal displacement, and improves the efficiency and accuracy of the dynamic monitoring of the subsidence of the mining area.

[0130] In an exemplary embodiment, the present application also provides a device for dynamic monitoring of subsidence in a mining area, including:

[0131] A data preprocessing module, configured to obtain the time-sequence data of UAV images and lidar point clouds observed in the same period in the mining area and perform data preprocessing to obtain DOM and DEM data of each period; where the same period refers to the same observation moment; each period refers to each observation moment;

[0132] A horizontal displacement extraction module, configured to extract the horizontal displacements in the east-west and north-south directions between any two-phase DOM data based on the DOM data of each period by using the normalized cross-correlation matching algorithm;

[0133] The horizontal displacement sequential adjustment module is used to dynamically update the historical east-west and north-south horizontal displacements and the newly added east-west and north-south horizontal displacements according to the sequential least squares adjustment principle for the newly added DOM data, and obtain the adjusted east-west and north-south horizontal displacements;

[0134] The plane position correction module is used to perform plane position correction on each period of DEM data based on the adjusted east-west and north-south horizontal displacements to obtain the corrected DEM data for each period;

[0135] The differential calculation module is used to perform differential calculation on the corrected DEM data for each period and the reference DEM data before deformation to obtain the time-series subsidence data of the mining area;

[0136] The mining area subsidence dynamic monitoring module is used to construct the time-series three-dimensional deformation data of the mining area based on the time-series subsidence data of the mining area and the adjusted east-west and north-south horizontal displacements to realize the dynamic monitoring of the mining area subsidence.

[0137] Compared with the existing methods, in this application, the time-series horizontal displacement data of the mining area is extracted from the DOM at times T0, T1,..., T (n-1) by the normalized cross-correlation matching algorithm, and for the newly added DOM at time T n the sequential least squares adjustment principle is adopted to dynamically update the historical horizontal displacement and the newly added horizontal displacement data, thus significantly improving the efficiency of the adjustment calculation and ensuring the efficient processing and accurate solution under large-scale observation data. During the adjustment process, the M-estimation method is introduced, and the influence of gross errors on the parameter estimation results is reduced by the weight selection iteration method, thereby effectively improving the accuracy and reliability of the horizontal displacement calculation. Further, based on the adjusted time-series horizontal displacement data, the plane position correction of each period of DEM after deformation is performed, so as to ensure that the plane positions of each ground feature point in each period of DEM before and after deformation are consistent. On this basis, the corrected DEM data for each period and the reference DEM before deformation are used for differential calculation to construct a time-series three-dimensional deformation model of the mining area subsidence basin, effectively weakening the influence of horizontal displacement on the calculation of vertical subsidence values and significantly improving the accuracy of subsidence monitoring.

[0138] Of course, the specific embodiments described in this application are merely illustrative of the inventive concept of this application. Those skilled in the art to which this application pertains can make various modifications or supplements to the described specific embodiments or use similar methods for substitution. For example, generating DOM can be replaced by other software; the improved progressive encryption triangular mesh filtering algorithm for removing non-ground points such as vegetation and buildings in the point cloud and then extracting ground points can be replaced by other algorithms; the algorithm for extracting the horizontal displacements in the east-west and north-south directions between any two periods of DOM can be replaced by other algorithms; the weight function of the weight selection and iteration method in sequential adjustment can be replaced by other weight functions; the inverse distance weighted interpolation algorithm in the process of correcting the deformed DEM can be replaced by other algorithms. The above substitution methods will not exceed the framework of the method proposed in this application, will not deviate from the inventive concept of this application, or exceed the scope protected by this application.

[0139] In an exemplary embodiment, this application further provides a computer device, which can be a server or a terminal. The computer device includes a processor, a memory, an input / output interface, and a communication interface. Among them, the processor, the memory, and the input / output interface are connected through a system bus, and the communication interface is connected to the system bus through the input / output interface. Among them, the processor of the computer device is used to provide computing and control capabilities. The memory of the computer device includes a non-volatile storage medium and an internal memory. The non-volatile storage medium stores an operating system, a computer program, and a database. The internal memory provides an environment for the operation of the operating system and the computer program in the non-volatile storage medium. The input / output interface of the computer device is used to exchange information between the processor and external devices. The communication interface of the computer device is used to communicate with external terminals through a network connection. The computer program, when executed by the processor, implements the described method for dynamic monitoring of mining area subsidence.

[0140] In an exemplary embodiment, this application further provides a computer-readable storage medium, on which a computer program is stored, and the computer program, when executed by the processor, implements the described method for dynamic monitoring of mining area subsidence.

[0141] In an exemplary embodiment, this application further provides a computer program product, including a computer program, and the computer program, when executed by the processor, implements the described method for dynamic monitoring of mining area subsidence.

[0142] Those of ordinary skill in the art can understand that all or part of the processes in the above-described embodiment methods can be completed by hardware related to computer program instructions. The computer program can be stored in a non-volatile computer-readable storage medium. When the computer program is executed, it can include the processes of the above-described method embodiments. Among them, any reference to a memory or other medium provided in the embodiments of the present application can include at least one of non-volatile and volatile memories. Non-volatile memory can include read-only memory (ROM), magnetic tape, floppy disk, flash memory, optical memory, high-density embedded non-volatile memory, resistive random access memory (ReRAM), magnetoresistive random access memory (MRAM), ferroelectric random access memory (FRAM), phase change memory (PCM), graphene memory, etc. Volatile memory can include random access memory (RAM) or external cache memory, etc. By way of illustration and not limitation, RAM can be in various forms, such as static random access memory (SRAM) or dynamic random access memory (DRAM), etc.

[0143] It should be noted that the information (including but not limited to user device information, user personal information, etc.) and data (including but not limited to data for analysis, stored data, displayed data, etc.) involved in this application are all information and data that have been authorized by the user or fully authorized by all parties. And the collection, use, and processing of relevant data need to comply with relevant regulations.

[0144] The technical features of the above embodiments can be combined arbitrarily. For the sake of brevity of description, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, it should be considered as the scope described in this specification.

[0145] In this article, specific examples are used to elaborate on the principles and implementation manners of the present application. The description of the above embodiments is only used to help understand the method and its core idea of the present application; at the same time, for those of ordinary skill in the art, according to the idea of the present application, there will be changes in the specific implementation manners and application scopes. In summary, the content of this specification should not be construed as a limitation to the present application.

Claims

1. A dynamic monitoring method for mining area subsidence, characterized in that, Including: Obtaining the time-series data of UAV images and LiDAR point clouds for synchronous observations in the mining area and performing data preprocessing to obtain DOM and DEM data for each period; Where "synchronous" refers to the same observation moment; the DOM and DEM data for each period refer to the DOM and DEM data at each observation moment; Based on the DOM data for each period, using the normalized cross-correlation matching algorithm to extract the horizontal displacements in the east-west and north-south directions between any two periods of DOM data; For the newly added DOM data, dynamically updating the historical horizontal displacements in the east-west and north-south directions and the newly added horizontal displacements in the east-west and north-south directions according to the principle of sequential least squares adjustment to obtain the adjusted horizontal displacements in the east-west and north-south directions; Based on the adjusted horizontal displacements in the east-west and north-south directions, performing plane position correction on the DEM data for each period to obtain the corrected DEM data for each period; Performing differential calculation on the corrected DEM data for each period and the reference DEM data before deformation to obtain the time-series subsidence data of the mining area; Constructing the time-series three-dimensional deformation data of the mining area based on the time-series subsidence data of the mining area and the adjusted horizontal displacements in the east-west and north-south directions to realize the dynamic monitoring of the subsidence in the mining area.

2. The dynamic monitoring method for mining area subsidence according to claim 1, characterized in that The obtaining of the time-series data of UAV images and LiDAR point clouds for synchronous observations in the mining area and performing data preprocessing to obtain DOM and DEM data for each period specifically includes: Using Pix4D software to process the UAV image data at each observation moment in the time-series data of UAV images to generate the corresponding DOM data for each period; For the LiDAR point cloud data at each observation moment in the time-series data of LiDAR point clouds, adopting an improved progressive densification triangulation filtering algorithm to remove the non-ground points in the LiDAR point cloud data, extract the ground points, and using the inverse distance weighted interpolation algorithm to interpolate the ground points to generate the corresponding DEM data.

3. The dynamic monitoring method for mining area subsidence according to claim 2, wherein, The extracting of the horizontal displacements in the east-west and north-south directions between any two periods of DOM data based on the DOM data for each period using the normalized cross-correlation matching algorithm specifically includes: Converting the selected two periods of DOM data into the corresponding grayscale images; Taking the grayscale image with an earlier observation moment as the reference image and the grayscale image with a later observation moment as the target image, and constructing the corresponding matching windows in the reference image and the target image; Calculating the similarity of the image blocks within the two matching windows through the normalized cross-correlation matching algorithm to obtain the normalized cross-correlation coefficient; Adopting the winner-takes-all strategy, taking the peak point of the normalized cross-correlation coefficient as the matching point, and the corresponding left-right disparity and up-down disparity are the relative displacements in the row and column directions respectively, and using the sinc interpolation function to interpolate the relative displacements in the row and column directions respectively to obtain the sub-pixel level relative displacements in the row and column directions; Combining the projection coordinate system parameters and the pixel resolution, converting the sub-pixel level relative displacements in the row and column directions into the horizontal displacements in the east-west and north-south directions in the projection coordinate system.

4. The dynamic monitoring method for mining area subsidence according to claim 3, wherein The dynamically updating of the historical horizontal displacements in the east-west and north-south directions and the newly added horizontal displacements in the east-west and north-south directions according to the principle of sequential least squares adjustment for the newly added DOM data to obtain the adjusted horizontal displacements in the east-west and north-south directions specifically includes: An error equation is established based on the relationship between the horizontal displacements in the east-west and north-south directions between any two periods of DOM data. According to the least squares principle, the corresponding normal equation is obtained from the error equation. Robust estimation is carried out based on the normal equation, and the historical horizontal displacements in the east-west and north-south directions are obtained based on the estimation results. When new DOM data is obtained, the new DOM data and the historical DOM data are jointly adjusted and calculated, and a sequential adjustment error equation is established. Combined with the sequential adjustment error equation, according to the Bayesian least squares criterion, a sequential adjustment normal equation is constructed. Robust estimation is carried out based on the sequential adjustment normal equation, and the adjusted horizontal displacements in the east-west and north-south directions are obtained based on the estimation results.

5. The dynamic monitoring method for mining area subsidence according to claim 4, wherein, Based on the adjusted horizontal displacements in the east-west and north-south directions, plane position correction is performed on each period of DEM data to obtain the corrected DEM data for each period, which specifically includes: For each period of DEM data, taking the plane positions of each feature point in the DEM data before deformation as the reference, and combining the adjusted horizontal displacements in the east-west and north-south directions, the plane positions of the corresponding feature points in each period of DEM data after deformation are corrected so that the plane positions of each corresponding feature point before and after deformation are in the same vertical direction. Using the inverse distance weighted interpolation algorithm, interpolation is performed on each period of DEM data after deformation at the plane positions of the corrected feature points to obtain the corrected elevation values of each feature point, forming the corrected DEM data for each period.

6. The dynamic monitoring method for mining area subsidence according to claim 5, characterized in that The differential calculation of the corrected DEM data for each period and the reference DEM data before deformation to obtain the time-series subsidence data of the mining area specifically includes: Directly subtracting the corrected DEM data for each period from the reference DEM data before deformation to obtain the DEM differential data for each period, and arranging them according to the time series of each period to form the time-series subsidence data of the mining area.

7. A dynamic monitoring device for mining area subsidence, characterized in that, It includes: A data preprocessing module for obtaining the time-series data of UAV images and lidar point clouds observed in the same period for the mining area and performing data preprocessing to obtain the DOM and DEM data for each period; where the same period refers to the same observation moment; each period refers to each observation moment. A horizontal displacement extraction module for extracting the horizontal displacements in the east-west and north-south directions between any two periods of DOM data based on the DOM data for each period using the normalized cross-correlation matching algorithm. A horizontal displacement sequential adjustment module for dynamically updating the historical horizontal displacements in the east-west and north-south directions and the new horizontal displacements in the east-west and north-south directions according to the sequential least squares adjustment principle for the new DOM data to obtain the adjusted horizontal displacements in the east-west and north-south directions. A plane position correction module for performing plane position correction on each period of DEM data based on the adjusted horizontal displacements in the east-west and north-south directions to obtain the corrected DEM data for each period. A differential calculation module for performing differential calculation on the corrected DEM data for each period and the reference DEM data before deformation to obtain the time-series subsidence data of the mining area. A mining area subsidence dynamic monitoring module for constructing the time-series three-dimensional deformation data of the mining area based on the time-series subsidence data of the mining area and the adjusted horizontal displacements in the east-west and north-south directions to realize the dynamic monitoring of the mining area subsidence.

8. A computer device, comprising: A memory, a processor, and a computer program stored on the memory and executable on the processor, wherein the processor executes the computer program to implement the mine subsidence dynamic monitoring method according to any one of claims 1 to 6.

9. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by the processor, it implements the mine subsidence dynamic monitoring method according to any one of claims 1 to 6.

10. A computer program product, comprising a computer program, characterized in that, When the computer program is executed by the processor, it implements the mine subsidence dynamic monitoring method according to any one of claims 1 to 6.

Images