Methods for Monitoring Fast and Large Gradient Subsidence in Coal Mining Areas using SAR Images: A Review

Ground subsidence caused by coal mining may lead to seriously damage to ecological environment system and economic and social development. SAR image interferometry technology has been widely used for monitoring coal mining subsidence. However, at present it is still a challenge for it to monitor fast and large gradient subsidence. In order to obtain fast and large gradient subsidence using SAR images, some scholars have resorted other methods and technology which doesn’t just depend on the interferograms, and they have obtained better results. In a summary, these methods and technology can be grouped into three main categories: (i) Combining InSAR with coal mining subsidence prediction models; (ii) Using SAR image amplitude information offset tracking technology; (iii) Combining InSAR with terrestrial 3D laser scanning technology. In this paper, the technology features, research progress, applications and limitations of these methods are analyzed; and some potential prospects of using SAR images to monitor fast and large subsidence in coal mining areas are discussed in the end.


I. INTRODUCTION
Underground coal mining has a history of thousands of years. Coal resource is still the most important source of energy at present. Underground coal mining may cause the rock layers and ground above the goaf to bend, move, subside, or even collapse, and then destroy ground facilities and ecology. China is the world's largest coal producer. In 2019, the annual coal production by China was 3692Mt, it was about half of the total world coal production. Since 2000, coal production by China has increased by 172.47% [1]. In China, coal is still one of the most important energy sources, and the most coal mined in China is for domestic consumption. Therefore, coal production will remain large for a long time [2] [3]. Large scale underground coal mining has caused serious ground subsidence in China. For example, some coal mining areas in Yulin subsided 3m within 11 days [5], the maximum subsidence rate reached 4m/yr in some high-intensity mining areas [6], severe surface subsidence could reach up to 6-12 meters [7]. Fast and large gradient ground subsidence may lead to seriously damage to farmland, roads, railways, buildings, underground pipeline facilities, water resources, ecological, environmental problems. About 4 hm 2 of land and 2.2 billion m 3 of groundwater resources are destroyed every year [4]. Fig. 1 shows some serious ground subsidence and damage in some coal mining areas. Serious subsidence land is unsuitable for agriculture or other economic and social activities, and inhospitable to wildlife [7]. Therefore, monitoring serious subsidence in coal mining areas has practical significance to subsequent coal mining, the safety of human property and life, land reclamation and ecological reconstruction in China and around the world.
The conventional coal mining large gradient subsidence monitoring methods are mainly topographic survey methods using Levels, total stations, GNSS receivers and so on instruments. These methods mainly monitor some sparse feature points of the subsidence area, and use these feature points to make a profile to show the subsidence. If we want to obtain the full subsidence basin of a coal mining area, we must measure a large number of points, and then draw the subsidence contour map to estimate the whole subsidence basin by interpolating and fitting. If we want to obtain the subsidence process of a subsidence basin, we must repeatedly measure these points. The methods and techniques have been in development for many years and have become relatively mature. These methods can achieve centimeter or even millimeter level monitoring accuracy whether small or large gradient ground subsidence. However, these conventional coal mining subsidence monitoring methods have some disadvantages: (i) The points where people can not arrive at and blind spots cannot be measured by Levels and total stations. (ii) Both temporal and spatial resolutions of these sparse feature points are very low. The distance between leveling points often exceeds 100m, the distance between GNSS points often exceeds several kilometers. It takes quite a long time to measure all feature points each time. (iii)It is time-consuming and laborious. Obtaining the full subsidence basin and its subsidence process is rather time-consuming and laborious. (iv)Interpolating and fitting will also bring additional errors. In a word, it is not economical and efficient for monitoring large-scale, high-frequency and continuous large gradient subsidence information in coal mining areas using the conventional topographic survey methods [8] [9] Compared to these conventional topographic survey methods, methods using Synthetic Aperture Radar (SAR) image interferometry technique, especially Differential Interferometric SAR(DInSAR), have many advantages for monitoring coal mining subsidence [10][11][12][13][14][15]: large area coverage, high spatial resolution, low cost, all day and all night, almost all-weather capabilities and the monitoring result is a region rather than a series of feature points, obtaining SAR images is getting easier, etc. Vertical subsidence is the main deformation in coal mining areas, DInSAR is sensitive to vertical deformation, therefore, DInSAR has become an extremely powerful tool for monitoring ground subsidence in coal mining areas and has been widely used [16][17][18][19][20][21]. However, these studies have mainly focused on small or slow deformation. Some scholars have found that only small or slow subsidence (e.g., several millimeters or centimeters between two acqusitions ) can be reliably detected and measured [22] [23]. When there is a fast and large subsidence (e.g., more than 100 millimeters between two acqusitions), in the center of subsidence basin, measurement can hardly be obtained from the SAR image interferograms [24][25][26][27]. For example, in a study, the maximum subsidence was 3000mm in a coal mine in China, but only the maximum 80mm was obtained by using TerraSAR-X images [28]. This is mainly due to the limitation of the maximum detectable deformation gradient(DDG) of DInSAR (Detailed in Section II). Some advanced time series InSAR methods have been proposed, such as SBAS-InSAR (Small Baseline Subset InSAR) [29][30], PS-InSAR (Persistent Scatter InSAR) [31], IPTA(interferometric point target analysis) [32], SqueeSAR [33], CRInSAR(corner reflector InSAR) [34], TCPInSAR (temporarily coherent point InSAR) [35], Ultrashort-Baseline TCP [36], and so on. These methods have been effective in decreasing the effects of temporal and spatial decoherence and atmospheric delay, especially in the urban areas with a large number of strong interference points and low deformation velocity. However, in large gradient subsidence coal mining areas, due to the limitation of the maximum DDG of interference phase, nonlinear deformation, few strong interference points, it was difficult to achieve ideal monitoring results by these methods [25] [26]. At present, it is still a challenge for using SAR image interferograms to monitor fast and large gradient coal mining subsidence.
In order to overcome above problems, and use SAR images to obtain a large number of high spatial resolution subsidence points, instead of rather than just some sparse feature points. Over the last twenty years, many scholars have proposed and experimented more than twenty different methods for monitoring large gradient subsidence in coal mining areas using SAR images, and some of them have obtained better results. In a summary, these methods can be grouped into three main categories: (i) combining InSAR with mining subsidence prediction models; (ii) SAR image amplitude information offset tracking technique; (iii) combining InSAR with terrestrial 3D laser scanning technique. This paper focused on providing a basic overview of these methods and applications, and summarized some existing limitations of these methods, prospected some future developments. We do not attempt to give a comprehensive review of the SAR technique for monitoring deformation. In this paper, SAR image means acquisitions by spaceborne SAR, the airborne SAR and ground-based SAR are not considered. And, the term 'small subsidence' means the subsidence basin is relatively gentle, not more than the detecting limit of InSAR; 'large gradient subsidence' means that the ground subsidence is very deep (As shown in Fig. 2).

II. THE PRINCIPLE OF DINSAR AND ITS LIMITATION FOR MONITORING FAST AND LARGE SUBSIDENCE
SAR systems emit microwaves and collect the returned amplitude and phase in the line of sight(LOS) direction from the ground targets. DInSAR is used to extract the deformation of the ground target in the LOS direction by calculating the returned phase difference of the same ground target acquired at different times. The basic principle of DInSAR is as follows [37] : Supposing there is a single pixel point P on the ground. The satellite sensor acquires a first SAR image and measures a returned phase of point P when it at the position M. Then there is a deformation, the point P moves to P1, and the satellite sensor acquires another SAR image and measures another returned phase of point P1 when it at the position S. As shown in Fig.3.
Where MP is the distance from sensor to P when the sensor at the position M, SP1 is the distance from sensor to P1 when the sensor at the position S;  is the radar wavelength; scatter is the target random scatter phase, is related to the character of target. It is assumed that the random scatter phase doesn't change, it is the same when the sensor at M and S.
By subtracting and adding the term SP/(/4),equation (4) becomes Where the first term is the phase difference at different times before deformation, it is the topographic phase component. The second term is the deformation phase component, related to the LOS deformation L shown in Figure 3. Theoretically we can simulate the topographic phase component by a Digital Elevation Model (DEM) or other ways and subtracte it from (5), the deformation L can be calculated by the interference phase int. More details on the DInSAR fundamentals can be found in many papers.

FIGURE 3. Principle sketch of DInSAR technology
The surface deformation in SAR interferograms is represented by a series of interference fringes, and each fringe corresponds to half of a radar wavelength deformation in the LOS direction [38]. If the deformation of a point between the two acquisitions is more than half a wavelength in the LOS direction, the deformation cannot be recovered from the SAR image interferograms [39]. Accordingly, Massonnet et al. first proposed the concept of deformation gradient, and the theoretical equation of the maximum detectable deformation gradient(DDG) by InSAR [24].
The maximum DDG obtained by (6) is a theoretical value. However, the detectable deformation gradient is affected by many factors: spatial and temporal decorrelation, Doppler centroid shift, atmospheric delay, noise and so on [39] [46].. Based on the coherence, Baran et al. (2005) constructed the function model between the maximum DDG and interferograms coherence of the ERS SAR image (as (7)). The results showed that the maximum DDG was 4×10 -5 for coherence equal to 0.32 with assuming the SAR image pixel size 20m×20m [43].
In order to reduce noise and achieve an approximate square pixel, multi-look processing is implemented. Multi-look  processing also degrades the resolution, which then greatly reduces the maximum DDG [41][42] [47]. Under the influence of various factors, in practical applications, the actual maximum DDG is far less than the theoretical maximum DDG. Experiments in a coal mining area used RadarSet-2 images showed that the DInSAR monitoring results were consistent with the leveling results when the vertical subsidence was less than 100mm, but when the vertical subsidence was more than 100mm, the DInSAR monitoring results were different from the leveling results and seriously inconsistent with the actual values [44]. Another experiment showed that 100mm deformation could make the L-band SAR images completely decoherent, and 30mm deformation could make the C-band images completely decoherent [48]. The limitation of maximum DDG is mainly due to phase interference. In order to avoid the limitation of obtaining fast and large gradient subsidence in coal mining areas by phase interference method, many scholars have resorted to other methods, such as amplitude offset tracking technique, or combining with coal mining subsidence prediction models or terrestrial 3D laser scanning point cloud data.

III. COMBINING WITH MINING SUBSIDENCE PREDICTION MODELS
A large number of measured data show that the process of ground subsidence caused by underground coal mining has certain laws. After mining for a period of time, the ground will begin to subside and move horizontally; the subsidence acceleration changes: 0→-amax→0→+amax→0; subsidence speed changes: 0→vmax→0; and the subsidence changes: 0→ Wmax. Based on these laws, many scholars proposed some coal mining subsidence prediction models to predict the development process of ground subsidence, including: probability integral model (PIM), Knothe model, hyperbolic model, logistic model and so on, which have been verified to be able to predict the coal mining surface subsidence fairly well. Some scholars have tried to combine SAR image interferometry technique with different types of coal mining subsidence prediction models to get large subsidence information. In this section, three main types of combining InSAR with coal mining subsidence prediction models are described and discussed in detail. Some of them are summarized in Table II.

A. InSAR AND PROBABILITY INTEGRAL MODEL
The probability integral model (PIM) is proposed by LIU B. & LIAO G in 1965 based on the stochastic medium theory [49]. With decades of development, it has become one of the most widely used methods to predict the ground subsidence in the coal mining areas in China. Its name comes from its movement and subsidence prediction formula in which there is a probability integral expression. For a ground point p(x,y) in the ground coordinate system xOy and any mined unit point B(u,v) in the coal seam coordinate system uO1v, according to the principle of probability integral model, the subsidence prediction formula is expressed as (8) [50] : Where ( , ) is the subsidence value of p(x,y), is the mined coal thickness in Normal direction, is the subsidence coefficient; is the dip angle of the coal seam, It is the propagation angle of mining influence, is the main is the tangent of the main influence angle, = − , = ( − ) , D3 and D1 are the length and width of coal seam, , , , are deviations of inflection point in different directions.
With the prediction parameters, PIM can obtain more accurate subsidences of any gradient. Since Wang et al. firstly proposed the method of combining the PIM and InSAR to dynamically predict the coal mining subsidence [51], many scholars have made further researches and applications on this method. Most of the researches focused on the following aspects: obtaining the subsidence values of some control points on the edge of subsidence basin by the DInSAR method; using DInSAR subsidence values to estimate and optimize the PIM parameters by some algorithm; deducing the large subsidence by PIM at last. In this process, estimating PIM parameters is the key problem.
Because of the advantages of Genetic algorithm (GA) in parameter estimating, it is the most used algorithm. The basic idea of estimating PIM parameters by GA is as follows [52]. Encoding all the PIM parameters into a GA chromosome strings. Evaluating the fitness of each individual in the population by a fitness function and eliminating the small fitness individuals. Using crossover, mutation and other operators to perform genetic operations on the most fitness individuals to generate the next generation population. Repeating the above steps until the optimal PIM parameter sequence is obtained (As shown in Fig. 4). Using this method, H. D. Fan et al. obtained the whole subsidence basin, and extracted maximum about 1000mm large subsidence at a coal mining goaf located in Huaibei (As shown in Fig. 5(a)) [52]. And the relative error of the maximum subsidence point was much better than the results of DInSAR. Other scholars obtained more than maximum 1500mm, even 4000mm large subsidence using SAR, GA and PIM method at other coal mines [53] [54]. VOLUME XX, 2017 9

FIGURE 4 Flowchart of estimating the PIM parameters by GA[52]
In most researches, the horizontal deformation was ignored and the DInSAR LOS deformation was directly converted into the vertical subsidence. The PIM parameters related to vertical subsidence were estimated based on the converted vertical deformation components. This may result in the improper conversion of the real horizontal deformation of the surface to the vertical deformation. Moreover, the parameters related to horizontal deformation of PIM could not be estimated, empirical values were used when using PIM to estimate ground deformation. Both would increase the error of the estimated vertical deformation by PIM.
In order to overcome above horizontal deformation problems, an InSAR-PIM method was proposed in [55]. Firstly, constructing the relationship between the LOS deformation and the full parameters of the PIM according to the principle that the deformation components in vertical, east-west and north-south predicted by PIM should be consistent with the deformation components monitored by InSAR (Supposing there is a ground point P(x, y) ) [55]: Where , , are the vertical subsidence, north-south deformation and west-east deformation respectively, is the incident angle, is the azimuth angle, and are azimuth angles from the north or the east positive direction along counterclockwise.
= [ , , , , , , ] is the geological and mining parameters set, = [ , , , , , , ] is the model parameters set. Then, an improved GA (GA with gross error elimination, GAGEE) were used to estimate and optimize the model parameters based on the LOS deformation observations of a large number of control points. Ultimately, the horizontal deformation and subsidence were obtained by PIM based on estimated model parameters above and geological parameters from the coal mining designing document and mining schedule [55]. The experiment with three ascending ALOS PALSAR images in Qianyingzi mining area indicated that the maximum subsidence up to 856mm and the maximum horizontal movement was about 294mm, and they were consistent with conventional topographic survey result(As shown in Fig. 5(b)) [55].
In critical and supercritical mining conditions, the PIM can predict the ground deformation well. But in subcritical mining condition, the PIM often overpredicts the maximum ground subsidence, horizontal movements [56] [57]. Aiming at this problem, Z. F. Yang et al. modified the conventional PIM with a simplified Boltzmann function, others were the same as above, called generalized InSAR-PIM (GInSAR-PIM) method [58]. The experiment in the same mining area above indicated that the maximum deformation was about 850mm, 320mm, and 240mm in the vertical, east, and north directions, and was consistent with other monitoring results [58].
Genetic algorithm has some instability, the results of the estimated parameters maybe different each time even if the input parameters are exactly the same, and the search speed are relatively slow in the later period [59] [60]. In view of the above problems, Pattern Search method was used to estimate the PIM parameters in the combination of InSAR and PIM. A maximum 1017 mm subsidence was obtained in a coal mining in Yanzhou by this method [61]. However, real coal mining experiments showed that there were some errors in the PIM parameters estimated by Pattern Search method, and the maximum error is up to 20%[62] [63]. Based on the relationship constructed by the InSAR-PIM method, the simulated annealing(SA) was used to estimate the PIM parameters [64]. The simulation and real data experiments showed that this method could accurately estimate the PIM parameters, and the deformation based on this method was consistent the field measurements [64]. However, In InSAR-PIM or GInSAR-PIM method, the PIM was simplified, a few PIM parameters were not involved, and assumed that s3 and s4 were equal. And the convergence performance of the Genetic Algorithm (GA) wasn't good when there were too many unknown parameters [65]. Then the Cultural Algorithm Random Particle Swarm Optimization (CA-rPSO) method was used to estimate full PIM parameters. The experiment in the Xuehu coal mining area showed that the results from the improved algorithm were consistent with field measurements, and the the maximum subsidence was about 748 mm. However, at the center of the subsidence basin, the maximum difference was up to 99 mm, and the absolute difference of about one third of the points was more than 60mm [65].
Xiaoqian et al used the Gauss function to fit D-InSAR monitoring value and PIM prediction value, reconstructed the subsidence characteristic curve of a longwall goaf in Bu'ertai mining area by combining DInSAR and PIM, the maximum subsidence value is 1485mm [66]. This method is applicable regardless of subcritical extraction or critical/supercritical extraction, but it has great limitations. It can only obtain the subsidence of points on the main section, but cannot obtain the horizontal movement. It cannot obtain the subsidence and horizontal movement of points on the non-main section.
The combination of the PIM and InSAR method greatly enhanced research progress of InSAR monitoring large gradient surface subsidence in coal mining areas, and it has become a research focus. However, in most researches, the horizontal deformation was ignored and the LOS deformation was directly converted into the vertical subsidence, which would affect the final subsidence accuracy. The PIM has too many unknown parameters, estimating PIM parameters needs more InSAR observations and long time. In order to estimate the PIM parameters, not only the DInSAR monitoring values were used, but also a small number of field measurements values near the maximum subsidence point and inflection point in the center of the subsidence basin were used [53]. In addition, all of these methods used the geography and mining condition parameters about the goaf obtained from the coal mining designing document and underground mining records. Can the PIM parameters be estimated only by InSAR monitoring data without given geographic and mining parameters or field measurement values? Further researches are needed.

B. InSAR AND LOGISTIC MODEL
Many studies show that the coal mining subsidence curve over time has an "S" shape [67]. The logistic model is a typical "S" shaped function and has been used to predict the coal mining subsidence [68]. The logistic model only has 3 unknown parameters, fewer parameters than that of the PIM. The logistic model function is: Where ( ) is the ground surface dynamic subsidence over time , when the subsidence first happens, = = and ( ) = ; is the possible maximum subsidence; is the model parameter; is the time influence coefficient.
Some scholars also researched how to combine the logistic model with InSAR to detect mining subsidence [69][70][71]. The key problem is also how to estimate the logistic model parameters and then deduce the large subsidence by the logistic model. GA and Levenberg Marquard(LM) were used to estimate and optimize the logistic model parameters serially. GA was mainly used for obtaining the initial values of the parameters, and LM was used to refine the GA results [70]. Using this method, a maximum 1300mm subsidence was obtained in a coal mining area in Yungang [70]. In order to more accurately estimate the logistic model parameters, the functional relationship between the DInSAR LOS deformation and the logistic mode full parameters was constructed (as (11)) [71].
Where , are the accumulative times from the starting time to , ( = − , = − ). Then, an algorithm based on the GA and the simplex algorithm (SA) was used to estimate these four unknown model parameters , , , and ∆ with at least four InSAR observations. An experiment was performed in the Datong coal mining area using this method. The results showed that this method could detect 1.26 m subsidence(As shown in Fig. 5(c)) [71]. In studies above, the horizontal movement was also ignored, the vertical subsidence was directly converted by the LOS deformation(dLOS/cos). And the accuracy of the estimated dynamic subsidence was very low only about 20mm. At present, related studies of combining InSAR and logistic model for monitoring coal mining subsidence are seldom, its applicability in different coal mining areas needs to be further tested.

C. InSAR AND EXPONENT KNOTHE MODEL
The exponent Knothe model is also a commonly used coal mining subsidence prediction model and has been used in the coal mining areas by some scholars [72] [73]. The exponent Knothe model only has three unknown parameters, much less than the PIM. The model function as follows.  [75]. They proposed two methods: (1) Using DInSAR to obtain the cumulative deformation of points and construct the relationship between the cumulative deformation and the model parameters (as (13)); (2) Using time series InSAR to obtain the relative deformation of points and construct the relationship between the relative deformation and the model parameters (as (14)). Then taking the least square method as the parameter estimation method, two sets of exponential Knothe model parameters were estimated by using some cumulative deformation values or some relative deformation values. In the end, they predicted the large gradient subsidence by the exponent Knothe model. Experiments show that the model parameters estimated by cumulative deformation were more accurate. Another experiment in Huainan coal mining area obtained the maximum 930 mm subsidence [74] [75]. As shown in Fig. 5 Where is the cumulative LOS deformation of point p at time , ∆ is the relative LOS deformation between time ti and tj, is the radar incidence angle, is the possible maximum subsidence; p, is the influence coefficient about the overlying rock stratum of the coal seam, is the fitting parameter.Although the subsidence curves obtained by these two methods were both consistent with the leveling data, and the maximum error were 20mm and 21mm. These methods had the same problems as other methods: the horizontal deformation was also ignored, and the LOS deformation was directly converted into the vertical subsidence. In the experiments, only six Leveling points were used to validate the fitting result and the prediction accuracy of exponential Knothe model, the number of feature points was too limited and more points were needed. At present, only L. CHEN et al. studied the large subsidence monitoring by combining InSAR and the exponential Knothe model in the coal mining area, and only used Radarsat-2 one type SAR images and experimented in one coal mining area [74] [75]. More types of SAR images in different coal mining areas should be tested in the future.

D. A SUMMARY
In summary, InSAR and PIM, logistic model, and exponent Knothe model, all can be used to estimate the large subsidence in the coal mining areas, the processing steps of these methods are similar (As shown in Fig. 6), and have achieved satisfactory results in practical applications. These methods greatly improved the gradient of InSAR deformation monitoring ability. However, there are still some shortcomings to be overcome.
• There is no doubt that there is horizontal movement on the surface during the ground deformation. However, these methods usually ignored the horizontal movement. Perhaps the horizontal deformation has been converted into the vertical deformation by mistake, which will increase the error of the estimated vertical deformation. Neglecting of the horizontal deformation also results in some parameters related to the horizontal deformation (such as the horizontal deformation coefficient) to be incorrectly estimated [15].
• These methods usually directly convert the LOS This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. deformation into vertical subsidence by means of over the cosine of the incident angle. However, the surface subsidence process in coal mining area is very complex and nonlinear [76]. This direct conversion may affect the accuracy of the estimated model parameter, and then affect the accuracy of estimated subsidence.
• The starting time of the model prediction is the day when the surface deformation starts. In many cases, the acquisition time of the first SAR image is not that day. Therefore, it is necessary to convert the two into the same day according to certain methods. This conversion will affect the accuracy of model parameter estimation, and then affect the estimation result of the subsidence.
• Compared to the conventional GNNS and leveling methods, the accuracy of the estimated subsidence results was still low and needs to be further improved.
• Compared to the logistic model and exponential Knothe model, the estimated results of PIM were more satisfactory. However, there are too many parameters to be estimated in the PIM, which requires a lot of calculation and time.
• The practical application of these methods in coal mining subsidence monitoring are still rare and they need to be tested with different types of SAR images in different coal mining areas.

IV. SAR IMAGE AMPLITUDE OFFSET TRACKING METHOD
In 1998, Gray et al. proposed an image registration method to estimate ice motion in both range and azimuth. Supposing that the orbits of the two observations are parallel, and the images are formed at zero Doppler points. The relationship between the surface deformation in ground range( ) and azimuth( ) directions and the pixel shifts in image slant range ( ) and image azimuth ( ) is as follows [77] : where is the baseline, is the baseline angle, is the radar look angle at the satellite, , are the terrain slopes in range and azimuth respectively (to the local horizontal), is the local incidence angle.
Later, many scholars applied this method to monitor the displacement due to earthquakes [78] [79], volcanos [80], landslides [81] [82], and the movement of glaciers [83] [84]. There are usually two kinds of offset tracking methods: (1) amplitude tracking method, (2) phase coherence tracking method [85] . The amplitude information or phase interference fringes of a pair of SAR images were analyzed to find out the best matching homogenous pixels of image pairs based on normalized cross-correlation maximization, then extracted the pixel offset between the corresponding pixels of the master and slave images in the azimuth and range direction. According to the pixel resolution, after removing the orbital and terrain offsets, the two-dimensional deformation of the ground along range and azimuth can be obtained. Phase coherence tracking method uses the phase information of the image, and the amplitude information is often ignored. Most of the coal mining areas are located in the suburbs or rural areas, and the ground surface is often covered by vegetations. Vegetations change every day. And the coal mining ground surface deformation range in horizontal direction is often small but the vertical subsidence gradient is large. Both vegetation's change and large gradient subsidence will lead to poor phase coherence or even incoherence of the image pairs. So, phase coherence tracking method often does not work well. The amplitude information offset tracking method does not need interferometric operation or phase unwrapping, and the high or low coherence of SAR image pairs has no effect on it. As long as there is obvious and comparable intensity information and the corresponding point pairs can be well matched, the deformation offset field can be obtained. Therefore, most of the existing researches on this technology are based on the amplitude information [86]. It overcomes some limitations of InSAR, such as relying on phase interference, insensitive to azimuth deformation and only one-dimensional deformation in LOS [83] [87]. Many scholars have applied the amplitude offset tracking to monitor coal mining large gradient subsidence. Some methods are described and discussed below in detail. Zhao et al. firstly applied this technology to monitor the large gradient subsidence in a coal mining area and obtained a maximum 4.5m vertical subsidence in Bulianta by using ALOS PALSAR data, and the accuracy was 1/25 pixels(about 18.8cm) [88]. Obviously, the accuracy is very low. If the subsidence is less than 20cm, the measurement results cannot be distinguished from measurement error, the small subsidence points information will lose.
In order to obtain both large and small subsidence information, some scholars proposed mixed methods: conventional DInSAR or SBAS-InSAR + conventional offset tracking method [42][5] [89] [90], SBAS-InSAR + SBASoffset-tracking method [40]. The small subsidence points were mornitored by conventional DInSAR or SBAS-InSAR, whereas those large subsidence points in the center of the subsidence basin were estimated by conventional offset tracking or SBAS-offset-tracking method. The SBAS-offsettracking method followed the basic strategys of SBAS-InSAR. It expanded offset tracking method from the conventional single pair of images to the time series measurement category. Its general process was as follows [91]. Grouped the SAR images into some small baseline sets according to the temporal or spatial baseline. Obtained the slant-range offset value ∆ and azimuth offset value ∆ of SAR image pairs by conventional offset tracking. Then estimated the average offset rate ̅ and ̅ between the adjacent SAR images by (16) according singular value decomposition. Finally, due to the nonlinearity and discontinuity of the coal mining subsidence in time and space, time dimension integration was used to get the final time series deformation.
where B is a SAR image coefficient matrix. Experiments of these mixed methods showed that the subsidence results were consistent with the field measurement data in the edge of subsidence basin, and the more reliable meter level large subsidences(about 4m) were obtained (As shown in Fig. 7(a)) [40].
However, above related mothods did not improve the monitoring accuracy obviously and the offset errors caused by topography, ionospheric anomalies and others were often ignored [40] [92]. Aiming at improving the low monitoring accuracy, some scholars constructed the offset error model to reduce different errors. After image registration, the globle offset of the same name pixel of the master and slave images mainly consists of many components [93] [94]. (17) where Odef, Oorbit, Otop, Oiono, Onois were the offset caused by deformation, orbit location and gesture difference, topographic relief, ionospheric, noise respectively.

Ototal=Odef+ Oorbit+ Otop+ Oiono+ Onoise
The orbit offset Oorbit was the most important error and that needs to be firstly removed. It was often considered as a low frequency signal, a systematic error, and could be described approximately by a polynomial model using no displacement points outside the deformation area in image pairs [95] [96].
Many scholars proposed different offset error models according to the characteristics of the interesting area and purposes. Such as, ignored Otop, Oiono ,Onois and constructed a bilinear polynomial to estimate the orbit offset errors (as (18)) [97]; considered the topographic relief would cause large error, added a topographic relief model to the range offest polynomial(as (19)) [98]; ignored Oiono and Onois, and constructed a quadratic polynomial model to remove both the orbit and topographic relief offset error(as (20)) [99]. c0c9, m0,m6 were the coefficients, r and a were the row and column numbers on the image in range and azimuth, H is the height of the point, H is the average height.

Oorbit+top=m0+m1×r+m2×a+m3×r×a +m4×r 2 +m5×a 2 + m6(H-H)
Huang et al. firstly used a locally adaptive template window size to improve registration accuracy, then used a quadratic polynomial (like the first six terms in (20) to remove the offset errors caused by including ionospheric effects, sensor attitude and partial topographic relief based on the stable points, and used another second-order polynomial(as (21)) to remove the residual topographic relief offset error by using external DEM [100]. The experimental results showed that the accuracy was improved by 3~4mm [100]. But this method required much time to find the optimal window, its calculation burden was too heavy, the efficiency was low.

Ores=b0+b1(H-H)+b2(H-H) 2 ( 21)
In order to improve the monitoring accuracy and obtain time-series 3D large deformation induced by coal mining, some scholars proposed the method of combining offest tracking with subsidence model. Yang et al. proposed single imaging geometry offest tracking (OT-SIG) method. Its process was roughly as follows [101]. Firstly, formed small temporal baseline and spatial baseline SAR amplitude pairs, obtained the LOS and azimuth deformations dLOS(i, j), dAzi(i, j) at a surface point (i, j) of a SAR amplitude pair in conventional offset tracking method. Then a priori model was used to construct the relationship (as (22)) between the offset deformation and the real surface 3D deformation in vertical, east, and north directions W(i, j), E(i, j), N(i, j). Finally, following a similar to SBAS-InSAR method to estimate the time-series 3D deformation. where ( , ) and ( , ) are the coefficients of the gradients of vertical subsidence, and are the spatial resolutions of offset tracking deformation map in east and north directions, is the mining depth, is the influence angle; is the horizontal movement coefficient, is the incidence angle, is the azimuth angle.
Using this method with 13 scenes about 1 m resolution spotlight TerraSAR-X amplitude images, the maximum more than 4000mm and 1400 mm deformation in vertical and north directions were obtained in Daliuta coal mining area (As shown in Fig. 7(b)(c)). The Root Mean Square error (RMSE) were 220mm and 110mm respectively [101].
Yang et al. assumed that the horizontal motions were proportional to the vertical subsidence and the 3D deformation was linear between two time-adjacent SAR images, and considered some model parameters, such as and , as constants [101]. It was obviously inconsistent with the practical subsidence characters in coal mining areas [102]. With consideration of the above problems, B. Chen et al. proposed an offset tracking measurements and clastic medium theory(OT-CMT) method [102]. This method constructed a three-dimensional deformation model of coal mining based on CMT theory and Knothe's theory, obtained the time series LOS deformation values of the ground points by offset tracking method. Then, the unknown parameters of the above model were estimated by these values. Ultimately, the 3D deformation was estimated by the model [102]. Using this method and the same SAR images as [101], B. Chen et al. obtained the maximum 4.43m subsidence and the maximum 1.4 m and 1.2 m horizontal movements in north and east directions respectively. The RMSEs were 74mm, 83mm and 93mm in the vertical, east and north directions [102]. Compared with [101], the accuracy was improved.
The results of both [101] and [102] can meet the coal mining deformation practical measuring requirements. But the applications of these methods are limited, because they are dependent on a specific subsidence model adapted to a specific coal mining area. The coal mining subsidence process is very complex, the applications in other coal mining subsidence cases should be further studied.
The amplitude information offset tracking method has obvious advantages in monitoring large gradient deformation. It does not use phase information. It also does not require phase unwrapping operations. It is also less affected by the temporal and spatial decorrelation problem. Although it is far inferior to InSAR in detecting small deformation, it can detect large deformations that exceed the maximum DDG of InSAR. Furthermore, it can obtain the deformation both in range and azimuth directions simultaneously rather than InSAR can only obtain one-dimensional deformation in LOS direction.
The offset tracking method makes up for the deficiency of InSAR method based on phase interferometry in monitoring large gradient deformation, and has been used by some scholars for monitoring large gradient subsidence in coal mining areas. However, the applications showed that the amplitude information offset tracking method still had some shortcomings to be overcome.
• The monitoring accuracy is still relatively low. It is directly proportional to the pixel resolution of the SAR image, and the theoretical accuracy can only reach 1/30 of the pixel resolution [103]. The highest resolution of a SAR image is 1 meter, therefore the theoretical maximum deformation accuracy of the offset tracking method can reach 3.33 cm. However, at present, the resolution of SAR images commonly used is less than 10 meters, that is to say, the deformation monitoring accuracy will be lower than 33.3cm using such SAR images to monitor large gradient deformation, which is far from enough for monitoring subsidence in coal mining areas.
• More accurate registration methods are needed. The core of the offset tracking method is to find the maximum value of the intensity normalized cross-correlation coefficient in search window. The current registration accuracy can only reach the sub-pixel level, and it is difficult to achieve 0.01 pixel registration accuracy, which makes it is difficult to extract more accurate surface deformation information [104], and it cannot accurately show the surface deformation [94] .
• Choosing an appropriate cross-correlation calculation window. Because the cross-correlation coefficient is dependent on search window size, the registration accuracy of the the corresponding pixels is sensitive to search window size. If the search window is too small, a large number of points will be mismatched due to the reduction of sampling information and noise, the matching reliability will be reduced [94] [105]. Theoretically, the larger the search window, the higher accuracy. But a larger template size results in greater computation time and deformation distortion [106]. The experiments show that the accuracy firstly increases with the increase of the search window, and then decreases with the increase of the search window. If the search window is too large, the deformation of the bottom of the subsidence basin will be compressed, and the larger the search window, the "compression" phenomenon becomes more serious [107]. How to choose a suitable cross-correlation window is one of the key factors to improve the monitoring accuracy of this offset tracking method.

V. COMBINING WITH TERRESTRIAL 3D LASER SCANNING
The working principle of terrestrial 3D laser scanning is similar to that of SAR: it actively mits electromagnetic waves, receives the reflected signals from the ground objects. The terrestrial 3D laser scanning can get high-density and highprecision point cloud. Under the condition of close observation (less than 100 meters), the distance accuracy between any two points and the position accuracy of a single point both can reach millimeter level[108] [109]. And the obtained point cloud data is a region rather than a series of feature points which can quickly and comprehensively show the surface morphology of the subsidence area. Many scholars have studied the use of terrestrial 3D laser scanning to monitor mining subsidence and have achieved satisfactory results [110][111][112]. The accuracy reached millimeter level [113] [114]. In view of the deficiency of InSAR in monitoring large gradient subsidence and the limited measurement range of terrestrial 3D laser scanning, some scholars have carried out relevant researches on the combination of InSAR and terrestrial 3D laser scanning to monitor large gradient mining subsidence. Generally speaking, the process is as follows (As shown in Fig. 8). InSAR methods were used to obtain the small subsidence around the subsidence basin, and generated the subsidence field, called InSAR subsidence field. The large gradient subsidence was obtained by terrestrial 3D laser scanning method, and generated the LiDAR subsidence field. Note that the 3D laser scanning data sets should cover the large gradient subsidence area completely and had the same time intervals as the SAR images. In the end, some methods were used to merge the two subsidence fields into a complete subsidence basin. For example, an inverse distance weighting method was proposed to merge the two subsidence fields.
The subsidence values of the central points in InSAR subsidence field were obtained by interpolation from the corresponding points in the LiDAR subsidence field. The farther the LiDAR point was from the InSAR subsidence field point, the smaller contribution to the subsidence value of the InSAR subsidence field point, and vice versa [115]. Kai et al. used the GNSS control points which were both contained in InSAR subsidence field and LiDAR subsidence field as the bridge to correct and merge the two monitoring results [116]. The experiments showed that the integrated subsidence fields were consistent with field measurement, and the maximum detectable subsidence reached more than 1250mm and 6000mm (The mining height is 8.6m) respectively in Baodian and Shangwan coal mining areas (As shown in Fig. 9(a), (b)). The maximum absolute error was 64 mm, and the mean absolute error is only 23.1 mm in Baodian [115] [116]. Its performance was satisfactory.
It can be seen that the combination of InSAR and 3D laser scanning can monitor the whole subsidence basin including central large gradient subsidence, and the accuracy is also very high. But the related researches are still very few, and there are still many problems to be solved.
• Using 3D laser scanning point cloud data to fill the "hole" of InSAR subsidene field. The key is the mosaic problem between the InSAR subsidence field and the 3D laser scanning point cloud subsidence field. There are no accurate and efficient methods to achieve this.
• SAR image acquisition time is not completely consistent with that of 3D laser scanning observation. Current researches and applications do not consider the time nonsynchronous problem.
• 3D laser scanning cannot obtain the horizontal movement information. What 3D laser scanning obtains is the digital surface model, it is easy for it to obtain the subsidence surface, but it is difficult to find the same feature points in different periods, it is difficult to obtain the horizontal movement information of the surface.
• Because of the limitation of measurement distance, in order to observe the whole subsidence basin, it is necessary to set up 3D laser scanner at different observation stations, and then to splice different point data sets together. Multi observation stations acquisition and point cloud data splice both will increase monitoring error.
• Data missing in some areas. For some reasons, such as 3D laser scanning light is blocked by terrain, objects and uneven distribution of observation stations. Some areas will not be scanned. In practice, interpolation is generally used to fill the missing data, which also leads to error from the actual situation.

VI. CONCLUSIONS AND FUTURE PROSPECTS
A. CONCLUSIONS InSAR techniques have many advantages for monitoring deformation. Since CARNEC et al (1996) initially used InSAR for monitoring surface subsidence caused by underground mining [117], InSAR has become an powerful tool for monitoring ground subsidence in coal mining areas and has been widely used. Although InSAR technology has limitations in terms of maximum detectable deformation gradient, in order to obtain the whole coal mining subsidence basin (including the large gradient subsidence in the middle area) and its subsidence process, scholars have proposed more than 20 different methods to monitor large gradient subsidence using SAR images (this number will certainly increase in the future). Some satisfactory results have been achieved in practical applications, and some of them can be improved further and used to monitor large gradient subsidence in coal mining areas in the future. In this paper we focus on summarizing the methods of monitoring large gradient subsidence in coal mining areas using SAR images. In general, these methods can be grouped into three categories: (i) Combining InSAR with coal mining subsidence prediction models; (ii) Using SAR image amplitude information offset tracking technology; (iii) Combining InSAR with terrestrial 3D laser scanning technology. All these methods could obtain large subsidence in coal mining areas, where maximum subsidence was more than 1000mm, even 6000mm. Among various methods proposed, researches related to combining InSAR with coal mining subsidence prediction models are in the majority, and the accuracy of the prediction results is also enough in practical applications. SAR image amplitude information offset tracking method does not use phase information, therefore it does not need phase unwrapping operation, and does not consider the influence of decoherence. It is simpler than InSAR. It has a strong ability to monitor large gradient deformation, and can obtain deformations both in the range and azimuth directions simultaneously. Combining InSAR and the offset tracking method can obtain both large and small subsidence information of a whole subsidence basin. The SBAS-offsettracking method expanded the offset tracking method from the conventional single pair images to time series images. The combination of InSAR and terrestrial 3D laser scanning technology explores a new idea of merging two different types of surface data. It can obtain relative high accuracy in the whole subsidence basin.
Despite the improvements that have been achieved in monitoring large gradient subsidence in coal mining areas using SAR images, there are still some shortcomings to be overcome and some room for more improvement.
1)The horizontal movement was often ignored or was difficult to obtain. Perhaps the horizontal deformation may have been converted into the vertical deformation by mistake.
2)The surface subsidence were complexity and nonlinearity, but the vertical subsidence was often obtained by directly dividing the InSAR LOS deformation by the cosine of the incident angle. Perhaps the result may not represent the true vertical subsidence.
3) Time nonsynchronous problem. The subsidence start time used in the subsidence prediction model, or 3D laser scanning field measurement time, was usually inconsistent with the acquisition time of the first SAR image. When they are combined, it is necessary to convert the two different times into the same day using to certain methods. This conversion may affect the result of the subsidence. 4) The current researches have only conducted experiments on a certain goaf, and extensive experiments in various different types of mining areas have not conducted. Whether they have wide applicability needs to be further verified. 5) Compared to the conventional GNNS and leveling methods, the accuracy of the estimated subsidence results is still low. Especially SAR image amplitude information offset tracking method, it can not meet the accuracy requirements in practical applications, and cannot be used alone for monitoring coal mining subsidence.

B. FUTURE PROSPECTS
After years of development, SAR has developed from low resolution, single polarization, single incident angle, single satellite and long revisit to high spatial resolution, high temporal resolution, full polarization, multi-mode, multiincident angles, multi satellite constellation, multi-channel and 3D imaging. The earth observation capability of SAR continues to improve. SAR has also become a high efficiency method for monitoring large gradient subsidence in coal mining areas, insteading of the conventional topographic survey methods which are time-consuming and laborious and can only monitor sparse points. If there is a widely applicable method for monitoring large gradient subsidence in coal mining areas using SAR images, it will greatly save time and effort. However, those proposed methods for monitoring large gradient subsidence in coal mining areas using SAR images are all tested in a certain coal mining area or a coal mining goaf. Although some satisfactory results have been achieved, none of them have been tested in other coal mining areas or coal mining goafs. Therefore, it is necessary to apply the existing methods to more different types of mining areas or goafs to verify whether they have universal applicability.
In addition, only one mode SAR images from a certain SAR sensor were used in the related experiments in each of the existing methods. There are many different SAR sensors currently in orbit, and there will certainly be more in the future. They fly in different directions and emit different frequency electromagnetic waves in different incident angles, observe the ground from different directions, and finally obtain the different LOS deformation. By combining the multiple different LOS deformations, we can more accurately obtain the ground deformation in the 3D direction, and overcome the limitation of obtaining 3D deformation from single LOS deformation. Therefore, it is very important to construct the relationship between multiple LOS deformation and the actual ground 3D deformation.
The more information about the subsidence area we get, the better we can understand the situation related to the subsidence. Even the same subsidence area, the images from different frequency SAR sensors may provide different information. The difference will be more between SAR sensors and optical remote sensing sensors. How to comprehensively utilize the information of the same subsidence area obtained by different sensors (including SAR sensors of different bands and optical remote sensing sensors) is also a practical research direction. And now, SAR data have become massive data, therefore, artificial intelligence technology, big data cloud processing technology, and information mining technology have increasingly become useful tools.
Both amplitude information offset tracking method and combining InSAR method are affected by SAR image resolution and registration accuracy. Although SAR image resolution and registration accuracy are constantly improving, the current situation cannot meet the actual needs of monitoring large gradient subsidence in coal mine areas. Higher resolution SAR image and higher precision SAR image registration method are expected in order to improve accuracy.
Phase incoherence and phase unwrapping are still the biggest difficulties for InSAR technology to monitor large gradient subsidence, and better methods are needed to solve these problems.