# Demonstration of Time-Series InSAR Processing in Beijing Using a Small Stack of Gaofen-3 Differential Interferograms

Authors: Jili Wang; Weidong Yu; Yunkai Deng; Robert Wang; Yingjie Wang; Heng Zhang; Mingjie Zheng
Journal: Journal of Sensors
 (2019)

---

## Abstract

More and more synthetic aperture radar (SAR) satellites in orbit provide abundant data for remote sensing applications. In August 2016, China launched a new Earth observation SAR satellite, Gaofen-3 (GF-3). In this paper, we utilize a small stack of GF-3 differential interferograms to map land subsidence in Beijing (China) using the time-series SAR interferometry (InSAR) technique. The small stack of differential interferograms is generated with 5 GF-3 SAR images from March 2017 to January 2018. Orbit errors are carefully addressed and removed during differential InSAR (DInSAR) processing. Truncated singular-value decomposition (TSVD) is applied to strengthen the robustness of deformation rate estimation. To validate the results of GF-3 data, an additional deformation measurement using 26 Sentinel-1B images from March 2017 to February 2018 is carried out using the persistent scatterer interferometry (PSI) technique. By implementing a cross-comparison, we find that the retrieved results from GF-3 images and Sentinel-1 images are spatially consistent. The standard deviation of vertical deformation rate differences between two data stacks is 11.24 mm/y in the study area. The results shown in this paper demonstrate the reasonable potential of GF-3 SAR images to monitor land subsidence.


---

## Body

## 1. Introduction

Synthetic aperture radar (SAR) has been widely used for Earth remote sensing for more than 30 years. Currently, more than 15 spaceborne SAR sensors are being operated and several new SAR systems are now under development or in planning [1, 2]. On November 19, 2012, China launched its first civil SAR satellite with an S-band imaging radar onboard, HuanJing-1C (HJ-1C) [3]. HJ-1C is the third satellite of China’s Environmental Protection and Disaster Monitoring Constellation. The overall objective of this constellation is to establish an operational Earth observing system for disaster monitoring and mitigation using remote sensing technology. On August 9, 2016, the Gaofen-3 (GF-3) high-resolution satellite was launched from the Taiyuan Satellite Launch Center. The Gaofen satellite series began in 2013 and are part of the China High-resolution Earth Observation System (CHEOS) which is aimed at forming an earth observation constellation consisting of 9 satellites to provide high-resolution data. The GF-3 satellite is China’s first C-band fully polarimetric SAR satellite [4, 5]. The satellite is designed by the China Academy of Space Technology (CAST). The SAR payload is designed and built by the Institute of Electronics, Chinese Academy of Sciences (IECAS). The GF-3 satellite has a designed lifespan of eight years flying in a highly inclined orbit. It has the ability of using its 12 imaging modes to scan swaths varying from 10 km to 650 km with a corresponding resolution varying from 1 m to 500 m [4]. The primary objective of GF-3 is to make multimode and high-resolution C-band data available for a wide spectrum of both scientific and civilization applications in such fields as ocean monitoring, disaster reduction, and water conservancy [6]. The GF-3 satellite was officially put into operation in January 2017. During the trial run in orbit, GF-3 has provided various users with close to 40,000 images taken from 150 million square kilometers of ocean and land.SAR interferometry (InSAR) has been widely used for investigating surface deformation phenomena over the last two decades [7]. More and more in-orbit SAR sensors provide abundant data for InSAR applications. With the accumulation of repeat-pass SAR images, various advanced time series InSAR (TS-InSAR) approaches have been proposed [7–10]. The TS-InSAR approaches require multiple SAR images to extract high-precision surface deformation measurements. The classical PSI (Permanent Scatterer Interferometry) approach has been proposed by Ferretti et al. [11, 12]. Persistent scatterers used in PSI suffer little from temporal and geometric decorrelation. The PSI method allows a displacement time series measurement with submillimeter accuracy [13]. The small baseline method (for example, Small BAseline Subset (SBAS) [14, 15]) is a multiple-master TS-InSAR approach. The SBAS approach imposes constraints on the maximum temporal and spatial separation between orbits during the interferometric pair combination. Coherent Pixel Technique (CPT) [16] and Interferometric Point Target Analysis (IPTA) [17] have been proposed later in 2003. Subsequently, Spatio-Temporal Unwrapping Network (STUN) [18], Stanford Method for Persistent Scatterers (StaMPS) [19, 20], SqueeSAR [21], Temporarily Coherent point InSAR (TCP-InSAR) [22], and other specific TS-InSAR approaches [7–10] have been proposed to extend the scope of application of TS-InSAR approaches.TS-InSAR has been used for the land subsidence measurements of many populated cities in China, for instance, Beijing [23–27], Tianjin [28], and Shanghai [29, 30]. Beijing, the capital city of China, has suffered from the groundwater-induced subsidence since the late 1950s [23]. Ng et al. studied the long-term displacements over the metropolitan area of Beijing City from 2003 to 2009 using 41 Envisat ASAR images and 24 ALOS-1 PALSAR images [26]. They carried out the three-dimension analysis to discriminate the vertical and horizontal deformation components. Their results pointed out that ground deformation over Beijing was mainly in the vertical direction ranging from -115 mm/y to 6 mm/y. Chen et al. exploited the SBAS method to investigate land subsidence in Beijing due to overextraction of groundwater from 2003 to 2010 using 41 Envisat ASAR images and 14 TerraSAR images [25]. They analyzed the relationships between land subsidence and geoinformation, for instance, groundwater level, active faults, accumulated soft soil thickness, and different aquifer types. Gao et al. investigated the surface deformation in Beijing Plain from 2003 to 2013 using Envisat and TerraSAR images [23]. They pointed out that land subsidence shows obvious seasonal characteristics as the groundwater level changes. The land subsidence time series lags several months behind the groundwater level change. Du et al. studied the ground deformation over the eastern Beijing city using 19 ALOS-1 PALSAR images (June 2007-January 2011), 24 Sentinel-1 images (June 2015-November 2016), and 9 ALOS-2 PALSAR images (September 2014-February 2017) [24]. The study pointed out an increasing trend in the rate and extension of land subsidence in most existing subsiding regions. The previous studies basically reached a consensus that the land subsidence over Beijing is caused by overextraction of groundwater. Therefore, it can be assumed that land subsidence in Beijing is mainly in the vertical direction.The data used in the previous studies with respect to the land subsidence in Beijing mainly focus on the time period before the year of 2017. In this paper, five GF-3 SAR images are utilized to investigate the surface deformation of the metropolitan area in Beijing from March 2017 to January 2018 for demonstrating the interferometric ability of GF-3 SAR data. An additional deformation rate measurement is carried out using the PSI method on 26 Sentinel-1B images from March 2017 to February 2018. We compared the deformation rate maps derived from the GF-3 data and the Sentinel-1 data. The deformation rate maps are consistent in terms of the spatial distribution. The results demonstrate the reasonable potential of GF-3 SAR to monitor surface deformation. Although there have been many studies showing surface deformation measurements in Beijing [23–27], the surface deformation results over the Beijing area using GF-3 data are shown for the first time.This paper is organized as follows. Section2 introduces the GF-3 data and the study area. The applied methodology is described in Section 3. Parameters used in the TS-InSAR processing followed by the experimental results and analysis in details are illustrated in Section 4. The discussions of the experimental results are provided in Section 5. The conclusions are presented in Section 6.


## 2. Data and the Study Area

The high-resolution SAR images can be applied for a multitude of scientific research ranging from geoscience to Earth system monitoring. GF-3 SAR images have been utilized for several applications, such as wind and wave retrieval [31] and vessel detection [32]. The SAR images enable full-time, all-weather ocean and land monitoring, disaster reduction, water conservancy, and meteorology. Moreover, GF-3 data has achieved the ability for InSAR applications [33, 34].In order to meet the accuracy of SAR imaging, the orbit restitution accuracy of the GF-3 satellite is better than 10 m (1σ), and the precise orbit accuracy is better than 20 cm (1σ) [4]. The inaccuracy of the orbit causes significant orbit errors in the InSAR processing. The inaccurate baseline would lead to the flat-earth phase errors, the topographic phase errors, and incomplete topographic contribution removal for D-InSAR [35]. Therefore, orbit errors should be carefully addressed and removed during DInSAR processing. Since the GF-3 satellite is not designed for InSAR applications, repeat-pass SAR images are few. Currently, due to few repeat-pass data and the orbit errors of the satellite platform, the ability to perform time series interferometric measurements using GF-3 data remains to be tested.Five fine stripmap I mode (FSM_I) GF-3 images are acquired in ascending mode, from March 2017 to January 2018 over the metropolitan area of Beijing (China). The map of the study area and the footprints of the GF-3 SAR images are shown in Figure1. The details of the acquired GF-3 data are listed in Table 1. The incidence angle at the center of the image scene is 34.9°. The FSM_I image has a resolution of 5 m and a swath of 50 km [4]. In the TS-InSAR processing, the image acquired on August 22, 2017, is selected as the reference image.Figure 1
Location of Beijing city. The footprints of five GF-3 SAR images are marked with rectangles.Table 1
Basic parameters of the used GF-3 images in the study.
No.
Imaging mode
Orbit ID
Swath width (km)
Spatial resolution (m)
Acquisition date
1
FSM_I
3351
50
5
2017/3/30
2
FSM_I
3770
50
5
2017/4/28
3
FSM_I
5441
50
5
2017/8/22
4
FSM_I
6695
50
5
2017/11/17
5
FSM_I
7530
50
5
2018/1/14Beijing, the capital of China, is located in the northern region of the North China Plain. The geographical coordinates of Beijing are 39°28′-41°05′N and 115°25′-117°35′E. The elevation of the urban area in Beijing is 40 to 60 m. The average annual precipitation is approximately 570 mm, concentrated in the summer season from June to September [27]. Due to the long-term groundwater withdrawal, land subsidence has been one of the most serious geological hazards in Beijing. Many TS-InSAR techniques and different satellite datasets such as Envisat ASAR, TerraSAR-X, RadarSat-2, ALOS-1, and Sentinel-1 have been applied to the surface deformation measurement of Beijing [23–27]. However, it is the first time that GF-3 data are used for the land subsidence measurement in Beijing.


## 3. Methodology

Workflow of a feasible time series InSAR processing chain based on approaches proposed in [14, 36] is shown in Figure 2. The processing chain is mainly divided into two stages, namely, DInSAR processing and time series processing. In DInSAR processing, SAR images are coregistered to the selected reference image. The DInSAR processing steps, i.e., interferogram generation, flat earth phase removal, topographic phase removal, interferogram filtering, phase unwrapping, and baseline refinement, are carried out on each small-baseline InSAR pair sequentially. The DInSAR processing obtains unwrapped differential interferograms from the input single look complex (SLC) SAR data. The primary objective of DInSAR processing is to keep the deformation phase and minimize the unwanted artefacts. In time series processing, a cubic deformation model in the typical SBAS-DInSAR approach [14] is implemented to separate the temporal low-pass deformation and residual topography. Truncated singular-value decomposition (TSVD) [37] is applied to retrieve the deformation rates robustly. Time series processing acquires deformation parameters from the unwrapped differential interferograms. In the following, a brief overview of the noteworthy processing steps is given.Figure 2
Diagram block of the TS-InSAR processing chain. The noteworthy processing steps are marked with red rectangles.


### 3.1. The Mitigation of Decorrelation and Orbit Errors in DInSAR Processing

The limited number of GF-3 images and satellite positioning errors need to be carefully addressed in DInSAR processing. The decorrelation and the orbit errors caused by the aforementioned two questions are identified as the main challenges in retrieving surface deformation using GF-3 images.In the generation of interferograms, range spectral shift filtering and azimuth common band filtering are applied to reduce the decorrelation [38]. These two filtering steps can keep the correlated phase contributions and avoid losing coherence. A complex multilook operation is exploited to reduce the impact of uncorrelated noise due to temporal decorrelation, baseline decorrelation, and volume decorrelation. Through the multilook operation, the speckle noise is reduced and the signal-to-noise ratio of interferogram is improved.The flat-earth phase is removed from the original interferogram. Inaccurate satellite positioning gives rise to time variant errors of the spatial baseline. As a result, there are residual fringes in the flattened interferograms. The residual fringe rates are used for baseline estimation [39, 40], and the flattened interferograms are reflattened. However, there are still few residual orbit errors in the interferogram. After the topography removal and spatial filtering, the differential interferograms are unwrapped. A polynomial-based method is used to perform accurate baseline re-estimation and residual orbit errors removal on unwrapped differential interfergrams [40]. Through the removal of the orbital residuals, the deformation rates can be estimated more accurately. The details will be presented in Section 4.1.


### 3.2. Truncated Singular-Value Decomposition

The critical linear system (Equation (1)) of the SBAS-DInSAR approach [36] is formed for each high coherent pixel. The mean phase velocities vector v can be estimated from the M redundant InSAR observations δϕobs with the rank-deficient design matrix A from the following formation:
(1)AM×N−1v=δϕobs,where N represents the number of SAR acquisitions and v represents the mean phase velocities between time-adjacent acquisitions. In the SBAS approach, the mean phase velocity vector is obtained from v=A‐1δϕobs, where A‐1 is the pseudoinverse of the design matrix A. This problem can be solved by applying singular-value decomposition (SVD) to the design matrix A as follows:
(2)A=UΣVT,Σ=diagσ1,σ2,…,σN−1,σ1≥σ2≥⋯≥σr>σr+1=⋯=σN−1=0,where the left and right singular matrices U∈RM×M and V∈RN−1×N−1 are orthogonal and r=rankA. The truncated singular-value decomposition (TSVD) of A is defined as the rank-k matrix
(3)Ak=UΣkVT=∑i=1kuiσiviT,Σk=diagσ1,…,σk,0,…,0∈RM×N−1,where ui and vi are the columns of the matrices U and V, respectively. Σk equals Σ with the smallest n-k singular values replaced by zeroes, and k<r. The pseudoinverse of the design matrix can be replaced by Ak−1 as follows:
(4)Ak−1=VΣk−1UT,Σk−1=diagσ1−1,…,σk−1,0,…,0∈RN−1×M.The neglected components have a large impact on the solution. By neglecting these components, the errors of solution are mitigated.


## 3.1. The Mitigation of Decorrelation and Orbit Errors in DInSAR Processing

The limited number of GF-3 images and satellite positioning errors need to be carefully addressed in DInSAR processing. The decorrelation and the orbit errors caused by the aforementioned two questions are identified as the main challenges in retrieving surface deformation using GF-3 images.In the generation of interferograms, range spectral shift filtering and azimuth common band filtering are applied to reduce the decorrelation [38]. These two filtering steps can keep the correlated phase contributions and avoid losing coherence. A complex multilook operation is exploited to reduce the impact of uncorrelated noise due to temporal decorrelation, baseline decorrelation, and volume decorrelation. Through the multilook operation, the speckle noise is reduced and the signal-to-noise ratio of interferogram is improved.The flat-earth phase is removed from the original interferogram. Inaccurate satellite positioning gives rise to time variant errors of the spatial baseline. As a result, there are residual fringes in the flattened interferograms. The residual fringe rates are used for baseline estimation [39, 40], and the flattened interferograms are reflattened. However, there are still few residual orbit errors in the interferogram. After the topography removal and spatial filtering, the differential interferograms are unwrapped. A polynomial-based method is used to perform accurate baseline re-estimation and residual orbit errors removal on unwrapped differential interfergrams [40]. Through the removal of the orbital residuals, the deformation rates can be estimated more accurately. The details will be presented in Section 4.1.


## 3.2. Truncated Singular-Value Decomposition

The critical linear system (Equation (1)) of the SBAS-DInSAR approach [36] is formed for each high coherent pixel. The mean phase velocities vector v can be estimated from the M redundant InSAR observations δϕobs with the rank-deficient design matrix A from the following formation:
(1)AM×N−1v=δϕobs,where N represents the number of SAR acquisitions and v represents the mean phase velocities between time-adjacent acquisitions. In the SBAS approach, the mean phase velocity vector is obtained from v=A‐1δϕobs, where A‐1 is the pseudoinverse of the design matrix A. This problem can be solved by applying singular-value decomposition (SVD) to the design matrix A as follows:
(2)A=UΣVT,Σ=diagσ1,σ2,…,σN−1,σ1≥σ2≥⋯≥σr>σr+1=⋯=σN−1=0,where the left and right singular matrices U∈RM×M and V∈RN−1×N−1 are orthogonal and r=rankA. The truncated singular-value decomposition (TSVD) of A is defined as the rank-k matrix
(3)Ak=UΣkVT=∑i=1kuiσiviT,Σk=diagσ1,…,σk,0,…,0∈RM×N−1,where ui and vi are the columns of the matrices U and V, respectively. Σk equals Σ with the smallest n-k singular values replaced by zeroes, and k<r. The pseudoinverse of the design matrix can be replaced by Ak−1 as follows:
(4)Ak−1=VΣk−1UT,Σk−1=diagσ1−1,…,σk−1,0,…,0∈RN−1×M.The neglected components have a large impact on the solution. By neglecting these components, the errors of solution are mitigated.


## 4. Results

### 4.1. DInSAR Processing of the GF-3 Data

First of all, four GF-3 images are coregistered to the selected reference image (20170822). To make full use of the available data stack, all ten interferograms were generated. In order to mitigate the decorrelation, range spectral shift filtering and azimuth common band filtering were applied during the generation of interferograms. A complex multilook operation with 5 looks in the slant range direction and 5 looks in the azimuth direction was exploited to further mitigate the uncorrelated phase noise. Next, the flat-earth phase was removed from the original interferogram (see Figure3(a)) according to the initial baseline. Due to the errors of the initial baseline, there were residual horizontal fringes in the flattened interferograms (see Figure 3(b)). Most of the study area has a relatively flat terrain. The residual orbit errors were removed by subtracting the residual fringes from the flattened interferogram, and the flattened interferogram was refined (see Figure 3(c)). Then, the topographic phase, which was generated using the Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM), was removed from the flattened interferograms. Generic Atmospheric Correction Online Service (GACOS) data [41] was used to remove the atmospheric artefacts from the differential interferogram. Finally, we obtained the differential interferograms (see Figure 3(d)) that contain the deformation phase, residual topographic phase, residual orbit errors, and partial atmospheric artefacts.Figure 3
The interferogram and differential interferogram of 20170330-20170428 interferometric pair. (a) The original interferogram, (b) the flattened interferogram, (c) the flattened interferogram after residual fringe removal, and (d) the wrapped differential interferogram after topographic phase removal.
(a)
(b)
(c)
(d)The coherence coefficient maps were generated using the differential interferograms and their corresponding intensity maps. Then, two InSAR pairs with poor coherence were deleted. The perpendicular baselines of these deleted InSAR pairs are more than 1500 m (1976.9 m and 2315.1 m, respectively). Long perpendicular baseline causes serious geometric decorrelation. Eight InSAR pairs were left for the time series processing (see Figure4). The details of the surplus InSAR pairs are listed in Table 2.Figure 4
Distribution of the interferometric InSAR pairs in the temporal/perpendicular baseline plane.Table 2
GF-3 InSAR pairs used in time series InSAR processing.
No.
Master
Slave
Bperp (meters)
Btemp (days)
Altitude of ambiguity (meters)
1
20170330
20170428
-136
29
-104.99
2
20170330
20170822
-1055.6
145
-13.53
3
20170330
20171117
921.3
232
15.50
4
20170330
20180114
1259.6
290
11.34
5
20170428
20170822
-919.5
116
-15.53
6
20170428
20171117
1057.3
203
13.50
7
20170428
20180114
1395.6
261
10.23
8
20171117
20180114
338.3
58
42.21The corresponding coherence coefficient maps are used to calculate the average coherence coefficient for each pixel. Pixels, whose average coherence coefficient is greater than 0.4 and average intensity is greater than 1000, were selected as candidate pixels of the time series processing. Goldstein filtering [42] was performed on each differential interferogram. To ensure a spatially smooth filtering effect, different window sizes (128, 64, and 32 pixels) were used for iterative filtering operations. The window size was reduced by half in each iteration. The phase unwrapping operation was guided by the generated candidate pixels’ mask using the minimum-cost flow (MCF) algorithm [43]. Because of inaccurate baseline estimation of GF-3, there were still residual orbit errors in the unwrapped interferogram. A polynomial surface was modeled and subtracted from the unwrapped differential interferogram to remove the residual orbit errors. At this point, the unwrapped differential interferogram, which mainly contains the deformation phase, the residual topography phase, the residual orbit errors, and the atmospheric artefacts, were obtained. Except for the deformation phase, the three unwanted artefacts have been minimized in the DInSAR processing as much as possible.


### 4.2. Deformation Parameter Inversion

The unwrapped differential interferograms are used for deformation parameter estimation. A cubic deformation model in the typical SBAS approach [14] was used to separate the temporal low-pass deformation and residual topography. To reduce the contributions from the unwrapping errors of the unwrapped differential interferometric phase, we neglected the smallest singular value in the deformation velocity retrieval according to (4). Here, the number k is chosen manually (k=3 in (3)). Thus, the solution is more stable using the TSVD algorithm.


### 4.3. Experimental Results

The deformation in the study area is assumed mainly in vertical direction, and the line-of-sight (LOS) deformation rates are directly back-projected into the vertical direction using a local incidence angle. The reference point was selected over a relatively stable region in the northeast corner of the Third Ring Road of Beijing. The vertical deformation velocity map of GF-3 data over the study area is shown in Figure5.Figure 5
Linear deformation velocity map of the study area in vertical direction using 5 GF-3 SAR images. Background image: Google Map. The reference point is highlighted by a white triangle. Three regions (a to e) marked with red circles are selected for cross-comparison with Sentinel-1 results.According to the deformation results of GF-3, the deformation rates are very small in the metropolitan area of Beijing. The majority of the subsidence rates in the metropolitan area of Beijing are between -2.0 and 2.0 cm/y, which agree with the measurement from previous studies [24, 26]. In region A (the metropolitan area near imperial palace) of Figure 5, the mean deformation velocity is -7.8 mm/y. The land subsidence of Chaoyang district in eastern Beijing is obvious. Two obvious subsidence regions (Cuigezhuang and Taipingzhuang) marked with red circles in Figure 5 (B and C regions) are observed. The mean deformation velocities of region B (Cuigezhuang) and region C (Taipingzhuang) are -54.2 and -52.3 mm/y, respectively. Since B and C regions are at the edge of the GF-3 data coverage, the estimated deformation velocities cannot represent the maximum deformation velocity of the subsidence bowls. Therefore, the results of Sentinel-1B using the PSI technique were acquired to evaluate the quality of the GF-3 results.Twenty-six Sentinel-1B images from March 2017 to February 2018 (see Table3) are used in PSI processing to get the deformation velocity map. The Sentinel-1B images are acquired in descending orbit (relative orbit 47). The image acquired on July 6, 2017, was selected as master image. The vertical deformation velocity map of Sentinel-1 data is shown in Figure 6. The reference point is the same as the reference point of GF-3 data. The results of GF-3 data and Sentinel-1 data are very similar to those of the previous literatures [25–27] in terms of the spatial distribution. There are several obvious subsidence areas in Chaoyang district and Changping district in Figure 6. Deformation rates in obvious subsidence areas reach or exceed -100 mm/y. From the comparison of the spatial deformation patterns in Figures 5 and 6, we found that the GF-3 results have similar spatial distribution with Sentinel-1 results.Table 3
Information of the used Sentinel-1 data.
No.
Date
Bperp (meters)
Btemp (days)
No.
Date
Bperp (meters)
Btemp (days)
1
20170308
16.3
-120
14
20170823
38
48
2
20170320
66.8
-108
15
20170916
80
72
3
20170401
27.1
-96
16
20170928
28.1
84
4
20170425
81.8
-72
17
20171010
11.7
96
5
20170507
71.1
-60
18
20171022
74.8
108
6
20170519
20.1
-48
19
20171103
107
120
7
20170531
129.1
-36
20
20171115
64.3
132
8
20170612
43.2
-24
21
20171127
11.3
144
9
20170624
37.7
-12
22
20171209
16.3
156
10
20170706
0
0
23
20171221
93.7
168
11
20170718
1.9
12
24
20180102
121.8
180
12
20170730
0.3
24
25
20180207
6.8
216
13
20170811
1.9
36
26
20180219
10.2
228Figure 6
Subsidence of Beijing (China) in the vertical direction mapped by Copernicus Sentinel-1 time series interferometry from March 2017 to February 2018. The reference point is highlighted by a white triangle.A cross-comparison of estimated deformation velocities between GF-3 and Sentinel-1 is investigated to validate the consistency in the study area (yellow rectangle in Figures5 and 6). The deformation velocities are resampled on a 60-meter geographic coordinate grid for a reasonable comparison. The cross-comparison results between two data stacks are shown in Figure 7. In Figure 7(a), when the absolute value of the deformation rate is small (the absolute value less than 40 mm/y), the points are distributed along the line y=x. When the absolute value of the deformation rate is large (the deformation rate less than -40 mm/y), the points are mainly distributed above the line y=x. It is indicated that GF-3 results are underestimated on the pixels whose deformation rate is large. In Figure 7(b), the mean of the deformation rate difference is -1.95 mm/y. The standard deviation of the differences is 11.24 mm/y. It is shown that the deformation rates detected by two data stacks show a relative good agreement.Figure 7
Cross-comparison results between two data stacks. (a) The cross-comparison of the deformation rate between Gaofen-3 and Sentinel-1 data stacks. (b) The histogram of the deformation rate difference for the study area.
(a)
(b)In Figure8, three small characteristic regions marked with red circles in Figure 5 are selected for comparison in details. Region A is located in the metropolitan area near the imperial palace. Region B is located in the subsidence area near Cuigezhuang. Region C is located in the subsidence area near Taipingzhuang. Regions B and C suffer from obvious land subsidence. It can be seen from Figure 8 that the spatial deformation patterns of GF-3 (Figures 8(a)–8(c)) are similar with those of Sentinel-1 (Figures 8(d)–8(f)). The mean time series deformation in three small characteristic regions (A and A′, B and B′, and C and C′) are plotted in Figures 8(h)–8(j), respectively. The time series deformation results of GF-3 and Sentinel-1 have a similar pattern, and both result have subsidence direction.Figure 8
Linear vertical deformation rate map comparison of three selected areas between GF-3 and Sentinel-1 results. GF-3-derived linear deformation rate map of (a) the metropolitan area near the imperial palace. (b) The subsidence area near Cuigezhuang. (c) The subsidence area near Taipingzhuang. Sentinel-1-derived linear deformation rate map of (d) the metropolitan area near the imperial palace. (e) The subsidence area near Cuigezhuang. (f) The subsidence area near Taipingzhuang. Mean time series deformation comparison between (h) regions A and A′, (i) regions B and B′, and (j) regions C and C′.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)We assume that the deformation rates using Sentinel-1 data and the PSI method are true values. Inverse distance-weighted spatial interpolation is performed on the deformation velocity results of Sentinel-1. The interpolated deformation velocities are then extracted to the pixels of GF-3 results for cross-comparison. As shown in Table4, the deformation rates of GF-3 and Sentinel-1 results are numerically approximate. In metropolitan areas, differences of the deformation velocity are small. There are 3036 analyzed pixels in region A used for the cross-comparison. The absolute difference of mean deformation velocities (DV) in region A is 2.18 mm/y, and the standard deviation of differential DV from the two results is 2.39 mm/y. The results in region A are highly consistent between two data stacks. In regions where land subsidence is obvious, the mean of DV differences in regions B and C is 18.1 mm/y and 16.8 mm/y, respectively. The standard deviations of DV differences between the derived results in regions B and C are 5.51 mm/y and 2.97 mm/y, respectively. The large mean of regions B and C shows that the results of GF-3 data in subsidence bowls are underestimated. A small amount of the GF-3 data stack leads to the bias in GF-3 deformation estimation. The absolute mean of deformation velocity difference of regions D and E is 19.6 mm/y and 19.7 mm/y, respectively. In these two regions, the surface deformation estimation of GF-3 data shows a bias. The bias is probably caused by the atmospheric phase.Table 4
Cross-comparison results between GF-3 with Sentinel-1 deformation maps in the selected regions.
Area
Number of points
Mean DV of GF-3 (mm/y)
Mean DV of S1 (mm/y)
Abs. mean of DV difference (mm/y)
Standard deviation of DV difference (mm/y)
A
3036
-7.8
-5.6
2.18
2.39
B
1709
-54.2
-72.3
18.1
5.51
C
1615
-52.3
-69.1
16.8
2.97
D
1050
-17.7
1.9
19.6
2.11
E
1063
-17.6
2.1
19.7
1.64
DV: deformation velocities.


## 4.1. DInSAR Processing of the GF-3 Data

First of all, four GF-3 images are coregistered to the selected reference image (20170822). To make full use of the available data stack, all ten interferograms were generated. In order to mitigate the decorrelation, range spectral shift filtering and azimuth common band filtering were applied during the generation of interferograms. A complex multilook operation with 5 looks in the slant range direction and 5 looks in the azimuth direction was exploited to further mitigate the uncorrelated phase noise. Next, the flat-earth phase was removed from the original interferogram (see Figure3(a)) according to the initial baseline. Due to the errors of the initial baseline, there were residual horizontal fringes in the flattened interferograms (see Figure 3(b)). Most of the study area has a relatively flat terrain. The residual orbit errors were removed by subtracting the residual fringes from the flattened interferogram, and the flattened interferogram was refined (see Figure 3(c)). Then, the topographic phase, which was generated using the Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM), was removed from the flattened interferograms. Generic Atmospheric Correction Online Service (GACOS) data [41] was used to remove the atmospheric artefacts from the differential interferogram. Finally, we obtained the differential interferograms (see Figure 3(d)) that contain the deformation phase, residual topographic phase, residual orbit errors, and partial atmospheric artefacts.Figure 3
The interferogram and differential interferogram of 20170330-20170428 interferometric pair. (a) The original interferogram, (b) the flattened interferogram, (c) the flattened interferogram after residual fringe removal, and (d) the wrapped differential interferogram after topographic phase removal.
(a)
(b)
(c)
(d)The coherence coefficient maps were generated using the differential interferograms and their corresponding intensity maps. Then, two InSAR pairs with poor coherence were deleted. The perpendicular baselines of these deleted InSAR pairs are more than 1500 m (1976.9 m and 2315.1 m, respectively). Long perpendicular baseline causes serious geometric decorrelation. Eight InSAR pairs were left for the time series processing (see Figure4). The details of the surplus InSAR pairs are listed in Table 2.Figure 4
Distribution of the interferometric InSAR pairs in the temporal/perpendicular baseline plane.Table 2
GF-3 InSAR pairs used in time series InSAR processing.
No.
Master
Slave
Bperp (meters)
Btemp (days)
Altitude of ambiguity (meters)
1
20170330
20170428
-136
29
-104.99
2
20170330
20170822
-1055.6
145
-13.53
3
20170330
20171117
921.3
232
15.50
4
20170330
20180114
1259.6
290
11.34
5
20170428
20170822
-919.5
116
-15.53
6
20170428
20171117
1057.3
203
13.50
7
20170428
20180114
1395.6
261
10.23
8
20171117
20180114
338.3
58
42.21The corresponding coherence coefficient maps are used to calculate the average coherence coefficient for each pixel. Pixels, whose average coherence coefficient is greater than 0.4 and average intensity is greater than 1000, were selected as candidate pixels of the time series processing. Goldstein filtering [42] was performed on each differential interferogram. To ensure a spatially smooth filtering effect, different window sizes (128, 64, and 32 pixels) were used for iterative filtering operations. The window size was reduced by half in each iteration. The phase unwrapping operation was guided by the generated candidate pixels’ mask using the minimum-cost flow (MCF) algorithm [43]. Because of inaccurate baseline estimation of GF-3, there were still residual orbit errors in the unwrapped interferogram. A polynomial surface was modeled and subtracted from the unwrapped differential interferogram to remove the residual orbit errors. At this point, the unwrapped differential interferogram, which mainly contains the deformation phase, the residual topography phase, the residual orbit errors, and the atmospheric artefacts, were obtained. Except for the deformation phase, the three unwanted artefacts have been minimized in the DInSAR processing as much as possible.


## 4.2. Deformation Parameter Inversion

The unwrapped differential interferograms are used for deformation parameter estimation. A cubic deformation model in the typical SBAS approach [14] was used to separate the temporal low-pass deformation and residual topography. To reduce the contributions from the unwrapping errors of the unwrapped differential interferometric phase, we neglected the smallest singular value in the deformation velocity retrieval according to (4). Here, the number k is chosen manually (k=3 in (3)). Thus, the solution is more stable using the TSVD algorithm.


## 4.3. Experimental Results

The deformation in the study area is assumed mainly in vertical direction, and the line-of-sight (LOS) deformation rates are directly back-projected into the vertical direction using a local incidence angle. The reference point was selected over a relatively stable region in the northeast corner of the Third Ring Road of Beijing. The vertical deformation velocity map of GF-3 data over the study area is shown in Figure5.Figure 5
Linear deformation velocity map of the study area in vertical direction using 5 GF-3 SAR images. Background image: Google Map. The reference point is highlighted by a white triangle. Three regions (a to e) marked with red circles are selected for cross-comparison with Sentinel-1 results.According to the deformation results of GF-3, the deformation rates are very small in the metropolitan area of Beijing. The majority of the subsidence rates in the metropolitan area of Beijing are between -2.0 and 2.0 cm/y, which agree with the measurement from previous studies [24, 26]. In region A (the metropolitan area near imperial palace) of Figure 5, the mean deformation velocity is -7.8 mm/y. The land subsidence of Chaoyang district in eastern Beijing is obvious. Two obvious subsidence regions (Cuigezhuang and Taipingzhuang) marked with red circles in Figure 5 (B and C regions) are observed. The mean deformation velocities of region B (Cuigezhuang) and region C (Taipingzhuang) are -54.2 and -52.3 mm/y, respectively. Since B and C regions are at the edge of the GF-3 data coverage, the estimated deformation velocities cannot represent the maximum deformation velocity of the subsidence bowls. Therefore, the results of Sentinel-1B using the PSI technique were acquired to evaluate the quality of the GF-3 results.Twenty-six Sentinel-1B images from March 2017 to February 2018 (see Table3) are used in PSI processing to get the deformation velocity map. The Sentinel-1B images are acquired in descending orbit (relative orbit 47). The image acquired on July 6, 2017, was selected as master image. The vertical deformation velocity map of Sentinel-1 data is shown in Figure 6. The reference point is the same as the reference point of GF-3 data. The results of GF-3 data and Sentinel-1 data are very similar to those of the previous literatures [25–27] in terms of the spatial distribution. There are several obvious subsidence areas in Chaoyang district and Changping district in Figure 6. Deformation rates in obvious subsidence areas reach or exceed -100 mm/y. From the comparison of the spatial deformation patterns in Figures 5 and 6, we found that the GF-3 results have similar spatial distribution with Sentinel-1 results.Table 3
Information of the used Sentinel-1 data.
No.
Date
Bperp (meters)
Btemp (days)
No.
Date
Bperp (meters)
Btemp (days)
1
20170308
16.3
-120
14
20170823
38
48
2
20170320
66.8
-108
15
20170916
80
72
3
20170401
27.1
-96
16
20170928
28.1
84
4
20170425
81.8
-72
17
20171010
11.7
96
5
20170507
71.1
-60
18
20171022
74.8
108
6
20170519
20.1
-48
19
20171103
107
120
7
20170531
129.1
-36
20
20171115
64.3
132
8
20170612
43.2
-24
21
20171127
11.3
144
9
20170624
37.7
-12
22
20171209
16.3
156
10
20170706
0
0
23
20171221
93.7
168
11
20170718
1.9
12
24
20180102
121.8
180
12
20170730
0.3
24
25
20180207
6.8
216
13
20170811
1.9
36
26
20180219
10.2
228Figure 6
Subsidence of Beijing (China) in the vertical direction mapped by Copernicus Sentinel-1 time series interferometry from March 2017 to February 2018. The reference point is highlighted by a white triangle.A cross-comparison of estimated deformation velocities between GF-3 and Sentinel-1 is investigated to validate the consistency in the study area (yellow rectangle in Figures5 and 6). The deformation velocities are resampled on a 60-meter geographic coordinate grid for a reasonable comparison. The cross-comparison results between two data stacks are shown in Figure 7. In Figure 7(a), when the absolute value of the deformation rate is small (the absolute value less than 40 mm/y), the points are distributed along the line y=x. When the absolute value of the deformation rate is large (the deformation rate less than -40 mm/y), the points are mainly distributed above the line y=x. It is indicated that GF-3 results are underestimated on the pixels whose deformation rate is large. In Figure 7(b), the mean of the deformation rate difference is -1.95 mm/y. The standard deviation of the differences is 11.24 mm/y. It is shown that the deformation rates detected by two data stacks show a relative good agreement.Figure 7
Cross-comparison results between two data stacks. (a) The cross-comparison of the deformation rate between Gaofen-3 and Sentinel-1 data stacks. (b) The histogram of the deformation rate difference for the study area.
(a)
(b)In Figure8, three small characteristic regions marked with red circles in Figure 5 are selected for comparison in details. Region A is located in the metropolitan area near the imperial palace. Region B is located in the subsidence area near Cuigezhuang. Region C is located in the subsidence area near Taipingzhuang. Regions B and C suffer from obvious land subsidence. It can be seen from Figure 8 that the spatial deformation patterns of GF-3 (Figures 8(a)–8(c)) are similar with those of Sentinel-1 (Figures 8(d)–8(f)). The mean time series deformation in three small characteristic regions (A and A′, B and B′, and C and C′) are plotted in Figures 8(h)–8(j), respectively. The time series deformation results of GF-3 and Sentinel-1 have a similar pattern, and both result have subsidence direction.Figure 8
Linear vertical deformation rate map comparison of three selected areas between GF-3 and Sentinel-1 results. GF-3-derived linear deformation rate map of (a) the metropolitan area near the imperial palace. (b) The subsidence area near Cuigezhuang. (c) The subsidence area near Taipingzhuang. Sentinel-1-derived linear deformation rate map of (d) the metropolitan area near the imperial palace. (e) The subsidence area near Cuigezhuang. (f) The subsidence area near Taipingzhuang. Mean time series deformation comparison between (h) regions A and A′, (i) regions B and B′, and (j) regions C and C′.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)We assume that the deformation rates using Sentinel-1 data and the PSI method are true values. Inverse distance-weighted spatial interpolation is performed on the deformation velocity results of Sentinel-1. The interpolated deformation velocities are then extracted to the pixels of GF-3 results for cross-comparison. As shown in Table4, the deformation rates of GF-3 and Sentinel-1 results are numerically approximate. In metropolitan areas, differences of the deformation velocity are small. There are 3036 analyzed pixels in region A used for the cross-comparison. The absolute difference of mean deformation velocities (DV) in region A is 2.18 mm/y, and the standard deviation of differential DV from the two results is 2.39 mm/y. The results in region A are highly consistent between two data stacks. In regions where land subsidence is obvious, the mean of DV differences in regions B and C is 18.1 mm/y and 16.8 mm/y, respectively. The standard deviations of DV differences between the derived results in regions B and C are 5.51 mm/y and 2.97 mm/y, respectively. The large mean of regions B and C shows that the results of GF-3 data in subsidence bowls are underestimated. A small amount of the GF-3 data stack leads to the bias in GF-3 deformation estimation. The absolute mean of deformation velocity difference of regions D and E is 19.6 mm/y and 19.7 mm/y, respectively. In these two regions, the surface deformation estimation of GF-3 data shows a bias. The bias is probably caused by the atmospheric phase.Table 4
Cross-comparison results between GF-3 with Sentinel-1 deformation maps in the selected regions.
Area
Number of points
Mean DV of GF-3 (mm/y)
Mean DV of S1 (mm/y)
Abs. mean of DV difference (mm/y)
Standard deviation of DV difference (mm/y)
A
3036
-7.8
-5.6
2.18
2.39
B
1709
-54.2
-72.3
18.1
5.51
C
1615
-52.3
-69.1
16.8
2.97
D
1050
-17.7
1.9
19.6
2.11
E
1063
-17.6
2.1
19.7
1.64
DV: deformation velocities.


## 5. Discussion

The spatial deformation maps derived from GF-3 images have a similar distribution pattern with Sentinel-1 results. However, we found that GF-3 results derived from the small stack of differential interferograms are underestimated in the obvious subsidence areas. The reasons are analyzed as follows.Firstly, the number of GF-3 SAR images is small. Therefore, the deformation inversion of GF-3 images is affected more seriously by atmospheric artefacts than that of Sentinel-1 images. In this paper, we used GACOS data to remove the atmospheric artefacts from the differential interferogram. But it cannot guarantee that all atmospheric artefacts are removed. These residual errors may bias the deformation estimation results. To mitigate the phase noise, the GF-3 differential interferograms are spatially filtered; the spatial filtering operation eliminates phase noise and also causes the spatial deformation component to be smoothed. This also degrades the deformation inversion of GF-3 results.Secondly, some seasonal deformation signals are observed in Sentinel-1 results. The time series deformation shows an uplift signal between July and September 2017. This may be because the groundwater level rises during the rainy season. GF-3 results do not demonstrate the uplift signal because of data limitation. In addition, the SVD algorithm in the SBAS method provides a minimum-norm solution; however, some nonlinear movements would be underestimated [16].Again, the GF-3 images and Sentinel-1 images are acquired in ascending and descending orbits, respectively. The results of the two datasets have to be converted to the vertical direction for comparison. The results of cross-comparison between data sets are usually worse than those of the comparison between individual data sets and in situ measurements. In [25], the standard deviation between TerraSAR-X and Envisat ASAR InSAR-derived subsidence rates is 7.48 mm/y. In [26], the standard deviation of the displacement rate difference on the common pixels between Envisat ASAR and ALOS-1 PALSAR is 13 mm/y. In [24], the derived RMSE of the displacement rate difference between Sentinel-1 and ALOS-1 PALSAR is 2.3 cm/y. Considering that the amount of GF-3 data we use is small, a standard deviation of 1.1 cm/y is acceptable. As can be seen in Figure 7(a), the large bias of deformation rate difference mainly occurs in the region where the absolute value of the deformation is large. Compared to the absolute value of deformation rate, the bias can be acceptable.Finally, large perpendicular baselines of GF-3 InSAR pairs have a negative impact on deformation estimation. Firstly, the large perpendicular baseline causes geometric decorrelation, which affects the accuracy of the deformation. In addition, the large perpendicular baseline amplifies the influence of the DEM errors. The residual topographic phase in DInSAR observation will also reduce the deformation estimation accuracy. The GF-3 satellite performs maneuvering to correct the orbit track. With the increase in the number of GF-3 images, the retrieved errors of the deformation rates could be reduced in the future.However, the differences in the estimated results are still acceptable compared with the value of the deformation velocities. With such a small amount of GF-3 SAR images, the precision of derived deformation rates is satisfactory. Moreover, the resolution of GF-3 data is higher than that of Sentinel-1 data, and the GF-3 deformation rate distribution is more detailed.


## 6. Conclusions

In this paper, surface deformation monitoring in Beijing city is investigated using five GF-3 SAR images. The corresponding InSAR processing using a small stack of GF-3 images is demonstrated. The capability of DInSAR interferometry using GF-3 data for detecting and monitoring ground deformation is validated. More importantly, it is the first time that China’s SAR images are used for time series InSAR analysis. Because the main task of a GF-3 satellite is ocean observation, repeat-pass data of land areas have a relatively long temporal and spatial baseline. Therefore, several issues need to be addressed carefully in InSAR processing.This paper describes the processing of the surface deformation measurement using GF-3 data. The main processing stages are summarized in the following:(a)
InSAR Pair Selection. In the urban area, GF-3 InSAR pairs with a long perpendicular baseline hold a moderate correlation. However, two interferograms with a perpendicular baseline above 1500 m are not used for deformation inversion because they are affected by an obvious decorrelation. Therefore, GF-3 InSAR pairs with a perpendicular baseline above 1500 m is not recommended to be used in surface deformation monitoring
(b)
Atmospheric Phase Removal. The differential interferometric phase of GF-3 data is affected by an atmospheric phase screen effect. In the DInSAR processing, GACOS data are used to partially remove the atmospheric phase
(c)
Robust Estimation. Truncated singular-value decomposition (TSVD) is applied to strengthen the robustness of deformation rate estimation. This method has been proved to be effective even if the dataset is smallThis research provides more recently information about the surface deformation in Beijing city using new GF-3 and Sentinel-1 SAR data. By cross-comparison between GF-3 and Sentinel-1 deformation maps, it shows that GF-3 data have the ability to monitor surface deformation. Due to the limitation of data amount, the overall deformation inversion accuracy is not satisfactory. However, this method provides a verification scheme for a small SAR dataset in the absence of in situ data.In the conclusion section, we discuss the impact on InSAR result difference of the data amount, atmospheric phase, seasonal deformation signals, satellite illumination direction, and long perpendicular baseline.



---

*Source: 4204580-2019-04-08.xml*