Monitoring and analysis of mining 3D deformation by multi-platform SAR images with the probability integral method

Abstract

Only one-dimensional (1D) deformation along the radar line of sight (LOS) can be obtained using differential interferometry synthetic aperture radar (D-InSAR), and D-InSAR observation is insensitive to deformation in the north direction. This study inferred three-dimensional (3D) deformation of a mining subsidence basin by combining the north-south deformation predicted by a probability integral method with the LOS deformation obtained by D-InSAR. The 15235 working face in Fengfeng mining area (Hebei Province, China) was used as the object of study. The north-south horizontal movement was predicted by the probability integral method according to the site’s geological and mining conditions. Then, the vertical and east-west deformation fields were solved by merging ascend-orbit RadarSAT-2, descend-orbit TerraSAR, and predicted north-south deformation based on a least squares method. Comparing with the leveling data, the results show that the vertical deformation accuracy of the experimental method is better than the inversed vertical deformation neglecting the horizontal deformation. Finally, the impact of the relationship between the azimuth of the working face and the SAR imaging geometry on the monitoring of the mining subsidence basin was analyzed. The results can be utilized in monitoring mining subsidence basins by single SAR image sources.

Keywords

D-InSAR / ascend-descend orbit data / subsidence prediction / probability integral method / 3D deformation

Introduction

Coal production is 3.75 billion tonnes in China within 2015, and production is growing at an annual rate of 0.8%. As a result of mining activities, the amount of land affected by subsidence will increase by 65,600 hectares by the year 2020. After underground coal seam exploitation, the overlying strata and the surface will move and be deformed (^{Deng et al., 2014}). This movement and deformation has a range of different effects on the buildings and structures located on the surface within the deformation zone. Therefore, accurate monitoring of the surface subsidence, horizontal movement, and vertical deformation at different locations are important to determine the scale of the surface subsidence, to evaluate the degree of structural damages, and to study the damage assessment standards of different buildings in mining areas (^{Fan, 2010}; ^{Liu et al., 2014}; ^{Fan et al., 2015}).

Differential interferometry synthetic aperture radar (D-InSAR) technology has the advantages of all-day, all-weather, and high spatial coverage. It has been widely used in the monitoring of volcanoes (^{Massonnet et al., 1995}; Henderson and Pritchard, 2013), earthquakes (^{Zebker et al., 1994}; ^{Fialko et al., 2001}), landslides (^{Refice et al., 2000}; ^{Liu et al., 2013}; ^{Motagh et al., 2013}), mining subsidence (Carnec and Delacourt, 2000; ^{Fan et al., 2017}), and other geophysical phenomena. To overcome the temporal and spatial decorrelation effects, some advanced D-InSAR techniques were proposed, such as Persistent Scatterer InSAR (^{Ferretti et al., 2000a}, ^b), Small Baseline Subset InSAR (Berardino et al., 2003), and Temporarily Coherent Point InSAR (^{Zhang et al., 2011}, ²⁰¹²). These technologies have been successfully used to monitor a wide range of surface deformation caused by urban groundwater exploitation (^{Dai et al., 2015}), to evaluate the active layer of soils and the permafrost beneath (^{Chen et al., 2013}), and to diagnose threats to heritage sites (^{Chen et al., 2017}).

The D-InSAR and advanced D-InSAR techniques can only obtain a one-dimensional (1D) deformation of target along the radar line of sight (LOS), and the LOS deformation is usually converted directly to vertical deformation without considering the horizontal movement. However, surface deformation caused by coal mining not only results in large vertical deformation but is often accompanied by large horizontal movement over a small range. When ignoring the effects of horizontal movement, it is difficult to guarantee the reliability of the results when the LOS deformation is transformed to the vertical deformation, and the shape of the acquired surface subsidence basin may also be distorted. Therefore, there is a need to utilize SAR data sources with different imaging geometries or by means of external GPS data, as well as new technologies such as Multiple Aperture InSAR (MAI) (Bechor and Zebker, 2006; ^{Yan et al., 2016}) and offset-tracking (^{Strozzi et al., 2002}) to accurately obtain the three-dimensional (3D) deformation of the surface subsidence basin. ^{Liu et al. (2012)} used a multi-satellite platform to establish a model of a 3D deformation field due to excessive pumping of ground water in Tianjin. ^{Zhu et al. (2014a}, ^b) analyzed the geometric relationship between the vertical, north, and east deformations and the radar LOS deformation using SAR data with three different imaging geometries. Then the vertical and east-west deformations of the mining area were solved using the least squares method. However, the use of multi-source SAR data does not guarantee the accurate extraction of 3D deformation. Due to the polar orbit of the satellite, the D-InSAR technique is not sensitive to the north-south deformation, and the monitoring result of the north-south horizontal movement is unreliable (^{Wright et al., 2004}). Moreover, many mines are located in remote mountainous areas and sufficient data sources cannot be guaranteed. ^{Chen et al. (2012)}, ^{Hu et al. (2013)}, and ^{Wang et al. (2016)} obtained 3D deformation fields of the study area using fusion GPS data and D-InSAR results. However, the characteristics of the non-linear horizontal movement of the mine are not known, therefore it is difficult to guarantee the precision of either the north-south deformation or east-west deformation obtained by GPS data interpolation. MAI (^{Jung et al., 2009}) and offset-tracking (^{Fialko et al., 2005}) techniques have been widely used to determine the 3D deformation caused by earthquakes and volcanoes. However, their accuracy is affected by the scattering characteristics of ground features, so these techniques cannot achieve better results in the complex environment of the mining area.

The probability integral method takes the normal distribution function as the influence function and uses the integral form to express the surface subsidence basin. The basis of this method is the stochastic medium theory. The stochastic medium theory was first used to study the movement of strata by the Polish scholar J. Litwiniszyn in 1957 (^{Litwiniszyn, 1957}). Chinese scholars Liu and Liao (1965) had further developed this technique to create the probability integral method. After decades of mining subsidence research in China, it has become the most mature and widely used prediction method in China (^{Deng et al., 2014}; ^{Fan et al., 2014}). Moreover, the fusion probability integral method and the results of D-InSAR have been well applied in monitoring mining 3D deformation (^{Li et al., 2015}; ^{Diao et al., 2016}), predicting the mining deformation (^{Yang et al., 2016}, ²⁰¹⁷), and large-gradient subsidence (^{Fan et al., 2014}). However, these methods only use a single SAR image source.

In this study, we first considered the characteristics of mining deformation according to the known geological and mining conditions. The surface movement parameters were obtained and the north-south deformation was predicted based on the probability integral method. Then, the vertical and the east-west deformations were acquired by combining the north-south deformation with the LOS deformation obtained from SAR data of ascend-descend orbits based on the least squares principle. The reliability of the vertical deformation was evaluated by a comparison with ground leveling data. Finally, the influence of the relationship between the azimuth of the working face and the imaging geometry of the SAR system on the monitoring of surface subsidence basin was analyzed.

Study area, data and methods

Study area and data

The study area is the 15235 working face in the Fengfeng mining area (Hebei province, China), which is surrounded by a large number of villages, industrial plants, and farmland. There is a circular railway passing through the top of the well field, which provides passenger and freight transport. The strike length and dip width of the working face are 935 m and 142 m, respectively. The average thickness of the coal seam is 5.9 m. The average depth is 740 m. The inclination of the coal seam is 13°. The actual mining thickness of the mine is 4.5 m. In the southwest there are three goafs. An overview of the study area is shown in Fig. 1, in which the arrow is the mining direction and the purple rectangle is the mining range of the working face during the study period. Ground leveling data were observed on December 19, 2015 and January 16, 2016, respectively.

The data source in this experiment are three descending TerraSAR images, collected from 30 December 2015 to 21 January 2016, and two ascending RadarSAT-2 images, collected from 24 December 2015 to 17 January 2016. The coverage of the two types of images are shown in Fig. 1. The parameters related to the two types of images are shown in Table 1. The temporal and perpendicular baselines of the interferometric pairs of the two types of images are shown in Table 2. The temporal and perpendicular baselines of the interferometric pairs are short, which enhance the interference coherence and mitigate residual topographic errors.

D-InSAR monitoring of 3D deformation

The D-InSAR observations only reflect the magnitude of the actual deformation in the LOS direction. If a ground object k has shifted during two imaging periods, then the interferometric phase can be expressed as:

(1){\varphi }_{\mathrm{int}}=\omega \{{\phi }_{\mathrm{los}}+{\phi }_{\mathrm{atm}}+{\phi }_{\mathrm{flat}}+{\phi }_{\mathrm{topo}}+{\phi }_{\mathrm{n}}\},

where \omega represents the wrapped algorithm; {\phi }_{\mathrm{los}}represents the LOS deformation phase; {\phi }_{\mathrm{atm}}represents the atmospheric phase; {\phi }_{\mathrm{flat}}and {\phi }_{\mathrm{topo}}represents the flat ground and topography phases, respectively; and {\phi }_{\mathrm{n}}represents the noise phase.

Through the use of the two-track difference method, {\phi }_{\mathrm{topo}} can be removed by simulating the topography phase from the external digital elevation model (DEM). {\phi }_{\mathrm{flat}}can be removed by baseline estimation. {\phi }_{\mathrm{n}} can be suppressed by adaptive filtering (^{Chen et al., 2014}). Finally, through the use of the minimum cost flow (MCF) unwrapping algorithm, the LOS deformation phase {\phi }_{\mathrm{los}} is acquired from the interferometric phase {\varphi }_{\mathrm{int}}. The LOS deformation {\phi }_{\mathrm{los}} can be obtained by Eq. (2):

(2)d_{\mathrm{los}}=-{\displaystyle \frac{\lambda }{\mathrm{4}\pi }}{\phi }_{\mathrm{los}}.

In order to obtain the 3D deformation field, it is necessary to use the SAR imaging geometry to determine the relationship between the north, east, up direction deformation and the LOS direction deformation. d_{\mathrm{n}}, d_{\mathrm{e}}, and d_{\mathrm{u}} are supposed to be the deformation decomposition of the ground target k in the direction of the coordinate axis of the horizontal coordinate system, the north, east, and up direction deformation is positive, and the d_{\mathrm{los}} point to the satellite direction is positive. The imaging geometry of Fig. 2 can be used to acquire Eq. (3):

(3)d_{\mathrm{los}}=d_{\mathrm{n}}\sin \alpha \sin \theta -d_{\mathrm{e}}\cos \alpha \sin \theta +d_{\mathrm{u}}\cos \theta ,

where \theta is the radar incidence angle at the center of the study region, and \alpha is the clockwise angle between the north direction and the satellite flight direction.

Theoretically, if three types of SAR data sources with different imaging geometries are known, the 3D deformation componentsd_{\mathrm{n}}, d_{\mathrm{e}}, and d_{\mathrm{u}} can be solved according to Eq. (3). However, due to the polar orbit of SAR satellites, the deformation sensitivity of the north-south direction is low (usually 0.06‒0.13) (^{Liu et al., 2012}). For some mines it is impossible to obtain sufficient SAR data sources; in the meanwhile, the north-south horizontal movement of the mining area cannot be ignored. Therefore, combined with the characteristics of mining deformation, the vertical and east-west deformations can be calculated accurately by combining the probability integral method (to predict the north-south horizontal deformation) with the LOS direction deformation obtained from SAR data of ascend-descend orbits.

Probability integral method

To simplify the derivation process, a horizontal rectangular working face is taken as an example. We use l and L to represent the length and width of the working face. The lower left corner of the working face is used as the origin of coordinates. The strike direction is the x-axis and the dip direction is the y-axis (as shown in Fig. 2) (^{Diao et al., 2016}). The surface deformation at any point on the surface is predicted as:

(4)\{\begin{array}{c}W{(x,y)}_k=W^0\left(x\right)W^0\left(y\right)/W_{\max }\\ U{(x,y,\phi )}_k=(U^0\left(x\right)U^0\left(y\right)\cos \phi +U^0\left(y\right)W^0\left(x\right)\sin \phi )/W_{\max }\end{array},

where W_{\max }=qm\cos {\alpha }_{\mathrm{1}} is the maximum subsidence value when mining is fully mined; {\alpha }_{\mathrm{1}} is the dip angle of the coal seam; m is the coal thickness; q is the subsidence factor. W{(x,y)}_kis the subsidence of point k with coordinates (x,  y); U{(x,y,\phi )}_kis the horizontal movement of point k with coordinates (x,  y) along the j direction; j is the horizontal angle that is rotated counterclockwise from the x axis positive direction to the surface movement direction.

Terms W^{\mathrm{0}}(x) and U^{\mathrm{0}}(x) are the subsidence and horizontal movement of those points on the strike main section whose abscissa is x when the coal seam is fully mined in the dip direction, and W^{\mathrm{0}}(y) and U^{\mathrm{0}}(y) are the same but for those points on the dip main section whose ordinate is y when the coal seam is fully mined in the strike direction:

(5)\{\begin{array}{l}{W}^{\mathrm{0}}(x)=W(x)-W(x-l)\\ W^{\mathrm{0}}(y)=W(y)-W(y-L)\\ U^{\mathrm{0}}(x)=U(x)-U(x-l)\\ U^{\mathrm{0}}(y)=U(y)-U(y-L)\end{array},

where, W(x) and U(x) are the predicted subsidence and horizontal movement, respectively, of the strike main section when the coal seam is fully mined in the dip direction, and semi-infinite mined in the strike direction:

(6)\begin{array}{l}W(x)={\displaystyle \frac{W\max }{\mathrm{2}}}\left[erf\left({\displaystyle \frac{\sqrt{\pi }}{r}}x\right)+\mathrm{1}\right]\\ U(x)=bW\max e^{-\pi \frac{x^{\mathrm{2}}}{r^{\mathrm{2}}}},\end{array}

where erf\left({\displaystyle \frac{\sqrt{\pi }}{r}}x\right) is a probability integral function; r={\displaystyle \raisebox{1ex}{$H$}\!\left/ \!\raisebox{-1ex}{$\tan \beta $}\right.} is the main influence radius; H is the depth of mining; is the tangent of the influence angle; b is the horizontal movement coefficient.

Solution strategy

The experimental data are obtained from TerraSAR and RadarSAT-2 sensor, and their imaging geometric parameters are shown in Table 1. The parameters given in Table 1 are incorporated into Eq. (3) to obtain Eq. (7). Due to the lack of available data sources for Eq. (7) the normal equation displays rank-deficiency, which shows that the north-south horizontal movement cannot be ignored when combined with the previous analysis. Therefore, the north-south deformation d_nPIM obtained by the probability integral method referred to in Section 2.3 was taken as a known quantity into Eq. (7) to obtain Eq. (8). Then the d_nPIM and the LOS direction deformation are combined to solve the vertical deformation and east-west deformation.

(7)\{\begin{array}{l}{d}_{\mathrm{los\_R}}=-\mathrm{0.109}{d}_{\mathrm{n}}+\mathrm{0.57}{d}_{\mathrm{e}}+\mathrm{0.814}{d}_{\mathrm{u}}\\ d_{\mathrm{los\_T}}=-\mathrm{0.111}{d}_{\mathrm{n}}-\mathrm{0.648}{d}_{\mathrm{e}}+\mathrm{0.754}{d}_{\mathrm{u}}\end{array},

where d_{\mathrm{los\_R}} and d_{\mathrm{los\_T}} are the LOS deformation obtained by RadarSAT-2 and TerraSAR data, respectively.

(8)\{\begin{array}{l}{d}_{\mathrm{los\_R}}+\mathrm{0.109}{d}_{\mathrm{nPIM}}=\mathrm{0.57}{d}_{\mathrm{e}}+\mathrm{0.814}{d}_{\mathrm{u}}\\ d_{\mathrm{los\_T}}+\mathrm{0.111}{d}_{\mathrm{nPIM}}=-\mathrm{0.648}{d}_{\mathrm{e}}+\mathrm{0.754}{d}_{\mathrm{u}}\end{array}.

Converting Eq. (8) to the form of the error equation:

(9)\mathit{V}=\mathit{A}\widehat{\mathit{X}}-\mathit{L},

where \mathit{A}=\left[\begin{array}{ll}\mathrm{0.57}& \mathrm{0.81}\\ -\mathrm{0.648}& \mathrm{0.754}\end{array}\right]; \widehat{\mathit{X}}={\left[\begin{array}{ll}{d}_e& d_u\end{array}\right]}^{\mathrm{T}}; \mathit{L}={\left[\begin{array}{ll}{d}_{los\_R}+\mathrm{0.109}{d}_{nPIM}& d_{los\_T}+\mathrm{0.111}{d}_{nPIM}\end{array}\right]}^{\mathrm{T}}; V is the residual vector of the observed values. The least squares can be solved for \widehat{\mathit{X}}:

(10)\widehat{\mathit{X}}={[{\mathit{A}}^{\mathrm{T}}\mathit{A}]}^{-\mathrm{1}}{\mathit{A}}^{\mathrm{T}}\mathit{L}.

According to linear algebra, the condition number of the matrix shows the sensitivity of the matrix operation to the error, the closer the condition number is to 1, the better the numerical stability. The condition number of the 2-norm of the coefficient matrix A is 1.28, which is close to 1. Therefore, a stable solution can be obtained.

Results and discussion

Experimental results and analysis

Through the use of the two-track D-InSAR technique, combined with the one arc-seconds Digital Elevation Model (DEM) obtained from the Shuttle Radar Topography Mission (SRTM), the LOS direction displacement of the three interferometric pairs, referred to in Table 2, are obtained. Because of the time span of the two different types of SAR, images are not uniform. Considering the short time span, the deformation velocity changed only a little, then the two sets of SAR data were normalized to the same time interval (20151224‒20160117) according to the average rate in the time period. To compare with the experimental vertical deformation at the same time, the vertical deformation of the two image sources was obtained using Eq. (3) when the horizontal movement is neglected. The results are shown in Figs. 3(a) and 3(b). After the LOS deformation was converted to the vertical deformation, the difference of the vertical deformation between C-band RadarSAT-2 and X-band TerraSAR data was compared. Figure 3(c) shows a map of the difference, while Figure 3(d) shows the corresponding histogram.

Figures 3(a) and (b) show that the shape of the surface subsidence basin and the maximum subsidence value obtained from RadarSAT-2 and TerraSAR data are very different when ignoring the horizontal movement. Because of the short wavelength and high resolution, TerraSAR data at the edge of the surface subsidence basin is better than the RadarSAT-2 data. Figure 3(c) shows that the two vertical deformations are quite different when ignoring the impact of horizontal movement.

Based on the geological and mining conditions of the working face, using the software of the mining subsidence prediction system (Wu and Zhou, 1999), the dynamic parameters of the working face could be obtained. The subsidence coefficient q= 0.4; the horizontal movement coefficient b= 0.25; the main effect angle tangent tanb= 1.79; and the maximum subsidence coefficient angle= 85°. The surface movement corresponding to the image acquisition time period is predicted by the dynamic prediction parameters based the probability integral method referred in Section 2.3, the results are shown in Fig. 4.

RadarSAT-2 images are created by multi-look operations with a pixel spacing of 5.82 m in azimuth and 4.58 m in ground range. TerraSAR images are single-look complex images with a pixel spacing of 1.89 m in azimuth and 2.07 m in ground range. The resolution of the result obtained by the probability integration method depends on the pixel space of SAR images. In order to unify the resolution of the three data sources, the grid size is set to 6 m× 6 m, and a Kriging interpolation of the semivariogram is used to ensure that the resolution of the results are consistent. Then, the vertical and east-west deformations are solved according to the solution strategy described in Section 2.4. Finally, the results are processed by low-pass filtering for noise reduction. The results are shown in Fig. 5.

As can be seen from the vertical deformation in Fig. 5(a), during the observation period the surface subsidence is between 0–35 mm, and the center of the surface subsidence basin is located directly above the goaf but not directly above the area that is being mined. This is because of the hysteresis of surface subsidence caused by mining underground working face. A regular elliptical shape did not appear because of the influence by the surrounding goaf, instead, the surface subsidence basin expanded toward the southwest. A small subsidence funnel appeared in the southeast of the surface subsidence basin. This may have been caused by movement of the overlying strata of due to the coal mining. Moreover, the activation range of overlying strata is larger than the size of the goaf. Therefore, when the adjacent working face is being mined, the overlying strata become more active and result in the faster subsidence velocity of the southeastern area. Comparing Figs. 3(a) and (b), it can be seen that the shape and edge of the surface subsidence basin are in accordance with the law of mining subsidence.

The east-west horizontal movement is between –15 and 13 mm, as shown in Fig. 5(b). With regard to the horizontal movement, the west side of the surface subsidence basin moves towards the east, the east side of the surface subsidence basin moves towards the west, and the horizontal movement value in the center of the surface subsidence basin is close to zero. The magnitude of the horizontal movement from the center to both sides first increases and then decreases, which is in line with the general law of mining subsidence. By comparison with Fig. 4(b), the east-west horizontal movement predicted by the probability integration method is between –18–18 mm. The two results agree well with the magnitude, and the position of the maximum horizontal movement is similar. Therefore, the results of both kinds of monitoring can be proven to be correct. In contrast to Fig. 3(c), it can be seen that the difference in the results of the vertical deformation calculated by the two single image sources are consistent with the trend distribution of the east-west horizontal movement result. It is apparent that neglecting the horizontal movement has a large influence on the vertical deformation.

To validate the accuracy of the vertical deformation obtained by the experimental method, ground leveling data of 39 observation stations in the study area are used for analysis. The specific locations of these observation stations are shown in Fig. 1. Because the observation time is not consistent with the D-InSAR monitoring time, the time interval is normalized by interpolation according to the subsidence velocity of each leveling point. The results of the comparison are shown in Fig. 6. It can be seen from Fig. 6 that the experimental results are consistent with the subsidence trend and magnitude of the ground leveling data, and the position of the maximum subsidence is consistent with the leveling data, which shows that the vertical deformation obtained by the experimental method is reliable.

To further illustrate the effect of horizontal movement on vertical subsidence, the results are compared with the vertical subsidence obtained by the single image source without considering horizontal movement. The results are shown in Fig. 6, where it can be seen that TerraSAR and RadarSAT-2 data are slightly worse than the experimental results. In addition, the surface subsidence basin obtained by TerraSAR data shifts to the northeast, and the surface subsidence basin obtained by RadarSAT-2 data shifts to the southwest. This can be observed from the locations of the ground observation stations given in Fig. 1. The reason is that the incidence direction of the two images is reversed, the results are shifted in the opposite direction when the horizontal deformation is ignored. However, the offset of RadarSAT-2 monitoring results is smaller than the TerraSAR results. This is because the incidence angle of the TerraSAR satellite is larger than that of the RadarSAT-2. This can be confirmed by the comparison of Figs. 3(a) and 3(b) with Fig. 5(a).

The accuracy of the three methods is quantified through the maximum deviation (MD), root mean square error (RMSE), and standard deviation (STD) with the results presented in Table 3.

From Table 3, the results of the quantitative analysis show that the accuracy of the experimental method is superior to the monitoring accuracy of the single image. Moreover, because of the small size and the flat terrain of the study area in conjunction with the short baselines of interferograms, the errors from residual topographical heights and atmospheric disturbances can be ignored. Therefore, the accuracy of the experimental method can reach the millimeter level. As can be seen from Eq. (3), when the effect of horizontal movement is ignored and the LOS direction deformation is a constant value, the vertical deformation increases with the increase of incidence angle. Therefore, the deviation between monitoring data and the actual subsidence should be greater, which is demonstrated by the accuracy of the quantitative analysis.

The relationship between the azimuth of the working face and SAR imaging geometry and its impact on monitor the surface subsidence basin

When studying the law of surface subsidence caused by mining, the 3D space problem of surface movement is usually divided into two plane problems: the main section along the strike and the main section along the dip. Therefore, in order to analyze the surface subsidence basin monitored by a single SAR image source, the main sections of the strike and dip are taken as the research objective. The influence of the relationship between the azimuth of the working face and the SAR imaging geometry on monitoring results is investigated and explained by taking the experimental data as an example.

The direction of horizontal movement of the main section is toward the goaf center and the horizontal movement of the goaf center is zero. Therefore, the goaf center O is taken as the demarcation point, and the strike and dip main sections are divided into four parts L, R, U, and D, as shown in Fig. 7. According to the azimuth of the working face β, the horizontal movement of the main section could be decomposed into the north and east directions, and Eq. (11) could be obtained by further developing Eq. (3). Directions to the north, east, and up are positive.

(11)\{\begin{array}{l}{d}_{\mathrm{los\_L}}=-U_y^{\mathrm{0}}\sin \beta \sin \alpha \sin \theta -U_y^{\mathrm{0}}\cos \beta \cos \alpha \sin \theta +d_{\mathrm{u\_L}}\cos \theta \\ d_{\mathrm{los\_R}}=U_y^{\mathrm{0}}\sin \beta \sin \alpha \sin \theta +U_y^{\mathrm{0}}\cos \beta \cos \alpha \sin \theta +d_{\mathrm{u\_R}}\cos \theta \\ d_{\mathrm{los\_U}}=-U_x^{\mathrm{0}}\cos \beta \sin \alpha \sin \theta +U_x^{\mathrm{0}}\sin \beta \cos \alpha \sin \theta +d_{\mathrm{u\_U}}\cos \theta \\ d_{\mathrm{los\_D}}=U_x^{\mathrm{0}}\cos \beta \sin \alpha \sin \theta -U_x^{\mathrm{0}}\sin \beta \cos \alpha \sin \theta +d_{\mathrm{u\_D}}\cos \theta \end{array},

(12)\{\begin{array}{l}{U}_x^{\mathrm{0}}=U(x)-U(x-l)\\ U_y^{\mathrm{0}}=U(y)-U(y-L)+W_y^{\mathrm{0}}\cot {\theta }_{\mathrm{1}}\end{array},

(13) U(x)=b{W}_{\max }{e}^{-\pi {\left(\frac{x}{r}\right)}^{\mathrm{2}}},U(y)=b{W}_{\max }{e}^{-\pi {\left(\frac{y}{r_{\mathrm{1}}}\right)}^{\mathrm{2}}},U(y-L)=b{W}_{\max }{e}^{-\pi {\left(\frac{y-L}{r_{\mathrm{2}}}\right)}^{\mathrm{2}}},

(14)W_y^{\mathrm{0}}={\displaystyle \frac{W_{\max }}{\mathrm{2}}}\left\{\left[erf\left(\sqrt{\pi }{\displaystyle \frac{y}{r_{\mathrm{1}}}}\right)+\mathrm{1}\right]-\left[erf\left(\sqrt{\pi }{\displaystyle \frac{y-L}{r_{\mathrm{2}}}}\right)\right]\right\},

where U_x^{\mathrm{0}}, U_y^{\mathrm{0}} is the horizontal movement of the strike and dip main section, respectively; l and L are the strike length and dip width of the working face, respectively; b is the horizontal movement coefficient; {\theta }_{\mathrm{1}} is the propagation angle of extraction. r, r₁, r₂ are the influence radii of the strike direction, the downhill direction and the uphill direction, respectively; W_{\max }=mq\cos {\alpha }_{\mathrm{1}} is the maximum subsidence value, {\alpha }_{\mathrm{1}} is the coal seam dip, when the {\alpha }_{\mathrm{1}}\le \mathrm{15}°, W_y^{\mathrm{0}}\cot {\theta }_{\mathrm{1}} can be ignored; d_{\mathrm{los\_L}}, d_{\mathrm{los\_R}}, d_{\mathrm{los\_U}}, d_{\mathrm{los\_D}} are the LOS direction deformations of the left, right, top, and bottom of the working face respectively; d_{\mathrm{u\_L}}, d_{\mathrm{u\_R}}, d_{\mathrm{u\_U}}, d_{\mathrm{u\_D}} are the vertical deformations of the left, right, top, and bottom of the working face respectively.

When the horizontal movement is neglected, the vertical deformation of the left d_{\mathrm{u\_L}}^{\text{'}}, right d_{\mathrm{u\_R}}^{\text{'}}, top d_{\mathrm{u\_U}}^{\text{'}}, and bottom d_{\mathrm{u\_D}}^{\text{'}} of the working face on the main section is shown as:

(15)\{\begin{array}{l}{d}_{\mathrm{los\_L}}=d_{\mathrm{u\_L}}^{\text{'}}\cos \theta \\ d_{\mathrm{l}\mathrm{os\_R}}=d_{\mathrm{u\_R}}^{\text{'}}\cos \theta \\ d_{\mathrm{los\_U}}=d_{\mathrm{u\_U}}^{\text{'}}\cos \theta \\ d_{\mathrm{los\_D}}=d_{\mathrm{u\_D}}^{\text{'}}\cos \theta \end{array}.

The vertical deformation error can be obtained from Eq. (11) and Eq. (15) when horizontal movement is neglected, \Delta d_{\mathrm{u}}=d_{\mathrm{u}}-d{\text{'}}_{\mathrm{u}}:

(16)\{\begin{array}{l}\Delta d_{\mathrm{u\_L}}\cos \theta =U_y^{\mathrm{0}}\sin \beta \sin \alpha \sin \theta +U_y^{\mathrm{0}}\cos \beta \cos \alpha \sin \theta \\ \Delta d_{\mathrm{u\_R}}\cos \theta =-U_y^{\mathrm{0}}\sin \beta \sin \alpha \sin \theta -U_y^{\mathrm{0}}\cos \beta \cos \alpha \sin \theta \\ \Delta d_{\mathrm{u\_U}}\cos \theta =U_x^{\mathrm{0}}\cos \beta \sin \alpha \sin \theta -U_x^{\mathrm{0}}\sin \beta \cos \alpha \sin \theta \\ \Delta d_{\mathrm{u\_D}}\cos \theta =-U_x^{\mathrm{0}}\cos \beta \sin \alpha \sin \theta +U_x^{\mathrm{0}}\sin \beta \cos \alpha \sin \theta \end{array}.

For the mining areas with known geological conditions, the horizontal movement on the main section can be predicted by the probability integral method; the vertical deformation on the main section monitored by a single source can be corrected, according to Eq. (16). For the mining areas whose azimuth of the working face is known but the geological conditions are unknown, the monitoring results can be qualitatively analyzed according to Eq. (16).

Taking the 15235 working face as an example, we can know from the geological conditions that the dip angle of the coal seam is 13° (west high), so the surface subsidence basin is an asymmetric ellipsoid biased towards the downhill direction (east direction), and the magnitude of horizontal movement in the downhill direction is large, which can be confirmed from the results of the prediction and monitoring in Fig. 4(b) and Fig. 5(b). Bringing the azimuth angle of the working face (50.9°), the incidence angle, and the heading angle of TerraSAR data into Eq. (16) yields results:

(17)\{\begin{array}{l}\Delta d_{\mathrm{u\_L}}=\mathrm{0.427}{U}_y^{\mathrm{0}}\\ \Delta d_{\mathrm{u\_R}}=-\mathrm{0.427}{U}_y^{\mathrm{0}}\\ \Delta d_{\mathrm{u\_U}}=-\mathrm{0.759}{U}_x^{\mathrm{0}}\\ \Delta d_{\mathrm{u\_D}}=\mathrm{0.759}{U}_x^{\mathrm{0}}\end{array}.

As can be seen from Eq. (17), when the effect of horizontal movement is neglected, in the left half of the dip main section, the results of TerraSAR monitoring become smaller and the results of the right half become larger. In the up half of the strike main section, the results of TerraSAR monitoring become larger and the results in the down half become smaller. This difference becomes more obvious when the incidence angle increases. Because the observation station and the strike main section of the working face are approximately parallel, as shown in Fig. 1, the leveling data could be used to verify the results of the analysis. As can be seen from Fig. 6, the results of TerraSAR monitoring are larger than the leveling data in the up half of the strike main section, and the monitoring results are smaller than the leveling data in the down half, which is consistent with the results of the analysis using Eq. (17).

Therefore, when monitoring the surface subsidence basin with a single kind of images and the vertical deformation is calculated without considering the horizontal movement of the ground surface, the reliability of the monitoring results can be analyzed qualitatively using the relationship between the azimuth of the working face and the imaging geometry of the SAR system. The monitoring results of the main section can be corrected using Eq. (16) with the known geological and mining conditions.

Conclusions

Through the method proposed in this paper, mining 3D deformation can be accurately obtained. By comparing experimental results with ground leveling data and monitoring the results from a single SAR data source, the reliability of the vertical deformation and the east-west horizontal movement is proven. The RMSE and STD of the experimental method are 2.2 mm and 1.5 mm, respectively, in this test case. This provides a feasible method for monitoring mining 3D deformation when lacking a SAR data source.

The influence of the relationship between the azimuth of the working face and the imaging geometry of the SAR system on the monitoring of mining subsidence basins is analyzed. The results of the analysis are consistent with the experimental and ground leveling data and can be used as a guide when analyzing the results of monitoring mining subsidence basin by a single SAR image source.

When monitoring mining deformation without considering the horizontal movement, the effect of the incidence angle on the monitoring results is very significant. The larger the incidence angle, the more obvious the distortions of the monitored surface subsidence basin. Therefore, the SAR data source should be carefully selected to monitor the mining subsidence basin.

References

[1] Bechor N B D, Zebker H A (2006). Measuring two-dimensional movements using a single InSAR pair. Geophys Res Lett, 33(16): L16311

[2] Berardino P, Fornaro G, Lanari R, Sansosti E (2003). A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms. IEEE Trans Geosci Remote Sens, 40(11): 2375–2383

[3] Carnec C, Delacourt C (2000). Three years of mining subsidence monitored by SAR interferometry, near Gardanne, France. J Appl Geophys, 43(1): 43–54

[4] Chen B Q, Deng K Z, Fan H D (2014). Combining D-InSAR and SVR for monitoring and prediction of mining subsidence. Journal of China University of Mining & Technology, 43(5): 880–886 (in Chinese)

[5] Chen F L, Guo H D, Ma P F, Lin H, Wang C, Ishwaran N, Hang P (2017). Radar interferometry offers new insights into threats to the Angkor site. Sci Adv, 3(3): e1601284

[6] Chen F L, Lin H, Zhou W, Hong T H, Wang G (2013). Surface deformation detected by ALOS PALSAR small baseline SAR interferometry over permafrost environment of Beiluhe section, Tibet Plateau, China. Remote Sens Environ, 138(2): 10–18

[7] Chen Q, Liu G X, Hu J C, Ding X L, Yang Y H (2012). Mapping ground 3-D displacement with GPS and PS-InSAR networking in the Pingtung area, southwestern Taiwan, China. Chin J Geophys, 55(10): 3248–3258

[8] Dai K R, Liu G X, Li Z H, Li T, Yu B, Wang X W, Singleton A (2015). Extracting vertical displacement rates in Shanghai (China) with multi-platform SAR images. Remote Sens, 7(8): 9542–9562

[9] Deng K Z, Tan Z X, Jiang Y, Dai H Y, Shi Y, Xu L J (2014). Deformation Monitoring and Subsidence Engineering. Xuzhou: China University of Mining and Technology Press (in Chinese)

[10] Diao X P, Wu K, Hu D H, Li L, Zhou D W (2016). Combining differential SAR interferometry and the probability integral method for three-dimensional deformation monitoring of mining areas. Int J Remote Sens, 37(21): 5196–5212

[11] Fan H D (2010).  Study on Several Key Algorithms of InSAR Technique and its Application. Dissertation for PhD degree. Xuzhou: China University of Mining and Technology (in Chinese)

[12] Fan H D, Gao X X, Yang J K, Deng K Z, Yu Y (2015). Monitoring mining subsidence using a combination of phase-stacking and offset-tracking methods. Remote Sens, 7(7): 9166–9183

[13] Fan H D, Gu W, Qin Y, Xue J Q, Chen B Q (2014). A model for extracting large deformation mining subsidence using D-InSAR technique and probability integral method. Trans Nonferrous Met Soc China, 24(4): 1242–1247

[14] Fan H D, Xu Q, Hu Z B, Du S (2017). Using temporarily coherent point interferometric synthetic aperture radar for land subsidence monitoring in a mining region of western China. J Appl Remote Sens, 11(2): 026003

[15] Ferretti A, Prati C, Rocca F (2000a). Analysis of permanent scatterers in SAR interferometry. In: Proceedings of IEEE 2000 International. Hawaii: Honolulu, 761–763

[16] Ferretti A, Prati C, Rocca F (2000b). Nonlinear subsidence rate estimation using permanent scatterers in differential SAR interferometry. IEEE Trans Geosci Remote Sens, 38(5): 2202–2212

[17] Fialko Y, Sandwell D, Simons M, Rosen P (2005). Three-dimensional deformation caused by the Bam, Iran, earthquake and the origin of shallow slip deficit. Nature, 435(7040): 295–299

[18] Fialko Y, Simons M, Agnew D (2001). The complete (3-d) surface displacement field in the epicentral area of the 1999 M_W 7.1 Hector mine earthquake, California, from space geodetic observations. Geophys Res Lett, 28(16): 3063–3066

[19] Henderson S T, Pritchard M E (2013). Decadal volcanic deformation in the Central Andes Volcanic Zone revealed by InSAR time series. Geochem Geophys Geosyst, 14(5): 1358–1374

[20] Hu J, Li Z W, Zhu J J, Ding X L, Wang C C, Feng G C, Sun Q (2013). Measuring three-dimensional surface displacements from combined InSAR and GPS data based on BFGS method. Chin J Geophys, 56(1): 117–126

[21] Jung H S, Won J S, Kim S W (2009). An improvement of the performance of multiple-aperture SAR interferometry (MAI). IEEE Trans Geosci Remote Sens, 47(8): 2859–2869

[22] Li Z W, Yang Z F, Zhu J J, Hu J, Wang Y J, Li P X, Chen G L (2015). Retrieving three-dimensional displacement fields of mining areas from a single InSAR pair. J Geod, 89(1): 17–32

[23] Litwiniszyn J (1957). The theories and model research of movements of ground masses. In: Proceedings of the European congress ground movement. UK: Leeds, 203–209

[24] Liu B C, Liao G H (1965). The Basic Law of Surface Movement in Coal Mine. Beijing: China Industry Press(in Chinese)

[25] Liu G X, Zhang R, Li T, Yu B, Li T, Jia H G, Nie Y J (2012). Extracting 3D ground deformation velocity field by multi-platform persistent scatterer SAR interferometry. Chin J Geophys, 55(8): 2598–2610 (in Chinese)

[26] Liu P, Li Z, Hoey T, Kincal C, Zhang J, Zeng Q, Muller J P (2013). Using advanced InSAR time series techniques to monitor landslide; movements in Badong of the Three Gorges region, China. Int J Appl Earth Obs Geoinf, 21(1): 253–264

[27] Liu X F, Deng K Z, Fan H D, Wang J T (2014). Study of old goaf residual deformation monitoring based on D-InSAR techniques. Journal of China Coal Society, 39(3): 467–472 (in Chinese)

[28] Massonnet D, Briole P, Arnaud A (1995). Deflation of Mount Etna monitored by spaceborne radar interferometry. Nature, 375(6532): 567–570

[29] Motagh M, Wetzel H U, Roessner S, Kaufmann H (2013). A TerraSAR-X InSAR study of landslides in southern Kyrgyzstan, Central Asia. Remote Sens Lett, 4(7): 657–666

[30] Refice A, Bovenga F, Wasowski J, Guerriero L (2000). Use of InSAR data for landslide monitoring: a case study from southern Italy. In: Proceedings of IEEE 2000 International. Hawaii: Honolulu, 2504–2506

[31] Strozzi T, Luckman A, Murray T, Wegmuller U, Werner C L (2002). Glacier motion estimation using SAR offset-tracking procedures. IEEE Trans Geosci Remote Sens, 40(11): 2384–2391

[32] Wang X, Zhang J, Zhang Q, Zhao C (2016). Inferring multi-dimensional deformation filed in Xi’an by combining InSAR of ascending and descending orbits with GPS data. Acta Geodaetica et Cartographica Sinica, 45(7): 810–817 (in Chinese)

[33] Wright T J, Parsons B E, Lu Z (2004). Toward mapping surface deformation in three dimensions using InSAR. Geophys Res Lett, 31(1): L01607

[34] Wu K, Zhou M (1999). Mining Subsidence Prediction System. Xuzhou: China University of Mining and Technology Press (In Chinese)

[35] Yan S Y, Ruan Z X, Liu G, Deng K Z, Lv M, Perski Z (2016). Deriving ice motion patterns in mountainous regions by integrating the intensity-based pixel-tracking and phase-based D-InSAR and MAI approaches: a case study of the Chongce glacier. Remote Sens, 8(7): 611

[36] Yang Z F, Li Z W, Zhu J J, Hu J, Wang Y J, Chen G L (2016). InSAR-based model parameter estimation of probability integral method and its application for predicting mining-induced horizontal and vertical displacements. IEEE Trans Geosci Remote Sens, 54(8): 4818–4832

[37] Yang Z F, Li Z W, Zhu J J, Preusse A, Yi H W, Wang Y J, Papst M (2017). An extension of the InSAR-based probability integral method and its application for predicting 3-D mining-induced displacements under different extraction conditions. IEEE Trans Geosci Remote Sens, 55(7): 3835–3845

[38] Zebker H A, Rosen P A, Goldstein R M, Gabriel A, Werner C L (1994). On the derivation of coseismic displacement fields using differential radar interferometry: the Landers earthquake. J Geophys Res Solid Earth, 99(B10): 19617–19634

[39] Zhang L, Ding X L, Lu Z (2011). Ground settlement monitoring based on temporarily coherent points between two SAR acquisitions. ISPRS J Photogramm Remote Sens, 66(1): 146–152

[40] Zhang L, Lu Z, Ding X L, Jung H S, Feng G, Lee C W (2012). Mapping ground surface deformation using temporarily coherent point SAR interferometry: application to Los Angeles basin. Remote Sens Environ, 117(1): 429–439

[41] Zhu C G, Deng K Z, Zhang J X, Zhang Y H, Fan H D, Zhang L Y (2014a). Three-dimensional deformation field detection based on multi-source SAR imagery in mining area. Journal of China Coal Society, 39(4): 673–678 (in Chinese)

[42] Zhu C G, Zhang J X, Deng K Z, Zhang Y H, Fan H D, Li P X (2014b). Three-dimensional displacement field of buildings detection from multi-source SAR imagery. Journal of China University of Mining and Technology, 43: 701–706 (in Chinese)

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag GmbH Germany, part of Springer Nature