Application of Gray-Markov Model to Land Subsidence Monitoring of a Mining Area

Land subsidence monitoring in mining areas is one of the main applications of surface deformation monitoring, which is of great significance for safety production. Using the IPTA (Interferometric Point Target Analysis) time-series InSAR (Interferometry Synthetic Aperture Radar) method, land subsidence data of the new exploration area in the Weizhou mining area were analyzed and compared with static GPS (Global Positioning System) monitoring data from 2017 to 2020. Gray-Markov model was established by combining the gray prediction model with the Markov model to predict the surface subsidence of the mining area. The results show that <xref rid="deqn1" ref-type="disp-formula">(1)</xref> InSAR data has high accuracy and application potential in prediction of long-term surface deformation in mining areas; <xref rid="deqn2" ref-type="disp-formula">(2)</xref> The Gray-Markov model can better reflect the volatility and practicality of subsidence data in mining areas; <xref rid="deqn3" ref-type="disp-formula">(3)</xref> The prediction results have high accuracy, and the Gray-Markov model can serve as an effective guide for long-term surface deformation monitoring and safety management.


I. INTRODUCTION
Coal has always occupied a high proportion in China's energy consumption. China is rich in coal resources, mainly distributed in Shanxi, Inner Mongolia and Shaanxi. The reserves of coal resources accounts for 33.8% of the world's total coal reserves, and the mining production ranks the first in the world [1]. While coal resources brings huge economic benefits, the land subsidence and impounded surface water will cause varying degrees of damage to buildings, farmland and roads. When the situation is serious, it will cause ecological and environmental problems, effect people's lives and hinder economic and social development. Therefore, accurate long-term dynamic monitoring of the mine is of great significance for the assessment and control of geological disaster and ecological reconstruction.
In recent years, increasing attention to safety and geological hazard assessment in mining areas has resulted in a greater need for land subsidence monitoring in these areas. The monitoring and prediction of surface deformation in mining areas play an important role in evaluating potential geological The associate editor coordinating the review of this manuscript and approving it for publication was Gerardo Di Martino . hazards and ensuring safety production. Traditional monitoring of mine surface deformation usually involves traditional geodetic techniques, such as leveling instruments and GPS (Global Positioning System), to observe the subsidence and regular or irregular horizontal movement along the two levelling lines [2]. Although the measurement accuracy of these methods is high, the monitoring range is small, the workload is large, the cost is high, the efficiency is low, the measuring points are easily destroyed, and continuous monitoring is not possible [2], [3]. Traditional methods usually only monitor the mining-induced ''linear'' deformation of a few working faces in the mining area, and the spatial and temporal resolution is low [3]- [5]. In fact, it is difficult for traditional ''linear'' deformation observation to reflect the spatiotemporal evolution of the three-dimensional deformation of the whole mining area. This limitation hindered the overall monitoring of the mining area and the theoretical research, thereby reducing the reliability of potential geological hazard prediction.
InSAR (Interferometry Synthetic Aperture Radar) technology has the advantages of all-weather use, low cost, large coverage, and high spatial resolution. One example is Sentinel-1 TOPS (Terrain Observation with Progressive Scans in azimuth) mode SAR data (open to all users) [6], [7].
Ideally, SAR data obtained using the TOPS mode can monitor about 40,000 km 2 for 6 days (with two satellites) with a return period of 12 days to produce surface deformation maps with a spatial resolution of about 20 m (including all mining areas in this range) [8]- [11]. Surface deformation maps with such high spatial and temporal resolution are almost impossible to obtain using traditional geodetic instruments [12].
There are mainly 4 types of subsidence prediction methods in the mining areas: (1) Theoretical research method. This method usually uses the elastic-plastic theory, according to the engineering conditions to calculate the surface subsidence. However, many factors involved in the modeling process. Some conditions need to be assumed as ideal and the results do not fully correspond to the actual subsidence.
(2) Mathematical method. The typical curve method usually collects a large number of measured data and processes it, draws the curve and finally obtains the land subsidence in a mining area. However, the application of this method is limited, and the data processing timeliness is not strong; The profile function method is modeled according to the measured data for the land subsidence calculation, and this method can not be used to analyze the subsidence mechanism. (3) Test model method. Before mining, the model of mining area is established according to a certain scale of reduction. Then we simulate the actual loading and mining process and force the model, find the changes of its force and displacement, infer the changes of land subsidence in the actual mining process, and then analyze the law of changes. But The process of establish model before mining is too complicated and each model can only be used in the specific area which is not universal. It's barely used at present. (4) The new prediction method mainly uses computer theory, nonlinear theory and mathematical theory, and achieved some results in practical application.
For the prediction of land subsidence in mining areas, the combination of the Gray model and the Markov model can effectively overcome the limitations of the traditional GM(1,1) (Single variable gray model of first order) and significantly improve the prediction accuracy [12]. The Gray-Markov model can play an important guiding role in land subsidence monitoring and safety management in mining areas [13], [14].
Studies have shown that the deformation monitoring data in mining areas has a certain gray characteristics, but most of these studies are based on observatory leveling or GPS measurements. It is independent of the subsequent evolution links, and cannot form a complete and systematic monitoring and prediction model [15]. This paper takes the new exploration area of Weizhou mining area as the research area, considering the dynamic nature of the actual surface subsidence, we collected InSAR technology and Gray-Markov model, proposed a new dynamic prediction method, processed the Sentinel-1 satellite data from 2017.01 to 2020.12 with IPTA (Interferometric Point Target Analysis) time series InSAR technology, and compared the result with the GPS data to verify the accuracy of the InSAR data. A combined model of Gray model and Markov model is established to predict the subsidence values of four typical subsidence points in the study area, which can provide reference for the dynamic monitoring, real-time prediction and subsequent treatment. Based on the subsidence results of 2017.01-2020.06, the surface subsidence of 2020.07-2020.12 is predicted every three months, and the predicted value is compared with the InSAR results. The accuracy of the model is analyzed to provide reference for long-term surface deformation monitoring and safety production in mining area.

II. RESEARCH METHOD
Subsidence from 2017.01.14 to 2020.12.24 in the new exploration area in the Weizhou mining area were obtained by processing time-series InSAR data, which were compared with existing static GPS data to verify the accuracy and reliability of the InSAR results. Based on the GM(1,1) model, the subsidence trend was predicted, and the subsidence values were optimized by the Markov model, which compensates for the low prediction accuracy of the Gray model and predicts the land subsidence in the mining area with high accuracy. The Gray-Markov model has good application potential for long-term land subsidence monitoring in mining areas.
The Gray-Markov model can effectively predict data in both the horizontal and vertical dimensions; the Gray model can accurately capture the trend of the data, and the Markov model can effectively solve the problem of data volatility and dispersion. The experimental data includes Sentinel-1 satellite image data, SRTM-30m DEM data, Sentinel-1 satellite precision orbit data, and static GPS field observation data. The study area is shown in figure 1.

III. IPTA TIMING InSAR PROCESSING
The main principle of coherent point target timing InSAR technology is to obtain the geometric characteristics of the target point and its related change information based on interference phases from satellite radar data [15]. The IPTA VOLUME 9, 2021 method can complete the extraction of high-precision deformation information in a large range, and in contrast to traditional InSAR technology, it is not seriously affected by incoherent and atmospheric delay factors [11], [16]. The phase information for the interference fringes of any two images includes (1) terrain phase information; (2) level flat-earth phase information; (3) the phase of radar line-of-sight (LOS) changes between two images; (4) the atmospheric delay difference; and (5) the image noise difference [17].
During IPTA processing, the interference phase can be decomposed into: where ϕ is the interference phase of the target point; ϕ top is the interference phase caused by topographic elevation, which can be calculated by adding external DEM data; ϕ f represents the flat ground phase, which can be calculated from the imaging geometric relationship; ϕ def represents the spatial vertical baseline of the directional alignment of the radar line-of-sight (LOS); ϕ atm represents the noise phase caused by the atmospheric delay; and ϕ n is the noise phase. Among the above variables, ϕ top and ϕ f can be removed by differential interference, resulting in the differential phase ϕ diff : where ϕ tope is the elevation error phase. The required surface deformation information can finally be extracted from the separation ϕ def .
In the above formula, B ⊥ indicates the spatial vertical baseline of differential interference image pairs; λ represents the wavelength of radar waves; R is the slant range from ground target points to radar sensors; and θ represents the radar incident angle.
where v represents the linear variable rate of the coherence point relative to the reference point, t represents the time baseline, and ϕ def_n represents the nonlinear deformation phase. The IPTA method uses the two-dimensional linear phase model, completes the separation and removal of various errors by constantly returning to the iterative scheme, and then obtains the linear deformation rate of the coherent point target [11]. During this iteration, the Delaunay triangulation network is constructed, and the phase unwrapping of the point target is realized by the MCF (Minimum Cost Flow) algorithm. The quality evaluation criteria for phase unwrapping and iterative regression fitting are the phase standard deviation and residual phase smoothness.

IV. GRAY SYSTEM GM(1,1) MODEL
Gray system prediction refers to the prediction of the eigenvalue change of a system's behavior, where the system contains both known and uncertain information; that is, changes are predicted in a certain threshold range and are related to a time series. The phenomena in the gray process are random and fluctuating but also orderly and bounded, so the data set has a potential law [16]. The aim of the gray system prediction model is to use this potential law to establish the model and then realize the prediction of the system.
The process of establishing the gray model is as follows: 1) Suppose that the original data sequence corresponding to time is: 2) In order to reduce the dynamic randomness of the data, the original data sequence is accumulated: 3) Construct the adjacent mean value sequence: The adjacent mean value y (1) (k) is the whitening of the background value. 4) Construct whitening differential equations: Use least square method to calculate: . . .
The tests of the GM(1,1) model, including the mean variance ratio (C) test and small probability error ( ) test,are performed to determine the merits of the model. The inspection criteria are shown in Table 1. If the model meets these two requirements, then it is considered to be qualified.
In the average variance ratio test, the original data variance is S 2 1 , and the variance of the residual series is S 2 2 : where x 0 (k) is the raw data sequence,x (0) is the mean value of the raw data sequence, ε (0) (k) is the residual sequence, andε (0) is the mean of residual sequence [17].

V. MARKOV PREDICTION MODEL
The Markov model was proposed by the Russian Mathematician Markov. It models a stochastic dynamic process called the Markov chain [18], [19]. The Markov transfer process is only directly related to the previous connected data but not to other data in the past, which is called ''post-invalidity'' [20]. The model expression is: The Markov model divides data into different state intervals (generally divided into 3 or 4 for both prediction accuracy and data complexity) [21], [22]. The optimal state is found by the step of the state transition matrix, and the future change is estimated [23].

A. STATE DIVISION
In this paper we take the ratio of the original InSAR subsidence data and the GM(1,1) model prediction data as a reference, and the data are divided into three intervals according to the ratio value expressed as S i ∈ [a i , b i ] , i = 1, 2, · · · n. The lower and upper limits of the interval are a i , b i , respectively.
The transfer probability from state S i to state S j by step k is expressed as: where k indicates the number of steps from S i to S j , P k ij expresses the transfer probability from state S i to state S j by step k, n i denotes the sample numbers of S i , and n k ij denotes sample numbers from state S i to state S j by step k [24].

B. CONSTRUCT THE STATE TRANSITION MATRIX
The state transition matrix composed of the transition probability P ij is:

C. DETERMINE THE OPTIMAL VALUE
The matrix P (k) is the state corresponding to the maximum value of the column vector sum in the state list of row vectors [25]; that is, Markov's ideal state. Its optimized forecast value is: where the Gray-Markov prediction model is the final optimized prediction value y (0) (i + 1), a i is the lower limit value, b i is the upper limit value, andx (0) (i + 1) is the forecast of the gray model [26].    and its technical mining conditions and market prospects are promising. We collected SLC data from 114 images captured by the Sentinel-1 satellite in an ascending orbit covering the study area. The time span is from 2017.01.14 to 2020.12.24. The DEM data is SRTM data, and the resolution is 30 m. Specific data parameters are shown in Table 2.

B. IPTA TIMING ANALYSIS
GAMMA software was used to process the satellite image data from the new exploration area in the Weizhou mining area for each year, and the results were compared with static GPS monitoring data for the same period and region. According to the InSAR time-series monitoring results, surface subsidence mainly occurs in the central and southern parts of the exploration area, and the northern subsidence is not obvious. Four typical subsidence points were selected for data analysis.   The annual subsidence results from the InSAR monitoring data are shown in figure 2: The distribution of existing GPS monitoring points and typical InSAR subsidence points in the study area is as shown in figure 3 (InSAR subsidence points are represented from I1 to I4,GPS monitoring points are G1 and G2): We have set the static GPS monitoring points (G1 and G2) near typical InSAR subsidence points (I1 and I3) in the monitoring area. Taking every three months as an observation period, we got the time-series InSAR surface subsidence results and static GPS monitoring results. The results of IPTA-InSAR in the period from 2017.01 to 2020.06 were compared with the data for adjacent static GPS monitoring points as shown in figure 4 and figure 5: The results show that although the InSAR data points do not completely overlap with the GPS monitoring points,    InSAR data has the same deformation trend and close numerical values with adjacent GPS monitoring data. Thus, the processing result is correct.

C. ESTABLISHMENT OF GRAY MODEL
Between 2017.01 and 2019.12.28, based on an observation period of every three months, the temporal InSAR surface subsidence results were obtained, and the original data sequence was constructed. By comparing the ratio of the actual value of InSAR temporal monitoring to the value predicted by the gray model, the reliability of the model was tested. Computer program was used for formulas (5)- (14) to realize the gray prediction process and obtain simulation prediction results, as shown in figure 6: The results of the model accuracy test are shown in Table 3. According to the test results, the accuracy of the gray model for Points I1, I2, and I4 meets the first standard in Table 1, and the accuracy of the gray model for typical subsidence Point I3 meets the second standard. The prediction results are therefore reliable, and the model can be optimized and used for prediction.

D. OPTIMIZATION OF THE GRAY-MARKOV MODEL
The GM(1,1) model can accurately capture data development trends, but its prediction accuracy is limited. The Markov model can be further optimized to enhance the accuracy of the volatility data prediction.  By using computer program, formulas (15)- (18) were programmed to realize the Markov random process, and the predicted results were obtained for the two models. The predicted results are as shown from figure 7 to figure 10: Table 4-7 show the predicted accuracy for subsidence data at each monitoring point after the optimization of the Gray-Markov model. We get the predicted accuracy by calculating the ratio of the true value to the predicted value. Apparently, predicted results of the Gray-Markov model is closer to the true value than gray model. The maximum improved accuracy is more than 24%. Majority of the accuracy of Gray-Markov model's results can reach 95%, some of them    The Gray-Markov model was used to predict the subsidence data in the second half of 2020 according to the original data sequence, and the values were compared with the InSAR monitoring results. The predicted results are in line with the trend of InSAR data; the accuracy is high, and the prediction data have high reliability . The results are shown  from table 8 to table 11. VII. CONCLUSION 1) Through the processing and prediction of InSAR data and comparison with adjacent GPS data, it has been shown that the InSAR data have high accuracy and application potential in the prediction of long-term surface deformation in mining areas.
2) After combining the traditional gray prediction model with the Markov model, the optimized prediction model is closer to objective reality. The prediction results show that the Gray-Markov model can better reflect the volatility and practicality of the subsidence data in the mining area than traditional gray model.
3) The Gray-Markov model was used to simulate and predict subsidence data from 2017 to 2020 in the Weizhou mining area. The prediction accuracy is high, so the Gray-Markov model can serve as an effective guide for long-term surface deformation monitoring and safety management in mining areas.
4) Although the Gray-Markov model can better predict the actual situation, it is not suitable for surface subsidence monitoring in all mining areas. The model has a certain applicability to the observation object, which is more suitable for the mine subsidence monitoring with a fast to slow tendency. He has been an Associate Professor with China University of Mining and Technology (Beijing), a master's supervisor, and the Director of Surveying and Mapping. He is the author of two books, more than 40 articles, and two inventions. He participated in the completion of two National Natural Science Fund projects and one of 1973 sub-projects. His research directions include GPS positioning and navigation, deformation disaster monitoring and data processing, 3S integration and application, 3-D laser measurement technology and its application, UAV disaster monitoring, and data processing. His awards include the second prize of Beijing Higher Education and Teaching Achievement Award, the second prize of Geographic Information Technology Progress Award, and the second prize of Coal Industry Education and Teaching Achievement Award.
CHENGXING GENG was born in November 1997. He received the bachelor's degree in engineering in surveying and mapping engineering from Liaoning Technical University, in 2019 and the M.Eng. degree in surveying and mapping engineering from China University of Mining and Technology (Beijing).
His research interests include deformation monitoring data processing and InSAR landslide deformation monitoring.
LING ZHANG received the M.S. degree in photogrammetry and remote sensing from China University of Mining and Technology (Beijing), in 2006, with a thesis on SAR interferometry application in land subsidence surveys. From 2006 to 2008, she worked as a Technical Support Engineer of RS and GIS data production. Since 2008, she has been with China Aero Geophysical Surveying and Remote Sensing Center for Natural Resources. Her research interest includes InSAR application/research on land deformation.
ZHENCHAO ZHANG was born in Heze, Shandong, China, in March 1994. He received the master's degree from China University of Mining and Technology (Beijing), in July 2021. His research interests include deformation monitoring data processing and InSAR landslide deformation monitoring. VOLUME 9, 2021