Method for monitoring surface subsidence large-deformation surface domain

Abstract

The invention relates to a surface subsidence large deformation area monitoring method. The method comprises the following steps: obtaining subsidence data of a research area in n-1 periods according to an SBAS-InSAR technology; using a dynamic settlement rate interpolation method to obtain a theoretical settlement value of a monitoring point position corresponding to each period of the SBAS-InSAR; according to the SBAS-InSAR deformation original result in each period and the theoretical settlement amount of the monitoring point position in the corresponding period, calculating to obtain a correction coefficient of the SBAS-InSAR deformation original result in each period at the monitoring point position in the corresponding period; and calculating the correction coefficient of the SBAS-InSAR deformation original result in each period at the monitoring point position in the corresponding period by adopting a Kriging interpolation method to obtain the correction coefficient of the SBAS-InSAR deformation original result in each period at any point position in the corresponding period, and finally obtaining a corrected mining area dynamic settlement result. According to the method, the problems of phase ambiguity and the like of the InSAR technology in large deformation monitoring are effectively solved, and the monitoring accuracy of the large-scale subsidence area of the mining area is improved.

Description

Technical Field

[0001] The invention belongs to the field of mining subsidence monitoring, and in particular relates to a method for monitoring large deformation areas of surface subsidence. Background Art

[0002] Coal mines occupy an important position in China's energy structure and are one of the pillars of the country's economic development. With the rapid development of the economy, the demand for coal continues to rise, and coal mining activities are becoming increasingly frequent. However, the coal mining process inevitably causes surface subsidence, which has a significant impact on the surrounding ecological environment and infrastructure. Ground subsidence may cause structural damage to buildings, road collapse, changes in groundwater levels, and water pollution, posing a potential threat to the lives and safety of surrounding residents. In addition, geological disasters caused by subsidence may also cause huge economic losses and social problems. Effective monitoring methods are urgently needed to prevent and control these risks.

[0003] Traditional mining subsidence monitoring methods mainly rely on ground monitoring points and manual observation. Although these methods provide relatively accurate subsidence data in local areas, due to the sparse distribution of monitoring points and the frequency limitation of manual observation, it is often difficult to fully and timely reflect the subsidence status of the entire mining area, resulting in a lag in the judgment of subsidence trends, which affects the effectiveness of emergency decision-making.

[0004] In recent years, with the continuous development of remote sensing technology, synthetic aperture radar interferometry (InSAR) technology has gradually become an important means of monitoring the overall subsidence of mining areas. InSAR technology uses radar images obtained by satellites to analyze the phase difference of ground reflection signals, and can efficiently and accurately obtain surface deformation information in a large area. Compared with traditional ground monitoring methods, InSAR has the advantages of wide area coverage, high spatial resolution and high temporal frequency, enabling it to capture the overall subsidence trend of mining areas in a timely manner. However, when faced with large deformations, the phase information of InSAR may experience "phase ambiguity", that is, the phase value cannot accurately reflect the true deformation of the surface. This is because InSAR technology relies on the measurement of phase difference. When the deformation exceeds the phase period of the radar wave, it is difficult to accurately calculate, resulting in distortion of the monitoring results. Summary of the invention

[0005] In order to solve the above problems, the present invention provides a large deformation area monitoring method for surface subsidence. This method combines the settlement data of known ground monitoring points with the SBAS-InSAR results through interpolation technology to generate more accurate mining area subsidence monitoring results. This method has the following advantages: first, it can realize dynamic and continuous monitoring of mining area subsidence and capture deformation changes in a timely manner; second, the SBAS-InSAR results are globally corrected through the monitoring data of ground monitoring points, thereby improving the accuracy of the monitoring results; finally, this method reduces the time and cost of manual monitoring, improves monitoring efficiency and response speed, and provides a scientific basis and strong support for safe production and environmental protection of coal mines. The application of this innovative technology will provide important guarantees for the sustainable development and safety management of the coal mining industry.

[0006] To achieve the above objectives, this application provides the following technical solutions:

[0007] A method for monitoring surface subsidence and large deformation, characterized by comprising the following steps:

[0008] S1) Obtain n SAR images during the mining period of the study area, and obtain the original SBAS-InSAR deformation results of the study area in n-1 periods based on the SBAS-InSAR technology;

[0009] S2) Collect the settlement data of known monitoring points in the study area, analyze and process the settlement data of the monitoring points using the dynamic settlement rate interpolation method, and obtain the theoretical settlement values ​​of the corresponding monitoring point positions at each period of SBAS-InSAR;

[0010] S3) calculating the settlement deviation of the monitoring point position at each period of SBAS-InSAR according to the original SBAS-InSAR deformation results at each period and the theoretical settlement amount of the monitoring point position at the corresponding period, and using it as the correction coefficient of the original SBAS-InSAR deformation results at each period of SBAS-InSAR at the monitoring point position at the corresponding period;

[0011] S4) using the Kriging interpolation method to solve the correction coefficient of the original SBAS-InSAR deformation result at the monitoring point position in each period to obtain the correction coefficient of the original SBAS-InSAR deformation result at any point position in the corresponding period;

[0012] S5) obtaining the corrected SBAS-InSAR deformation result of any point according to the original SBAS-InSAR deformation result of each period and the correction coefficient of any point position in the corresponding period, thereby obtaining the corrected dynamic settlement result of the mining area.

[0013] Furthermore, the step S2) specifically includes:

[0014] Obtain the location information of known monitoring points in the study area and the settlement monitoring value S at each time node t;

[0015] Calculate the sedimentation rate of adjacent monitoring rounds at the monitoring point. The calculation expression is:

[0016]

[0017] Where i represents the monitoring rounds of each monitoring point, t i represents the time node corresponding to the i-th round of monitoring, Indicates that the monitoring point is at t i+1 -t i The sedimentation rate over time, It means that at t = t i The settlement monitoring value of the monitoring point at that time;

[0018] According to the dynamic sedimentation rate of each monitoring round, the time point T at which n long-time series SAR image data are acquired is calculated. x When , the theoretical settlement value of each monitoring point is calculated as:

[0019]

[0020] 1 of them <x≤n,t i ≤T x <t i+1 , Indicates time node t = T x The theoretical settlement value corresponding to the monitoring point at that time.

[0021] Furthermore, the step S3) specifically includes:

[0022] Calculate the theoretical settlement amount of the corresponding monitoring points of SBAS-InSAR in each period The calculation expression is:

[0023]

[0024] According to the original SBAS-InSAR deformation results at different periods and the theoretical settlement of the monitoring point in the corresponding period Calculate the settlement deviation of the monitoring point position at each period of SBAS-InSAR, and the calculation expression is:

[0025]

[0026] in It represents the correction coefficient of the original SBAS-InSAR deformation results at the monitoring point position in each period.

[0027] Furthermore, the step S4) specifically includes:

[0028] The correction coefficients of the original SBAS-InSAR deformation results at the corresponding monitoring point positions in each period are analyzed. Assuming that there are k monitoring points in the study area, the calculation expression of the correction coefficient of any point in the study area in each period is:

[0029]

[0030] in It means t = T x The correction coefficient for any position (x,y) at the time, θ j represents the weight coefficient, It means t = T x The correction coefficient corresponding to the monitoring point j at time .

[0031] Furthermore, the step S4) specifically includes:

[0032] Calculate the weight coefficients, assuming that the correction coefficients at any point (x, y) in the study area have the same expected value c and variance σ 2 , that is, it holds true for any point:

[0033]

[0034] Calculate the distance D and the corresponding semi-variance γ(D) between the known monitoring points. The distance is obtained based on the coordinate value:

[0035]

[0036] Where D 1,2 represents the straight-line distance between monitoring point 1 and monitoring point 2, Δx and Δy represent the differences in the horizontal and vertical coordinates of the two monitoring points respectively;

[0037] The semivariogram function is defined as:

[0038]

[0039] in represents the semivariance between monitoring point 1 and monitoring point 2;

[0040] All the obtained (D,γ) are fitted with the spherical model, and the function is:

[0041]

[0042] Where C0 represents the nugget constant, C0+C represents the base value, C represents the arch height, and a represents the range;

[0043] For any point in the study area, firstly calculate the distance D and the corresponding semivariance C from all known monitoring points according to the above method;

[0044] The goal of the Kriging interpolation method is to find a set of weight coefficients that can satisfy the minimum difference between the estimated value and the true value of the correction coefficient at any point in the study area, that is:

[0045]

[0046] in, represents the estimated value of the correction coefficient at any point in the study area, ε(x i ,y i ) represents the true value of the correction coefficient at any point in the study area;

[0047] At the same time, it satisfies the unbiased estimation conditions:

[0048]

[0049] By constructing a short covariance matrix and a covariance vector, combined with unbiased constraints, a linear equation system is established and the optimal weight coefficient vector can be obtained by solving it.

[0050] Furthermore, the step S5) specifically includes:

[0051] The original SBAS-InSAR deformation results of the study area in n-1 periods are corrected, and the calculation expression is:

[0052]

[0053] in Represents the corrected SBAS-InSAR deformation results.

[0054] Beneficial effects:

[0055] (1) The present invention uses the settlement data of known ground monitoring points and SBAS-InSAR results to construct a correction model and correct the SBAS-InSAR results of each period, which effectively overcomes the problems such as "phase ambiguity" existing in InSAR technology in large deformation monitoring and improves the monitoring accuracy of large-scale settlement areas in mining areas.

[0056] (2) The present invention adopts the dynamic settlement rate interpolation method to further study the local settlement law of the monitoring point, and can interpolate the settlement rate of the monitoring point in different time periods to accurately reflect the dynamic changes of the settlement in the mining area.

[0057] (3) The Kriging interpolation method in the present invention assigns different weights to each monitoring point to minimize the interpolation error. This method ensures that the interpolation results are both highly accurate and avoid systematic deviations, making the correction model and SBAS-InSAR results more reliable. It can effectively capture the spatial variation trend of the subsidence data, so that the method can show different subsidence laws in different areas of the mining area, adapt to the different effects of geological conditions and mining activities on ground deformation, thereby improving the interpolation accuracy. BRIEF DESCRIPTION OF THE DRAWINGS

[0058] Figure 1 It is a flow chart of the surface subsidence large deformation area monitoring method of the present invention;

[0059] Figure 2 This is a flow chart of the Kriging interpolation method of the present invention;

[0060] Figure 3 This is a diagram of the experimental results of the study area of ​​the present invention. DETAILED DESCRIPTION

[0061] The technical solution provided by the present invention is more clearly described below in conjunction with the embodiments and drawings, but the scope of protection claimed by the present invention is not limited to the following embodiments.

[0062] In order to verify the effectiveness of the method of the present invention, we selected the Zhenchengdi mining area as the research area and took the correction process of the last SBAS result data (i.e., total settlement) as an example. The present invention is further described in detail with reference to the accompanying drawings, and the specific implementation steps are as follows:

[0063] (1) Data preparation: The mining in the study area started on March 25, 2016 and ended in April 2017. We obtained 35 SAR data from March 2016 to September 2017 in the study area, and collected monitoring data from known monitoring points in the mining area.

[0064] (2) Perform SBAS-InSAR processing on the acquired multi-temporal SAR image data: First, select the study area and crop the SAR image to generate a short baseline connection map and select the super master image.

[0065] Interference processing is performed to eliminate terrain phase errors, weaken the influence of noise, and obtain high-precision ground deformation information. Through track refinement and secondary inversion calculation, atmospheric effects and systematic errors are eliminated to obtain the surface deformation results of the mining area time series, such as Figure 3 As shown in the figure, by processing the 35 SAR data collected in the study area, the maximum total sedimentation value in the SBAS results obtained is 169.6 mm.

[0066] (3) The sedimentation data of each monitoring point were analyzed. First, the sedimentation rate of the monitoring point in different time periods was calculated. The rate was calculated using the sedimentation data of two periods. The dynamic rate interpolation method was used to generate the sedimentation value of the monitoring point at the corresponding time of the SAR image. The maximum sedimentation value of the monitoring point was 1410.3 mm, indicating that the SBAS monitoring results are unreliable in areas with severe sedimentation.

[0067] The monitoring point and SBAS results are spatially aligned to obtain the settlement value corresponding to the SBAS at the monitoring point, and the correction coefficient at the monitoring point is calculated according to the correction model.

[0068] (4) Use the Kriging interpolation method to calculate the correction coefficients for the entire study area at each period, such as Figure 2 As shown: Calculate the distance D between each monitoring point and the corresponding semi-variance γ(D);

[0069] The spherical model is used to fit the relationship between the distance D and the corresponding semivariance γ(D) at the time node, and the distance and semivariance between each point in the study area and the monitoring point are calculated based on the trained spherical model;

[0070] After solving the objective function, the optimal weight coefficient θ of the point position can be obtained j , and find the correction coefficient;

[0071] According to the above process, the corresponding correction coefficients are solved for each point in each study area to obtain the correction coefficients of the study area, such as Figure 3 As shown;

[0072] The same process is used to construct the correction model for data from different periods, and the correction coefficients for any point in the study area are solved.

[0073] (5) Superimpose the obtained correction coefficient and SBAS result, such as Figure 3 As shown in the figure, the corrected SBAS results can more comprehensively and accurately represent the overall settlement results of the mining area. The SBAS settlement data of each period are corrected according to the correction coefficients of different periods to obtain more accurate dynamic monitoring results of the study area.

[0074] (6) In order to verify the accuracy of the corrected results, we use a monitoring point that is not involved in the calculation ( Figure 3 The measured data of five time nodes and the corrected SBAS results were selected for verification. It can be seen from Table 1 that the maximum error between the corrected SABS results and the measured values ​​(both are cumulative settlement) is 4.8 mm, and the root mean square error is 3.79 mm, which verifies the effectiveness of the method of the present invention.

[0075] Table 1 Measured values ​​and calibration results

[0076] time 2016-10-02 2016-11-19 2017-02-11 2017-03-19 2017-05-06 Measured value 1182.3mm 1283.2mm 1300.8mm 1303.2mm 1307.2mm Corrected SBAS results 1180.6mm 1286.4mm 1296.0mm 1301.5mm 1302.8mm

[0077] The above-mentioned embodiments only express several implementation methods of the present invention, and the description is relatively specific and detailed, but it cannot be understood as limiting the scope of the invention patent. It should be pointed out that for ordinary technicians in this field, several modifications and improvements can be made without departing from the concept of the present invention, which all belong to the protection scope of the present invention.

Claims

1. A method for monitoring surface subsidence and large deformation, characterized in that: The steps include: S1) Obtain n SAR images during the mining period of the study area, and obtain the original SBAS-InSAR deformation results of the study area in n-1 periods based on the SBAS-InSAR technology; S2) Collect the settlement data of known monitoring points in the study area, analyze and process the settlement data of the monitoring points using the dynamic settlement rate interpolation method, and obtain the theoretical settlement values ​​of the corresponding monitoring point positions at each period of SBAS-InSAR; S3) calculating the settlement deviation of the monitoring point position at each period of SBAS-InSAR according to the original SBAS-InSAR deformation results at each period and the theoretical settlement amount of the monitoring point position at the corresponding period, and using it as the correction coefficient of the original SBAS-InSAR deformation results at each period of SBAS-InSAR at the monitoring point position at the corresponding period; S4) using the Kriging interpolation method to solve the correction coefficient of the original SBAS-InSAR deformation result at the monitoring point position in each period to obtain the correction coefficient of the original SBAS-InSAR deformation result at any point position in the corresponding period; S5) obtaining the corrected SBAS-InSAR deformation result of any point according to the original SBAS-InSAR deformation result of each period and the correction coefficient of any point position in the corresponding period, thereby obtaining the corrected dynamic settlement result of the mining area.

2. The method for monitoring surface subsidence and large deformation according to claim 1 is characterized in that: The step S2) specifically includes: Obtain the location information of known monitoring points in the study area and the settlement monitoring value S at each time node t; Calculate the sedimentation rate of adjacent monitoring rounds at the monitoring point. The calculation expression is: Where i represents the monitoring rounds of each monitoring point, t i represents the time node corresponding to the i-th round of monitoring, Indicates that the monitoring point is at t i+1 -t i The sedimentation rate over time, It means that at t = t i The settlement monitoring value of the monitoring point at that time; According to the dynamic sedimentation rate of each monitoring round, the time point T at which n long-time series SAR image data are acquired is calculated. x When , the theoretical settlement value of each monitoring point is calculated as: 1 of them <x≤n,t i ≤T x <t i+1 , Indicates time node t = T x The theoretical settlement value corresponding to the monitoring point at that time.

3. The method for monitoring surface subsidence and large deformation according to claim 2 is characterized in that: The step S3) specifically includes: Calculate the theoretical settlement amount of the corresponding monitoring points of SBAS-InSAR in each period The calculation expression is: According to the original SBAS-InSAR deformation results at different periods and the theoretical settlement of the monitoring point in the corresponding period Calculate the settlement deviation of the monitoring point position at each period of SBAS-InSAR, and the calculation expression is: in It represents the correction coefficient of the original SBAS-InSAR deformation results at the monitoring point position in each period.

4. The method for monitoring surface subsidence and large deformation according to claim 3 is characterized in that: The step S4) specifically includes: The correction coefficients of the original SBAS-InSAR deformation results at the corresponding monitoring point positions in each period are analyzed. Assuming that there are k monitoring points in the study area, the calculation expression of the correction coefficient of any point in the study area in each period is: in It means t = T x The correction coefficient for any position (x,y) at the time, θ j represents the weight coefficient, It means t = T x The correction coefficient corresponding to the monitoring point j at time .

5. The method for monitoring surface subsidence and large deformation according to claim 4 is characterized in that: The step S4) specifically includes: Calculate the weight coefficients, assuming that the correction coefficients at any point (x, y) in the study area have the same expected value c and variance σ 2 , that is, it holds true for any point: Calculate the distance D and the corresponding semi-variance γ(D) between the known monitoring points. The distance is obtained based on the coordinate value: Where D 1,2 represents the straight-line distance between monitoring point 1 and monitoring point 2, Δx and Δy represent the differences in the horizontal and vertical coordinates of the two monitoring points respectively; The semivariogram function is defined as: in represents the semivariance between monitoring point 1 and monitoring point 2; All the obtained (D,γ) are fitted with the spherical model, and the function is: Where C0 represents the nugget constant, C0+C represents the base value, C represents the arch height, and a represents the range; For any point in the study area, firstly calculate the distance D from all known monitoring points and the corresponding semivariance γ(D) according to the above method; The goal of the Kriging interpolation method is to find a set of weight coefficients that can satisfy the minimum difference between the estimated value and the true value of the correction coefficient at any point in the study area, that is: in, represents the estimated value of the correction coefficient at any point in the study area, ε(x i ,y i ) represents the true value of the correction coefficient at any point in the study area; At the same time, it satisfies the unbiased estimation condition: By constructing a short covariance matrix and a covariance vector, combined with unbiased constraints, a linear equation system is established and the optimal weight coefficient vector can be obtained by solving it.

6. The method for monitoring surface subsidence and large deformation according to claim 4 is characterized in that: The step S5) specifically includes: The original SBAS-InSAR deformation results of the study area in n-1 periods are corrected, and the calculation expression is: in Represents the corrected SBAS-InSAR deformation results.

Images