Sequential InSAR Time Series Deformation Monitoring of Land Subsidence and Rebound in **’an, China

Abstract

Interferometric synthetic aperture radar (InSAR) time series deformation monitoring plays an important role in revealing historical displacement of the Earth’s surface. **’an, China, has suffered from severe land subsidence along with ground fissure development since the 1960s, which has threatened and will continue to threaten the stability of urban artificial constructions. In addition, some local areas in **’an suffered from uplifting for some specific period. Time series deformation derived from multi-temporal InSAR techniques makes it possible to obtain the temporal evolution of land subsidence and rebound in **’an. In this paper, we used the sequential InSAR time series estimation method to map the ground subsidence and rebound in **’an with Sentinel-1A data during 2015 to 2019, allowing estimation of surface deformation dynamically and quickly. From 20 June 2015 to 17 July 2019, two areas subsided continuously (Sanyaocun-Fengqiyuan and Qujiang New District), while **’an City Wall area uplifted with a maximum deformation rate of 12 mm/year. Furthermore, Yuhuazhai subsided from 20 June 2015 to 14 October 2018, and rebound occurred from 14 October 2018 to 17 July 2019, which can be explained as the response to artificial water injection. In the process of artificial water injection, the rebound pattern can be further divided into immediate elastic recovery deformation and time-dependent visco-elastic recovery deformation.

1. Introduction

The interferometric synthetic aperture radar (InSAR) is a remote sensing technique, which has been commonly used in the investigation of large-scale ground deformation. Land subsidence in urban areas has been investigated by the InSAR technique in Las Vegas, USA [1], Houston–Galveston, USA [2], Mexico City, Mexico [3], northeast Iran [4], West Thessaly Basin, Greece [5], Pisa urban area, Italy [6], Rome metropolitan area, Italy [7], Bei**g [8], Tian** [9], Taiyuan [10] and Datong, China [11].

** of groundwater was issued, and the deformation rate began to decrease [17]. When an aquifer water level rises during artificial water injection, the rebound can be divided into short-term elastic recovery and time-dependent visco-elastic recovery [18]. Previous study revealed the uplift phenomena at Jixiangcun (point D in Figure 4) in ** and sequential estimation of deformation parameters.

2.1. Selection of Coherent Pixels

To mitigate the effects of decorrelation and retrieve large-gradient deformation, the small baseline subset (SBAS) InSAR method was proposed based on the interferograms with short spatial and temporal baselines [20]. In this paper, we use the temporal coherence to select coherent pixels [21,22], which is defined in Equation (1) for one generic pixel x:

$${\gamma }_x=\frac{1}{N}\left|{\displaystyle \sum _{i=1}^N\exp \{\sqrt{-1}({\psi }_{x,i}-{\tilde{\psi }}_{x,i}-\Delta {\stackrel{⌢}{\varphi }}_{\theta }^u{}_{,x,i})\}}\right|$$

(1)

where N is the number of interferograms, $\psi$ represents the flattened and topographically corrected interferogram, $\tilde{\psi }$ represents the spatially correlated phase component, and ${\stackrel{⌢}{\varphi }}_{\theta }^u$ represents the spatially uncorrelated phase component (look-angle error phase). A more detailed introduction is provided in [21,22].

Followed by the phase unwrap**, the 3D phase unwrap** method was employed to mitigate the closed-loop discontinuities error in two-dimensional (2D) phase unwrap** [23]. It was used to explore the spatial and temporal relationships within the multi-interferograms, i.e., involving the computation of two Delaunay triangulations, which are usually referred to as “temporal” and “spatial” triangulations, respectively [24,25].

2.2. Sequential InSAR Time Series Estimation

After the atmospheric phase, orbital and digital elevation model (DEM) errors were removed from the interferograms, and we estimated the time series deformation phases by using the following function model:

$$\underset{{\mathbf{A}}_{}}{\underbrace{\left[\begin{array}{cccc}-1& 1& 0& \\ -1& 0& 1& \\ ⋮& ⋮& ⋮& \\ 0& 0& 0& \end{array}\begin{array}{cccc}\dots & 0& 0& 0\\ \cdots & 0& 0& 0\\ ⋮& ⋮& ⋮& ⋮\\ \cdots & 0& -1& 1\end{array}\right]}}\underset{{\mathit{X}}_{}}{\underbrace{\left[\begin{array}{c}{\phi }_1\\ {\phi }_2\\ ⋮\\ {\phi }_N\end{array}\right]}}=\underset{\mathit{L}}{\underbrace{\left[\begin{array}{c}un{w}_1\\ un{w}_2\\ ⋮\\ un{w}_M\end{array}\right]}}$$

(2)

where ${\phi }_i(i=1,\cdots ,N)$ denotes the deformation phase at the different synthetic aperture radar (SAR) acquisition date and $N$ is the number of SAR images. Note the deformation at the first SAR acquisition date is set to zero, i.e., ${\phi }_0=0$. To estimate the deformation time series, the archived SAR data are modeled as:

$$\begin{array}{l}{\mathit{V}}_1={\mathbf{A}}_1{\mathit{X}}^{(1)}-{\mathit{L}}_1\begin{array}{cc},& {\mathbf{P}}_1\end{array}\\ {\mathit{X}}^{(1)}={({\mathbf{A}}_1{}^{\mathrm{T}}{\mathbf{P}}_1{\mathbf{A}}_1)}^{-1}{\mathbf{A}}_1{}^{\mathrm{T}}{\mathbf{P}}_1{\mathit{L}}_1\\ {\mathbf{Q}}_{X^{(1)}}={({\mathbf{A}}_1{}^{\mathrm{T}}{\mathbf{P}}_1{\mathbf{A}}_1)}^{-1}\end{array}$$

(3)

where ${\mathit{L}}_1$ is archived SAR data with design matrix ${\mathbf{A}}_1$ and weight matrix ${\mathbf{P}}_1$. ${\mathit{X}}^{(1)}$ indicates the first estimation of parameter $\mathit{X}$, and ${\mathbf{Q}}_{X^{(1)}}$ is its cofactor matrix. The superscript T stands for the transposition of a matrix.

When we acquire a new SAR image, unlike conventional SBAS InSAR, to estimate deformation time series for all SAR images again we use the sequential estimation to update dynamically the deformation time series by only considering the unwrapped interferograms related to the new SAR image. Assuming the new measurements ${\mathit{L}}_2$ are the unwrapped interferograms related to the $(N+2)-\mathrm{th}$ new SAR acquisition, the weight matrix is ${\mathbf{P}}_2$, the design matrixes are ${\mathbf{A}}_2$ and $\mathit{B}$, and parameters are $\mathit{X}$ and $\mathit{Y}$, we can write its observational equation as follows:

$${\mathit{V}}_2=\left[\begin{array}{cc}{\mathbf{A}}_2& \mathit{B}\end{array}\right]\left[\begin{array}{c}{\mathit{X}}^{(2)}\\ \mathit{Y}\end{array}\right]-{\mathit{L}}_2\begin{array}{cc},& {\mathbf{P}}_2\end{array}$$

(4)

According to the principle of least squares (LS) Bayesian estimation [26], it holds that:

$${\mathit{V}}_2^{\mathrm{T}}{\mathbf{P}}_2{\mathit{V}}_2+{({\mathit{X}}^{(2)}-{\mathit{X}}^{(1)})}^{\mathrm{T}}{\mathbf{Q}}_{X^{(1)}}^{-1}({\mathit{X}}^{(2)}-{\mathit{X}}^{(1)})=\min$$

(5)

Then, we can deduce Equation (6) through Equations (3), (4) and (5) [26,27] as follows:

$$\begin{array}{l}\left[\begin{array}{c}{\mathit{X}}^{(2)}\\ {\mathit{Y}}_{}\end{array}\right]=\left[\begin{array}{c}{\mathit{X}}^{(1)}+{\mathbf{J}}_x({\overline{\mathit{L}}}_2-{\mathbf{B}}_{}{\mathit{Y}}_{})\\ {({\mathbf{B}}^{\mathrm{T}}{\mathbf{Q}}_J{}^{-1}\mathbf{B})}^{-1}{\mathbf{B}}^{\mathrm{T}}{\mathbf{Q}}_J{}^{-1}({\mathit{L}}_2-{\mathbf{A}}_2{\mathit{X}}^{(1)})\end{array}\right]\\ {\mathbf{Q}}_{[X^{(2)};Y]}=\left[\begin{array}{cc}{\mathbf{Q}}_{X^{(2)}}& {\mathbf{Q}}_{X^{(2)},Y}\\ {\mathbf{Q}}^{\mathrm{T}}{}_{X^{(2)},Y}& {\mathbf{Q}}_Y\end{array}\right]\\ {\mathbf{Q}}_{X^{(2)}}={\mathbf{Q}}_{X^{(1)}}-{\mathbf{J}}_x{\mathbf{A}}_2{\mathbf{Q}}_{X^{(1)}}+{\mathbf{J}}_x\mathbf{B}{\mathbf{Q}}_Y{\mathbf{B}}^{\mathrm{T}}{\mathbf{J}}_x^{\mathrm{T}}\\ {\mathbf{Q}}_{X^{(2)},Y}=-{\mathbf{J}}_x{\mathbf{B}}_{}{\mathbf{Q}}_{Y,Y}\\ {\mathbf{Q}}_Y={({\mathbf{B}}^{\mathrm{T}}{\mathbf{Q}}_J{}^{-1}\mathbf{B})}^{-1}\\ {\mathbf{J}}_x={\mathbf{Q}}_{X^{(1)}}{\mathbf{A}}_2{}^{\mathrm{T}}{\mathbf{Q}}_J{}^{-1}\\ {\mathbf{Q}}_J={\mathbf{P}}_2^{-1}+{\mathbf{A}}_2{\mathbf{Q}}_{X^{(1)}}{\mathbf{A}}_2{}^{\mathrm{T}}\end{array}$$

(6)

where $[X^{(2)};Y]$ is the updated deformation time series, in which Y is the cumulative deformation at the new SAR acquisition date, and ${\mathbf{Q}}_{[X^{(2)};Y]}$ represents their cofactor matrixes. ${\mathbf{J}}_x$ is the gain matrix, in which ${\mathbf{Q}}_J$ is the updated cofactor matrix with the new SAR image. Therefore, we can update the deformation parameters as quickly as possible, once a new SAR image is presented. For a more detailed discussion of sequential estimation of SBAS-InSAR dynamic deformation parameter methods, readers can refer to [28].

3. Study Area and Data

3.1. Study Area

** wells, a pressure difference between the aquifer pressure in pum** wells and the surrounding aquifer occurred, which drove the water in the surrounding soil to move toward pum** wells [31]. Under this situation, land subsidence accelerated and a subsided funnel was formed in the aquifer system. Over-exploitation of groundwater leads to a decrease of groundwater level. Although the confined aquifer is elastic, the continuous over-exploitation of groundwater leads to irrecoverable confined aquifer deformation, which further leads to land subsidence [18].

In order to study the spatiotemporal characteristics of subsidence and ground fissures, both global positioning system (GPS) and InSAR observations between 2005 and 2007 were employed [32,33]. Then, Envisat, advanced land observation satellite (ALOS) and TerraSAR SAR datasets were also used to investigate the two-dimensional deformation in **, DEM error correction, and atmospheric and orbital error correction, followed by the inversion of deformation parameters and their cofactor matrixes.

To update the time series deformation on the new SAR acquisition date, there are usually two ways to generate interferograms among newly received SAR images and the archived SAR images: single-link configuration to unwrap interferograms in the spatial domain, as shown in Figure 3A, and network-link configuration, shown in Figure 3B, where interferograms can be unwrapped in both spatial and time domains, i.e., 3D phase unwrap** [23,24,25]. We used the latter method to update the deformation time series.

The selected coherent pixels in the first group of SAR data were used to extract the phase for the interferograms and connect to the newly received SAR images. After 3D phase unwrap**, followed by the correction of DEM, atmospheric and orbital errors, the deformation rate and time series were updated by sequential estimation.

4. Results

4.1. Deformation Rate Map

Figure 4 presents the annual deformation rate map in the vertical direction over the main ** well to avoid pollution of groundwater. The result shows that Yuhuazhai continuously subsided from 20 June 2015 to 14 October 2018, then rebounded from 14 October 2018 to 17 July 2019, owing to the artificial water injection operation.

Owing to the high water quality of the deep confined aquifer, it was used as the main source of residential water supply. However, long-term over-exploitation of groundwater led to the decline of the groundwater level. The cumulative decline of the water level was nearly 70 m in Yuhuazhai areas [29]. Approximately 85% of the pum** wells were located to the south of the ground fissure F4 [30], where serious land subsidence occurred; the existence of ground fissure F4 hindered the flow of groundwater and limited the expansion of the surface deformation. For the phreatic aquifer head, the south of ground fissure F4 was 13.0–13.52 m in depth, while the north of ground fissure F4 was 14.15–16.48 m in depth [30]. For the confined aquifer head, the south of the ground fissure F4 was 103–106.4 m in depth, while the north of the ground fissure F4 was 49.57 m in depth [30]. Therefore, during the water injection, the groundwater to the south of the ground fissure recovered rapidly. After long-term land subsidence of Yuhuazhai (i.e., 1.8–2 m land subsidence from 1996 to 2016), a rapid rebound with a magnitude of about 130 mm occurred from 14 October 2018 to 1 December 2018, as a recoverable elastic deformation. Then, the time-dependent rebound of visco-elastic deformation continued from 1 December 2018 to 17 July 2019, which will likely continue for some time.

Moreover, there are 45 deep wells in the whole Yuhuazhai area, i.e., approximately 1 well per 1.7 km² [30]. We show the optical image of the largest rebound center in Figure 12, where we visually identify seven pum** wells. Figure 13 shows three images and photos corresponding to pum** wells 1, 2, and 3 of Figure 12.

6. Conclusions

In the process of urbanization, over-exploitation of groundwater in the ** well of Yuhuazhai around October 2018. The rebound pattern comprises two stages: the elastic deformation and visco-elastic deformation. The former can recover immediately, while the latter is time-dependent. The complex surface deformation in **’an reflects the changes in the aquifer system. Therefore, the control of groundwater balance can alleviate surface deformation.