Time Series InSAR Three-Dimensional Displacement Inversion Model of Coal Mining Areas Based on Symmetrical Features of Mining Subsidence

: The three-dimensional (3-D) displacements of mining areas is the basis of studying the mining subsidence law and obtaining surface movement parameters. The traditional multi-temporal interferometry synthetic aperture radar (InSAR) technology can only obtain the surface deformation in line-of-sight (LOS) direction, even if some methods can obtain the 3-D displacements of mining area based on InSAR. However, it has high data requirements for data types, which are not conducive to the inversion of 3-D displacements. In this paper, the symmetry of the surface basin caused by mining subsidence under different mining degrees is analyzed. According to the basic symmetrical features of mining subsidence—that the surface vertical displacement and horizontal displacement in near horizontal coal seam is symmetrical with respect to the main section of the basin, combined with time series InSAR technology (TS-InSAR)—a novel method for retrieving the 3-D displacement results from a single-geometry InSAR dataset based on symmetrical features (hereafter referred to as the SGI-SF method) is proposed. The SGI-SF method ﬁrst generates multi-temporal observations of LOS displacement from a single-geometry InSAR dataset, and then transforms them into multitemporal observations of 3-D displacement datasets according to symmetrical features. There is no necessity to obtain the surface movement parameters from the measured data to calculate 3-D displacement ﬁelds. Finally, the time series of 3-D displacements are estimated from multi-temporal 3-D displacements using the singular value decomposition (SVD) method. Nine descending Sentinel-1 images from the Yulin mining area of China are used to test the proposed SGI-SF method. The results show that the average root mean square errors (RMSE) in the vertical and horizontal direction of the three-dimensional deformations are approximately 9.28 mm and 13.10 mm, respectively, which are much smaller than mining-induced displacements and can provide support for deformation monitoring in mining areas.


Introduction
Most advanced monitoring technologies are used to monitor the surface subsidence and movement in the mining area of the working panel regularly and periodically and then to understand the development trend and dynamic changes in real time, thus preventing upcoming disasters and reducing losses [1][2][3][4][5]. Various mining movement parameters can be obtained by obtaining the in-situ surface deformation of observation stations in mining areas [2]. The commonly used methods of monitoring surface displacement are traditional geodetic methods, such as leveling surveying, total station observation, Global Navigation Satellite System (GNSS) measurement [6][7][8], conventional photogrammetry based on airborne imagery [9] and airborne laser (Light Detection and Ranging or LiDAR) [10]. However, these methods are usually time-consuming, climate-sensitive, expensive and labor-intensive and cannot obtain the mining-induced surface deformation of the whole basin [11].
In the past, most of the surface deformation monitoring results based on InSAR of coal mining areas are 1-D LOS deformation. Wright and Stow obtained LOS deformation maps over the Selby coalfield in the United Kingdom from European remote sensing (ERS) satellite SAR images in 1999 and confirmed the subsidence can be measured by D-InSAR to an accuracy of a few centimetres [30]. The application of TS-InSAR in coal mining areas can better reveal the temporal evolution law of mining deformation and SBAS-InSAR performs better than PS-InSAR in deformation monitoring of mining area [31][32][33][34]. Zhang and Wang et al. proposed an algorithm that modified time-series InSAR technique combining persistent scatterers (PSs) and distributed scatterers (DSs) for surface deformation monitoring in coal mining areas, which increased the density of measurement points [35,36]. In addition, the SAR datasets with longer wavelengths (e.g., L band) have a greater ability to measure mining-induced displacement than those with shorter wavelengths (e.g., X and C bands) [37,38].
The LOS displacement result is a superposition of the surface point 3-D displacements in the up-down, north-south and west-east directions [39]. To retrieve full 3-D displacement fields, various methods based on InSAR technology have been developed, such as the integration of measurements from multiple sensors by multitrack D-InSAR [40][41][42], D-InSAR combined with offset tracking [43,44], D-InSAR combined with multi-aperture interferometry [45,46], D-InSAR combined with GPS [47], the method of establishing physical relationship between horizontal displacement and vertical deformation and (hereafter referred to as the PRHV method) [48], and retrieved three-dimensional displacement fields of mining areas from a single InSAR pair (SIP) [49,50]. The PRHV method and SIP-based method are both inversed the 3D deformation by establishing the gradient relationship between subsidence and horizontal movement. The PRHV method is based on the hypothesis proposed by Kratzsch that the horizontal displacement is proportional to the tilt; the SIP-based method is based on the stochastic medium model [51] proposed by Litwiniszyn. The SIP-based method has the problem of error propagation: when calculating the deformation along the north-south and west-east directions, it is affected by the subsidence error, which leads to the error of horizontal displacement becoming larger [50]. Additionally, the SIP-based method needs surface movement parameters, such as a horizontal movement constant and major influence angle, that are obtained from the observation data. The dynamic distribution of vertical displacement and horizontal displacement characteristics are very complicated in the process of mining and the existing mining subsidence theory has not yet established a functional model for vertical displacement with horizontal displacement under dynamic mining condition. Therefore, the 3-D displacement inversion method based on the deterministic function relationship between subsidence and horizontal displacement has some limitations. In this paper, we analyze the symmetry of mining subsidence from the theory of mining subsidence. According to the basic characteristics of the symmetry of the surface subsidence and horizontal displacement in a near-horizontal coal seam, combined with time series InSAR technology, a novel method for retrieving the 3-D displacement results from a single-geometry InSAR dataset based on symmetrical features is proposed. As surface subsidence lags behind underground mining, it is necessary to determine the center of surface movement basin at each interference period of image dataset. Based on the characteristics about the displacement of basin center on near horizontal coal seam, we propose a novel method for automatically searching the center of a basin. In the inversion of 3-D deformation, the proposed method does not need the prior surface movement parameters obtained from the measured data, so it has a wide range of applicability. The three-dimensional deformation field retrieved by SGI-SF method is verified by the sentinel-1 data and measured data of Yulin mining area and the result show the effectiveness of the method.

Symmetry Analysis of Mining Areas Based on Mining Theory
According to the existing theory of mining subsidence, under normal horizontal coal seam mining conditions, because the dip angle of the coal seam is close to zero, the geometric characteristics and physical properties of the overlying rock formation and the surface are all isotropic. The overlying strata and surface subsidence basins caused by mining are symmetrically distributed, and their symmetry axes are major cross section along dip direction and along strike direction, which go through the center of the mined-out areas and are parallel and perpendicular to the direction of the working panel, respectively [52]. The subsidence and horizontal displacement of any point on the surface of the subsidence basin have a symmetrical point about the major cross section along dip direction and along strike direction under a stable surface [53][54][55][56]. Therefore, only half of the observation line is set on one side of the open cutting or stopping line of the working panel, and the characteristics of the whole surface of the moving basin are described by the observed values on the main section of half of the moving basin [52].
According to the basic principle of mining subsidence in horizontal coal seams, the distribution of subsidence and horizontal displacement in the main section under different mining levels is shown in Figure 1 [57,58]. In Figure 1, under the conditions of non-full mining and super-full mining, the subsidence and horizontal displacements in the main section of the movement surface are all symmetrically distributed around the surface point O (the point of maximum subsidence above the center of the basin). The horizontal displacement direction at all points in the surface movement basin points to point O, which corresponds to the center of the mined-out areas, as shown in Figure 2. The major cross section along dip direction Z 1 Z 2 and along strike direction Q 1 Q 2 intersect at the O point, and the surface movement basins are divided into four regions: I, II, III and IV. It can be seen from the horizontal movement curve in Figure 1 and the subsidence isoline in Figure 2 that, (1) any two symmetrical points P 1 and P 3 in symmetrical regions I and III have the same subsidence, and the horizontal displacements are equal in value and opposite in direction, as are the two symmetrical points P 2 and P 4 in symmetrical regions II and IV and, (2) any two symmetrical points P 1 and P 2 in symmetrical regions I and II have the same subsidence and same horizontal displacement along strike direction, and the horizontal displacements along the dip direction are equal in value and opposite in direction, as are the two symmetrical points P 3 and P 4 in symmetrical regions III and IV [52].

Geometric Relationship between the Radar LOS Direction and 3-D Direction
The displacements of the surface point along the LOS direction are projected from the east-west, north-south and up-down 3-D displacements of the movement surface, and the geometric relation between them is shown in Figure 3 [49].
As shown in Figure 3, θ is the incident angle of the radar, and αh is the heading angle of the platform, which is defined as the angle between the flight direction of the satellite and the north direction (the clockwise direction is positive). The relationship between the radar LOS displacements and the 3-D displacements of the surface points can be described by the spatial geometric relationship of Figure 3, which is shown as Equation (1).
As shown in Equation (1), the west-east, north-south and up-down components of the real ground displacements are UN, UE and W, respectively. InSAR-derived deformation is along the radar line-of-sight (LOS) direction, which is the joint contribution of actual 3-D ground surface displacements in accordance with the imaging geometries of SAR sensors [59,60].

Geometric Relationship between the Radar LOS Direction and 3-D Direction
The displacements of the surface point along the LOS direction are projected from the east-west, north-south and up-down 3-D displacements of the movement surface, and the geometric relation between them is shown in Figure 3 [49].
As shown in Figure 3, θ is the incident angle of the radar, and αh is the heading angle of the platform, which is defined as the angle between the flight direction of the satellite and the north direction (the clockwise direction is positive). The relationship between the radar LOS displacements and the 3-D displacements of the surface points can be described by the spatial geometric relationship of Figure 3, which is shown as Equation (1).
As shown in Equation (1), the west-east, north-south and up-down components of the real ground displacements are UN, UE and W, respectively. InSAR-derived deformation is along the radar line-of-sight (LOS) direction, which is the joint contribution of actual 3-D ground surface displacements in accordance with the imaging geometries of SAR sensors [59,60].

Geometric Relationship between the Radar LOS Direction and 3-D Direction
The displacements of the surface point along the LOS direction are projected from the east-west, north-south and up-down 3-D displacements of the movement surface, and the geometric relation between them is shown in Figure 3   As shown in Figure 3, θ is the incident angle of the radar, and α h is the heading angle of the platform, which is defined as the angle between the flight direction of the satellite and the north direction (the clockwise direction is positive). The relationship between the radar LOS displacements and the 3-D displacements of the surface points can be described by the spatial geometric relationship of Figure 3, which is shown as Equation (1).
As shown in Equation (1), the west-east, north-south and up-down components of the real ground displacements are U N , U E and W, respectively. InSAR-derived deformation is along the radar line-of-sight (LOS) direction, which is the joint contribution of actual 3-D ground surface displacements in accordance with the imaging geometries of SAR sensors [59,60].

The SGI-SF Inversion Method at Pixel Level
As the dip angle of the near-horizontal coal seam is close to zero, the subsidence and horizontal displacement in the movement surface state are always symmetrical with the major cross section along strike direction. As shown in Figure 4, the main strike section Z 1 Z 2 is always taken as the symmetry axis of the movement and deformation of the surface point during the advance of the working panel. The subsidence values of P 1 and P 2 and the horizontal displacement value along the strike direction are equal, while the horizontal displacement value along the dip direction is equal and the direction is opposite. In Figure 4, the horizontal displacement direction of P 1 points to the dynamic surface movement basin with center O', and the angle between this direction and the north direction is ω1. Due to the lag in the surface subsidence, the center O' of the mining movement basin is not located in the center of the dynamic mined-out areas and generally tends towards the side of the open cut. The angle between the advancing direction of the working panel and the north direction is the azimuth of the advancing direction of the working panel, and the angle between the horizontal displacement direction of symmetrical points P 1 and P 2 and Z 1 Z 2 is (ω 1 − α).

The SGI-SF Inversion Method at Pixel Level
As the dip angle of the near-horizontal coal seam is close to zero, the subsidence and horizontal displacement in the movement surface state are always symmetrical with the major cross section along strike direction. As shown in Figure 4, the main strike section Z1Z2 is always taken as the symmetry axis of the movement and deformation of the surface point during the advance of the working panel. The subsidence values of P1 and P2 and the horizontal displacement value along the strike direction are equal, while the horizontal displacement value along the dip direction is equal and the direction is opposite. In Figure 4, the horizontal displacement direction of P1 points to the dynamic surface movement basin with center O', and the angle between this direction and the north direction is ω1. Due to the lag in the surface subsidence, the center O' of the mining movement basin is not located in the center of the dynamic mined-out areas and generally tends towards the side of the open cut. The angle between the advancing direction of the working panel and the north direction is the azimuth of the advancing direction of the working panel, and the angle between the horizontal displacement direction of symmetrical points P1 and P2 and Z1Z2 is (ω1-α).  The horizontal displacements of the strike and dip are U T and U I, respectively, and the horizontal displacement components U N and U E are along the north-south and east-west directions; then, the conversion relationship between the two displacement components is as follows: According to the symmetrical characteristics of the moving surface in near-horizontal coal seam mining, the following relations can be obtained: Let the LOS displacement of symmetry points P 1 and P 2 in Figure 4 be d LOS1 and d LOS2 , respectively. Refer to Equation (1) can obtain the relationship between the LOS displacement and the vertical, south-north, east-west displacement W 1 , W 2 , U N1 , U N2 , U E1 and U E2 . According to Equations (1)-(3), the relations between the LOS displacement and the vertical, dip direction, strike direction displacement W 1 , W 2 , U I1 , U I2 , U T1 and U T2 can be obtained: In Equation (4), B = −sin θ cos (αh − 3/2π), C = −sin θ sin (αh − 3/2π). The horizontal displacement components U I1 and U I2 and the horizontal displacement components of strike U T1 and U T2 can be obtained by the difference of the two expressions in Equation (4): Adding the two equations in Equation (4), and then substituting Equation (6), the sinking components W 1 and W 2 can be obtained as follows: Substituting Equations (5) and (6) into Equation (2) obtains the displacement of points P 1 and P 2 in the north-south and east-west directions: Equations (7) and (8) are the inversion models of the three-dimensional displacement from InSAR LOS deformation. The intermediate variables B and C in the formula are the same as in Equation (4).

Deriving Time Series 3-D Displacement Fields
Assuming that there are N+1 SAR images acquired from a single imaging geometry and relative to the same mining areas, we define these acquisition times as t = [t 0 , t 1 . . . t N ]. M small-baseline interferograms are formed from the N+1 SAR images by setting a given threshold for the spatiotemporal baselines. M interference pairs are processed by differential interferometry, and M interference phases are obtained. The minimum cost flow method (MFC) [61] is used to unwrap M interference phases, M unwrapped phases are obtained and the unwrapped phase is expressed as follows: According to Equations (7) and (8), the relationship between the phase and vertical deformation and horizontal deformation can be obtained: Superscripts 1 and 2 represent point 1 and corresponding point 2 in Figure 4, respectively. Other parameters have the same meaning as in Section 2.3, and i represents the i-th interference pair I = [1, 2, 3 . . . . . . M]. We get the vector of the mean subsidence velocity and deformation velocity along north-south and west-east directions between each two time-adjacent SAR acquisitions, and thus, we can obtain the following system for each pixel: where F is a coefficient matrix that depends on the formed small-baseline InSAR pairs. To easily and effectively combine all the available small-baseline interferograms, the minimumnorm criterion and the singular value decomposition (SVD) method are used to obtain the 3-D displacement velocity between each pair of time-adjacent SAR acquisitions from M small-baseline interferograms. Then, the time series 3-D displacements of the mining area can be estimated by the following: After obtaining the three-dimensional deformation, the movement basin can be analyzed and the surface movement parameters can be inversed. In addition, the displacements along the up-down, west-east and north-south directions can convert to displacements along the vertical and horizontal directions.

Determining the Center of the Dynamic Surface Movement Basin
To obtain the azimuth of horizontal displacement ω 1 when the working panel is advancing, it is necessary to determine the dynamic center O' of the surface movement basin. It takes time for the influence of underground mining to propagate from underground to the surface, which makes the surface movement basin center O' lag behind the corresponding central position of the mined-out areas, and the specific lag position is difficult to determine due to different geological mining conditions. According to the basic law of mining subsidence, the dynamic movement basin center O' point of the near horizontal coal seam has following characteristics [10]: 1.
The subsidence value of the center point O' in the moving basin should be the maximum value in the whole major cross section along strike direction; and According to the above characteristics, based on the superposition calculation of the D-InSAR LOS displacement on the major cross line along strike direction Z 1 Z 2 , the moving basin center O' is examined. In the geometric relationship between the LOS displacement and subsidence and the dip horizontal displacement and strike horizontal displacement, since the horizontal displacement perpendicular to the major cross line along the strike direction is zero, the three-dimensional displacement decomposition formula of any point of the major cross line along strike direction can be obtained from Equation (1): The d LOSA and d LOSB in the above formula are superposed and differentiated, respectively, and the subsidence and horizontal displacement values of pixel 1 and pixel 2 are calculated and recorded as the first group displacement (W 1 , U T1 ).

•
The displacement values of pixel 2 and pixel 3 are selected and assigned to d LOSA and d LOSB, respectively. The subsidence and horizontal displacement values of pixel 2 and pixel 3 are calculated from Equations (7) and (8)  calculated and recorded as the first group displacement (W1, UT1).

•
The displacement values of pixel 2 and pixel 3 are selected and assigned to dLOSA and dLOSB, respectively. The subsidence and horizontal displacement values of pixel 2 and pixel 3 are calculated from Equations (7) and (8) and recorded as the second group displacement (W2, UT2).
• This process continues until the subsidence and horizontal displacement values of pixels N-1 and N are calculated and recorded as a group of N-1 data, forming an array of displacement values (W1, UT1), (W2, UT2).... (WN-1, UT(N-1)). The calculation flow is shown in Figure 5.  The "flat bottom" phenomenon appears in the center of the moving basin after the local surface subsidence reaches super-full mining. According to the above search method, multiple groups of pixel regions satisfying the conditions are obtained. In this case, the horizontal displacement in the half basin on one side of the cutting hole points to the first group of qualified pixel positions, while the horizontal displacement in the half basin on the other side of the advancing boundary points to the last group of qualified pixel positions.

Study Area
The Yulin coal mining area (marked by the red dashed line in Figure 6), which is located northwest of the city of Yulin, Shanxi Province, China, was selected to test the SGI-SF method.

Study Area
The Yulin coal mining area (marked by the red dashed line in Figure 6), which is located northwest of the city of Yulin, Shanxi Province, China, was selected to test the SGI-SF method.
The minerals in the geological strata above the coal seam are mainly quartz and feldspar. The surface of the working panel is covered by thick eolian sand and vegetation. To mitigate the potentially adverse effects caused by the interpolation error of the external digital elevation model (DEM) and the atmospheric artifact phase, we chose a study area with a relatively flat terrain (see the black point in Figure 6a).

SAR Data Selection and Data Processing
Sentinel 1 contains C-band data, and the wavelength is 5.6 cm. It is difficult to monitor the large deformation gradient in mining areas due to serious decorrelation phenomena. The dramatic seasonal variation in the land cover can give rise to severe decorrelation in interferograms. To ensure the reliability of the processing results, the winter with a high coherence is selected as the experimental period. We select 9 images from 26 November 2018 to 7 April 2019 for processing. The footprint of these SAR images is marked by the red rectangle in Figure 6a.
To select suitable thresholds for the perpendicular and temporal baselines of the InSAR pairs, we set 500 m and 48 days as the thresholds of the perpendicular and tem-  Figure 6b. The mining length of the 301 working panels is 1078 m, the average mining width is 241 m, the average mining thickness is 5.3 m and the coal seam dip angle is 0.5 • .

SAR Data Selection and Data Processing
Sentinel 1 contains C-band data, and the wavelength is 5.6 cm. It is difficult to monitor the large deformation gradient in mining areas due to serious decorrelation phenomena. The dramatic seasonal variation in the land cover can give rise to severe decorrelation in interferograms. To ensure the reliability of the processing results, the winter with a high coherence is selected as the experimental period. We select 9 images from 26 November 2018 to 7 April 2019 for processing. The footprint of these SAR images is marked by the red rectangle in Figure 6a.
To select suitable thresholds for the perpendicular and temporal baselines of the InSAR pairs, we set 500 m and 48 days as the thresholds of the perpendicular and temporal baselines. After eliminating the interferograms with severe decorrelation, we obtain 18 available InSAR pairs, and their perpendicular and temporal baselines are plotted in Figure 7. We use SRTM DEM and precise orbit data to preprocess the sentinel data and then generate the single look complex (SLC) data after registration. The SRTM DEM is applied to remove the topographic phase component of the interferograms and generate 18 pairs of interference phase, and the temporal and perpendicular baselines are shown in Figure  8. In this procedure, a multi-look operation, of 2:8 pixels in the range and azimuth directions, is carried out. To further suppress the noise in the interferograms, Goldstein filtering is performed. The minimum cost flow algorithm [61,64] is applied to unwrap the differential interferograms with the same reference point. Finally, the SVD method is used to obtain the time series LOS deformation. Before the 3-D displacements are generated, the basic parameters of the working panel and the basic information of the Sentinel 1 satellite parameters need to be determined. We only need to determine the mining azimuth angle and the spatial position of the major cross section along dip direction of the working panel and do not need to determine the mining height, horizontal motion constant, mining depth, or major influence angle of the working panel. When the major cross section along dip direction is uncertain, the maximum subsidence point of the major cross section along dip direction and strike direction can be evaluated by twice basin center search methods to determine the basin center. The incidence angle and heading angle of the SAR data are 36.67° and 67° and 347.66°, respectively. After acquiring the values of these parameters, we apply the SGI-SF method to process the 8 geocoded LOS deformation maps to generate the multi-temporal observations of vertical subsidence.

Estimation of Time Series 3-D Displacements
For each pixel in the study area, the eight unknown 3-D displacements between timeadjacent SAR acquisitions needed to be estimated by 18 multi-temporal observation of phase and transformed into multi-temporal observation of 3-D displacement datasets according. The time series displacements along the up-down, west-east and north-south directions in the study area with respect to the reference acquisition of 26 November 2018 are estimated by Equation (20). For analysis, displacements along the up-down, west-east Subsequently, the two-pass D-InSAR procedure [62] is utilized to process the 18 available InSAR pairs. In this procedure, the 1 arc-second Shuttle Radar Topography Mission (SRTM) DEM is applied to remove the topographic phase component of the interferograms. Additionally, the least squares-based filter [63] and time filter is utilized to suppress the noise in the interferograms.
We use SRTM DEM and precise orbit data to preprocess the sentinel data and then generate the single look complex (SLC) data after registration. The SRTM DEM is applied to remove the topographic phase component of the interferograms and generate 18 pairs of interference phase, and the temporal and perpendicular baselines are shown in Figure 8. In this procedure, a multi-look operation, of 2:8 pixels in the range and azimuth directions, is carried out. To further suppress the noise in the interferograms, Goldstein filtering is performed. The minimum cost flow algorithm [61,64] is applied to unwrap the differential interferograms with the same reference point. Finally, the SVD method is used to obtain the time series LOS deformation. tions are converted to displacements along the up-down, strike and dip directions. The time series displacement results along the up-down, strike and dip directions are shown in Figures 8-10, respectively. With the advancement of underground mining (see the white arrow in Figure 6b), the subsidence basin expanded extensively in the following 5 months. The maximum subsidence value reached −0.30 m, and the maximum absolute horizontal displacement in the strike and dip directions increased to 0.15 m and 0.08 m, respectively.  Before the 3-D displacements are generated, the basic parameters of the working panel and the basic information of the Sentinel 1 satellite parameters need to be determined. We only need to determine the mining azimuth angle and the spatial position of the major cross section along dip direction of the working panel and do not need to determine the mining height, horizontal motion constant, mining depth, or major influence angle of the working panel. When the major cross section along dip direction is uncertain, the maximum subsidence point of the major cross section along dip direction and strike direction can be evaluated by twice basin center search methods to determine the basin center. The incidence angle and heading angle of the SAR data are 36.67 • and 67 • and 347.66 • , respectively. After acquiring the values of these parameters, we apply the SGI-SF method to process the 8 geocoded LOS deformation maps to generate the multi-temporal observations of vertical subsidence.

Estimation of Time Series 3-D Displacements
For each pixel in the study area, the eight unknown 3-D displacements between time-adjacent SAR acquisitions needed to be estimated by 18 multi-temporal observation of phase and transformed into multi-temporal observation of 3-D displacement datasets according. The time series displacements along the up-down, west-east and north-south directions in the study area with respect to the reference acquisition of 26 November 2018 are estimated by Equation (20). For analysis, displacements along the up-down, west-east and north-south directions are converted to displacements along the up-down, strike and dip directions. The time series displacement results along the up-down, strike and dip directions are shown in Figures 8-10, respectively.
With the advancement of underground mining (see the red arrow in Figure 6a), the subsidence of the basin expanded extensively in the following 5 months. The maximum subsidence value reached −0.30 m, and the maximum absolute horizontal displacement in the strike and dip directions increased to 0.15 m and 0.08 m, respectively.
For each pixel in the study area, the eight unknown 3-D displacements between timeadjacent SAR acquisitions needed to be estimated by 18 SIP-derived multi-temporal 3-D displacements. The time series displacements along the up-down, west-east and northsouth directions in the study area with respect to the reference acquisition of 26 November 2018 were estimated by Equation (20).

Verification
The in-situ GPS measurements of the mining-induced vertical subsidence and horizontal motions at 45 surface observation points along the dip observation line (see line AA' in Figure 6) from 10 April 2018 to 17 February 2019 are used as validation data. In these measurements, the observations of the horizontal motions are generated by synthe-

Verification
The in-situ GPS measurements of the mining-induced vertical subsidence and horizontal motions at 45 surface observation points along the dip observation line (see line AA' in Figure 6) from 10 April 2018 to 17 February 2019 are used as validation data. In these measurements, the observations of the horizontal motions are generated by synthe-

Verification
The in-situ GPS measurements of the mining-induced vertical subsidence and horizontal motions at 45 surface observation points along the dip observation line (see line AA' in Figure 6) from 10 April 2018 to 17 February 2019 are used as validation data. In these measurements, the observations of the horizontal motions are generated by synthesizing the 2-D horizontal motions in the east and north directions (or strike and dip directions). The 3-D displacement results on 7 April 2019 are compared with the verification data, and the two time periods are more consistent, which allows us to quantitatively evaluate the accuracies of the estimated 3-D displacements. We then synthesize the accumulated 2-D horizontal motions in the east and north directions. Finally, we compare the vertical subsidence and horizontal motions with the in-situ GPS measurements, and the results are plotted in Figure 11.  The comparison shows that the trend of the deformation reflected by the proposed method is basically consistent with the in-situ GPS observation value. The average root mean square error (RMSE) of the estimated subsidence over all of the observation points is approximately 9.28 mm. The RMSE of the estimated horizontal motions over all of the observation points is approximately 13.10 mm. These RMSEs equate to 3.6% and 6.3% of the maximum SGI-SF estimated 3-D displacements in the vertical (0.305 m) and horizontal (0.213 m) directions, respectively. This result indicates that the SGI-SF estimated time series 3-D displacements are reliable. As InSAR is not sensitive to north-south deformation, the estimation error of horizontal movement is larger than that of subsidence.

Discussion
As we know, the PRHV method proposed by S. Samieie-Esfahany et al. and SIPbased method proposed by Z.W. Li et al. are both inversed the 3-D deformation by establishing the gradient relationship between subsidence and horizontal movement, but the latter method has better performance than the PRHV method because it assumes a simpler gradient relationship. Different from the two methods, this paper proposes SGI-SF method based on the analysis of symmetry characteristics of surface movement caused by coal mining, which has never been used in 3-D deformation inversion of coal mining area. This method has the following advantages: (1) there is no error propagation problem, and there is no need to calculate horizontal displacement by settlement; (2) in the calculation, it is not necessary to know the surface movement parameters obtained from the measured data, but only need to know the location of mined-out areas, which is easier to obtain than Figure 11. Comparison between the in-situ GPS measurements (orange line) and the displacements estimated by the proposed method (blue line) over the field observation points (i.e., profile AA' in Figure 6) in the vertical (a) and horizontal (b) directions. The gray line is the difference between the GPS and InSAR measurements at the monitoring points.
The comparison shows that the trend of the deformation reflected by the proposed method is basically consistent with the in-situ GPS observation value. The average root mean square error (RMSE) of the estimated subsidence over all of the observation points is approximately 9.28 mm. The RMSE of the estimated horizontal motions over all of the observation points is approximately 13.10 mm. These RMSEs equate to 3.6% and 6.3% of the maximum SGI-SF estimated 3-D displacements in the vertical (0.305 m) and horizontal (0.213 m) directions, respectively. This result indicates that the SGI-SF estimated time series 3-D displacements are reliable. As InSAR is not sensitive to north-south deformation, the estimation error of horizontal movement is larger than that of subsidence.

Discussion
As we know, the PRHV method proposed by S. Samieie-Esfahany et al. and SIP-based method proposed by Z.W. Li et al. are both inversed the 3-D deformation by establishing the gradient relationship between subsidence and horizontal movement, but the latter method has better performance than the PRHV method because it assumes a simpler gradient relationship. Different from the two methods, this paper proposes SGI-SF method based on the analysis of symmetry characteristics of surface movement caused by coal mining, which has never been used in 3-D deformation inversion of coal mining area. This method has the following advantages: (1) there is no error propagation problem, and there is no need to calculate horizontal displacement by settlement; (2) in the calculation, it is not necessary to know the surface movement parameters obtained from the measured data, but only need to know the location of mined-out areas, which is easier to obtain than measured data. Therefore, the proposed method has a wider range of application scenarios than other methods.
Existing mining subsidence research confirmed that in the case of near-horizontal coal seam longwall face mining with a small topographic relief, the surface movement basin does have obvious symmetrical morphological characteristics. In the process of advancing the working panel, the dynamic subsidence and horizontal displacement of the surface also have symmetrical characteristics with major cross section along strike direction. In the most ideal state, the surface movement basin above the working panel must have the symmetry characteristics mentioned above, but the measurement results are not completely symmetrical due to the influence of complex geological conditions on the surface movement and deformation. Under complex geological conditions (such as large faults, folds, aquifers or mountain surface conditions), the movement and deformation trends of surface basins are covered, showing unique trends that vary from place to place, and the research in this area is still insufficient, which is a difficult problem in the field of mining damage research. Several complex geologies influence the law of surface subsidence and are listed as follows: (1) the influence of faults on mining subsidence; (2) the characteristics of surface movement and deformation in mountainous terrain; (3) the rich water layer greatly influences the surface movement and deformation; (4) the influence of the mining method and roof management method; and (5) the influence of the uneven nature of the rock stratum [65]. In reality, it is almost impossible to find a place without faults, dip angles, water layers, flat terrain and coal mines, so in most cases, the ground movement of coal mining is affected by complex geological conditions [66].

Conclusions
In this paper, we proposed a new model of time series 3-D displacement inversion from a single-geometry InSAR dataset according to the symmetrical characteristics of subsidence and horizontal displacement in the surface basin of near the horizontal coal seam. We selected nine descending Sentinel-1 images over the Yulin mining area of China from 26 November 2018 to 7 April 2019 to test the proposed SGI-SF method. The results show that the maximum time series 3-D displacements in the time period of SAR data acquisition reach 0.305, 0.151 and 0.82 m in the vertical, strike, and dip directions, respectively. In addition, the accuracy of the SGI-SF method was evaluated by comparing time series 3-D displacements results with in-site GPS or leveling measurements. The comparison result shows that the average root mean square error (RMSE) of the estimated subsidence and horizontal motions over all of the observation points was approximately 9.28 mm and 13.10 mm, respectively. These RMSEs equate to 3.6% and 6.3% of the maximum SGI-SF-estimated 3-D displacements in the vertical (0.305 m) and horizontal (0.213 m) directions, respectively. This result indicates that the SGI-SF-estimated time series 3-D displacements are reliable. Testing of the model with practical examples reveals that the 3-D displacements estimated by the model are effective. In the inversion of 3-D deformation, the proposed method does not need the prior surface movement parameters obtained from the measured data, so it has a wide range of application scenarios.  Data Availability Statement: Sentinel-1 data is provided by the European Space Agency (ESA) and available from the Alaska Satellite Facility (ASF) (https://vertex.daac.asf.alaska.edu). SRTM DEM data is available at https://srtm.csi.cgiar.org/srtmdata/.