Deriving Dynamic Subsidence of Coal Mining Areas Using InSAR and Logistic Model

Abstract

The seasonal variation of land cover and the large deformation gradients in coal mining areas often give rise to severe temporal and geometrical decorrelation in interferometric synthetic aperture radar (InSAR) interferograms. Consequently, it is common that the available InSAR pairs do not cover the entire time period of SAR acquisitions, i.e., temporal gaps exist in the multi-temporal InSAR observations. In this case, it is very difficult to accurately estimate mining-induced dynamic subsidence using the traditional time-series InSAR techniques. In this investigation, we employ a logistic model which has been widely applied to describe mining-related dynamic subsidence, to bridge the temporal gaps in multi-temporal InSAR observations. More specifically, we first construct a functional relationship between the InSAR observations and the logistic model, and we then develop a method to estimate the model parameters of the logistic model from the InSAR observations with temporal gaps. Having obtained these model parameters, the dynamic subsidence can be estimated with the logistic model. Simulated and real data experiments in the Datong coal mining area, China, were carried out in this study, in order to test the proposed method. The results show that the maximum subsidence in the Datong coal mining area reached about 1.26 m between 1 July 2007 and 28 February 2009, and the accuracy of the estimated dynamic subsidence is about 0.017 m. Compared with the linear and cubic polynomial models of the traditional time-series InSAR techniques, the accuracy of dynamic subsidence derived by the logistic model is increased by about 50.0% and 45.2%, respectively.

1. Introduction

Mapping ground surface dynamic subsidence caused by underground extraction is essential to assess mining-related geo-hazards and understand the dynamics of mining subsidence. The potential of time-series interferometric synthetic aperture radar (TS-InSAR) techniques, such as the persistent scatterer interferometry and the small baseline subset technique, have been investigated for measuring ground surface time-series deformation caused by mining activities [1,2,3,4,5,6,7]. However, these studies have mainly focused on abandoned coal mines, salt mines, iron ore areas, etc., where slow deformation (e.g., on the order of several or dozens of millimeters per year) happens, on the basis of a linear deformation model (i.e., mean velocity). If these traditional TS-InSAR techniques with a linear deformation model are applied to estimate dynamic subsidence associated with coal mining, large errors may occur. This is because coal mining-related dynamic subsidence is generally characterized by high speed (e.g., dozens of centimeters per month) and high nonlinearity [8,9,10,11,12].

To overcome the above-mentioned limit of the TS-InSAR techniques, a strategy which first estimates the mean deformation between two time-adjacent SAR acquisitions and then integrates them to determine the time-series deformation has been proposed [13,14,15]. It is an alternative strategy for deriving ground surface time-series deformation with a highly nonlinear behavior. This is because the core idea behind this strategy approximates the highly nonlinear time-series deformation as quickly as possible by accumulating the differential displacements between two time-adjacent SAR acquisitions. Theoretically, this approximation becomes more accurate if the time interval of the time-adjacent SAR acquisitions is reduced [16].

However, this strategy also has some limitations. For instance, we need to estimate the N differential displacements between the two time-adjacent SAR acquisitions while N + 1 SAR acquisitions are obtained over the mining area. Therefore, the many unknown parameters may lead to over-parameterization, resulting in instability of the solutions of the differential displacements [17]. Furthermore, the errors in any one of the estimated displacements between two time-adjacent SAR acquisitions would be propagated to the following time-series deformation by the accumulation operation. Most importantly, the InSAR observations have to cover the entire time period of SAR acquisitions, i.e., no temporal gaps occur in the InSAR observations. Due to the lack of observations in the temporal gaps, the displacements between two time-adjacent SAR acquisitions are nearly impossible to estimate accurately. Hence, this strategy will fail if the available InSAR observations contain temporal gaps.

It is common that temporal gaps exist in multi-temporal InSAR observations in coal mining areas. The temporal gaps occur for the following reasons. The main limitation in the detection of ground surface deformation using the InSAR technique is the loss of coherence (i.e., decorrelation) [18]. For an InSAR pair with severe decorrelation, the speckle noise phase usually dominates the interferometric phase of the InSAR pair [18,19]. This significantly hampers the accurate detection of ground surface deformation. However, many factors, such as the dramatic temporal variation of the land cover and the large baselines, may give rise to severe geometrical and temporal decorrelation in interferograms [19,20,21]. The decorrelation phenomenon is more common in coal mining areas due to their non-urban land coverage and large deformation gradients (we illustrate this issue in Section 2). As a consequence, it is commonplace that only InSAR pairs with temporal gaps are available to estimate ground surface time-series deformation.

In this work, we propose a method for deriving mining-induced dynamic subsidence using InSAR by bridging the temporal gaps in multi-temporal InSAR observations with a logistic model. The logistic model has been widely applied to describe the temporal behavior of mining-induced subsidence at a single ground surface point [12,22,23]. Only three unknown model parameters require estimating if the logistic model is used to describe or predict mining-related dynamic subsidence. Thus, the proposed method first constructs a functional relationship between the InSAR observations and the model parameters of the logistic model. An algorithm based on the genetic algorithm and the simplex algorithm is then used for estimating these three model parameters from the InSAR observations with temporal gaps. Finally, the dynamic subsidence in the mining area is estimated using the logistic model.

In this paper, we first analyze the factors causing temporal gaps in multi-temporal InSAR observations, and the effect of these gaps on dynamic subsidence estimation with the traditional TS-InSAR techniques. Subsequently, the proposed method is described in Section 3. The simulated and real data experiments are described in Section 4 and Section 5, respectively. Finally, a comparison between the dynamic subsidence estimated by the logistic model, the linear model, and the cubic polynomial model of the traditional TS-InSAR techniques is provided in Section 6.

2. Factors Causing Temporal Gaps in InSAR Observations in Coal Mining Areas

Severe temporal and geometrical decorrelation in interferometric radar echoes is the main factor causing the temporal gaps in InSAR observations. Generally, it is common that severe temporal and geometric decorrelation occurs in coal mining areas. In this section, we analyze the reasons behind this issue.

2.1. Variation of Non-Urban Land Cover

Most coal mining areas are located in non-urban areas covered by farming land and/or vegetation. Therefore, the variation of the land cover of mining areas in the physical distribution of the elementary scatterers are more dramatic than that of urban areas, resulting in more severe temporal decorrelation in interferograms than that of urban areas [19,21]. Additionally, anthropogenic changes, such as farming activities and construction works, also worsen the temporal decorrelation.

2.2. Large Deformation Gradients

The phase gradient between two spatially adjacent pixels of an interferogram should be less than π for the sake of estimating the integer ambiguity number from the wrapped phase [24]. Hence, decorrelation or phase discontinuity would occur once the surface deformation gradient exceeds the maximum detectable one [25,26]. However, ground surface displacements associated with underground coal mining are usually characterized by large deformation gradients and high speeds [11,27]. For example, in some coal mining sites in China, deformation of dozens or even hundreds of centimeters per month can occur within a small area (e.g., 1 km by 1 km) [9,28,29,30]. Such large deformation gradients in coal mining areas cause severe decorrelation or phase discontinuity in interferograms.

2.3. The Large Spatial-Temporal Baselines of InSAR Pairs

The large spatial-temporal baselines of InSAR pairs are an important factor causing decorrelation. Firstly, the large spatial baseline leads to a difference in the local incidence angle between the two SAR acquisitions [20], thereby decreasing the interferometric coherence and even triggering decorrelation (more than the critical baseline). Additionally, the long temporal baseline generally means that more severe seasonal variation takes place in the surface cover, causing more severe temporal decorrelation [18].

2.4. Limited SAR Data

The wavelengths of most spaceborne SAR satellite systems (e.g., ERS-1/2, Envisat, TerraSAR-X, Radarsat-2) are C-band or X-band, whereas only a few ones (e.g., JERS and PALSAR-1/2) feature an L-band instrument. Generally, SAR sensors with the longer L-band wavelength often perform better in measuring coal mining induced subsidence [31]. Conversely, severe decorrelation may occurs in the interferograms generated by SAR acquisitions with a shorter wavelength (e.g., C- and X-band), because of the fast and large deformation in coal mining areas [25]. This significantly reduces the options for InSAR pairs in the monitoring of mining-related deformation, thereby increasing the possibility of temporal gaps occurring in the multi-temporal InSAR observations.

3. Methodology

3.1. Logistic Model

The dynamic subsidence of a ground surface point caused by underground coal extraction generally develops like an S-shaped curve, i.e., it starts slowly (initial phase), then increases rapidly (main phase), and finally levels off (residual phase) [8,27]. The logistic model is a typical S-shaped function and has been widely utilized to describe the dynamic subsidence of a single surface point caused by underground mining, tunneling, compression of a building’s foundation, etc. [12,23,32], and it is expressed as:

$$W(t)=\frac{W_0}{1+a{e}^{-bt}}$$

(1)

where $W(t)$ is the ground surface dynamic subsidence at time $t$ with respect to the reference at the starting time $t_0=0$ (i.e., $W(t_0=0)\equiv 0$). $W_0$ denotes the possible maximum subsidence. $a$ and $b$ represent the shape parameters of the logistic curve. $W_0$, $a$, and $b$ are the parameters of the logistic model.

3.2. Bridging the Temporal Gaps in InSAR Observations with the Logistic Model

Only three model parameters (i.e., $W_0$, $a$, and $b$) require estimating (see Equation (1)) if the logistic model is applied to describe or predict the mining-related dynamic subsidence at a single surface point. Therefore, it is possible to estimate these parameters if at least three effective observations of subsidence are obtained, even though some temporal gaps exist in the observations. This is the core idea behind the proposed method. To this end, we first construct a functional relationship between the InSAR observations and the model parameters of the logistic model. Subsequently, an algorithm based on the genetic algorithm (GA) and the simplex algorithm (SA) is used to estimate these model parameters from the InSAR observations with temporal gaps. Having obtained these model parameters, the dynamic subsidence can be estimated based on the logistic model.

3.2.1. InSAR Observations

For an differential InSAR pair generated by two SAR images acquired on $t_{\mathrm{A}}$ and $t_{\mathrm{B}}$, its unwrapped phase $\delta {\varphi }_{\mathrm{obs}}$ can be roughly expressed as follows [18,33]:

$$\delta {\varphi }_{\mathrm{obs}}\approx \delta {\varphi }_{\mathrm{def}}+\delta {\varphi }_{\mathrm{topo}}+\delta {\varphi }_{\mathrm{orbit}}+\delta {\varphi }_{\mathrm{atm}}+\delta {\varphi }_{\mathrm{noise}}$$

(2)

where $\delta {\varphi }_{\mathrm{def}}$, $\delta {\varphi }_{\mathrm{topo}}$, $\delta {\varphi }_{\mathrm{orbit}}$, $\delta {\varphi }_{\mathrm{atm}}$, and $\delta {\varphi }_{\mathrm{noise}}$ are the phase components associated with the ground deformation, the topographic errors of the external digital elevation model (DEM) data, the orbit inaccuracies, the differential atmospheric artifacts, and the noise, respectively.

In Equation (2), the phase component caused by the orbit inaccuracies $\delta {\varphi }_{\mathrm{orbit}}$ and the long-wavelength component of differential atmospheric phase $\delta {\varphi }_{\mathrm{atm}}$ can be mitigated significantly using a biquadratic polynomial model [34,35]. Furthermore, the phase noise $\delta {\varphi }_{\mathrm{noise}}$ also can be reduced effectively using phase filtering [36,37,38]. After the mitigation of these three phase components, the unwrapped phase $\delta {\varphi }_{\mathrm{obs}}$ can be rewritten as [21]:

$$\delta {\varphi }_{\mathrm{obs}}\approx \delta {\varphi }_{\mathrm{def}}+\delta {\varphi }_{\mathrm{topo}}+\delta {\varphi }_{\mathrm{residual}}\approx \frac{4\pi }{\lambda }\left[d_{LOS}\left(t_{\mathrm{B}}\right)-d_{LOS}\left(t_{\mathrm{A}}\right)\right]+\frac{4\pi }{\lambda }\frac{B_{\perp }\Delta h}{r\sin \theta }+\delta {\varphi }_{\mathrm{residual}}$$

(3)

where $d_{LOS}\left(t_{\mathrm{B}}\right)$ and $d_{LOS}\left(t_{\mathrm{A}}\right)$ are the time-series displacements along the radial line-of-sight (LOS) direction at times $t_{\mathrm{A}}$ and $t_{\mathrm{B}}$ with respect to the starting time $t_0$ (i.e., $d_{LOS}\left(t_0\right)\equiv 0$); $\lambda$, $r$, $\theta$, and $B_{\perp }$ are the radar wavelength, distance of sensor to target, radar incidence angle, and perpendicular baseline of the two SAR acquisitions, respectively; $\Delta h$ denotes the topographic error of the external DEM; $\delta {\varphi }_{\mathrm{residual}}$ is the residual phase caused by the orbital errors, the noise, and the atmospheric artifacts. Transforming the unwrapped phase $\delta {\varphi }_{\mathrm{obs}}$ into the LOS change $d_{\mathrm{obs}}$, we obtain:

$$d_{\mathrm{obs}}=\delta {\varphi }_{\mathrm{obs}}\times \frac{\lambda }{4\pi }=d_{LOS}\left(t_{\mathrm{B}}\right)-d_{LOS}\left(t_{\mathrm{A}}\right)+\frac{B_{\perp }\Delta h}{r\sin \theta }+\delta {\varphi }_{\mathrm{residual}}\cdot \frac{\lambda }{4\pi }$$

(4)

Generally, it is inevitable that the InSAR-derived LOS change $d_{\mathrm{obs}}$ contains the residual phase $\delta {\varphi }_{\mathrm{residual}}$ caused by orbit errors, noise, and atmospheric artifacts, although some mitigation strategies (e.g., filtering and fitting) were applied in this study. However, mining-related deformation is characterized by large magnitude and it usually dominates the InSAR-derived LOS change after performing these mitigation strategies. Hence, the influence of these residual phases on the following estimation of the logistic model should be limited. Thus, the last term (i.e., $\delta {\varphi }_{\mathrm{residual}}$) can be safely ignored in this study.

3.2.2. Functional Relationship between InSAR Observations and the Logistic Model

As can be seen from Equation (4), the InSAR observation mainly contains two components, i.e., differential LOS deformation $d_{LOS}\left(t_{\mathrm{B}}\right)-d_{LOS}\left(t_{\mathrm{A}}\right)$ between the two SAR acquisitions and the topographic error $B_{\perp }\Delta h/\left(r\sin \theta \right)$. Since the LOS deformation is the joint contribution of the actual three-dimensional (3-D) ground surface displacements [9,21], and the logistic model can only describe the mining-induced vertical subsidence, it is difficult to directly estimate the model parameters of the logistic model from the InSAR observations without some assumptions.

Due to the small incidence angles of the current SAR sensors, the vertical subsidence actually dominates the contribution of the 3-D displacements to LOS deformation [21,39]. Therefore, with the assumption that the mining-induced deformation in the vertical direction in near field (close to the centre of deformation basin) is larger than that in the horizontal direction [9,11], we transform the LOS deformation directly to the vertical subsidence, i.e., $W=d_{LOS}/\cos \theta$ [12]. In fact, this is a common strategy when the mining-induced 3-D displacements are difficult to derive from the available InSAR data [15]. Expressing the vertical subsidence with the logistic model, we have:

$$d_{LOS}\left(t\right)=\frac{W_0}{1+a{e}^{-bt}}\times \cos \theta$$

(5)

Substituting Equation (5) into Equation (4), we construct the functional relationship between the InSAR observation $d_{\mathrm{obs}}$ and the logistic model:

$$d_{\mathrm{obs}}==\left[\frac{W_0}{1+a{e}^{-b{t}_{\mathrm{B}}^{\prime }}}-\frac{W_0}{1+a{e}^{-b{t}_{\mathrm{A}}^{\prime }}}\right]\times \cos \theta +\frac{B_{\perp }}{r\sin \theta }\times \Delta h$$

(6)

where $t_A^{\prime }$ and $t_B^{\prime }$ are the cumulative times of $t_A$ and $t_B$ with respect to the initial time $t_0$, i.e., $t_A^{\prime }=t_A-t_0$ and $t_B^{\prime }=t_B-t_0$.

3.2.3. Model Parameter Estimation

We assume that n+1 SAR acquisitions acquired on the ordered times $(t_0,t_1,\cdots ,t_n)$ cover the mining area, and m available InSAR pairs with temporal gaps are generated due to the decorrelation factors (see Section 2.1). For an arbitrary pixel $\left(i,j\right)$ with coherence stability, m multi-temporal InSAR observations with temporal gaps can be derived from the m available InSAR pairs. In this case, m equations can be obtained based on the constructed functional relationship between the InSAR observations and the logistic model. Considering that only four unknown parameters, i.e., $W_0$, $a$, $b$, and $\Delta h$, require estimating at this pixel, it is possible to solve them if the number of InSAR observations m is over four.

Due to the highly nonlinear behavior of Equation (6), a nonlinear searching method based on the GA and the SA (hereafter referred to as GA + SA) is developed in this study for estimating the unknowns in Equation (6). Firstly, we obtain globally optimal solutions using the GA procedure based on the InSAR observations with temporal gaps. In fact, the GA-derived solutions contain large random errors due to the large number of random operations in the GA [40]. Consequently, we refine these solutions using the SA procedure with the initial starts of the GA-derived solutions, in view of the good local optimization and sensitivity to the initial starts of the SA.

3.2.4. Dynamic subsidence Estimation with the Logistic Model

Having obtained the model parameters of the logistic model at pixel $(i,j)$, its dynamic subsidence at the times of SAR acquisitions $(t_0,t_1,\cdots ,t_n)$ can be estimated with:

$$W(t_k,i,j)=\frac{W_0(i,j)}{1+a(i,j)e^{-b(i,j)\left(t_k-t_0\right)}}k=0,1,\cdots ,n.$$

(7)

Similarly, the dynamic subsidence at the other pixels with coherence stability can be estimated on a pixel-by-pixel basis. Figure 1 shows the block diagram of the proposed method.

4. Simulated Experiment

As stated in Section 3.2, the proposed method is performed on a pixel-by-pixel basis. For the sake of simplicity, we just focused on a single surface point in the mining area and tested the proposed method by the use of a simulated experiment.

4.1. Simulation of InSAR Observations

We take the in situ leveling measurements of dynamic subsidence caused by a longwall working panel of a coal mining area in China (see the black triangles in Figure 2) as an example. Fitting the in situ measurements of subsidence using the logistic model and the GA + SA, we obtained the following parameters for the logistic model: $W_0=-1.10m$, $a=3998$, and $b=0.041$. As can be seen in Figure 2, the dynamic subsidence fitted by the logistic model (see the green line) is in good agreement with the in situ measurements, with an STD of 0.014 m. Afterwards, we selected this surface point as the source of the simulated experiment in this study, and we simulated the InSAR observations obtained by five common satellites with SAR sensors, i.e., ALOS1, Envisat, COSMO-SkyMed, ALOS2, and TerraSAR-X (their revisit cycles are listed in

Table 1. Revisit cycles of the five common satellites with SAR sensors, and the number of their simulated InSAR observations.

Satellites	Revisit Cycle (Days)	Number of InSAR Observations
ALOS1	46	10
EnviSAT	35	14
COSMO-SkyMed	16	31
ALOS2	14	35
TerraSAR	11	45

), based on the logistic model and its estimated model parameters at this surface point.

We use the simulation of the ALOS1 observations as an example in order to demonstrate the simulation procedure. We first calculated the dynamic subsidence at this point using the logistic model and its estimated model parameters, and the results are shown in Figure 2 (green line). Assuming that the time of the first ALOS SAR image is $t_0=0$, then the times of the other 10 ALOS1 SAR images in the following 500 days (i.e., $t_0$~$t_{10}$, see the blue diamonds in Figure 2) can be obtained according to its revisit cycle (i.e., 46 days). Having obtained the times of the SAR acquisitions, the dynamic subsidence at these times can be calculated by the logistic model and its model parameters (see Equation (7)). For the sake of simplicity, we just generated 10 InSAR pairs using two time-adjacent SAR acquisitions. Consequently, the differences in the dynamic subsidence of two time-adjacent SAR acquisitions give 10 error-free InSAR observations. Finally, we added Gaussian noise with a mean value of zero and an STD of 2 cm to the 10 InSAR observations, to simulate the inaccuracies of the InSAR observations. At this stage, 10 simulated InSAR observations with errors were obtained.

4.2. Simulation of Temporal Gaps

Theoretically, to analyze the influence of the temporal gaps on the dynamic subsidence derived by the proposed method, we should simulate different numbers of temporal gaps at different locations. For instance, 45 TerraSAR InSAR observations are available (see Table 1), whereas only three model parameters of the logistic model require estimation. Therefore, theoretically, no more than 42 temporal gaps should be allowed to exist in any one of the TerraSAR InSAR observations. In this case, about $3.52\times {10}^{13}$ possible combinations (calculated by $C_{45}^1+C_{45}^2+\cdots +C_{45}^{42}$, where $C_m^k$ means choosing $k$ combinatorial subsets from $m$ numbers) need to be tested. However, such a strategy is very time-consuming with the GA + SA procedure.

In fact, there are five feature epochs at the logistic curve [41,42], i.e., the subsidence start time T1, the first extreme value of acceleration T2, the maximum velocity T3, the second extreme value of acceleration T4, and the stability start time T5 (e.g., the red circles in Figure 2). The observations at these five feature times (termed feature observations) play an essential role in the accurate estimation of the model parameters of the logistic model. Generally, the model parameters can be accurately estimated if the subsidence observations cover all five feature times [42]. Hence, it should be sufficient to simulate five groups of InSAR observations with temporal gaps that miss one to all of the five feature observations, in order to analyze the influence of the temporal gaps on the dynamic subsidence derived by the proposed method.

We also take the ALOS1 data as an example to show the simulation procedure of the five groups of InSAR observations with temporal gaps. Firstly, we simulated the first group of ALOS1 InSAR observations with a single temporal gap that misses one of the five feature times (i.e., T1−T5, see the red circles in Figure 2). For example, the InSAR observation generated by the SAR images $t_3$ and $t_4$ will be eliminated if the temporal gap does not cover the feature time T2 (see Figure 2). Likewise, the other four groups of ALOS1 InSAR observations with 2−5 temporal gaps that miss 2−5 feature observations were simulated. The numbers of InSAR observations with different temporal gaps for the five common satellites are shown in Figure 3a.

4.3. Influence of the Temporal Gaps on the Derived Dynamic Subsidence

Having obtained the InSAR observations with different numbers of temporal gaps (i.e., 1−5), we can derive the dynamic subsidence from these observations using the proposed method and then analyze the influence of the temporal gaps on the derived dynamic subsidence. We again take the ALOS1 InSAR observations as an example to demonstrate this procedure. Firstly, the proposed method was applied to derive the dynamic subsidence based on the first group of InSAR observations with a single temporal gap, and we then calculated the mean STD of the differences between the simulated dynamic subsidence and those derived by the proposed method. Subsequently, the same procedure was carried out on the other four groups of InSAR observations, yielding four corresponding mean STDs.

Figure 3b shows the mean STDs of the dynamic subsidence derived from the InSAR observations of the five common satellites. As can be seen from this figure, the mean STDs of the derived dynamic subsidence from the ALOS1 (see the red bars in Figure 3b) and Envisat InSAR observations (see the blue bars in Figure 3b) are about 2.1 cm and 1.2 cm, respectively, in the case of one or two temporal gaps. Such accuracies amount to 1.9% and 1.1% of the maximum subsidence at this point (110 cm). However, the mean STDs reach as high as 41 cm and 34 cm, respectively, when the number of temporal gaps increases from three to five, which amounts to 37.2% and 30.9% of the maximum subsidence. In other words, the proposed method nearly fails when the ALOS1 and Envisat InSAR observations miss all five feature observations. This is mainly because fewer effective InSAR observations can be used to constrain the logistic model [42]. We again take the ALOS 1 observations as an example. Only 10 ALOS1 InSAR observations are available due to its long revisit cycle (46 days). Once InSAR observations at the five feature times (see the red circles in Figure 2) are unavailable, the main subsidence phase (roughly from t₃ to t₇) cannot be constrained by the InSAR observations. In this case, the estimated model parameters may contain large errors.

However, due to the short revisit cycles of the COSMO-SkyMed, ALOS2, and TerraSAR satellites, more InSAR observations are available than for ALOS1 and Envisat (see Figure 3a). Therefore, the mean STD of the dynamic subsidence estimated from these three sets of satellite observations is only about 0.028 m when the number of temporal gaps reaches five (see Figure 3b). This result suggests that more InSAR observations can help reconstruct the dynamic subsidence of coal mining areas by the proposed method.

5. Real Data Experiment

5.1. Study Area and SAR Acquisitions

Datong coal mining area (marked by the blue dashed line in Figure 4), which is located to the northwest of the city of Datong (the black circle in Figure 4), China, was selected to test the proposed method. The land surface of this mining area is mainly covered by farming land and deciduous trees (e.g., poplar and Sophora japonica Linn). In this case, the dramatic seasonal variation of the land cover can give rise to severe decorrelation in interferograms. Furthermore, the surface displacements in this coal mining area are often characterized by large deformation gradients and high speeds (e.g., 0.6 m in 46 days at the peak), due to its shallow mining depth (about 230−300 m) and large mining height (about 3–8 m). Such large deformation gradients also trigger severe decorrelation in interferograms.

Given the fact that L-band wavelength SAR data perform much better than C-band and X-band data in monitoring mining-related deformation [31], 12 ascending L-band wavelength ALOS PALSAR images spanning from 1 July 2007 to 28 February 2009 were selected to monitor the mining-induced deformation. The footprints of the SAR acquisitions and their acquisition times are shown in Figure 4a (red rectangle) and Figure 5 (blue circles), respectively. From 01 July 2007 to 28 February 2009, three working panels (entitled WP1, WP2, and WP3, respectively) were operated at a mean mining depth of around 250 m, with a mean advancing velocity of about 3 m/day. The geographic locations of these three working panels are shown in Figure 4b.

5.2. SAR Acquisitions and Data Processing

To select suitable thresholds for the perpendicular and temporal baselines of the InSAR pairs, we generated 30 differential interferograms whose perpendicular baselines and temporal separations were less than 3000 m and 250 days. By the visual analysis of the generated differential interferograms, we find that severe decorrelation occurs if the perpendicular and temporal baselines of the InSAR pairs are more than 1800 m and 138 days. Therefore, we set 1800 m and 138 days as the thresholds of the perpendicular and temporal baselines. After eliminating the interferograms with severe decorrelation, we obtained 17 available InSAR pairs, and their perpendicular and temporal baselines are plotted in Figure 5.

Figure 5 demonstrates that the available InSAR pairs cannot cover the entire time period of SAR acquisitions, and a temporal gap (marked by the light yellow box) occurs from 18 May to 3 July 2008. In addition to the large deformation gradient and the dramatic variation of the land cover in this mining area, the large perpendicular baselines between the SAR image acquired on 18 May 2008 and the later ones (from −3232 to −6025 m) may be the main factor causing this temporal gap. For the InSAR pairs with a temporal gap, the traditional TS-InSAR techniques cannot accurately estimate the dynamic subsidence in the entire period of SAR acquisitions. However, it is possible to derive the dynamic subsidence using the proposed method.

To generate InSAR observations from the 17 InSAR pairs, we first coregistered the other 11 SAR acquisitions with the one acquired on 16 February 2008. Subsequently, the two-pass D-InSAR procedure [43] was utilized to process the 17 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. Furthermore, a quadratic polynomial model is used to mitigate the phase ramps associated with possible orbit inaccuracies and the long wavelength part of atmospheric artifact screening. Additionally, the least squares-based filter [38] is utilized to suppress the noise in interferograms.

5.3. Deriving Dynamic Subsidence from InSAR Observations

5.3.1. Determining the Model Parameters of the Logistic Model and the Topographic Errors

We first identified the pixels with coherence stability using the coherence threshold of 0.3, i.e., the available InSAR pairs with a coherence of over 0.3. The pixels whose deformations were insignificant (less than 1 cm) over the entire time period of SAR acquisitions were then masked, in order to reduce the time consumption of the parameter estimation. For an arbitrary pixel with coherence stability, 17 multi-temporal InSAR observations were available and only four parameters needed to be determined (i.e., $W_0$, $a$, $b$, and $\Delta h$). Hence it was feasible to solve these parameters based on the developed GA + SA from the 17 InSAR observations, although they contained a temporal gap.

More specifically, setting the population size, crossover fraction, and maximum generations of the GA as 100, 0.65, and 700, respectively, we then estimated the solutions of the four parameters at each pixel with coherence stability using the GA procedure. Afterwards, the SA was applied to refine the solutions of the four parameters by using the GA-derived solutions as initial starts, and the results are shown in Figure 6.

As can be seen in Figure 6a, the possible maximum subsidence $W_0$ is about 1.26 m in this mining area from 01 July 2007 to 28 February 2009. Additionally, the estimated topographic errors vary from −16.3 to 14.7 m, which are in good agreement with the absolute vertical height error of the SRTM DEM data (about 16 m) [44]. This indicates that the estimated topographic errors in this study are reliable.

5.3.2. Estimating Dynamic subsidence with the Logistic Model

Having determined the model parameters of the logistic model at those pixels with coherence stability, the dynamic subsidence induced by underground extraction at the times of the SAR acquisitions over these pixels can be estimated using Equation (7), and the results are shown in Figure 7. As can be seen from Figure 7a,b, the mining activities mainly operated in working panels WP1 and WP3 (see the top and bottom rectangles in Figure 7) between 1 July and 1 October 2007. In this period, the maximum ground surface subsidence associated with the mining of working panel WP1 was about 0.15 m, due to its small advancing distance of about 110 m according to in-site measurements. On the other hand, the maximum ground surface subsidence induced by the extraction of working panel WP3 reached 1.06 m, because of its larger advancing distance of about 600 m according to in-site measurements.

Extraction in working panel WP2 started at the end of September 2007. It advanced only about 80 m up until 1 January 2008, and the resulting subsidence was insignificant (about 0.05 m, see Figure 7c,d). This is expected because the advancing distance of WP2 in this time period amounted to just 0.32 times the mining depth H = 250 m. According to [11], underground mining has an insignificant influence on the ground surface when the advancing distance is less than a threshold (generally H/6–H/3). However, the subsidence basin expanded gradually and the maximum deformation increased dramatically to 1.26 m (see Figure 7e−h) with the extraction of WP2 in the following five months (i.e., by 3 July 2008). Thereafter, the subsidence basin continued to expand but the maximum subsidence remained nearly stable with the extraction of WP2 by 3 January 2009. This is because the advancing distance of WP2 (about 3.5H) reached the supercritical size (generally over 0.9H to 2.3H) after 3 July 2008. These features are in good agreement with the theoretical temporal-spatial development of ground subsidence caused by longwall mining [11,27].

It is noted that an area emerged (marked by the cyan rectangle in Figure 7h) in the subsidence basin where the maximum subsidence (about 0.96 m) was obviously less than the surrounding area (about 1.25 m). Usually, surface subsidence in the central portion of a subsidence basin caused by supercritical longwall mining is approximately equal when there is no significant variation in the geomining condition [11,27]. Since the geomining condition in the concerned mining area is smooth, it is believed that the subsidence in this area marked by the cyan rectangle in Figure 7h was underestimated. Interestingly, this area with underestimated maximum subsidence covers the mined-out area in the time period of the temporal gap (i.e., between 18 May and 3 July 2008, see Figure 7g,h), and some feature times of the logistic model there cannot be covered by the InSAR observations. According to the analysis in Section 4.3, missing feature observations greatly affect the accuracy of the estimated dynamic subsidence. Therefore, this may be the main reason for the underestimation. How to address this issue will be our future research topic.

5.4. Accuracy Validation of the Dynamic Subsidence

Four ascending ALOS PALSAR images from a different orbit (Frame: 790, Path: 453) acquired on 15 December 2007, 30 January, 16 September, and 17 December 2008 were used to validate the accuracy of the estimated dynamic subsidence. Firstly, we formed two InSAR pairs from the first two and the last two SAR acquisitions, respectively, and we then generated the two LOS deformation maps by processing these two InSAR pairs with the two-pass D-InSAR procedure (see Section 5.2). Subsequently, the vertical subsidence during the time periods of the two InSAR pairs was obtained by means of ignoring the horizontal components of the LOS deformation (see Section 3.2.2), and the results are shown in Figure 8a,b. Meanwhile, the synchronous vertical subsidence corresponding to the two InSAR pairs at those pixels with coherence stability was calculated using the logistic model and its determined parameters (i.e., Equation (7)), and the results are shown in Figure 8c,d. Figure 8e,f represents the differences between the InSAR-derived result and the logistic model derived subsidence in the two time periods, respectively. As can be seen from Figure 8, the logistic model interpolated subsidence in the two time periods is in good agreement with the InSAR-derived subsidence.

We first masked those pixels whose maximum subsidence was less than 0.02 m (see the area outside the blue dashed circles in Figure 8c,d) before quantitatively validating the accuracy of the derived dynamic subsidence. This is because these pixels were insignificantly affected by the underground extraction in this period, and the logistic model may be unsuitable to model the temporal behavior of the dynamic subsidence at these pixels. Subsequently, we calculated the difference between the logistic model derived subsidence and the InSAR-measured subsidence at the remaining pixels (the area inside the blue dashed circles in Figure 8c,d). The STDs of the differences from 15 December 2007 to 30 January 2008 and from 16 September to 17 December 2008 are about 0.020 and 0.014 m, and the mean STD in the two time periods is about 0.017 m. These results indicate that the dynamic subsidence estimated by the proposed method is reliable.

However, it should be pointed out that maximum discrepancies of −0.132 and 0.051 m occur in the two periods, respectively. Interestingly, these points with large discrepancies are almost all situated at the pixels affected by the extraction of the multiple working panels WP1, WP2, and WP3 (see the red ellipses in Figure 8e,f). Thus, it is believed that the extraction of the multiple working panels may be responsible for these large discrepancies, and we clarify this issue in Section 6.2.

6. Discussions

6.1. Comparison Between the Logistic, Linear, and Cubic Polynomial Models

There are two popular models in the traditional TS-InSAR techniques, i.e., linear (i.e., mean velocity) and cubic polynomial models [33,45]. When there is no prior knowledge of the temporal ground surface deformation, these two models are usually used to estimate time-series deformation with temporally linear behavior and with low nonlinearity, respectively. Theoretically, these two models could also be applied to derive dynamic subsidence from InSAR observations with temporal gaps. Admittedly, it is expected that large errors would arise in the dynamic subsidence estimated by these two models, because the mining-related dynamic subsidence is usually highly nonlinear. However, how large the errors may occur is unknown when these two models are utilized to derive dynamic subsidence from InSAR observations with temporal gaps. Hence, it is necessary to carry out the comparisons of the dynamic subsidence estimated by the logistic, linear, and cubic polynomial models.

For this purpose, we first estimated the dynamic subsidence using the linear and cubic polynomial models with a least squares method (further details can be found in [33,46]) in this coal mining area from the InSAR observations with temporal gaps (see Section 5.2). Subsequently, we calculated the fitting residuals of the linear model and the cubic polynomial model, and the results are shown in Figure 9a,b. Meanwhile, the fitting residuals of the logistic model are also shown in Figure 9c for the sake of comparison.

As can be seen from Figure 9, the fitting residuals of the linear model (see Figure 9a) are primarily in the range of 0.01 to 0.24 m, with a mean value of 0.097 m. Additionally, the fitting residuals of the cubic polynomial model (see Figure 9b) mainly vary from 0.01 to 0.22 m, with a mean of 0.084 m. Such large residuals indicate that the linear and cubic polynomial models cannot accurately describe the mining-related dynamic subsidence, although the latter performs a little bit better than the former. However, the fitting residuals of the logistic model (see Figure 9c) mainly vary between 0.002 and 0.04 m, with a mean value of 0.021 m. The mean fitting residual of the logistic model reduces by about 78.4% and 75.0% when compared to the mean fitting residuals of the linear (i.e., 0.097 m) and cubic polynomial models (i.e., 0.084 m), respectively. This suggests that the logistic model performs much better than the linear and cubic polynomial models in describing mining-induced dynamic subsidence.

Afterwards, we evaluated the accuracies of the mining-related dynamic subsidence estimated by the linear and cubic polynomial models based on the four ALOS PALSAR images mentioned in Section 5.4. The results show that the mean STDs of the linear and cubic polynomial models are about 0.034 and 0.031 m. Such accuracies amount to 18.4% and 16.8% of the mean maximum subsidence in the two periods (about 0.185 m, see Figure 8b,e). In addition, the accuracies of the logistic model derived dynamic subsidence (with a mean STD of about 0.017 m) are roughly 50.0% and 45.2% greater than those of the linear and cubic polynomial models.

To further intuitively demonstrate the performances of these three models in describing mining-induced dynamic subsidence, we fitted the in situ leveling measurements of the ground surface point mentioned in Section 4.1 using the linear and the cubic polynomial models. The results are shown in Figure 10. Meanwhile, the fitted results of the logistic model are also added into Figure 10 for the sake of comparison.

As can be seen from Figure 10, large discrepancies occur between the in situ measurements and the ones fitted by the linear (green diamonds) and the cubic polynomial models (asterisks), with the STDs of 0.21 and 0.104 m. On the other hand, Figure 10 also shows that the fitted results of the logistic model agree very well with the in situ measurements, with the STD of about 0.014 m. This again indicates that the logistic model is capable of describing mining-related dynamic subsidence at a single surface point.

6.2. Influence of the Multiple Working Panel Extractions

As mentioned in Section 5.4, the multiple working panel extractions of WP1, WP2, and WP3 may have given rise to large discrepancies between the InSAR observations of ground surface subsidence and the logistic model estimated synchronous results. To clarify this issue, we selected four ground surface points from the areas affected by the extractions of WP1, WP2, and WP3 (denoted as P1 to P4 and marked by the red points in Figure 8e,f) as examples. Figure 11 shows the comparison between the InSAR observations of the differential subsidence at the four selected points (see the blue circles) and the logistic model estimated synchronous results (see the red diamonds).

The surface points P1 and P2 were first affected by the extraction of WP1 from 1 July to 1 October 2007. Afterwards, they were influenced by the extraction of WP2 (see Figure 4). In this case, points P1 and P2 actually witnessed two main subsidence phases [8]. This can be validated from Figure 11a,b. From 1 July to 1 October 2007 (i.e., the first three InSAR observations in Figure 11a,b), points P1 and P2 subsided rapidly (i.e., the main phase of subsidence) before the extraction of working panel WP1 ceased. Afterwards, these two points then experienced the initial phase, the main phase, and the residual phase again, due to the extraction of working panel WP2 after 1 October 2007. Similarly, points P3 and P4 also experienced two main subsidence phases caused by the extraction of WP3 and WP2 from 1 July 2007 to 1 January 2008 (see the first four InSAR observations in Figure 11c,d) and from 16 February 2008 to 3 January 2009 (see the 11th to 16th InSAR observations in Figure 11c,d), respectively.

The logistic model can only describe the S-shaped behavior of dynamic subsidence, i.e., only one main subsidence phase. However, two or more main phases occur in the dynamic subsidence of a ground surface point affected by multi working panel extractions. If the logistic model is used to describe dynamic subsidence with two or more main phases, only one can be fitted well, and large discrepancies will occur in other main phases. This can also be seen in Figure 11. Therefore, the dynamic subsidence caused by the extraction of the multiple working panels cannot be accurately estimated by the proposed method. This prior knowledge of the dynamic subsidence caused by multi working panel extractions, in return, could guide us to choose a suitable model for describing it in our future study.

7. Conclusions

The logistic model is aimed at bridging the temporal gaps in multi-temporal InSAR observations, thereby deriving mining-induced dynamic subsidence from InSAR observations with temporal gaps. A simulated experiment was first carried out in this study. The results show that the proposed method can obtain accurate mining-related dynamic subsidence in the case of a few missing feature observations (e.g., two or three). Furthermore, more InSAR observations can mitigate the influence of the temporal gaps on the derived dynamic subsidence. Subsequently, a real data experiment in the Datong coal mining area, China, was performed. The results show that the maximum subsidence reached 1.26 m from 1 July 2007 to 28 February 2009, and the accuracy of the dynamic subsidence estimated by the proposed method is about 0.017 m. Compared with the linear and cubic polynomial models of the traditional TS-InSAR techniques, the accuracy of the logistic model derived dynamic subsidence is 50.0% and 45.2% greater than those of the linear and cubic polynomial models, respectively.

However, the temporal gaps in the InSAR observations greatly affect the accuracies of the estimated dynamic subsidence. This issue will be our future research topic. In addition, the operation of multiple working panel extractions gives rise to large discrepancies in the dynamic subsidence estimated by the proposed method. This is mainly because the logistic model is only capable of describing dynamic subsidence with one main subsidence phase. How to overcome this limitation will also be addressed in our future study.

We are aware that the proposed method was only tested in a coal mining area with a single temporal gap, although a simulated experiment was also carried out in this study. Hence, it would certainly be worthwhile to test the proposed method in mining areas where the InSAR observations contain more temporal gaps, in order to further assess the feasibility of the proposed method.