A mine subsidence dynamic monitoring method, device, equipment, medium and product

Abstract

The application discloses a kind of mining area subsidence dynamic monitoring method, device, equipment, medium and product, involve mining area exploitation subsidence monitoring field.First, the same period observation unmanned aerial vehicle image and laser radar point cloud time series data are carried out data preprocessing, obtain each period DOM and DEM data;Based on each period DOM data, the horizontal displacement between any two period DOM data in east-west and north-south direction is extracted using normalized cross-correlation matching algorithm;For new DOM data, according to the sequential least square adjustment principle, the horizontal displacement in east-west and north-south direction of history and new is dynamically updated, to obtain the horizontal displacement in east-west and north-south direction after adjustment;Based on the horizontal displacement in east-west and north-south direction after adjustment, the plane position of each period DEM data is corrected;The corrected each period DEM data and the reference DEM data before deformation are differentially calculated, to obtain the time series subsidence data of mining area, realize the dynamic monitoring of mining area subsidence, can improve the precision and efficiency of mining area subsidence monitoring.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of mining area subsidence monitoring, and in particular to a mining area subsidence dynamic monitoring method, device, equipment, medium and product. BACKGROUND

[0002] Surface subsidence, building destruction, surface cracks and landslides caused by coal mining and other geological disasters seriously threaten the ecological environment and the safety of residents in the mining area. Therefore, obtaining high-precision three-dimensional deformation data of the mining area is of great significance to protect the ecological environment and the safety of the residential area in the mining area. 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 the disadvantages of large workload, high cost and harsh operating conditions. In comparison, unmanned aerial 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, the existing unmanned aerial photogrammetry and LiDAR technology mainly use DEM (Digital Elevation Model) difference method to construct the subsidence basin, which only reflects the ground elevation change under the same geographical spatial coordinates and ignores the influence of horizontal displacement on vertical deformation. SUMMARY

[0003] In view of the problems pointed out in the background section, the present application provides a mining area subsidence dynamic monitoring method, device, equipment, medium and product, which considers the influence of horizontal displacement on vertical deformation in the subsidence value calculation process, thereby improving the mining area subsidence monitoring precision and efficiency.

[0004] To achieve the above-mentioned purpose, the present application provides the following solutions.

[0005] In a first aspect, the present application provides a mining area subsidence dynamic monitoring method, comprising:

[0006] Obtaining the unmanned aerial image and laser radar point cloud time series data for the same period observation of the mining area and performing data preprocessing to obtain the DOM and DEM data of each period; wherein the same period refers to the same observation time; the DOM and DEM data of each period refer to the DOM and DEM data of each observation time;

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

[0008] For the newly added DOM data, the horizontal displacements in the east-west and north-south directions of the history and the newly added east-west and north-south directions are dynamically updated according to the sequential least square adjustment principle, 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 of each period DEM data is corrected to obtain the corrected DEM data of each period;

[0010] The corrected DEM data of 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 mining area subsidence.

[0012] Optionally, the unmanned aerial vehicle image and the laser radar point cloud time series data of the same period of the mining area are obtained and preprocessed to obtain the DOM and DEM data of each period, specifically including:

[0013] The Pix4D software is used to process the unmanned aerial vehicle image data of each observation time in the unmanned aerial vehicle image time series data to generate the corresponding DOM data of each period;

[0014] For the laser radar point cloud data of each observation time in the laser radar point cloud time series data, the improved progressive encryption triangulation network filtering algorithm is used to remove the non-ground points in the laser radar point cloud data, the ground points are extracted, and the inverse distance weighted interpolation algorithm is used to interpolate the ground points to generate the corresponding DEM data.

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

[0016] The selected any two periods of DOM data are converted into corresponding gray images;

[0017] The gray image of the earlier observation time is taken as the reference image, and the gray image of the later observation time is taken as the target image, and the corresponding matching window is constructed in the reference image and the target image;

[0018] The similarity of the image blocks in the two matching windows is calculated by the normalized cross-correlation matching algorithm to obtain the normalized cross-correlation coefficient;

[0019] The peak point of the normalized cross-correlation coefficient is taken as a matching point by using the winner-takes-all strategy, and the corresponding left and right parallax and the up and down parallax are the relative displacements in the row and column directions, respectively, and the relative displacements in the row and column directions are interpolated by using a sinc interpolation function, respectively, to obtain the relative displacements in the row and column directions at a sub-pixel level.

[0020] The relative displacements in the row and column directions at a sub-pixel level are converted into horizontal displacements in the east-west and north-south directions in the projection coordinate system in combination with the projection coordinate system parameters and the pixel resolution.

[0021] Optionally, 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 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, and specifically comprising:

[0022] 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;

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

[0024] Robust estimation is performed based on the normal equation, and the historical horizontal displacements in the east-west and north-south directions are obtained based on the estimation result;

[0025] When the new DOM data is obtained, the new DOM data and the historical DOM data are jointly adjusted to establish a sequential adjustment error equation;

[0026] In combination with the sequential adjustment error equation, a sequential adjustment normal equation is constructed according to the Bayesian least squares criterion;

[0027] Robust estimation is performed 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 result.

[0028] Optionally, the horizontal displacements in the east-west and north-south directions after adjustment are used to correct the plane positions of each period of DEM data to obtain the corrected each period of DEM data, specifically comprising:

[0029] For each period of DEM data, the plane positions of each feature point in the pre-deformation DEM data are taken as a reference, and the plane positions of the same feature points in each period of post-deformation DEM data are corrected in combination with the horizontal displacements in the east-west and north-south directions after adjustment, so that the plane positions of each same feature point before and after deformation are in the same vertical direction;

[0030] An inverse distance weighted interpolation algorithm is used to interpolate each period of post-deformation DEM data at the plane positions of the corrected feature points to obtain the corrected elevation values of each feature point, thereby constructing each period of corrected DEM data.

[0031] Optionally, the corrected DEM data of each period is differentially calculated with the reference DEM data before deformation to obtain time-series subsidence data of the mining area, specifically comprising:

[0032] The corrected DEM data of each period is directly subtracted from the reference DEM data before deformation to obtain DEM differential data of each period, which is arranged in time series to form time-series subsidence data of the mining area.

[0033] In a second aspect, the present application provides a mining area subsidence dynamic monitoring device, comprising:

[0034] A data preprocessing module is configured to obtain unmanned aerial vehicle images and laser radar point cloud time-series data observed in the same period of the mining area and perform data preprocessing to obtain DOM and DEM data of each period; wherein the same period refers to the same observation time; and each period refers to each observation time.

[0035] A horizontal displacement extraction module is configured to extract horizontal displacement 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 a normalized cross-correlation matching algorithm.

[0036] A horizontal displacement sequential adjustment module is configured to dynamically update the horizontal displacement in the east-west and north-south directions of the history and the newly added east-west and north-south directions according to the principle of sequential least squares adjustment to obtain adjusted horizontal displacement in the east-west and north-south directions for the newly added DOM data.

[0037] A plane position correction module is configured to correct the plane position of the DEM data of each period based on the adjusted horizontal displacement in the east-west and north-south directions to obtain corrected DEM data of each period.

[0038] A differential calculation module is configured to differentially calculate the corrected DEM data of each period with the reference DEM data before deformation to obtain time-series subsidence data of the mining area.

[0039] A mining area subsidence dynamic monitoring module is configured to construct 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 displacement in the east-west and north-south directions to realize dynamic monitoring of the mining area subsidence.

[0040] In a third aspect, the present application provides a computer device, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the computer program to realize the mining area subsidence dynamic monitoring method.

[0041] In a fourth aspect, the present application provides a computer readable storage medium having a computer program stored thereon, wherein the computer program is executed by a processor to realize the mining area subsidence dynamic monitoring method.

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

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

[0044] In the method, device, equipment, medium and product for monitoring the dynamic subsidence of a mining area provided in the present application, on one hand, the influence of the horizontal displacement in the east-west and south-north directions on the vertical deformation is considered in the calculation process of the time-series subsidence data of the mining area, the horizontal displacement is used to correct the plane position of the DEM data, and the dynamic monitoring of the subsidence of the mining area is performed based on the corrected DEM data, which can significantly improve the subsidence monitoring accuracy of the mining area; on the other hand, the efficient and dynamic updating of the historical and newly added horizontal displacement in the east-west and south-north directions can be realized for the historical and newly added DOM data, thereby improving the subsidence monitoring efficiency of the mining area. BRIEF DESCRIPTION OF DRAWINGS

[0045] In order to more clearly illustrate the technical solutions in the embodiments of the present application or the prior art, the drawings needed in the embodiments will be briefly introduced as follows. Obviously, the drawings in the following description are only some embodiments of the present application, and other drawings can be obtained by those skilled in the art without any creative effort.

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

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

[0048] Figure 3 FIG. 3 is a schematic diagram of the process of correcting the plane position of the DEM data. DETAILED DESCRIPTION

[0049] The technical solutions in the embodiments of the present application will be described clearly and completely with reference to the drawings in the embodiments of the present application. Obviously, the described embodiments are only some of the embodiments of the present application, but not all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without any creative effort are within the scope of protection of the present application.

[0050] The present application provides a method, device, equipment, medium and product for monitoring the dynamic subsidence of a mining area, which aims to establish a dynamic monitoring system for the subsidence of a mining area based on the time-series DOM / DEM data of an unmanned aerial vehicle and considering the horizontal displacement, and the influence of the horizontal displacement on the vertical deformation is considered in the calculation process of the subsidence value, thereby improving the subsidence monitoring accuracy and efficiency of the mining area.

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

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

[0053] Step 1: Obtain the unmanned aerial image and laser radar point cloud time series data for the same period observation of the mining area and perform data preprocessing to obtain the DOM and DEM data of each period.

[0054] Figure 2 The main process of monitoring the dynamic subsidence of a mining area using the method of the present application is shown. Referring to Figure 2 , the present application uses unmanned aerial photogrammetry and LiDAR technology to observe the target mining area and obtain the unmanned aerial image and laser radar point cloud time series data for the same period observation. Among them, the unmanned aerial image time series data includes unmanned aerial image data of T0, T1, …, T (n-1) n observation moments (also referred to as moments). Each observation moment is also referred to as each period. Correspondingly, the laser radar point cloud time series data also includes laser radar point cloud data of T0, T1, …, T (n-1) n observation moments. Wherein n is usually an integer greater than 3. For example, T (n-1) The unmanned aerial image data observed at T (n-1) The laser radar point cloud data observed at T (n-1) The purpose of the data preprocessing stage is to generate DOM (Digital Orthophoto Map) and DEM data of T0, T1, …, T (n-1) n moments based on the unmanned aerial image and laser radar point cloud time series data.

[0055] On the one hand, for the unmanned aerial image time series data, the Pix4D software is used to process the unmanned aerial image data of T0, T1, …, T (n-1) n observation moments in the unmanned aerial image time series data to generate the corresponding DOM data of each period, that is, to generate T0, T1, …, T (n-1) n period DOM data.

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

[0057] Further, the generated T0, T1, …, T (n-1) The DOM data and DEM data (hereinafter also directly referred to as DOM and DEM for the convenience of description) at the observation time are geographically registered to ensure that the coordinate system and spatial resolution of each period of DOM and DEM are consistent. The generated DOM and DEM usually do not completely cover the same area, and the common area of each period of DOM and DEM is cropped to obtain each period of DOM and DEM data with consistent spatial coverage and consistent planar position of each same named ground feature point.

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

[0059] Based on T0, T1, …, T (n-1) n observation times, the normalized cross-correlation matching algorithm is used to extract the east-west and north-south horizontal displacement between any two observation times (for example, 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) n-2 and n-1 represent two observation times. Referring to Figure 2 , the step 2 specifically includes the following steps 2.1 to 2.5.

[0060] Step 2.1: Convert the selected any two periods of DOM data into corresponding gray scale 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 gray image; R, G, and B are respectively the red, green, and blue channel pixel values of the DOM.

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

[0065] Select the gray image corresponding to the DOM at an earlier acquisition time as the reference image I1, and the gray image corresponding to the DOM at a later acquisition time as the target image I2. In the reference image I1, a matching window with a size of h x h is constructed with the position of the to-be-matched pixel as the center. In the target image I2, a matching window with the same size is constructed with the position of the target pixel as the center.

[0066] Step 2.3: Perform similarity calculation on the image blocks in the two matching windows by the normalized cross-correlation matching algorithm to obtain the normalized cross-correlation coefficient.

[0067] Perform similarity calculation on the image blocks in the two matching windows by the normalized cross-correlation matching algorithm, and the 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 of the to-be-matched pixel (x, y) in the reference image I1, and the mean value of all pixels in the matching window is W p is the matching window centered on the to-be-matched pixel. I2(x+u, y+v) is the pixel value of the center pixel (i.e., the target pixel) of the matching window in the target image I2, and the mean value of the corresponding matching window is (x+u, y+v) is the corresponding target pixel coordinate in the target image I2. u and v are the relative displacements in the row and column directions in the traversal process, respectively. The peak point of the normalized cross-correlation function is taken as the row and column direction relative displacements u' and v' of the matching point.

[0070] Step 2.4: Take the peak point of the normalized cross-correlation coefficient as the matching point using the winner-takes-all strategy, and the corresponding left and right disparities and the top and bottom disparities are the relative displacements in the row and column directions, respectively. The sinc interpolation function is used to perform interpolation processing on the relative displacements in the row and column directions to obtain the relative displacements in the row and column directions at the sub-pixel level.

[0071] The peak point of the normalized cross-correlation coefficient NCC(u, v) is taken as the matching point by using the winner-takes-all strategy, and the corresponding left and right disparities and the up and down disparities 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 value of the corresponding normalized cross-correlation function often locates at the sub-pixel position. In order to accurately locate the peak value, the sinc interpolation function is used to perform interpolation processing on the relative displacements u, v in the row and column directions, respectively, so as to realize the sub-pixel level disparity estimation. The relative displacement in the row direction at the sub-pixel level is denoted as u', and the relative displacement in the column direction at the sub-pixel level is denoted as v'.

[0072] Step 2.5: Combined with the projection coordinate system parameters and the pixel resolution, the relative displacements in the row and column directions at the sub-pixel level are converted into the horizontal displacements in the east-west and north-south directions under the projection coordinate system.

[0073] Combined with the projection coordinate system parameters and the pixel resolution, the relative displacements in the row and column directions at the sub-pixel level calculated in pixels are converted into the horizontal displacements in the east-west and north-south directions under 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, the horizontal displacement in the east-west direction between the n-2 and n-1 periods of DOM data is denoted as The horizontal displacement in the north-south direction between the n-2 and n-1 periods of DOM data is denoted as

[0076] Step 3: For the newly added DOM data, 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 according to the sequential least squares adjustment principle to obtain the adjusted horizontal displacements in the east-west and north-south directions.

[0077] As mentioned earlier, the present application has pre-collected the DOM data at T0, T1, …, T (n-1) These n observation time points are also called historical DOM data. Then for the DOM data at the subsequent n observation time point T n , the present application calls it the newly added DOM data. Correspondingly, and are called the historical horizontal displacements in the east-west and north-south directions, respectively; and The horizontal displacements in the east and north directions are called as new east and north displacements. Next, the historical east and north displacements and the new east and north displacements are dynamically updated according to the sequential least square adjustment principle, and the M-estimation method is introduced to reduce the influence of gross errors on the parameter estimation results, so as to improve the accuracy and reliability of the parameter solution. See Figure 2 The step 3 specifically includes the following steps 3.1 to 3.6.

[0078] Step 3.1: Based on the relationship between the east and north horizontal displacements between any two periods of DOM data, an error equation is established.

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

[0080]

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

[0082]

[0083] In the formula, Δx0 and Δx1 represent the historical east or north horizontal displacement parameter correction values, Δx0 and Δx1 are elements in the historical east or north horizontal displacement parameter correction value matrix Δx, and Δx0 and Δx1 represent the historical east or north horizontal displacement parameter correction values between the 0th and (n-1)th periods of DOM.

[0084] is a coefficient matrix of (n(n-1) / 2)*(n-1), which is related to the east or north horizontal displacement observations extracted from any two periods of DOM. Each row of the matrix is generated according to the following rules. For the east or north horizontal displacement extracted from the i,j period of DOM, it satisfies 0≤i When i=0, the jth column element of the corresponding row in the matrix ​​​​​​The i-th column element of the corresponding row in the matrix is -1, the j-th column element is 1, and the rest of the elements are 0. For example, when i = 0, j = 1, the corresponding correction number is Matrix The first column element of the corresponding row in the matrix is 1, and the rest of the elements are 0. When i = 0, j = 2, the corresponding correction number is Matrix The second column element of the corresponding row in the matrix is 1, and the rest of the elements are 0. When i = 1, j = 2, the corresponding correction number is Matrix The first column element of the corresponding row in the matrix is -1, and the second column element is 1, and the rest of the elements are 0. And so on.

[0085] For a better understanding, the meaning of each element in the entire error equation (5) is illustrated below. For example, assume that the east-west direction horizontal displacement between the 3rd period DOM and the 4th period DOM is extracted as Because of the existence of errors, the observed values need to be adjusted. The east-west direction horizontal displacement after adjustment is Then According to the principle of indirect adjustment, the east-west direction horizontal displacement calculated at the observation time T 0,1 ,T 0,2 ,T 0,3 ,…,T 0,(n-1) is taken as the parameter When adjusting, the approximate value D' EW of the parameter is generally taken, for example, The approximate value of is When adjusting the north-south direction horizontal displacement, the principle is the same as above.

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

[0087] According to the least squares principle, the normal equation of formula (4), (5) is as follows:

[0088]

[0089] wherein, is the covariance matrix of historical east-west or north-south direction horizontal displacement, is the weight matrix of historical east-west or north-south direction horizontal displacement observation values.

[0090] Step 3.3: Robust estimation is carried out based on the normal equation, and historical east-west and north-south direction horizontal displacement is obtained based on the estimation result.

[0091] In the process of robust estimation, M-estimation method is introduced to reduce the influence of gross errors on the estimation of the horizontal displacement parameters in the east or north direction, and the historical horizontal displacement parameters in the east or north direction are solved by the selected weight iteration method. The iteration formula is as follows:

[0092]

[0093] In the formula, w are the historical horizontal displacement parameter correction matrix in the east or north direction, the historical horizontal displacement correction number matrix in the east or north direction, the historical horizontal displacement correlation factor matrix in the east or north direction, and the equivalent weight matrix of the historical horizontal displacement observation in the east or north direction respectively in the kth iteration calculation. w(k-1) is the weight factor matrix in the (k-1)th iteration calculation.

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

[0095]

[0096] In the formula, w mm (k) is the weight factor in the (m, m) position of the weight factor matrix w(k) obtained in the kth iteration calculation; that is, w mm (k) is the weight factor in the (m, m) position of the weight factor matrix w(k) obtained in the kth iteration calculation; that is, w is the standardized residual, v m (k) is the corresponding residual, σ m (k) is the corresponding standard deviation. c is a constant, generally selected as 2.0-2.5.

[0097] In the iteration process, the difference between the historical horizontal displacement parameters in the east or north direction obtained by continuous two times of calculation satisfies the preset iteration convergence criterion, and the iteration process is terminated. The historical horizontal displacement parameters in the east or north direction obtained by iteration calculation are substituted into formula (4), and the correction number matrix can be solved. Since the elements in the correction number matrix are all less than 1, the horizontal displacement calculated in step 2 is added to the correction number calculated in step 3.3, and the more accurate horizontal displacement, that is, the adjusted historical horizontal displacement in the east or north direction is obtained. In the formula, Δx The corresponding horizontal displacement in the east-west or north-south direction after adjustment is also called the historical horizontal displacement in the east-west or north-south direction after adjustment since it is based on historical DOM data.

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

[0099] When the new DOM data is obtained, the method in step 2 is used to extract the east-west or north-south horizontal displacement observation value between the new DOM data and the previous DOM data (referred to as historical DOM data), and the new observation value is combined with the existing observation value to perform adjustment calculation. 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 displacement; 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 new east-west or north-south horizontal displacement and the approximate value of the historical east-west or north-south horizontal displacement. is the matrix composed of the estimated value of the parameter correction of the sequential adjustment. Among them, is the estimated value matrix of the new horizontal displacement parameter correction value based on the sequential adjustment of the new horizontal displacement. is the estimated value matrix of the parameter correction value of the historical horizontal displacement updated when the historical horizontal displacement and the new horizontal displacement are combined to perform sequential adjustment after the adjustment (after k iterations) after the addition of the new horizontal displacement data. The parameter subscript EW / NS refers to the east-west or north-south horizontal displacement, which is not specifically distinguished in some parameter meaning explanations. When performing east-west or north-south horizontal displacement adjustment, the corresponding east-west or north-south calculation value is substituted.

[0102] Step 3.5: According to the Bayesian least squares criterion, the sequential adjustment equation is constructed based on the sequential adjustment error equation.

[0103] According to the Bayesian least squares criterion, the objective function shown in formula (10) is constructed based on the sequential adjustment error equation of formula (9), and the corresponding sequential adjustment equation is further constructed (similar to the construction method of equation (6)).

[0104]

[0105] Among them, min represents the minimum value on the right side of the equal sign. The weight matrix of the newly added east-west or north-south horizontal displacement. is the parameter correction value matrix of the historical horizontal displacement obtained after k iterations of adjustment only using the historical horizontal displacement. is the parameter correction value matrix of the historical horizontal displacement obtained after k iterations of adjustment only using the historical horizontal displacement.

[0106] Step 3.6: Robust estimation based on the sequential adjustment equation, and the east-west and north-south horizontal displacement after adjustment based on the estimation results.

[0107] Combined with the error equation, according to the Bayesian least squares criterion, the objective function is constructed to obtain the estimated value of the parameter correction value of the sequential adjustment and the cofactor matrix, so as to realize the dynamic update of the east-west and north-south horizontal displacement after adjustment. Specifically, by solving formula (9) and (10), the estimated value matrix of the parameter correction value of the sequential adjustment can be obtained as follows:

[0108]

[0109] After obtaining the initial sequential solution, in order to resist the influence of gross errors in the newly added east-west or north-south horizontal displacement observation value on parameter estimation, the M-estimation method is further introduced to effectively suppress the influence of gross errors on parameter estimation. The weight iteration method is used to solve the correction value matrix of the newly added east-west or north-south horizontal displacement parameter and the equivalent weight matrix of the newly added east-west or north-south horizontal displacement is determined by iteration Finally, according to the matrix inversion theorem, the recursive calculation formula of the sequential solution is derived as follows:

[0110]

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

[0112]

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

[0114]

[0115] Among them,

[0116]

[0117] At this point, by applying formula (12)~(15), the dynamic update of the east-west and north-south horizontal displacement after adjustment can be realized. This process is based on the recursive principle of sequential adjustment, and the estimated value matrix of the parameter correction value is updated in real time after each new observation data Ensure the global consistency and dynamic optimality of the solution, and put the updated parameter correction value estimate into equation (9) to obtain the correction matrix of the new and historical east or north horizontal displacement Correct the extracted east or north horizontal displacement to obtain more accurate east or north horizontal displacement. Finally, the adjusted east and north horizontal displacement of the new and historical is collectively referred to as the adjusted east and north horizontal displacement.

[0118] Next time the horizontal displacement data is added, the horizontal displacement data updated this time will be included as historical horizontal displacement data. To perform sequential adjustment on the next added horizontal displacement, with each added horizontal displacement data, the estimated value matrix of the historical horizontal displacement parameter correction value calculated by equation (12) And the estimated value matrix of the new horizontal displacement parameter correction value will be dynamically updated.

[0119] As the mining workface in the mining area continues to advance, the flight period of the unmanned aerial vehicle and the amount of observation data gradually increases. These large amounts of observation data provide valuable opportunities to reveal the deformation law of the mining area, but at the same time, they also pose new challenges to quickly and efficiently obtain more accurate deformation data. With the increase in the frequency of unmanned aerial vehicle data collection, 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. The traditional overall adjustment calculation method faces a large amount of historical horizontal displacement data and newly added horizontal displacement data, and the calculation amount increases dramatically, resulting in a significant decrease in the efficiency of overall adjustment, which in turn affects the real-time dynamic update of horizontal displacement data. However, by using the adjustment method proposed in step 3 of the present application, efficient dynamic updating of east and north horizontal displacement data can be achieved.

[0120] Step 4: Based on the adjusted east and north horizontal displacement, the plane position of each period of DEM data is corrected to obtain the corrected DEM data of each period.

[0121] The present application uses the adjusted horizontal displacement data to correct the plane position of each period of DEM after deformation, so when performing data preprocessing in step 1, each period of DOM and DEM data needs to be cropped to maintain the same coverage. The DEM data at the T0 observation time is referred to as the pre-deformation DEM data, and the subsequent T1, T2, …, T nThe DEM data collected at these n observation times are called the post-deformation DEM data. For each period of the post-deformation DEM data, using the planar positions of all feature points in the pre-deformation DEM data as a benchmark, and combining the adjusted east-west and north-south horizontal displacements, the planar positions of the corresponding feature points in each period of the post-deformation DEM data are corrected to ensure consistency with the planar positions of the corresponding feature points in the pre-deformation DEM data. The ultimate goal of this step is to ensure that the planar positions of all corresponding feature points before and after deformation are in the same vertical direction.

[0122] Furthermore, an inverse distance weighted interpolation algorithm is used to interpolate the deformed DEM data for each period at the corrected planar location of the ground features to obtain the corrected elevation values ​​of each ground feature, thus forming the corrected DEM data for each period.

[0123] The specific correction process is as follows: Figure 3 As shown, Figure 3 Points A and A′ represent the same ground feature before and after the horizontal displacement, respectively. EW and D NS These represent the east-west and north-south horizontal displacements from point A to point A′, respectively. To accurately determine the vertical displacement (elevation change) of point A, point A′ is first corrected to its position in the deformed DEM using the adjusted east-west and north-south horizontal displacement information. Next, using the planar position of point A as a reference, and combining the adjusted east-west and north-south horizontal displacements, an inverse distance weighted interpolation algorithm is used to assign the elevation value stored in the A′ pixel to the A pixel, ensuring that the elevation values ​​of the pixels containing point A before and after deformation are in the same vertical direction. Other feature points are processed using the same method, constructing post-deformation DEMs corresponding to each feature point in the pre-deformation DEM, ensuring that the corrected post-deformation DEMs are in the same vertical direction as the corresponding feature points in the pre-deformation DEM.

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

[0125] During the calibration process, the baseline DEM data before deformation usually refers to the DEM data observed and acquired at time T0, while the DEM data after deformation is usually T1, T2, ..., T n DEM data acquired through continuous monitoring. Figure 3 Taking point A as an example, suppose that at time T1, A undergoes a horizontal displacement and moves to the horizontal position of point A′, where D EW and D NSThe east-west and north-south horizontal displacement is corrected from A' to A in the T0-T1 period A point, and the corrected position of the A point in the T1 time is obtained by using the east-west and north-south horizontal displacement after the adjustment in step 3. The height of the A' point in the T1 time is assigned to the position corresponding to the A point in the T1 time in the corrected DEM, so as to obtain the differential data of the DEM in the subsequent two periods.

[0126] Specifically, the correction process in step 4 has ensured that the same named ground points before and after deformation are in the same vertical direction. At this time, the differential data of the DEM in each period can be obtained by directly subtracting the corrected DEM data from the reference DEM data before deformation. Then, the time series subsidence data of the mining area can be obtained by arranging the differential data of the DEM in each period according to the time sequence.

[0127] Step 6: Based on the time series subsidence data of the mining area and the east-west and north-south horizontal displacement after adjustment, the time series three-dimensional deformation data of the mining area is constructed, and the dynamic monitoring of the subsidence of the mining area is realized.

[0128] With the reference of the reference DEM before deformation, the differential data of the corrected DEM and the reference DEM is obtained, and the time series subsidence data of the mining area after deformation is obtained. The east-west and north-south horizontal displacement after adjustment can be arranged according to the time sequence, and the time series data of the horizontal displacement after adjustment can be obtained. The time series three-dimensional deformation data of the mining area can be obtained by combining the time series subsidence data of the mining area with the time series data of the horizontal displacement after adjustment. Based on the time series three-dimensional deformation data of the mining area, the time series three-dimensional deformation model of the subsidence basin of the mining area is constructed, and the dynamic monitoring of the subsidence of the mining area is realized.

[0129] The method solves the problem that the influence of the horizontal displacement on the calculation of the subsidence value is not considered in the existing subsidence basin construction technology, avoids the error problem caused by the horizontal displacement in the calculation of the subsidence value in the area with large terrain undulations, and realizes the efficient dynamic updating of the time series observation data of the horizontal displacement, thereby improving the dynamic monitoring efficiency and accuracy of the subsidence of the mining area.

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

[0131] A data preprocessing module is configured to obtain the unmanned aerial vehicle image and the laser radar point cloud time series data observed in the same period of the mining area and perform data preprocessing to obtain the DOM and DEM data in each period. The same period refers to the same observation time, and each period refers to each observation time.

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

[0133] horizontal displacement in the east and south directions and the newly added horizontal displacement in the east and south directions according to the sequential least square adjustment principle, to obtain the horizontal displacement in the east and south directions after adjustment;

[0134] The plane position correction module is configured to correct the plane positions of the DEM data of each period based on the horizontal displacement in the east and south directions after adjustment, to obtain the corrected DEM data of each period.

[0135] The difference calculation module is configured to perform difference calculation on the corrected DEM data of each period and the reference DEM data before deformation, to obtain time-series subsidence data of the mining area.

[0136] The mining area subsidence dynamic monitoring module is configured to construct time-series three-dimensional deformation data of the mining area based on the time-series subsidence data of the mining area and the horizontal displacement in the east and south directions after adjustment, to realize dynamic monitoring of the subsidence of the mining area.

[0137] Compared with the existing method, the normalized cross-correlation matching algorithm is used to extract the horizontal displacement time-series data of the mining area from T0, T1, …, Tn in the present application. (n-1) The sequential least square adjustment principle is used to dynamically update the historical horizontal displacement and the newly added horizontal displacement data, thereby significantly improving the efficiency of adjustment calculation and ensuring efficient processing and accurate calculation under large-scale observation data. n In the adjustment process, the M-estimation method is introduced, and the influence of gross errors on the parameter estimation results is reduced through the selected weight iteration method, thereby effectively improving the accuracy and reliability of the horizontal displacement calculation. Further, the plane position correction module is configured to correct the plane positions of the DEM after deformation based on the horizontal displacement time-series data after adjustment, thereby ensuring that the plane positions of each feature point in the DEM before and after deformation are consistent. On this basis, the difference calculation module is configured to perform difference calculation on the corrected DEM of each period and the reference DEM before deformation, to construct a time-series three-dimensional deformation model of the subsidence basin of the mining area, thereby effectively reducing 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 the present application are only illustrative of the inventive concept of the present application. Those skilled in the art of the present application can make various modifications or supplements to the described specific embodiments or replace them with similar ways. For example, the DOM can be generated by other software; the improved progressive TIN filtering algorithm for extracting ground points by removing vegetation, buildings and other non-ground points in the point cloud can be replaced by other algorithms; the algorithm for extracting east-west and north-south horizontal displacement between any two periods of DOM can be replaced by other algorithms; the weight function of the weight function selected in the sequential adjustment can be replaced by other weight functions; the inverse distance weighted interpolation algorithm in the post-deformation DEM correction process can be replaced by other algorithms. The above replacement methods do not exceed the framework of the method of the present application, do not deviate from the inventive concept of the present application, or exceed the scope of protection of the present application.

[0139] In an exemplary embodiment, the present application also 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. Wherein 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. Wherein 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 operating system and the computer program in the non-volatile storage medium to run. 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 network connection. The computer program is executed by the processor to realize the mining area subsidence dynamic monitoring method.

[0140] In an exemplary embodiment, the present application also provides a computer readable storage medium having a computer program stored thereon, wherein the computer program is executed by a processor to realize the mining area subsidence dynamic monitoring method.

[0141] In an exemplary embodiment, the present application also provides a computer program product comprising a computer program, wherein the computer program is executed by a processor to realize the mining area subsidence dynamic monitoring method.

[0142] Those skilled in the art can understand that all or part of the processes in the above-mentioned embodiment methods can be completed by computer program instruction related hardware. The computer program can be stored in a non-volatile computer readable storage medium. When the computer program is executed, the processes of the above-mentioned embodiment methods can be included. Any reference to memory or other medium in the embodiments provided by the present application can include at least one of non-volatile and volatile memory. Non-volatile memory can include read only memory (ROM), magnetic tape, floppy disk, flash memory, optical storage, high-density embedded non-volatile memory, resistive memory (ReRAM), magnetoresistive random access memory (MRAM), ferroelectric memory (FRAM), phase change memory (PCM), graphene memory, etc. Volatile memory can include random access memory (RAM) or external cache memory, etc. As an illustration but 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 the present application are all information and data authorized by the user or authorized by all parties, and the collection, use and processing of related data need to comply with relevant regulations.

[0144] The technical features of the above embodiments can be combined arbitrarily. To make the description concise, all possible combinations of the technical features in the above embodiments are not described, but as long as the combinations of the technical features do not exist contradictory, they should be considered as the scope of the present application.

[0145] The principles and implementation modes of the present application are described by applying specific examples in this paper. The above description of the embodiments is only used to help understand the method and its core idea of the present application. For those skilled in the art, according to the idea of the present application, the specific implementation mode and application range can be changed. In conclusion, the content of the present application should not be understood as a limitation.

Claims

1. A method of dynamic monitoring of mining subsidence, characterized in that, The method comprises the following steps: acquiring unmanned aerial vehicle image and laser radar point cloud time series data for the same period observation of a mining area and performing data preprocessing to obtain DOM and DEM data of each period; wherein the same period refers to the same observation time; the DOM and DEM data of each period refer to the DOM and DEM data of each observation time; based on the DOM data of each period, using a normalized cross-correlation matching algorithm to extract the east-west and north-south horizontal displacement between any two periods of DOM data; for the new DOM data, according to the sequential least squares adjustment principle, dynamically updating the historical east-west and north-south horizontal displacement and the new east-west and north-south horizontal displacement to obtain the adjusted east-west and north-south horizontal displacement; The T0, T1, …, Tn-2, Tn-1 are collected in advance (n-1) The DOM data of the n observation moments are called historical DOM data; the DOM data of the subsequent n observation moment Tn+1 is called new DOM data n The historical east-west and north-south horizontal displacements are called historical east-west and north-south horizontal displacements The new east-west and north-south horizontal displacements are called new east-west and north-south horizontal displacements (n-2),(n-1) n-2 and n-1 represent the two observation moments​​​ the step of dynamically updating the historical east-west and north-south horizontal displacement and the new east-west and north-south horizontal displacement according to the sequential least squares adjustment principle to obtain the adjusted east-west and north-south horizontal displacement for the new DOM data specifically comprises: based on the relationship between the east-west and north-south horizontal displacement between any two periods of DOM data, an error equation is established; Based on the acquired T 0,1 ,T 0,2 ,…,T (n-2),(n-1) The relationship between the east-west and north-south horizontal displacement observations is established as follows: the extended form of the error equation is as follows: wherein is a matrix of historical east or north horizontal displacement parameter correction values, wherein the element represents the historical east or north horizontal displacement parameter correction value between the 0th and the n-1th DOMs; is a matrix of the difference between historical east or north horizontal displacement and the east or north horizontal displacement approximation, wherein the element represents the difference between the historical east or north horizontal displacement and the east or north horizontal displacement approximation between the n-2th and the n-1th DOMs; is a matrix of correction numbers of historical east or north horizontal displacement, wherein the element represents the correction number of historical east or north horizontal displacement between the n-2th and the n-1th DOMs; It is a coefficient matrix of (n(n-1) / 2)*(n-1), which is related to the east-west or north-south horizontal displacement observations extracted from any two periods of DOM; the matrix The generation rules for each row are as follows: For the east-west or north-south horizontal displacement extracted by the DOM in the i and j periods, 0 ≤ i < j ≤ (n-1); when i = 0, the matrix... The element in the j-th column of the corresponding row is 1, and all other 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 rest are 0; according to the least squares principle, the corresponding normal equation is obtained from the error equation; according to the least squares principle, the normal equation is obtained from formulas (4) and (5) as follows: wherein, is a matrix of covariants of the historical east or north horizontal displacement, is a weight matrix of the historical east or north horizontal displacement observations; robust estimation is performed based on the normal equation, and the historical east-west and north-south horizontal displacement is obtained based on the estimation result; in the robust estimation process, the M-estimation method is introduced, and the correction value of the historical east-west or north-south horizontal displacement parameter is solved through the selected weight iteration method, and the iteration calculation formula is as follows: wherein respectively, are the historical east or north horizontal displacement parameter correction value matrix, the historical east or north horizontal displacement correction number matrix, the historical east or north horizontal displacement cofactor matrix, and the equivalent weight matrix of the historical east or north horizontal displacement observation value of the kth iteration calculation; w(k-1) is the weight factor matrix of the (k-1)th iteration calculation; The iteration process is terminated when the difference between the historical east-west or north-south horizontal displacement parameter correction values obtained by two consecutive calculations satisfies the preset iteration convergence criterion; and the historical east-west or north-south horizontal displacement parameter correction value obtained by the iteration calculation is taken as the final historical east-west or north-south horizontal displacement parameter correction value. The correction number matrix is obtained by substituting formula (4) Since the element in the matrix The adjusted historical east-west or north-south horizontal displacement is obtained by adding the calculated horizontal displacement and the calculated correction number wherein, is the east-west or north-south horizontal displacement between the n-2 period DOM and the n-1 period DOM; is the corresponding adjusted east-west or north-south horizontal displacement, and since it is directed at historical DOM data, it is also referred to as the adjusted historical east-west or north-south horizontal displacement. when the new DOM data is acquired, the new DOM data and the historical DOM data are jointly adjusted to establish a sequential adjustment error equation; combined with the sequential adjustment error equation, according to the Bayesian least squares criterion, a sequential adjustment normal equation is constructed; robust estimation is performed based on the sequential adjustment normal equation, and the adjusted east-west and north-south horizontal displacement is obtained based on the estimation result; based on the adjusted east-west and north-south horizontal displacement, the plane position of the DEM data of each period is corrected to obtain the corrected DEM data of each period; the corrected DEM data of each period is differentially calculated with the reference DEM data before deformation to obtain the time series subsidence 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 displacement, the time series three-dimensional deformation data of the mining area is constructed to realize dynamic monitoring of the subsidence of the mining area.

2. The method of dynamic monitoring of mining area subsidence according to claim 1, characterised by, the step of acquiring unmanned aerial vehicle image and laser radar point cloud time series data for the same period observation of a mining area and performing data preprocessing to obtain DOM and DEM data of each period specifically comprises: using Pix4D software to process the unmanned aerial vehicle image data of each observation time in the unmanned aerial vehicle image time series data to generate the corresponding DOM data of each period; for the laser radar point cloud data of each observation time in the laser radar point cloud time series data, using an improved progressive encryption triangulation network filtering algorithm to remove the non-ground points in the laser radar point cloud data, extracting the ground points, and using an inverse distance weighted interpolation algorithm to interpolate the ground points to generate the corresponding DEM data.

3. The method of dynamic monitoring of mining area subsidence according to claim 2, characterised by, The horizontal displacement in the east-west and north-south directions between any two periods of DOM data is extracted based on the DOM data of each period by using a normalized cross-correlation matching algorithm, and specifically includes: The selected DOM data of any two periods is converted into corresponding gray-scale images; The gray-scale image of an earlier observation time is taken as a reference image, and the gray-scale image of a later observation time is taken as a target image, and a corresponding matching window is constructed in the reference image and the target image; Similarity calculation is performed on the image blocks in the two matching windows by using a normalized cross-correlation matching algorithm to obtain a normalized cross-correlation coefficient; A winner-takes-all strategy is adopted, the peak point of the normalized cross-correlation coefficient is taken as a matching point, and the corresponding left and right parallaxes and the up and down parallaxes are the relative displacements in the row and column directions, respectively, and a sinc interpolation function is used to perform interpolation processing on the relative displacements in the row and column directions to obtain the relative displacements in the row and column directions at a sub-pixel level; The relative displacements in the row and column directions at a sub-pixel level are converted into horizontal displacements in the east-west and north-south directions in a projection coordinate system in combination with the projection coordinate system parameters and the pixel resolution.

4. The method of dynamic monitoring of mining area subsidence according to claim 3, characterised by, The horizontal displacements in the east-west and north-south directions after adjustment are used to correct the plane positions of the DEM data of each period to obtain the corrected DEM data of each period, and specifically includes: For each DEM data, the plane positions of the same feature points in the DEM data after deformation are corrected based on the plane positions of the feature points in the DEM data before deformation and in combination with the horizontal displacements in the east-west and north-south directions after adjustment, so that the plane positions of the same feature points before and after deformation are in the same vertical direction; An inverse distance weighted interpolation algorithm is used to interpolate the DEM data after deformation at the plane positions of the corrected feature points to obtain the corrected elevation values of the feature points and form the corrected DEM data of each period.

5. The method of dynamic monitoring of mining area subsidence according to claim 4, characterised by, The corrected DEM data of each period and the reference DEM data before deformation are subjected to difference calculation to obtain time-series subsidence data of the mining area, and specifically includes: The corrected DEM data of each period and the reference DEM data before deformation are directly subtracted to obtain DEM difference data of each period, which is arranged in time series to form the time-series subsidence data of the mining area.

6. A device for dynamic monitoring of mining subsidence, characterized in that, The mining area subsidence dynamic monitoring method in any one of claims 1 to 5 is implemented; and the mining area subsidence dynamic monitoring device includes: A data preprocessing module is configured to obtain unmanned aerial vehicle image and laser radar point cloud time-series data observed at the same period of the mining area and perform data preprocessing to obtain DEM and DOM data of each period; wherein the same period refers to the same observation time; and each period refers to each observation time; A horizontal displacement extraction module is configured to extract horizontal displacement 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 a normalized cross-correlation matching algorithm; A horizontal displacement sequential adjustment module is configured to dynamically update the horizontal displacement in the east-west and north-south directions before and after adjustment based on the principle of sequential least squares adjustment for new DOM data and new horizontal displacement in the east-west and north-south directions to obtain the horizontal displacement in the east-west and north-south directions after adjustment. A plane position correction module is configured to correct the plane positions of the DEM data of each period based on the horizontal displacements in the east-west and north-south directions after the adjustment, to obtain the corrected DEM data of each period; A difference calculation module is configured to calculate the differences between the corrected DEM data of each period and the reference DEM data before the deformation, to obtain the time-series subsidence data of the mining area; A mining area subsidence dynamic monitoring module is 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 horizontal displacements in the east-west and north-south directions after the adjustment, to realize the dynamic monitoring of the subsidence of the mining area.

7. A computer device comprising: A memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that the processor executes the computer program to implement the mining area subsidence dynamic monitoring method of any one of claims 1 to 5.

8. A computer-readable storage medium having stored thereon a computer program, characterized in that, The computer program is executed by the processor to implement the mining area subsidence dynamic monitoring method of any one of claims 1 to 5.

9. A computer program product comprising a computer program, characterized in that, The computer program is executed by the processor to implement the mining area subsidence dynamic monitoring method of any one of claims 1 to 5.

Images