Revealing Land Surface Deformation Over the Yineng Backfilling Mining Area, China, by Integrating Distributed Scatterer SAR Interferometry and a Mining Subsidence Model

Monitoring the surface deformation of filling coal mines regularly and understanding their spatiotemporal evolution characteristics for mining management and disaster warning is greatly significant. However, there is a lack of research on extracting spatiotemporal evolution characteristics of large-scale and high-resolution surface deformation in backfill mining areas and on evaluating the filling effectiveness of subsidence restraint. In this study, we took the Yineng Coal Mine in the Shandong Province of China and the surrounding area of the coal mine as the study area. The advanced distributed scatterer interferometric synthetic aperture radar (DS InSAR) technique was adopted for time-series analysis. The probability integral method (PIM) model for backfilling mining and an arctangent time function were integrated with DS InSAR to overcome the sparsity of InSAR observation points due to the temporal decorrelation caused by vegetation coverage. The results show that the proposed integration strategy is helpful in improving the number of effective monitoring points and obtaining complete spatiotemporal information of the surface deformation of the working face. The whole study area has eight key deformation zones. All six working faces in the Yineng initial minery display different degrees of deformation during the study period. The comparison between the deformation results of backfill mining and the simulated deformation of nonbackfill mining in the CG1312 working face shows that backfilling mining technology effectively reduces the range (22.00%) and magnitude (61.50%) of surface subsidence and significantly reduces the potential threat to surface buildings.


I. INTRODUCTION
C OAL resources provide energy security for rapid social and economic development in China [1]. On account of decades of large-scale development and utilization of coal resources, a series of ecological and environmental problems, such as land subsidence, ground infrastructure damage, and the pollution of groundwater resources have arisen. More severely, the aforementioned issues could be coupled with each other and thereby seriously damage the environment around the coal mine and create potential safety hazards for people's lives and property. Backfilling mining is an effective subsidence-reducing technology that fills mined-out areas with some noncombustible materials to prevent roof caving [2]. According to the survey data, 51 mining companies involving 75 coal mines in China adopted backfilling mining and successfully reduced subsidence until August 2018. The study area-the Yineng Coal Mine, Shandong, China (see Fig. 1)-has adopted backfilling mining in 2014. Compared with the total caving method, the effectiveness of subsidence reduction in backfilling mining reaches 60% according to the ground monitoring data. Nevertheless, the magnitude of surface subsidence (approximately 500 mm between January 2016 and July 2020) still cannot be ignored. Any acceleration of the surface subsidence would result in possible hazards. Therefore, it is necessary to monitor the surface subsidence in backfilling mining areas.
Many geodetic surveying techniques including ground-based and space-based methods have been used to monitor surface deformation induced by mining activity. Ground-based methods, such as leveling, global navigation satellite system methods, and geophysical investigation techniques, can provide highprecision measurements with dense temporal sampling [3], [4]. However, these measurements are sparse in space, and thus, it is extremely difficult to adequately characterize the spatial scale and pattern of ground deformation. The advent of the spaceborne interferometric synthetic aperture radar (InSAR) technique has brought unprecedented possibilities and convenience for observing land surface deformation with both high accuracy and unparalleled spatiotemporal resolution. In particular, the successful operation of the Sentinel-1 C-Band SAR satellite (from the Copernicus initiative of the European Space Agency) has escalated InSAR from a few, limited practices to plentiful, various applications in earthquakes [5], [6], landslides [7], [8], volcanic activities [9], [10], urban deformation [11], [12], and mining activities [13], [14]. However, despite its success in measuring surface subsidence associated with mining activities, InSAR is still limited in two ways [4], [15], [16], [17]: 1) the magnitude of surface subsidence usually reaches the submeter or meter level in a relatively short period, thereby introducing phase unwrapping error or even decorrelation; 2) most coal mines are situated in nonurban areas with temporally unstable scattering properties, consequently, obtaining sufficient coherent targets during interference tends to be difficult.
To overcome the first limitation, some studies have integrated InSAR and pixel offset tracking (POT) techniques to obtain the large-scale deformation of coal mines [13], [18], [19]. The deformation must be at least one-third of the pixel size that can be acquired by POT (approximately 80 cm in line-of-sight (LOS) of the satellite direction) [20]. However, the maximum accumulated subsidence monitored (approximately 160 mm in the vertical direction) was far less than the minimum theoretical observation value of POT over the period 10 January 2019-3 July 2020 according to the leveling data in the Yineng Coal Mine. Furthermore, backfilling mining tends to be able to largely reduce the surface subsidence in the Yineng Coal Mine (by approximately 60% according to incomplete statistics). Thus, it is theoretically more suitable for InSAR to measure the deformation in backfilling mines than unbackfilling mines under equal conditions. However, few studies have been carried out on this topic, and adding to the literature on this topic is one of the motivations of this study. Compared with addressing the first limitation (large deformation gradient), addressing the second issue (temporal decorrelation due to land cover) was a more severe challenge to obtain comprehensive information of deformation over the study area. Long-wavelength electromagnetic waves (L-bands) are indeed more able than C or X bands to retain the temporal coherence. However, applications of L-band InSAR in coal mines always depend on the availability of data. Persistent scatterer interferometry techniques, such as permanent scatterers InSAR (PSInSAR) [21] and the Stanford method for persistent scatterers (StaMPS) [22], [23], [24], were developed to reduce spatial and temporal decorrelation by identifying persistent scatterers (PSs). However, the density of PS pixels is typically limited in nonurban areas where most coal mines are located. Then, the concept of distributed scatterers (DSs) was introduced to improve the density of coherent measurements. For instance, small baselines subset (SBAS) InSAR methods [24], [25] extract DS pixels by limiting spatiotemporal baselines. SqueeSAR [26], in contrast to SBAS methods, provides another way to derive DS pixels by using a typical two-step strategy -statistically homogeneous pixel (SHP) identification and the phase optimization algorithm. Many enhanced SHP identification algorithms and phase optimization algorithms [27], [28], [29], [30], [31], [32], [33] were subsequently proposed to improve the performance of DS InSAR time-series analysis. In engineering applications, the probability integral method (PIM) is usually used for mining subsidence prediction. The PIM prediction results are the final cumulative deformation, including the deformation during the mining period and residual deformation, and can supplement the lack of InSAR monitoring results in the center of the subsidence basin. However, most existing studies use the InSAR monitoring results covering the mining period as the final cumulative deformation; however, by ignoring residual deformation, this approach introduces uncertainty into the results [34], [35], [36].
Therefore, the main objectives of this study are 1) to comprehensively characterize the spatial and temporal evolution of surface deformation in the Yineng Coal Mine, 2) to explore the applicability of InSAR for mining subsidence monitoring in backfilling coal mines, and 3) to further discuss the effectiveness of backfilling mining on subsidence restraint. To fulfill these goals, we implemented a systematic methodology (see Section III) based on the specific conditions of the study area. We demonstrated that the capability of deriving sufficient coherent points for traditional SBAS InSAR was unsatisfactory. Thus, we adopted a DS InSAR time-series analysis strategy to perform a time-series analysis using 26 C-band Sentinel-1 SAR images acquired from January 2019 to July 2020. The strategy combines the SBAS method implemented in StaMPS [24], the hypothesis test of confidence interval (HTCI) SHP identification algorithm proposed by [30], a spatial adaptive filtering and the eigen-decomposition-based maximum-likelihood-estimator of the interferometric phase (EMI) proposed by Jiang and Guarnieri [31]. Then, the InSAR deformation time series were integrated with a mining subsidence model specifically for backfilling coal mines [37] by using the PIM [38] to invert the surface deformation in the areas where even the DS InSAR time-series analysis strategy failed to obtain valid measurements. The time inconsistency between the results from InSAR and the PIM was considered by using an arctangent time function. The whole methodology allowed us to obtain accurate and detailed spatiotemporal surface deformation results which are elaborated in Section IV. On this basis, the effectiveness of backfilling mining on subsidence restraint in the Yineng Coal Mine and the limitations and prospects of research methods are discussed in Section V. Finally, Section VI concludes this article.

A. Study Area
The Yineng Coal Mine is located in northern Jining City, Shandong Province, in eastern China (see Fig. 1). The mine is adjacent to the Yiqiao Coal Mine in the northwest, the Tangyang Coal Mine in the southwest, the Luxi Coal Mine in the south, and the Xinyi Coal Mine in the southeast. The geographical ranges of the mine field are 35°36 30 N-35°40 00 N and 116°32 45 E-116°37 30 E, covering a total area of approximately 30 km 2 [see Fig. 1(a)]. The area is situated in the alluvial plain of the lower reaches of the Yellow River. The terrain is flat and tends to be slightly higher in the northeast and relatively lower in the southwest, with an elevation of 41-48 above sea level [see Fig. 1 Most of the ground surface is covered by farmland and village construction [see Fig. 1(c)]. The mine field is a fully concealed Carboniferous and Permian coal-bearing stratum of North China type, within which the Shanxi Formation and Taiyuan Formation are the main coal-bearing strata. According to the location of four major faults, namely, the Sunshidian fault (120-170 m), the Yuncheng fault (10-180 m), DF4, and DF25 (100-270 m), the Yineng Coal Mine is spatially divided into four mines. The initial minery, which is located in the north, is the main study area [see Fig. 1(a), (d)]. There are many minor faults, several villages (Wang, Xisunwu, Dongsunwu, Hu-Gong, Zhao-Lv, etc.,) and six mechanized working faces (CG1301, CG1302, CG1303, CG1304, CG1306, and CG1312) that were mined between 2014 and 2020 in the initial minery [see Fig. 1(d)]. The detailed parameters of the working faces, including the mining period, length and width of the working face, average mining depth, mining thickness, and inclination angle of the coal seam are presented in Table I. Unlike the neighboring coal mines, the Yineng Coal Mine adopts "Mining-filling-retention" coordinated backfilling mining technology based on ultrahigh water  I  DETAILED PARAMETERS OF THE WORKING FACES IN THE YINENG COAL MINE filling materials. The filling material not only absorbs the energy accumulated in the overburden but also restrains the deformation and slows down the release of energy accumulated in the overburden.

B. Datasets
Taking into account the study objectives, the actual mining situation of the Yineng Coal Mine, and the availability of SAR and leveling data, 26 ascending C-band Sentinel-1A SLC SAR images spanning from January 2019 to July 2020 provided by the European Space Agency were used to perform the InSAR time-series analysis to obtain the spatiotemporally detailed surface deformation in the study area. These images were captured in vertical-vertical (VV) polarization, interferometric wide (IW) swath mode, and TOPS acquisition mode. The track number was 142, and the resolutions were approximately 2.3 m and 13.9 m in the slant range and azimuth, respectively. The central incidence angle over the study area was 39.2°. The heading angle of the radar is −10.178°. Detailed information regarding the processed SAR images is given in Table II. The corresponding precise orbit data from the European Space Agency were collected to correct the orbital errors. The corresponding shuttle radar topography mission (SRTM) digital elevation model (DEM) with a resolution of 30 m was downloaded from the United States Geological Survey to georeference the interferograms and to remove the topographic phases. In addition, 88 leveling stations [indicated by the gray triangles in Fig. 1(d)] with a time resolution of two months are evenly arranged at intervals of 20 m on the main road of the villages in the initial minery of the Yineng Coal Mine. In total, 13 epochs of leveling data during the study period were used for accuracy validation and model constraints. The boundary map of the coal mine, the mining engineering map, and other mining area data were also obtained.

III. METHODOLOGY
Although the SBAS InSAR method was designed to improve the correlations of the interferograms by possibly minimizing the spatial and temporal baselines, according to experimental results [see Fig. 2(a)], the density of the coherent points selected by the SBAS method implemented in StaMPS [24] could not meet the needs of this study. There was always a tradeoff between the amount and the reliability of the detected coherent points. For instance, only a few points (127 points over working faces, Fig. 2(a), area B) were selected without compromising quality partly because of the relatively week ability to detect DS points Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply. and partly because of the specific land cover (cropland) over the main working faces. Therefore, we adopted a two-step systematic methodology (see Fig. 3) to maximally increase the spatial density of valid measurements and derive comprehensive information of deformation over the study area. First, we performed a DS InSAR time-series analysis strategy to achieve a first-tier densification of coherent points and generate an accurate deformation time series (see Section III-A). Second, the derived InSAR deformation results were integrated with a mining subsidence model specifically for backfilling mines by using the PIM to retrieve the displacement information in the decorrelated areas and thus to reach the second-tier densification of measurements (see Section III-B). The integration of DS InSAR and PIM allowed us to get access to comprehensive information of deformation in the study area.

A. DS InSAR Processing
The DS InSAR time-series analysis strategy combines the SBAS method implemented in StaMPS [24], the HTCI SHP identification algorithm [30], spatial adaptive filtering, and EMI [31]. Thus, the time-series processing consists of five main steps: generating a subset of small baseline (SB) interferograms (see Section III-A1), identifying SHPs (see Section III-A2), estimating the optimal phase (see Section III-A3), extracting coherent pixels (see Section III-A4), and retrieving InSAR deformation time series (see Section III-A5).
1) Generating a Subset of SB Interferograms: The processing in this section was conducted using GAMMA software provided by GAMMA Remote Sensing AG. The image acquired on 18 November 2019 was selected as the supermaster image by analyzing the overall correlation coefficient. All the secondary images were strictly coregistered and resampled to the supermaster image. Then, the SLC images were deramped for the azimuth spectrum ramp before interferogram formation. No multilooking was applied to retain the highest spatial resolution. The maximum perpendicular and temporal absolute baselines were set to 200 m and 40 days, respectively, to construct the subset of SB interferograms, under the premise of connecting all images in the spatiotemporal baseline network (see Fig. 4). The 30 m SRTM DEM was used to simulate and subtract the topographic contribution. A subset of 67 SB interferograms was eventually obtained.
2) Identifying SHPs: The SHP identification algorithm determines the similarity between neighboring pixels and center pixels by statistical inference. The algorithm is the basis of the DS InSAR technique, and the algorithm's estimation accuracy directly affects the precision of subsequent parameter estimation [27], [39]. In this article, we used the HTCI algorithm [29] that combines the advantages of both the likelihood ratio test (LRT) [40] and the fast statistically homogeneous pixel selection (FaSHPS) [27]. The algorithm takes the unbiased estimation of the LRT as the truth value of the reference pixel. At a given confidence level, if a temporal averaged intensity value of the neighboring pixel falls within the interval of truth value construction, it is considered a homogeneous sample of reference pixel, which simultaneously guarantees the speed and accuracy of the SHP identification. In the progress, the search window and the significance level are the main parameters. After many adjustments, we set them as 15 × 15 and 0.05, respectively. Finally, the matrix Ω init containing the count of the homogeneous pixels of every pixel was obtained.
3) Estimating the Optimal Phase: A DS resolution cell consists of multiple independent scatterers leading to temporal geometrical decorrelation. Therefore, the estimation of systematic phase series, which is referred to as phase optimization or phase linking according to pioneering studies, after SHP identification is essential to improve the signal-to-noise ratio of the interferometric phase and increase the density of valid measurements. In this article, we adopted the EMI that considers simultaneously both the estimation and computational efficiency as the core algorithm in our phase optimization strategy.
According to [26], most PSs have fewer than 20 SHPs. Therefore, to preserve PS information, the pixels with more than 20 SHPs were selected to form a set of DS candidates Ω DSC for optimal phase estimation. The pixels with fewer than 20 SHPs were regarded as PS candidates Ω PSC and classified as pending.
Phase noise introduced by scene-dependent or system-related decorrelation can present challenges for phase estimation. Thus, unlike traditional EMI-based phase optimization methods, we applied spatial adaptive filtering [41] to the SB interferometric phase of DS candidates Ω DSC before phase estimation to improve fringe visibility and avoid potentially blurring or destroying phase information.
The complex coherence matrix (CCM) that describes the statistical characteristics of DSs is the basis of all phase optimization algorithms. The CCMĈ could be regarded as the Hadamard product of two matrices Φ andΓ reflecting the interferometric phase and coherence, respectively. To improve the estimation efficiency of EMI, we utilized the filtered interferometric phase and its corresponding coherence. On account of the biased estimation of coherence, the second kind statistics based on the Mellin transform [42] was applied to correct the above estimated coherence. Next, the optimized phase can be estimated by maximizing the absolute value of the logarithm of the joint-probability distribution function of the DS candidates. In the EMI, the solution is the eigenvector corresponding to the minimum eigenvalue of the Hadamard productΓ −1 •Ĉ, which is the final optimized phase series corresponding to SAR acquisitions for the reference DS candidate pixel.

4) Extracting Coherent Pixels:
Determining whether a DS candidate is reserved as a DS pixel relied on the quality of the candidate's estimated optimal phase. Hence, the goodness of fit between the original interferometric phase and the optimized phase was used as the evaluation index for optimized phase quality. The goodness of fit was calculated as where γ TC represents the goodness-of-fit index, known as the temporal coherence [26], [43], which was further used for DS pixel selection; ϕ m,n is the original interferometric phase between the mth and nth acquisitions; and θ m and θ n are the reconstructed optimal phases. Only those DS candidates with γ TC values greater than 0.5 were regarded as high-quality DS pixels. Then, the PS candidate pixels (with less than 20 SHPs, obtained in Section III-A3) were further screened by a threshold of amplitude dispersion (0.6). Finally, the high-quality DS pixels and the screened PS pixels were combined to form a coherent pixel set (Ω CP ) (the first-tier densification of measurements) for retrieving of InSAR deformation time series in the following section. As shown in Fig. 2(b), the number of measurement points has significantly increased (approximately 200%) compared with the number of measurement points shown in Fig. 2(a).

5) Retrieving InSAR Deformation Time Series:
The whole set of coherent pixels extracted from Section III-A4 was input into the StaMPS package [24] for further analysis by using the SBAS strategy. The atmospheric phase delays were corrected using the Generic Atmospheric Correction Online Service for In-SAR (GACOS) products developed based on numerical weather models [44]. Ultimately, the InSAR deformation time series corresponding to each SAR acquisition and the deformation rate for each coherent pixel were inverted from the space of corrected SB interferograms.

B. Integrating DS InSAR and the PIM for Backfilling Mining
Compared with the SBAS method, DS InSAR notably improved the density of coherence points (see Fig. 2, areas A and C) and the coherence [see Fig. 2(c) and (d)] in town and village areas. However, the improvement in the cropland area where the original mean coherences are approximately 0.1 is not evident. The optimized mean coherence is still mostly less than 0.3 [see Fig. 2(e)], resulting in the lack of coherence points above the CG1312 working face (see Fig. 2(b), area B). Therefore, we performed second-tier processing by integrating the derived InSAR deformation results with a mining subsidence prediction model for backfilling mining by using the PIM to obtain more detailed deformation information on the study area. In this section, the basic theory of the traditional PIM, the PIM for backfilling mining, estimating model parameters by using the DS InSAR results, and retrieving the final deformation time series are elaborated.
1) PIM for Backfilling Mining: Although many subsidence prediction methods have been developed, the PIM is the most widely used and recommended mining subsidence prediction method in China's "Three Under" (under buildings, railways, and water bodies) mining regulations [45]. The PIM connects the surface deformation with coal seam information by using the stochastic medium theory model (SMTM), which was introduced for the first time in research on strata movement in the 1960 s [38]. The SMTM assumes the rock strata above the working face to be a homogeneous discontinuous medium, and the movement of rock strata can be regarded as a linear superposition of the movement caused by infinite tiny mining units. This superposition can be accomplished by the integral of the probability distribution density curve.
As shown in Fig. 5(a) and (b), under full mining conditions, the subsidence at ground surface point A(x, y) is caused by exploiting any unit B(u, v) in an inclined coal seam; this subsidence is expressed as y) is the ground subsidence at point A(x, y); exploiting unit B(u, v) is the cause of W u (x, y); r is the major influence radius of mining unit B(u, v); r = H/ tan β ; H and tan β represent the mining depth of mining unit B(u, v) and the tangent of the major influential angle, respectively; and θ is the propagation angle of mining influence.
According to the superposition principle, at ground point A(x, y), the total subsidence induced by the whole underground working zone can be written as where the maximum subsidence of full extraction is W 0 = Mq cos α; M , q, and α represent the normal mining thickness, the subsidence coefficient, and the coal seam inclination, respectively; l = D 1 − s 3 ; L = (D 2 − s 2 ) cos α; D 1 and D 2 are the strike and dip lengths of the working face, respectively; and s 1 , s 2 , s 3 , and s 4 are the left, upper, right, and lower deviation of the reflection points, respectively. The aforementioned is traditionally adopted to predict the deformation induced by the caving mining method. Although the characteristics of the subsidence basin are similar between backfilling mining and caving mining methods, the strata movement mechanism of the former differs from that of the latter; consequently, the PIM parameters have different values and meanings [37]. For instance, backfilling mining reduces the actual mining thickness of a coal seam by backfilling specific materials; relying on this technique, the concept of the equivalent mining height was proposed, expressed as [37], [46], where σ is the compression ratio and is set to 90% according to the existing mine data. Then, considering M e and (2) and (3), the total subsidence at ground point A(x, y) in the PIM model for backfilling mining can be expressed as Correspondingly, the horizontal displacement (U ϕ (x, y)) of point A(x, y) along a certain direction ϕ can be determined based on its relationship with the subsidence as follows: where I ϕ (x, y) represents the tilt along direction ϕ, which is the first spatial derivative of the subsidence W (x, y) [47], and b is the coefficient of the horizontal displacement.

2) Model Constraints by Using DS InSAR Results:
To use DS InSAR results to constrain the PIM model for backfilling mining, two inconsistencies between these results and the PIM model must be resolved.
One inconsistency that must be resolved is that the displacements derived from InSAR are along the (LOS direction, while the displacements predicted by the PIM are vertical and horizontal displacements. Here, we adopted the following equation [48] to construct the relationship between the two types of deformations: , which is obtained from (4), is the predicted subsidence of pixel (i, j); U NS (i, j) and U EW (i, j), which are obtained from (5), are the predicted displacements of pixel (i, j) in the north and east directions, respectively; θ inc denotes the radar incidence angle; and α azi represents the radar heading angle. The other inconsistency that must be resolved is that the predicted displacements are the final cumulative deformation, while the InSAR-derived deformation time series are within one stage in the entire deformation. To make them unified, we fitted the DS InSAR deformation time series with an arctangent time function [see (7)] [49]. Compared with the classic Knothe time function [50], the arctangent time function is more suitable for this study according to the temporal behavior of the deformation of monitoring points (see Fig. 6). Once the function reached the best fit using the least-squares estimation, the optimal parameters and the S-shaped curve for each monitoring point were obtained. Then, the final cumulative deformation in the LOS direction where δ is the accumulative displacement value (cm), t is the monitoring time (days), p 0 represents the speed at which the deformation curve approaches the asymptote, p 1 represents the size of the asymptote numerical, and p 2 corresponds to the inflection point of time (day) on the deformation curve; this inflection point represents the time of maximum settling velocity at the monitoring point. p 1 , p 2 , and p 3 vary with the change in the monitoring data. The root-mean-square error (RMSE) between D LOS and D LOS in the overlaying area was used to evaluate the fitting degree and, thus, constrain the PIM model.

3) Estimation of Model Parameters:
In (4) and (5), parameters M , H, D 1 , D 2 , α, and s 1 -s 4 are constant over the continuous working face CG1312. To estimate the remaining unknown parameters θ, tan β, q, and b, an iterative approach was used as follows. 1) Step 1: Initial values and value ranges of the parameters are set by comprehensively considering the mining scheme, deformation characteristics, historical ground observation data, and previous studies on backfilling mining [37], [51].

2)
Step 2: The cumulative subsidence (W ) and horizontal deformation in the north and east directions (U NS and U EW ) are predicted according to Section III-B1. Then, the D LOS predicted deformation in the LOS direction is calculated using (6).

3)
Step 3: The final cumulative deformation in the LOS direction D LOS is simulated based on (7) and the RMSE, say ε, between D LOS and D LOS in the overlay area is generated.

4)
Step 4: E is compared with a given threshold ε T . If ε > ε T , the parameters are adjusted, and Steps 2 and 3 are repeated; if ε ≤ ε T , the parameters are outputted. With the aid of the estimated parameters, the optimal final cumulative subsidence and horizontal displacements and then the final cumulative displacement and the displacement time series corresponding to SAR acquisitions in LOS directions over working face CG1312 were generated according to (4)- (7). Finally, the DS InSAR-derived deformation time series and the PIM inverted deformation time series were fused according to the following rules: where D PIM represents the deformation inverted by the PIM for backfilling mining with DS InSAR constraints; D InSAR represents the deformation derived from DS InSAR processing; x is a pixel; and Ω CP is the set of coherent points extracted in Section III-A4.

IV. RESULTS
Based on the above strategy, the annual surface deformation rate covering the Yineng Coal Mine and its surrounding coal mines and the complete deformation time series of the CG1312 working face were obtained. In this section, on the basis of the accuracy validation and uncertainty analysis of the results (see Section IV-A), we comprehensively analyze the spatiotemporal evolution characteristics of the surface deformation of every coal mine (see Section IV-B), the initial minery of the Yineng Coal Mine (see Section IV-C), and the CG1312 working face in the initial minery of the Yineng Coal Mine (see Section IV-D).
A. Uncertainty Analysis of the Results 1) Accuracy Analysis of the DS InSAR Results: An internal coincidence accuracy analysis was carried out by comparing and analyzing the consistency of the deformation rate results obtained by the DS InSAR and SBAS methods. With 5 m as the buffer radius, 32 552 homonymy points of the two results were selected; these points were distributed mainly in cities and villages and a few in cropland. As shown in Fig. 7(a), the Pearson correlation coefficient is 0.77. Most of the SBAS monitoring values are slightly larger than those of DS InSAR. As the histogram of the deformation rate difference [see Fig. 7(b)] shows, the deformation rate difference between the DS InSAR and SBAS results statistically meets the normal distribution with a mean value of 5.54 mm/a and a standard deviation of 4.02 mm/a. There are 29 037 homonymy points with absolute difference values lower than 10 mm/a, accounting for 89.2% of the total, indicating that the DS InSAR and SBAS results have good consistency.
To verify the external coincidence accuracy of the DS InSAR results, they were compared with data from 20 ground leveling stations; the data were acquired at seven adjacent epochs (corresponding to the blue image points in Fig. 4). The leveling stations are distributed mainly along the northeast border of Xisunwu, Hu-Gong, and Zhao-Lv Villages [locations are indicated by the cyan triangles in Fig. 1(d)]. Spatiotemporal and geometric consistency are prerequisites for the comparative analysis of different data sources. Therefore, a nondeformation area with a known elevation south of Wenshang County [indicated by the red cross in Fig. 1(d)] was selected to replace the origin of Chinese leveling (Qingdao geodetic origin) as the reference point of surface deformation to ensure the spatial consistency of the data. The date 2 January 2019 was taken as the starting time of the leveling data to approximately ensure the time consistency of the data. Assuming that the north-south and east-west horizontal deformation of the mining area was the product of the maximum subsidence value and the horizontal movement coefficient (which, referring to a priori knowledge, is 0.3), the vertical leveling data were converted into LOS-directed deformation according to (6), ensuring the geometric consistency of the two sets of data. Fig. 7(c) represents the correlation between the DS InSAR results and the corresponding leveling data. The Pearson correlation coefficient reaches 0.9 and the RMSE is 10.1 mm, indicating that the two datasets are consistent. Fig. 7(d) represents the correlation between the SBAS results and the corresponding leveling data. The Pearson correlation coefficient reaches 0.8 and the RMSE is 14.1 mm. Compared with the SBAS results, it shows that the DS InSAR results are more accurate. In addition, these results indirectly suggest that the deformation model and horizontal movement coefficient used in this study are suitable in this mining area. In summary, the DS InSAR results are accurate and reliable.
2) Accuracy Analysis of the Integrated Results: The accuracy of the complete surface deformation time series of the CG1312 working face obtained by integrating DS InSAR and the PIM was evaluated in two ways in this section. First, the conformity between the integrated deformation time series and the leveling data in the working face [gray triangles in Fig. 1(d)] was validated. Second, the residuals between the integrated deformation time series and DS InSAR deformation time series were statistically analyzed in the overlapping area (Hu-Gong village). Fig. 8(a1)-(e1) shows the correlations between the five integrated time-series results and the leveling data recorded at adjacent times. Fig. 8(f1) illustrates the spatiotemporal mean RMSEs of the two above datasets at the leveling station positions. Overall, the data values are basically within double standard deviation confidence intervals. Over time, as the deformation increased, the RMSE increased slightly, ranging from 3.93 to 14.67 mm. The spatiotemporal mean RMSE fluctuated in the range of 0.5-24.8 mm. Since the cumulative deformation at the first time point is too small, the Pearson correlation coefficient at this time point is relatively low, while the correlation coefficients at other time points are all high. Therefore, the above two indicators comprehensively suggest that the integrated results of DS InSAR and the PIM accorded well with the leveling data. Fig. 8(a2)-(e2) and (f2) shows the spatial pattern and statistical boxline distribution of residuals between the integrated deformation time series of the CG1312 working face and the corresponding DS InSAR deformation time series in the overlapping area (Hu-Gong village), respectively. The residual values increased slightly over time and were mostly concentrated between −10 and 10 mm. The statistics shows a normal distribution with a mean value of approximately 0 mm. Therefore, the accuracy of the integrated results under the constraints of the DS InSAR results was demonstrated.
Combining both of the afore-mentioned two approaches to accuracy analysis demonstrates that the integrated deformation time series of the CG1312 working face can comprehensively and accurately reveal the spatiotemporal behavior of the surface deformation above the working face; this finding is helpful in further identifying the deformation characteristics of the working face and discussing the effect of backfilling mining. Fig. 9 shows the DS InSAR surface deformation rate for the entire monitoring area. The decorrelation was evident in the high vegetation coverage area. Although the coherent points were maximally extracted by the advanced DS InSAR method, there were still obvious "holes" with missing deformation information. On the whole, there are several deformation zones caused by coal mining activities in the Yineng Coal Mine and its surrounding mining areas. Among these zones, the deformation zone of the Yineng Coal Mine is located in the initial minery, and the maximum deformation rate is approximately −83.60 mm/a. A more detailed analysis of the spatiotemporal evolution characteristics of surface deformation in the Yineng Coal Mine is presented in Section IV-C. There is a deformation zone distributed mainly near Cao village in the south of the Yiqiao Coal Mine, which is located in the northwest of the Yineng Coal Mine [DZ1 in Fig. 9(a) and (b)]. The maximum deformation rate is approximately −48.40 mm/a, and the observed area affected by deformation is about approximately 0.45 km 2 . The Tangyang Coal Mine is located in the southwest of the Yineng Coal Mine. There are two main deformation zones. One is located in the north and connected to the deformation zone within the Yiqiao Coal Mine [DZ2 in Fig. 9(a) and (c)]. The maximum deformation rate is approximately −23.00 mm/a and the observed area affected by deformation is approximately 1.84 km 2 ; the other is located in Xitangyang village in the south of the Tangyang Coal Mine [DZ3 in Fig. 9(a)

C. Spatiotemporal Characteristic of Surface Deformation in the Yineng Coal Mine
This section focuses on analyzing the spatiotemporal evolution characteristics of surface deformation in the initial minery of the Yineng Coal Mine during the study period. As shown in Fig. 10(a), the obtained monitoring points are distributed mainly in the village areas and partially cover the working faces CG1301, CG1302, CG1303, CG1304, and CG1306. Xisunwu, Dongsunwu, Wang, and Hu-Gong villages are all affected by surface deformation induced by mining activities. Except for the CG1302, CG1304, and CG1312 working faces, all other working faces were mined during the study period (10 January 2019-3 July 2020); hence, the corresponding observed deformation is the residual deformation after mining. To analyze the spatial pattern of the accumulated surface deformation above each working face in detail, a profile along the strike direction of each working face was established, and the accumulated deformation corresponding to the deformation monitoring points within a 10 m buffer area was extracted, as shown in Fig. 10(b).
A part of the CG1301 and a part of the CG1302 working faces cross Dongsunwu village. The surface deformation rate gradually increases along the northeast-southwest direction,   consistent with the mining direction, and the maximum value of the surface deformation is −76.40 mm/a, occurring at the southwest corner of Dongsunwu village above the CG1302 working face [see Fig. 10(a)]. The last part of the CG1302 working face was being mined during the study period, but InSAR failed to obtain a coherent point in that area due to vegetation coverage. The deformation observed above the working faces, corresponding to areas that have been mined during the study period, is thus residual deformation. By analyzing the profiles of the two working faces [see Fig. 10(b)], the cumulative deformation of the two working faces displayed a similar trend of variation. The farther away from points A/B, the smaller the cumulative deformation. The mean deformation rate of monitoring point P1 located in the mining area of the CG1302 working face in 2018 is approximately −39.50 mm/a. The mining of the CG1301 working face was completed in October 2015, but the results showed that there was still an average land subsidence of −48.60 mm within 4-5 years after the mining. It is partly attributed to the on-going compaction of overlying strata and partly due to the influence of the mining of the CG1302 working face.
The CG1304 working face and part of the CG1303 working face cross Wang village. The surface deformation rate gradually increases along the northeast-southwest direction, which is consistent with the mining direction. The maximum value is −69.30 mm/a, and the deformation is located at the southwest corner of Wang village above the CG1303 working face [see Fig. 10(a)]. As shown in Fig. 10(b), the cumulative deformation of the two working faces presents a similar trend. The farther away from points C/D, the smaller the cumulative deformation. Notably, the mining of the CG1303 working face was completed in April 2017, but the residual deformation rate during the study period was even greater than that of the CG1304 working face that was being mined possibly because the CG1303 working face has a longer strike length (see Table I). The deformation rate of monitoring point P2 above the CG1303 working face is approximately −51.20 mm/a. The deformation rate slightly slowed down and tended to stabilize around July 2020.
Part of the CG1306 working face crosses the eastern edge of Xisunwu village. The deformation rate of Xisunwu village decreases along the southeast-northwest direction due to the backfilling mining of the working face, and the maximum residual deformation rate point is located in the southeast corner of Xisunwu village [see Fig. 10(a)] with a value of −83.60 mm/a. The profile of the working face [see Fig. 10(b)] shows that the farther away from point E, the greater the accumulated deformation. The maximum residual deformation reaches approximately −120.00 mm. According to the basic data of the mining area, monitoring point P3 above the CG1306 working face [see Fig. 10(a)] is within the mining area in 2017. The deformation time series (see Fig. 10(a), P3) shows that the mean deformation rate during the entire study period is approximately −54.00 mm/a. Specifically, the surface deformation showed a piecewise linear trend. The surface deformation rate slowed down around April 2019 and then did not slow down until July 2020; this finding indicates that large residual deformation may still exist within 2-3 years after mining. Therefore, the significance of regularly monitoring the residual deformation over the backfilling working faces cannot be ignored.
The CG1312 working face was in the mining state during the study period, and the cropland was covered directly above it, resulting in almost no reliable InSAR deformation information. However, the results show that Hu-Gong village, approximately 200 m northwest of the working face, was obviously affected by mining activity [see Fig. 10(a)]. The deformation rate of Hu-Gong village is correlated with its distance from the working face and decreases along the southeast-northwest direction, with a maximum deformation rate of approximately −61.00 mm/a.
The extracted time series of the monitoring points in the northeast and southwest of Hu-Gong village (see Fig. 10(a), P4 and P5) show that the two areas have different deformation behaviors. Specifically, P4 began to subside before P5 and was basically stable in July 2020, while P5 still declined; these results accord with the time course of coal mining advancement. The detailed analysis of the integrated inversion results of the surface deformation above the CG1312 working face is elaborated in Section IV-D.

D. Spatiotemporal Characteristic of the Complete Deformation Field in the CG1312 Working Face
By using the method described in Section III-B, the complete deformation fields along the LOS direction during the study period were obtained.
Five of the deformation time-series results and their corresponding optimal inversion parameters are shown in Fig. 11 and Table III, respectively. The scale and magnitude of the surface deformation constantly expanded with continuous mining and filling in the CG1312 working face. The maximum deformation along the LOS is −138 mm, and the affected area expanded nearly 1 km away from the working face by 2.5 km 2 .
To further analyze the temporal evolution of the surface deformation above the CG1312 working face, a profile FF' along the strike direction passing through the center of the working face and a profile GG' along the dip direction passing through Hu-Gong and Zhao-Lv villages were established [blue dotted line in Fig. 11(e)]. Fig. 12 is an overlay of the time-series profiles of FF' and GG'. The profiles of FF' show that each curve is funnel-shaped and generally symmetrical. The deformation decreases from the center of the funnel to both sides [see Fig. 12(a)]. The maximum deformation and influence radius of each curve gradually increased with time, and the center of the subsidence funnel gradually moved along the mining direction while lagging behind the mining stop line; these results accord with the mining subsidence law. Specifically, the maximum deformation reached −13 mm in the fourth month of mining (14 August 2019), and then with the continuous mining of the working face, the growth rate first increased and then decreased taking 5 January 2020 as the inflection point [see Fig. 12(c)]. The maximum deformation value reached −135 mm in the end. In particular, the deformation near leveling station N15 [as indicated by the gray triangle in Fig. 10 and the gray dotted line in Fig. 12(a)] on the profile was small, reaching a temporary steady state during the period from April 2020 to July 2020. According to the line graph of the maximum deformation and strike mining length [see Fig. 12(c)], the maximum deformation is positively correlated with the mining length of the working face, and the rate of the maximum deformation is also highly consistent with the mining rate.
In the profiles of GG' [see Fig. 12(b)], each curve is funnelshaped, and the deformation value decreases from the center of the funnel to both sides and is symmetrically distributed around the center 775 m away from point G [as indicated by the gray dotted line in Fig. 12(b)]. The deformation center deviates from the working face centerline (the strike section line) by  Fig. 12(d)], as the mining gradually approached profile GG', whose excavation time was approximately October 2019 according to the mining and excavation time data, the deformation and deformation rate continued to increase. As the mining gradually moved away from profile GG', although the deformation at the profile still increased, the rate of increase gradually decelerated.

A. Effectiveness of Backfilling Mining for Subsidence Restraint in the Yineng Coal Mine
The purpose of backfilling mining in coal mines is to reduce pressure and subsidence by using backfill materials and to prevent personnel and property losses caused by overlying rock subsidence and violent ground movement. In current research and engineering practice, the effectiveness of backfilling technology for subsidence reduction is generally evaluated by comparing the maximum subsidence value under backfilling and total caving mining conditions. For instance, Feng et al. [52] applied the ultrahigh water material filling technology to the Taoyi Coal Mine and used the equivalent mining height theory and the PIM method to calculate the ratio of the difference between the filling thickness and the subsidence and the mining thickness under the condition of the backfilling method. The calculated filling ratio is 86.5%, which accords with China's national first-class standard. Other scholars comprehensively evaluate the filling effectiveness according to a variety of factors. For instance, Wu et al. [53] established a backfilling effectiveness analysis model (including the main control factors of the cooperative bearing structure of stope surrounding rock, backfilling properties, mining geological conditions, and mining and backfilling technology) based on the hierarchical extension fuzzy evaluation method. The analytical model was successfully applied to the Xinjulong Coal Mine in Shandong Province, China.
In this section, by using the predicted final vertical cumulative displacement above the CG1312 working face in Section IV-D [see Fig. 13(a)] and the corresponding optimal parameters (see Table III), the predicted vertical cumulative displacement of total caving mining under the same conditions was obtained by using the PIM model. Then, the difference in surface subsidence caused by backfilling and total caving mining was calculated [see Fig. 13(b) and (c)]. The effectiveness of backfilling at the CG1312 working face of the Yineng Coal Mine was analyzed from the perspective of the spatial pattern and magnitude variation. Fig. 13 shows that the area of the deformation zone under the backfilling mining condition is approximately 1.25 km 2 , and the maximum cumulative deformation is −213 mm, while the area of the deformation zone under the total caving mining condition is approximately 1.60 km 2 , and the maximum cumulative deformation is −553 mm. By comparison, under backfilling mining conditions, the deformation influence range is reduced by 22%, and the cumulative maximum deformation is reduced by 61.50%. The mean deformation is reduced by −96.50 mm, and the mean deformation reduction ratio is 17.50%. Under total caving mining conditions, the area with a deformation greater than −200 mm accounted for 26.30% of the total displacement area, and the area with displacement greater than −500 mm accounted for 1.40%. In particular, the mean displacement was reduced from −75.30 to −30.10 mm, and the mean inclination was reduced from 0.57 to 0.22 mm/m with the reduction ratio of 61.40%. The maximum tilt is reduced from 1.50 to 0.60 mm/m due to the backfilling mining just above Hu-Gong village. Generally, in the stability evaluation of buildings in the mining area, when the inclination reaches 3.00 mm/m, the inclination constitutes the first-level damage of brick-concrete buildings [54]. The results show that the predicted maximum deformation without backfilling has not yet reached the damage threshold, but further deformation may still be caused by later mining and adjacent working face mining. Backfilling mining reduces the maximum inclination to less than 1 mm/m, thereby significantly reducing potential threats to the stability of village buildings.

B. Interpretation of the Deformation Inconsistency Revealed by DS InSAR and the PIM
Due to the influence of temporal decorrelation, the DS InSAR technique used in this study can obtain only limited monitoring points, but the introduction of a mining subsidence model makes it possible to obtain complete surface deformation information of the working face, which is of great significance to understand the law of surface movement caused by coal mining and the early warning of potential geological disasters. However, the PIM model assumes that the rock formation above the working face is a homogenous medium and does not consider the complex geological environment that may exist in reality, e.g., the distribution of faults. Therefore, the spatial distribution of the obtained deformation results is relatively uniform. The InSAR observation reflects the surface deformation characteristics under real and objective conditions, and the spatial distribution of the deformation is relatively complex; these results also explain the residual errors between the integrated results and InSAR/leveling data (see Fig. 8). Therefore, it is often challenging to use InSAR results as the boundary constraints on the PIM.
There are four major faults and several minor faults in the Yineng Coal Mine. The CG1312 working face is located between the two major faults: the Sunshidian fault and the Yuncheng fault. The two villages adjacent to the working face, namely Hu-Gong and Zhao-Lv, are located on both sides of the Sunshidian fault with a drop of 120-170 m. Hu-Gong village, located on the west side of the fault, is on the same side as the working face, while Zhao-Lv village is located on the east side of the Sunshidian fault. Visual interpretation of the DS InSAR results [see Fig. 10(a)] shows that although the distances from the working face to Hu-Gong and Zhao-Lv villages are not much different, there is obvious deformation in Hu-Gong (ranging from −20 to −60 mm/a), but almost no deformation occurred in Zhao-Lv village. The cumulative deformations of leveling stations D34 and D35 near Zhao-Lv village were −9 and 0 mm, respectively, which were consistent with the DS InSAR results. After excluding the cause of deformation accuracy, the author considers that this phenomenon may be due to the influence of the Sunshidian fault. Therefore, the deformation of Zhao-Lv village was not considered when using the DS InSAR results to constrain the PIM model for inversion, and only the deformation results of Hu-Gong village were used. Although the reliability of the integrated results has been validated (see Section IV-A2), there was a difference between the integrated results and the In-SAR/leveling measurements in Zhao-Lv village. The integrated deformations of D34 and D35 were −30 and −22 mm, respectively; these deformations are nearly 20 mm from the real deformation. In addition, the InSAR monitoring points in Hu-Gong village fit well with the time-series profiles [see Fig. 12(b)], but there is a large difference between the in Zhao-Lv village [as indicated by the red dotted circle in Fig. 12(b)]. Specifically, almost no deformation is observed by InSAR in Zhao-Lv village, but the integrated results show that the deformation values in this area increase with time. The above discussion indicates that the Sunshidian fault plays a certain role in controlling the surface deformation characteristics of Zhao-Lv village. The fault blocks the propagation of the stress changes from the overlying strata above the working face. However, the influence of complex geological conditions, such as faults, was not considered in the PIM model; thus, the inverted results would be different from the actual situation. In addition, according to the spatiotemporal mean RMSEs between the leveling data and the integrated results [see Fig. 8(f1)], the spatial heterogeneity of their magnitudes has a certain correlation with the distribution of major and minor faults. The local RMSE value in the fault area is relatively large.
The discussion in this section also indirectly demonstrates that both the InSAR and theoretical models have their own advantages and limitations in investigating the law and characteristics of surface movement in mining areas. The combination of the two can improve the monitoring effectiveness but a series key issues still need to be resolved in future studies.

C. Integration Methods for Monitoring Surface Deformation in Backfilling Mining Areas
Compared with general mining areas, InSAR surface deformation monitoring in backfill mining areas is also affected by temporal decorrelation, which however is mainly caused by surface vegetation cover, rather than by excessive subsidence gradient. It is precisely because the surface deformation gradient of backfilling mining area is relatively small; the accuracy of monitoring results is required to be higher. Therefore, a single monitoring method is difficult to meet the requirements of surface deformation monitoring in backfilling mining areas. In this study, we adopted the method of integration of DS InSAR and the PIM model, which not only realized high-precision monitoring but also achieved the purpose of spatial densification of observations. Besides, with the iterative update of data sources, methods, and models, more integrated ideas are provided for surface deformation monitoring in backfilling mining areas. 1) Although Sentinel-1 data have good availability and continuity, the ability of C-band microwave signals to penetrate vegetation to obtain surface scattering information is limited, thus the integration potential of long-wavelength spaceborne SAR data in mining areas with high vegetation coverage can be explored; these data include the L-band data of China's LT-1 and Japan's ALOS-2 satellites. In addition, the deformation monitoring can also be performed in combination with high-resolution optical remote sensing images, such as multiple-pairwise image correlation toolbox for processing OPTical images developed by Provost et al. [55]. 2) This study found that local geological faults would affect the spatial distribution of surface deformation in mining areas, which was not taken into account by the traditional PIM model. Therefore, the PIM model can further be enhanced and constrained according to specific conditions to improve the applicability of the model. 3) POT is not affected by temporal decorrelation of inteferometric signals, and some scholars have integrated it with InSAR [18] or the PIM model [56] to monitor the surface deformation in mining areas. However, its monitoring accuracy is often restricted by the spatial resolution of SAR images and it is more suitable for extracting large gradient deformation information. Therefore, its applicability in backfilling mining area needs to be investigated and discussed.

VI. CONCLUSION
In this study, DS InSAR was integrated with the PIM model and an arctangent time function to jointly obtain the surface deformation field in and around the Yineng Coal Mine. First, DS InSAR was used to perform a detailed time series analysis on 26 C-band Sentinel-1A SAR images from 10 January 2019 to 3 July 2020. Then, the PIM model specific for backfilling mining and an arctangent time function was integrated with DS InSAR results. The spatial and temporal evolution characteristics of surface deformation are analyzed in detail, and the subsidence reduction effect of backfill mining technology is investigated for the first time with the aid of a remote sensing method. The effectiveness and limitations of the proposed strategy are discussed, and the surface deformation monitoring by integration methods in backfilling mining areas is prospected. The following conclusions can be drawn from this study.
1) The proposed integration strategy of the DS InSAR and PIM model is helpful in obtaining the complete spatiotemporal information of the surface deformation of the working face. 2) Eight key deformation zones are induced by mining in the whole study area. Except for the Yineng initial minery, the observed area affected by deformation is approximately 14.26 km 2 , and the maximum deformation rate is −48.40 mm/a. All six working faces in the Yineng initial minery display different degrees of surface deformation during the study period. There was still an average residual deformation of −48.60 mm within 4-5 years after the mining of the CG1301 working face; this finding indicates that regular monitoring of the residual deformation of the backfilling mining area is essential.
3) The maximum deformation in the CG1312 working face along the LOS is −138 mm, and the affected area is approximately 2.5 km 2 . Backfilling mining technology reduced the surface subsidence range of the working face by 22.00%, and the accumulated maximum subsidence decreased by 61.50%; this decrease significantly reduced the potential threat to surface buildings. 4) The distribution of faults in the initial minery of the Yineng Coal Mine increases the complexity of the spatial pattern of surface deformation. The local spatial heterogeneity can be captured by InSAR and leveling but was not considered by the mining subsidence model, which indicates that the integration of the two means is of great significance. This research can provide useful insights for overcoming the limitation of InSAR applications on mining subsidence monitoring and provide an essential reference for disaster prevention and overall planning in backfilling mining areas. He is currently an Associate Professor with the School of Environment Science and Spatial Informatics, China University of Mining and Technology. His research interests include the fields of interferometric phase filtering, phase unwrapping, and synthetic aperture radar interferometry signal processing and applications. He is currently a Postdoctoral Research Fellow with the School of Environment Science and Spatial Informatics, China University of Mining and Technology. His current research interests include the fields of interferometric synthetic aperture radar phase filtering, phase unwrapping, multitemporal In-SAR data processing, and deformation measurement.

Shijin
Si Chen received the B.E. degree in geographical information science from the China University of Mining and Technology, Xuzhou, China, in 2020. She is a graduate student majoring in photogrammetry and remote sensing with the China University of Mining and Technology.
Her research interest is synthetic aperture radar interferometry.
Guangli Guo received the Ph.D. degree in geodesy and surveying engineering from the China University of Mining and Technology (CUMT), Xuzhou, China, in 1999.
He was a Senior Visiting Scholar at the University of Nottingham, Nottingham, U.K., in 2017. Since 2001, he has been a Professor with CUMT. His research interests include mining subsidence control, stability and construction suitability evaluation in goaf, deformation monitoring, and data processing. He was an Associate Professor at the School of Mines, CUMT, in 2015. He is currently a Professor with CUMT. His research interest is safe and efficient longwall mining and overlying strata movement control technology, with particular emphasis on its application for underground coal mine.
Dongsheng Zhao received the bachelor's degree in geodesy and geomatics engineering from Wuhan University, Wuhan, China, in 2014, and the Ph.D. degree in science and technology of surveying and mapping from the University of Nottingham, Nottingham, U.K., in 2019.
He is currently a Lecturer with the China University of Mining and Technology, Beijing, China. His research interest focusses on GNSS remote sensing, with particular emphasis on its application for ionosphere monitoring and modeling.
Kefei Zhang received the Ph.D. degree in science and technology of surveying and mapping from Curtin University, Perth, WA, Australia, in 1998.
He is currently a Chair Professor and the Director of the Institute of Resources and Environment Technologies with the China University of Mining and Technology (CUMT), Beijing, China. He is also an honorary Professor with RMIT University, Melbourne VIC, Australia, where he is the Founding Director of the Satellite Positioning for Atmosphere, Climate and Environment (SPACE) Research Centre. He has more than 30 years of research experience in geodesy, satellite positioning, and geospatial sciences in general. He has authored about 400 peer-reviewed publications in these fields and is a coinventor of 30 patents, and attracted in excess of 50 million dollars in funding from the Australian Research Council, the National Science Foundation of China, as well as funds from other national and international governments and industries. He is an Australian pioneer in cutting-edge technologies for smart tracking, GNSS atmospheric sounding for weather and climate, and space environment management. His satellite-to-satellite tracking frontier research has led to over 10-h weather forecast improvement and the successful integration of the GPS radio occultation data into the Australian weather forecasting system in 2012. He currently leads a team of more than 40 researchers at CUMT. He is a Fellow of the International Association of Geodesy. His current research interests are primarily in algorithm development and innovative applications of satellite technologies for high-accuracy positioning, atmospheric studies, space situational awareness, and space resource utilization. Her research interest focusses on InSAR technology for surface deformation monitoring.