Mining Subsidence Monitoring With Modified Time-Series SAR Interferometry Method Based on the Multi-Level Processing Strategy

The supply of coal is related to the stability of social development in China. Coal mining faces are generally located under farmland and wasteland to avoid destroying the livelihood and property of residents, which seriously affects the selection of permanent scatterer (PS) with time-series interferometric synthetic aperture radar (TS-InSAR) method. The advanced TS-InSAR (ATS-InSAR) method combining PS and distributed scatterer (DS) monitoring modules is a useful tool to increase measurement points. However, the number of DS and the accuracy of measurement are different with different thresholds of temporal coherence. On the basis of ATS-InSAR, we propose a modified method applying the multi-level processing strategy to obtain more reliable deformation information in this paper. 16 sentinel-1A images are used to monitor mining subsidence in Peixian, China. The results of three different methods are cross-verified. Meanwhile, reliable DS pixels are identified by using the thresholding of both temporal coherence and Pearson correlation coefficient through the hierarchical processing. The results show that the deformation of the modified method reveals a large subsidence with the maximum rate of −563 mm/yr. The number of measurement points selected by the modified method is about 6.6 times that of the TS-InSAR method, and 1.3 times that of the ATS-InSAR method. The modified strategy can extract a great number of reliable pixels and reduce error propagation to ensure measurement accuracy. This research offers information to relevant departments for risk management purpose.


I. INTRODUCTION
Coal is the most abundant and widely distributed conventional energy in the world. In recent years, the rapid development of electric power, building materials, chemical industry, and other industries has led to a substantial increase in the demand for coal in China [1]. However, coal mining is easily to cause land subsidence, which brings geo-hazards and structural damage to people. Synthetic aperture radar (SAR) system is a remote-sensing imaging system working in microwave band, which can continuously observe the Earth's The associate editor coordinating the review of this manuscript and approving it for publication was Cheng Hu . surface over a large-scale area with a fixed revisit period [2]. Many domestic and foreign scholars have applied time-series SAR interferometry (TS-InSAR) techniques to monitor mining subsidence [3]- [10].
TS-InSAR techniques that are aimed to detect point-like targets with high reflectivity and slight influence of temporal and geometric decoherence, generally include two types, permanent scatterer (PS) interferometry with single master image and small baseline subset (SBAS) interferometry with multiple master images [11]- [18]. In addition, a new method named Multiple-master Coherent Target Small-Baseline InSAR (MCTSB-InSAR) was presented in [19] and [20], which integrated the advantages of and PS-InSAR and VOLUME 9, 2021 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ SBAS-InSAR and proved to be a useful tool for measuring ground displacement. Nevertheless, it is difficult for these methods to obtain enough measurement points (MP) and deformation information of the mining area. Advanced time-series interferometric synthetic aperture radar (ATS-InSAR) method refers to those TS-InSAR methods that combine the PS and distributed scatterer (DS) monitoring modules [21], [22]. As we all know, the SqueeSAR algorithm [23] was put forward to select MP correspond to image pixels belonging to the available dataset of no less than moderate coherence. Deformation estimation with this method is the result of processing the selected DS jointly with the PS. Since then, SqueeSAR has gradually become a research hot topic in the field of ATS-InSAR. Most notably, the MP with the number of homogeneous pixels greater than a pre-set value are regarded as the DS candidates, and a certain threshold of the temporal coherence κ is utilized to extract reliable DS afterwards. Two key steps are statistically homogeneous pixel (SHP) identification and phase optimization.
Previous studies showed that fast statistically homogeneous pixel selection (FaSHPS) algorithm [24], [25], transforming the hypothesis testing problem from significant difference judgment to confidence interval estimation, had the advantage of high computational efficiency compared with the hypothesis test. Besides, the interferometric phases reconstructed by applying the phase triangulation (PT) algorithm produce less noise than the spatially filtered interferometric phases [23], [26]. Due to the larger weight on phases with higher temporal coherence, the coherence-weighted PT algorithm is more robust to decorrelation than the equal-weighted PT algorithm [26]. Therefore, FaSHPS and coherence-weighted PT algorithm have an advantage in phase optimization.
It is noteworthy that the number of DS and the value of deformation are different with different temporal coherence thresholds. To be specific, the larger the value of κ is, the smaller the number of DS is, which cannot provide more details for deformation. The smaller the value of κ is, the larger the number of DS is, which has a negative impact on the measurement accuracy. The focus of our research is how to select as many DS pixels as possible from moderate coherence pixels without loss of measurement accuracy.
To increase spatial density of deformation measurements, we propose a modified ATS-InSAR method based on the multi-level processing strategy. Due to the Pearson correlation coefficient (PCC) working well for measuring the data correlation in time and space [27], PCC is used to ensure the reliability of pixels with moderate coherence. To reduce error propagation, the hierarchical strategy is utilized to estimate deformation at reliable pixels.
In this paper, we test the modified ATS-InSAR method for subsidence monitoring over Peixian in China using 16 Sentinel-1A images. Furthermore, the monitoring results derived from different methods (TS-InSAR, ATS-InSAR and modified ATS-InSAR method) are applied for comparison, verification and analysis.

II. STUDY AREA AND DATASETS A. STUDY AREA
The study area is located in Peixian, Xuzhou City, China (see Figure 1). The climate is hot and rainy in summer, and the water resources are abundant. Eight coal mines of Peixian are underground covered coalfields, one of which named ZSL (Zhangshuanglou) mine. The land cover classification of study area mainly includes water bodies (like Dasha River, Anguo Wetland), vegetation (like farmland, grassland), bare land, buildings and structures. The location of mining area is around the working face (NO.92606 and NO.74101).  Table 1. After multi-looks and fine registration, 800 × 800 pixels were cropped. Additionally, 30 m resolution ALOS DEM data was used to generate differential interferogram through removing the topographic phase contribution.

III. METHODOLOGY
Time-series SAR images with size N are used for fine registration, and M interferograms with free combination are generated according to the pre-set perpendicular and temporal baseline. The modified strategy is implemented by applying hierarchical processing, which has the advantage that high-quality MP cannot be affected by low-quality MP. For the selection of reliable PS, the method is adopted refer to [19], [20]. The amplitude dispersion index and mean coherence are utilized to extract reliable PS. This section addresses DS selection and phase optimization, as well as deformation estimation at reliable pixels.

A. SHP IDENTIFICATION
The theoretical basis of SHP identification is that homogeneous pixels have similar statistical characteristics in time series. Assuming that there are a large number of distributed targets in a pixel of a SAR image, the amplitude A(p) for pixel p has a Rayleigh distribution with mean u(p) and the variance Var(A(p)) for single-loook image [24]. On the basis of the central limit theorem, the average amplitudeĀ(p) of SAR images with size N is supposed to follow a normal distribution with mean u(p) and the variance Var(A(p))/N .
The FaSHPS algorithm converts the problem of hypothesis test from the significant difference judgment to the confidence interval estimation under this assumption. As long as the confidence interval is determined, the similarity between pixels can be judged by simple logical operations. Based on the definition of the standard normal distribution and the connection between the population mean and variance, the confidence interval ofĀ(p) is then given by: where P{·} stands for the probability, and Z 1−α/2 is the 1−α/2 percentile of the standard normal probability density function. We find that the pixels falling into the interval of above formula are homogeneous with the pixel p when the significance level α is given. Furthermore, pixels with SHP number greater than 20 are regarded as the DS candidates [21], [23].

B. PHASE OPTIMIZATION
The coherence-weighted PT algorithm based on coherence matrixes is adopted for phase optimization [26]. The reason why we work on coherence matrixes rather than covariance matrixes is that it can be beneficial to compensate for possible backscattered power unbalances among all the images [23]. The coherence matrix T is calculated as follows: where E[·] represents the expectation. Y = [y 1 , y 2 , · · · y N ] T is the normalized complex observation of N SAR images, and H denotes the conjugate transpose. stands for the set of homogeneous pixels containing N P adjacent pixels with similar backscattering properties. The matrix T can be further expressed as: whereγ m,n and e jφ m,n respectively represent the estimated coherence and interference phase between the mth image and nth image. • indicates the Hadamard product. |T |, the modulus of T , is an N × N matrix with elementγ m,n , and is an N × N matrix with element e jφ m,n . During a series of mathematics deduction, the coherence-weighted phase is derived as follows: The ''goodness of fit'' measure κ, firstly introduced in [23], has been used as quality indicator for the optimal phase θ. Indeed, κ can be considered as the temporal coherence of DS [26], [28]. It is an important parameter for DS selection. The original phase is compared with optimal phase based on the following formula.
where Re(·) represents the real part of complex data. The larger the value of κ is, the higher the coherence and the signal-to-noise ratio are. Conversely, pixels with smaller value of κ are considered more serious affected by speckle noise. Pixels with temporal coherence higher than a certain threshold are extracted, and the corresponding phase values of the original SAR images are substituted with their optimized values in the process of phase optimization.

2) PEARSON CORRELATION COEFFICIENT
Pearson correlation coefficient (PCC) is the degree of correlation between the variables. PCC can be used to reflect the level of spatially uncorrelated noises and evaluate the quality of phase values in the time domain [27], [29]. For two pixels i and j, the PCC is expressed as: where ϕ s i and ϕ s j are the phase of i and j in the sth differential interferometry.φ i andφ j are the mean differential phase VOLUME 9, 2021 of i and j in interferograms, respectively. The value of ρ is between 0 and 1. The larger the ρ is, the stronger the correlation between i and j is. Thus, we used the PCC index to measure the relative phase quality of pixels. The connection between the variation of PCC and the correlation of phase data was discussed in [27], as shown in Table 2. It can be known that the pixel of i is strong correlative to the pixel of j when 0.5 ≤ ρ i,j ≤ 1.

3) MULTI-LEVEL PROCESSING
The analysis of temporal coherence can result in the success to identify a useful DS with the temporal steadiness of radar reflectivity [23]. As mentioned in the first section, the number of DS varies with the temporal coherence threshold. All the DS candidate pixels are further treated through the multi-levels of processing to reduce the possibility of error propagation. Pixels with 0 ≤ κ < 0.5 are considered as unreliable pixels on account of the low coherence, which can reduce the accuracy of the deformation estimation to a great extent. There is no denying that pixels with large value of κ possess satisfactory signal-to-noise ratio, which is our basis to regard the pixels with 0.7 ≤ κ ≤ 1 as reliable pixels. Moreover, pixels with 0.5 ≤ κ < 0.7 are defined as MP of moderate temporal coherence. To reduce the memory cost, these pixels are divided into two classes at equal intervals instead of multiple classes in this work. In view of this, we classify all the DS candidate pixels into four groups, namely pixels with temporal coherence values of [0.7, 1], [0.6, 0.7), [0.5, 0.6), [0, 0.5).
It is necessary to discard unreliable pixels with 0 ≤ κ < 0.5 before the process of triangulated network construction. Due to the time-consuming calculation, the spatiotemporal correlation analysis with the PCC threshold is only applied to select pixels with 0.5 ≤ κ < 0.7 in this work. The PCC values of [0.5, 1] mean the strong correlation between the phase of two adjacent pixels [27]. Therefore, we set the PCC threshold to 0.5 to ensure the reliability of pixels with moderate temporal coherence. Once the reliable DS pixels are selected, the phase values of the original SAR images are substituted with their optimal values. The DS pixels with estimated optimal phases are equivalent to quasi-PS pixels [21].

D. ESTIMATING DEFORMATION RATES AT RELIABLE PIXELS
The reliable DS and PS pixels are processed using the MCTSB-InSAR algorithm for the deformation of each measurement point. The local Delaunay triangulation is constructed to connect the selected reliable pixels for phase analysis. The weighted least square algorithm is applied to calculate the absolute value of the DEM error and linear deformation rate for each reliable pixel. After determining the linear deformation, the nonlinear deformation is estimated by temporal and spatial filtering of residual phase [19]. The deformation estimation of PS pixels is called TS-InSAR here. The process of extracting the reliable DS pixels from the temporal coherence obtained by the coherence-weighted PT algorithm, jointly with the PS for the deformation estimation, is defined as ATS-InSAR. For better understanding, Figure 2 shows the flowchart of the hierarchical processing with modified ATS-InSAR method.  Figure 3(a-b) is the average amplitude and classification map in the study area, which is only the basis for subsequent analysis. The area of buildings and structures accounts for a small percentage (approximately 10%), but the area proportion of vegetation coverage is the largest (nearly 53%) in the study area. Figure 3(c) is SHP number when α = 5%. The value of each pixel represents the number of SHP identified in a 15 × 15 search window. SHP number is 0, which demonstrates that there are no homogeneous pixels in the search window perhaps because of isolated PS points or noise points. SHP number is 224, which indicates that all the pixels in the search window are homogeneous. Figure 3(d) is the phase optimization quality evaluation map, which illustrates the temporal coherence value of each pixel.

A. SHP IDENTIFICATION AND PHASE OPTIMIZATION RESULTS
We find that SHP number of pixels is related to the type of land cover. SHP number of buildings and structures with  strong scattering is relatively few. Those targets whose neighboring pixels share similar reflectivity values in farmland imply a larger number of SHP. To reduce the error of SHP identification, the central pixel of search window would be considered as DS candidate pixel in our experiment if the SHP number is larger than 20. Besides, the temporal coherence value in the area of buildings and structures is large (close to 1), but it is just the opposite in the water area (less than 0.5). The proportion of pixels with temporal coherence greater than 0.7 (0.7 ≤ κ ≤ 1) is about 33%, and more than half of the pixels in the study area exhibit temporal coherence greater than 0.5 (0.5 ≤ κ ≤ 1). It makes the modified strategy to increase the possibility of the high-density MP.

B. DEFORMATION RATES ANALYSIS
A small baseline set of multiple-master images is generated under the design of temporal baseline smaller than 120 days and perpendicular baseline smaller than 150 m for images, respectively. 32 interferograms with less noise and clearer fringes are selected afterwards, as shown in Figure 4.
Many PS points are concentrated in areas where buildings and structures gather. There are few PS points in vegetation areas such as farmland and grassland, but no PS points in water area such as Dashahe River and Anguo Wetland. Only 30,146 points are shown in Figure 5(a) with the maximum VOLUME 9, 2021  deformation rate of −542 mm/yr. MP are sparse and unable to describe the deformation in detail. Figure 5(b) is the deformation result calculated by ATS-InSAR method. Reliable DS pixels are extracted from DS candidates with 0.7 ≤ κ ≤ 1. Figure 5(c) is the deformation result obtained by the modified ATS-InSAR method with multi-level processing strategy. As mentioned in the third section, DS candidate pixels with 0 ≤ κ < 0.5 are discarded. Three classes of pixels with temporal coherence values of [0.7, 1], [0.6, 0.7), [0.5, 0.6), respectively, are processed on a group-by-group basis, and the PCC threshold is set to 0.5.
The location of deformation in Figure 5(c) are consistent with those in Figure5(a-b). 150,528 points are visible in Figure 5(b) with the maximum rate of −555 mm/yr. There are 197,781 points in Figure 5(c) with the maximum rate value of −563 mm/yr. The density of MP in Figure 5(c) is about 6.6 times that in Figure 5(a), and 1.3 times that in Figure 5(b). More MP exert a beneficial impact on constructing a dense Delaunay triangulation for phase unwrapping, thereby reducing the influence of noise.
Most notably, the deformation caused by coal mining reaches a large magnitude in a relatively short period of time, which leads to the formation of water bodies in the south of the NO.74101 working face. The radar scattering echo signals of water bodies are so weak that it is difficult to show stable characteristics in repeated observations of SAR images. Thus, the coherence values in water area are relatively low and no reliable pixels are selected.
There are two obvious subsidence basins (ROI1 and ROI2) in the study area. The deformation of the area around the working face (NO.92606 and NO.74101) is larger than that of other areas. It is known that the coal faces of ZSL mine are located under farmland and wasteland. The echo signals of vegetation change greatly with time due to the scattering intensity greatly affected by seasonal changes and human factors (such as irrigation and farming) and the temporal coherence of some pixels near the working faces is less than 0.5, which has a very negative impact on the selection of reliable pixels in the subsidence basin. As a result, the deformation information is incomplete in the area near the working face. Even if that's the case, though, the number of MP in Figure 5(c) is more than that of in Figure 5(a-b).
From the results, the modified ATS-InSAR method with multi-level processing strategy has considerable potential and offers significant advantages for dealing with non-urban areas with few MP, which is effective in high density of deformation measurements. The modified strategy can extract more MP than ATS-InSAR method, especially in the north and east of the study area.

A. VERIFICATION OF MONITORING RESULTS
During the monitoring period, no leveling data are acquired. In order to compare the differences of monitoring results derived from different methods (TS-InSAR, ATS-InSAR and modified ATS-InSAR method), the correlation coefficient r is calculated from the deformation rates of all the completely coincident points. The absolute value of r (ranging from 0 to 1) measures the degree of close correlation between two methods [30], [31]. Assuming that v and u represent the deformation rates obtained by two different methods, their correlation coefficient r can be expressed as: where Cov(·) and D(·) represent the covariance and variance, respectively. The value of r is 0.903 in Figure 6(a), which indicates that Figure 5(b) is relatively close to Figure 5(a). We find that the correlation coefficient between Figure 5(a) and Figure 5(c) is VOLUME 9, 2021 0.872, which can be accepted in the experiment. That is to say, the monitoring results of modified ATS-InSAR method are in consistency with those of TS-InSAR and ATS-InSAR method.

B. DAMAGE CAUSED BY COAL MINING
To obtain the damage information caused by coal mining, the contour map and dangerous assessment are carried out based on the modified ATS-InSAR method. Referring to the existing literatures [32], [33], the grade of land subsidence dangerous assessment indexes is different in different areas and actual situation. In this study, we consider the area with deformation rate less than −10 mm/yr as dangerous area and then classify these dangerous areas into four groups. The contour map and area statistics are shown in Figure 7 and Table 3, respectively.
The total subsidence area (less than −10mm/yr) of ROI1 and ROI2 affected by mining activities accounts for approximately 2.47 km 2 and 3.62 km 2 , respectively. The dangerous area is situated in the subsidence basin and the very high dangerous area is close to the central area of the working face. Heavily populated areas, such as the factory near the NO.92606 working face and Wangtang village located to the north of the NO.74101 working face, have been seriously damaged, which brings hidden dangers to the safety of personal life and property. Farmland is the basis for the survival of residents, and its damage will bring inconvenience to economic construction. In view of this, the mining industry and related department should pay special attentions to deal with this matter.  In our experiment, 36,322 points of the ZSL mine within the study area are obtained by modified ATS-InSAR method. Although the study area we cropped do not include the western part of ZSL mine (see Figure 1), the comparison of our results with previous results is not affected due to the less contribution of the western part for mining subsidence of the whole mine. It proves that modified strategy has an advantage in the selection of MP.
This research represents the first investigation of mining subsidence over the ZSL mine using time-series SAR interferometry combined with PS and DS. There is detailed information for dangerous assessment and area statistics of the mining area. The potential danger of surface subsidence is analyzed during the monitoring period. All the information is useful to make the plan for land reclamation and ecological restoration.

VI. CONCLUSION
In this paper, the modified strategy is proposed to improve the ability of ATS-InSAR stability monitoring in an area of moderate coherence. Two important parameters, the temporal coherence and PCC, are applied to identify the reliable DS pixels. All the image pixels are classified into four groups by thresholding the temporal coherence values and processed group-by-group applying the multi-level processing strategy. 16 images of Sentinel-1A data are used to detect mining subsidence. The results of different methods (TS-InSAR, ATS-InSAR and modified ATS-InSAR) are cross-verified according to the deformation rates of all the completely coincident points. In addition, the dangerous assessment and area statistics are carried out for coal mine.
The results show that the deformation of the area around the working face is larger than that of other areas, and the location of subsidence obtained by the modified ATS-InSAR method are consistent with that obtained by the other two methods. The deformation rate map of the modified method derived from all the reliable pixels reveals a large subsidence with the maximum rate of −563 mm/yr. The correlation analysis of the completely coincident points indicates that the deformation rates estimated by the modified method are close to that calculated by the other two methods. Our experiments demonstrate that the modified ATS-InSAR method can increase the spatial density of deformation measurement with accuracy ensured.
Combined with the available mining information and monitoring results, we find that the very high dangerous area is located in the central area of the working face. These results provide the basis for taking corresponding prevention actions. In the future, we plan to use long time-series SAR images to enhance the real-time monitoring of mining subsidence. Further analysis will be conducted on these subsiding regions to offer more deformation information to relevant management departments.

AUTHOR CONTRIBUTIONS
Qian He and Yong Luo prepared data. Qian He designed and performed the experiments and wrote original draft. Yonghong Zhang and Hongan Wu provided technique support for MCTSB modeling. Yonghong Zhang and Hongan Wu gave comments and suggestions on the manuscript. Qian He revised the paper.
QIAN HE received the M.S. degree in geodesy and survey engineering from the China University of Mining and Technology, Xuzhou. She is currently pursuing the Ph.D. degree in photogrammetry and remote sensing. Her current research interests include InSAR technology and application, remote sensing image processing, and land subsidence monitoring.
YONGHONG ZHANG received the Ph.D. degree in photogrammetry and remote sensing from Wuhan University, China. He is currently a Research Fellow with the Institute of Photogrammetry and Remote Sensing, Chinese Academy of Surveying and Mapping, Beijing. His research interests include time-series SAR interferometry and change detection and disaster monitoring of remote sensing images.
HONGAN WU received the Ph.D. degree in geographic information system from the Chinese Academy of Sciences, Beijing. He is currently an Associate Research Fellow with the Institute of Photogrammetry and Remote Sensing, Chinese Academy of Surveying and Mapping, Beijing. His research interests include time-series SAR interferometry, SAR image processing, and information extraction.
YONG LUO received the M.S. degree in surveying and mapping from Liaoning Technical University, Fuxin, China. His research interests include geophysical simulation and InSAR technology and application.