# Estimating Actual 2D Ground Deformations Induced by Underground Activities with Cross-Heading InSAR Measurements

Authors: Xiaoying Yu; Jun Hu; Qian Sun
Journal: Journal of Sensors
 (2017)

---

## Abstract

InSAR can only monitor relative ground deformations with respect to a reference area. In order to obtain actual deformations, GCPs or stable area is required in the study area, which, however, may be unavailable in the investigating of geohazards associated with underground activities (i.e., groundwater pumping, underground mining, and oil/gas exploitation). We propose a novel approach to estimate actual 2D deformations based on the InSAR relative LOS measurements acquired from cross-heading datasets. The errors induced by the arbitrary selection of reference areas can thus be avoided. The performance of the proposed approach is validated by a series of simulations. By providing the ascending and descending measurements with errors of 2 and 1.5 mm/year STDs, respectively, the RMSEs are 2.1 and 2.6 mm/year for the estimated vertical and east deformations, respectively. A case study is carried out in Cangzhou, China, for estimating the actual 2D ground deformations associated with groundwater pumping. By integrating ALOS ascending and ENVISAT descending datasets acquired between 2007 and 2010, we found that the Cangzhou area experienced ground subsidence of up to 23.4 mm/year in the suburbs but ground uplift of up to 20.9 mm/year in the urban area, both of which are accompanied by considerable lateral deformations.


---

## Body

## 1. Introduction

Recently, Interferometric Synthetic Aperture Radar (InSAR) has been used as a routine tool in monitoring ground deformations associated with geohazards such as earthquake [1], volcano eruption [2], glacier movement [3], and landslide [4]. With the appearance of a number of multitemporal InSAR (MT-InSAR) algorithms [5–9], the inherent errors (i.e., decorrelation noises, topographic residuals, orbital ramps, and atmospheric artifacts) in the differential InSAR (D-InSAR) measurements can be well suppressed by analyzing a time series of SAR datasets. This allows us to map slow and long-term ground deformations caused by groundwater pumping [10], oil and gas exploitations [11], underground mining [12], consolidation of compressible deposits [13], permafrost freezing and thawing [14], and so forth.However, there are still two problems impeding the application of InSAR in the estimation of ground deformations. First, only the projection on the line-of-sight (LOS) direction of the three-dimensional (3D) ground deformations can be reflected in D-InSAR or MT-InSAR measurements, because of the polar orbit and side imaging mode of the current SAR satellites [15]. Second, relative ground deformations with respect to a reference point can only be achieved by InSAR, since the deformation measurements are reestablished from the interferometric phase differences between the adjacent pixels by conducting a phase unwrapping operation [16].During the past decade or so, many efforts have been made to solve the first problem. The 3D ground deformation can now be determined by integrating the InSAR-derived LOS measurements with offset-tracking or multiaperture InSAR-derived azimuth measurements from cross-heading orbits (i.e., ascending and descending tracks) [17–22]. In addition, 3D measurements provided by sparse GPS sites can also be employed to aid the LOS measurements for mapping spatial continuous 3D deformations [23–25].However, until now no much attention has been paid to the second issue. It is acknowledged that relative deformations can be calibrated by ground control points (GCPs) with known actual deformations [26]. But GCPs are not always available. More importantly, except for corner reflectors, it is quite difficult to distinguish the GCPs in the SAR image [27]. In areas without GCPs, the common method is to assume a relatively stable region as the reference area [16]. For instance, in researches on earthquake or volcano eruption, far field areas are generally assumed as the reference area in order to determine the actual ground deformation in the near field [17–20]. Obviously, this method is not suitable in monitoring large-scale ground deformations, which exceed the scope of a single frame SAR image. Therefore, it is of great importance to develop a method for estimating actual deformations that is not limited by the aforementioned assumptions, especially for the InSAR measurements acquired from current high-resolution SAR data with relatively small coverage.In this paper, a novel approach is proposed for estimating actual two-dimensional (2D) ground deformations based on cross-heading InSAR relative deformation measurements. This approach is designed to monitor long-term deformations associated with underground activities (e.g., groundwater pumping, underground mining, and oil/gas exploitation). With this algorithm, we can determine the actual deformation velocities in both vertical and east directions using only InSAR LOS measurements acquired from ascending and descending tracks, which means GCPs or the assumed stable area is not required. The proposed approach is first assessed by the simulated experiments and then applied to map the actual 2D ground deformation velocities in Cangzhou, China, which are caused by groundwater pumping.


## 2. Methodology

Ground deformations are the direct consequences of geohazards, but InSAR only can measure the relative LOS deformation at pointxi with respect to a reference point:(1)Drexi=Dabxi-K,where Drexiis the InSAR-derived relative LOS deformation measurement at point xi; Dabxi is the actual LOS deformation at point xi; K is the actual offset, that is, the actual LOS deformation at the reference point. On the other hand, the InSAR-derived relative LOS deformation measurement is the projection of 3D relative deformations onto the LOS direction [28]:(2)Drexi=Suxi,Sexi,Snxi·Durexi,Derexi,DnrexiT,where Durexi, Derexi, and Dnrexi are the relative deformations at point xiin the vertical, east, and north directions, respectively; Suxi, Sexi, and Snxi are the projection coefficients of the InSAR-derived LOS measurement at point xi in the vertical, east, and north directions, respectively:(3)Suxi=cos⁡θxi,Sexi=-cos⁡αxi·sin⁡θxi,Snxi=sin⁡αxi·sin⁡θxi,where θxi and αxi are the incidence and azimuth angles (clockwise from the north) of point xi, respectively.By providing the InSAR-derived LOS relative measurements from ascending and descending tracks, (2) can be rewritten as(4)DasrexiDdesrexi=SasuxiSasexiSasnxiSdesuxiSdesexiSdesnxi·DurexiDerexiDnrexi+ηasrexiηdesrexi,where the subscripts as and des indicate ascending and descending tracks, respectively; ηasrexi and ηdesrexi are the errors of the ascending and descending InSAR-derived LOS relative measurements at the point xi, respectively. It should be noted that the spatial-temporal basis of the ascending and descending measurements should be unified. Therefore, the ascending and descending deformation measurements are transformed into the average velocity values by dividing their corresponding time intervals. The 3D deformations thus become 3D deformation velocities. In addition, the ascending and descending measurements are geocoded into the same lattice.Equation (4) cannot be solved since the unknowns are more than the observations. Even though there are more observations, under the current SAR imaging geometry, the 3D deformations cannot be accurately determined by the InSAR-derived LOS measurements, because of their insensitivities to the north component [15]. A feasible method is to estimate the vertical and east deformation components by neglecting the contribution of north deformation component in the InSAR-derived LOS measurements [28]. So (4) can be rewritten as(5)DasrexiDdesrexi=SasuxiSasexiSdesuxiSdesexi·DurexiDerexi+ηasrexiηdesrexi.The relative deformations in the vertical and east directions can then be estimated by conducting a least squares inversion on a pixel-by-pixel basis.The actual horizontal deformations induced by underground activities (e.g., groundwater pumping, underground mining, and oil/gas exploitation) are generally symmetrically distributed with respect to the center of the deformation funnel, and the deformation rates vary gradually from the center to the margin [29]. This is an ideal but applicable assumption in the investigating of underground activities, even though faults and layering are included in the surrounding earth medium [30, 31]. Therefore, the absolute value of the actual vertical deformation at the funnel center is the maximum, while the actual east deformation is close to zero. For the InSAR-derived relative deformations, the absolute value of relative vertical deformation at the funnel center is also the maximum, despite the arbitrarily selected reference point for the InSAR measurements. Based on this, we can identify the location of the center of the largest deformation funnel in the field. Taking that location as the reference area, we then calibrate the ascending and descending InSAR LOS measurements. Following (5), the relative deformations in the vertical and east directions are reestimated from the calibrated InSAR measurements.At this step, the relative east deformations should symmetrically distribute with respect to the center of the deformation funnel and are very small near the center of the deformation funnel. Based on (1) and (5), we select the points with negligible east relative deformations (i.e., Derexi≤ε) to construct a joint model:(6)Dasrey1⋮DasreyNDdesrey1⋮DdesreyN=Sasuy1-10⋱⋮⋮SasuyN-10Sdesuy10-1⋱⋮⋮SdesuyN0-1·Duaby1⋮DuabyNKasKdes+ηasrey1⋮ηasreyNηdesrey1⋮ηdesreyN,where ε is the east deformation threshold, usually very small (e.g., 1 mm/year); N represents the number of selected points; Kas and Kdes are the offsets for the calibrated ascending and descending measurements, respectively. It can be observed that the number of observations (i.e., 2N) is always larger than the number of unknowns (i.e., N+2). A sparse least squares algorithm [32] can be applied to resolve (6) and thus obtain the actual 2D deformations at the selected points as well as the offsets for the calibrated ascending and descending InSAR LOS measurements.After correcting the calibrated LOS deformations with the offsets, the actual vertical and east deformations of the whole field can finally be estimated by applying the least squares inversion to the following model:(7)DasabxiDdesabxi=Dasrexi+KasDdesrexi+Kdes=SasuxiSasexiSdesuxiSdesexi·DuabxiDeabxi+ηasabxiηdesabxi.The flowchart of the approach is shown in Figure1. It should be highlighted that the ascending and descending InSAR-derived relative LOS deformation measurements can be obtained by different sensors. Therefore, the proposed approach can be easily applied using the SAR datasets provided by various satellites. However, the ascending and descending SAR datasets are not acquired simultaneously over the same area, which affect the application of the proposed method. Therefore, only the average actual velocity can be estimated by the proposed method, which might be insufficient to describe the ground deformations dominated by nonlinear behavior.Figure 1
Flowchart of the proposed approach.


## 3. Simulated Experiments

A simulation experiment is first conducted to evaluate the performance of the proposed approach. Actual 3D ground deformations caused by underground exploitation are simulated over a 400 × 450 grid, with a grid size of 10 m. The actual deformation velocities fields in the vertical, east, and north directions are shown in Figures2(a), 2(b), and 2(c), respectively. In order to make the simulation realistic, an uplift funnel and a subsidence funnel are generated in the field. The InSAR-derived relative LOS measurements from ascending and descending tracks are then obtained using (2). The SAR system parameters (i.e., the incidence and azimuth angles) are adopted from the ALOS PALSAR ascending and ENVISAT ASAR descending data over Cangzhou, China, which will also be employed in the next section. As shown in Figure 3, zero mean additive Gaussian noises with 2 and 1.5 mm/year standard deviations (STDs) are added to the ascending and descending measurements, respectively. The reference points are arbitrarily selected for the ascending and descending InSAR relative LOS measurements (see the triangles in Figure 3). Obviously, the InSAR relative measurements cannot describe actual ground deformations and even induce misinterpretation of the geohazards due to the arbitrarily selected reference point.Figure 2
Simulated actual deformation velocities fields in the vertical (a), east (b), and north (c) directions. Units: cm/year.
(a)
(b)
(c)Figure 3
Simulated relative LOS measurement from ascending track (a) and descending track (b). Units: cm/year. The triangles represent the arbitrarily selected reference points.
(a)
(b)The identified center of the largest deformation funnel is shown by black dots in Figure4. We select 16125 points with the reestimated relative east deformations smaller than the deformation threshold (1 mm/year in this study). These points are then used in (6), from which the offsets are estimated for the calibrated ascending and descending LOS measurements. Figures 4(a) and 4(b) show the vertical and east deformation velocity fields, respectively, which are estimated from the offsets-corrected ascending and descending LOS measurements. The two components both agree well with the simulated ones. The differences between the estimated and the simulated actual deformations in the vertical and east directions are exhibited in Figures 4(c) and 4(d), respectively. It seems that the differences are dominated by north-south trending ramps and Gaussian noises. This is somewhat expected since neglecting the actual north deformations will inevitably lead to some errors to the 2D estimations. In addition, the errors will be propagated from the InSAR measurements to the offsets of the calibrated InSAR measurements, which will be a constant basis for the 2D deformation estimations. In order to quantitatively assess the proposed approach, the root mean square errors (RMSEs) of the differences are calculated. The RMSEs of the actual vertical and east velocities are 2.1 and 2.6 mm/year, respectively, both of which are comparable to the added InSAR observation noises. This demonstrates that the proposed approach can obtain accurate actual ground deformations in both the vertical and east directions.Figure 4
(a, b) Estimated actual deformation velocities from the proposed approach. Black dots represent the identified location of the center of the largest deformation funnel. (c, d) Differences between the estimated and simulated actual deformation velocities. (a, c) The vertical component. (b, d) The east component. Units: cm/year.
(a)
(b)
(c)
(d)In order to further assess the performance of the proposed approach in the presence of noises, we conduct a series of simulated experiments providing the InSAR relative LOS measurements with different levels of noises. With respect to the ascending measurements, additive Gaussian noises with 0–10 mm/year STDs are added, respectively. Meanwhile, 75% of the noises of the ascending measurements are added to the corresponding descending measurements. As shown in Figure5, the RMSEs of the estimations are positively correlated to the STDs of the InSAR relative LOS measurements for both of the vertical and the east components. The east component increases faster than the vertical component, indicating that the east estimation is more sensitive to the InSAR noises. If the noises of the ascending InSAR measurements are up to 10 mm/year, the accuracies of the vertical and east estimations are higher than 8 and 13 mm/year, respectively. The RMSEs of the vertical component is nearly twice of that of the east component when no noise is added to the InSAR relative LOS measurements. This might be induced by the neglecting of the north deformation component in the estimation of the other two components. The results reveal that the vertical estimation can better resist the noises of the InSAR relative LOS measurements but is more sensitive to model misfit errors than the east estimation.Figure 5
RMSEs of the vertical and east deformation velocity estimations versus the noises of the InSAR relative LOS measurements.


## 4. Case Study in Cangzhou, China

The proposed approach is applied in the investigation of the actual ground deformations in Cangzhou, China. As shown in Figure6, Cangzhou is located in the eastern part of Hebei province, close to the Pohai Sea. Deep groundwater has long been the main source of the industrial, agricultural, and domestic water in Cangzhou [33]. Because of the overexploitation of deep groundwater, Cangzhou has experienced severe ground subsidence since the 1970s [34, 35], which has caused a lot of the environmental problems and geohazards. In order to control and solve the problems, a series of measures (e.g., shutting down wells, pumping-limit, and seawater irrigation) have been taken since 1998 [35].Figure 6
Optical map of the study area. Boxes indicate the coverage of the ALOS PALSAR ascending and ENVISAT ASAR descending datasets. The inset map shows the location of the study area in China.In total, 21 ALOS PALSAR images acquired in ascending track and 16 ENVISAT ASAR images acquired in descending track are collected in the study. The spatial coverage of the two SAR datasets is outlined in Figure6. The time intervals of the SAR images are nearly four years, that is, from January 17, 2007, to October 28, 2010, for ALOS PALSAR ascending dataset and from February 23, 2007, to October 15, 2010, for ENVISAT ASAR descending dataset. As shown in Tables 1 and 2, 44 and 32 small-baseline interferograms are generated from the 21 PALSAR ascending images and 16 ASAR descending images, respectively. A well-developed MT-InSAR algorithm, that is, Temporarily Coherent Point InSAR (TCP-InSAR) [9, 36, 37], is used to process the two tracks of datasets. The phase unwrapping errors and atmospheric observations can be well suppressed by the outlier detector and local Delaunay triangulation, respectively, both of which are the uniqueness of the TCP-InSAR algorithm. Note that external DEM data is not required to remove the topographic contribution since the terrain of Cangzhou is quite flat. As shown in Figure 7, the relative LOS deformation velocity fields are estimated from the PALSAR ascending and ASAR descending tracks, respectively, by selecting an arbitrary point as the reference area (see the black triangles in Figure 7). Several large-scale funnels are formed in both of the relative LOS deformation velocity fields (Figure 7). Since almost the whole investigated area is affected by the ground deformation, it is impossible to find a stable area as the reference area. However, because of the arbitrarily selected reference points, the intervals of the relative deformation fields acquired from ascending and descending tracks are quite different, probably resulting in a misinterpretation for the geohazards in Cangzhou.Table 1
Basic information of the InSAR pairs for the case study in Cangzhou, China, corresponding to ALOS PALSAR ascending dataset.
Number
Master
Slave
Perpendicular baseline (m)
Temporal baseline (days)
1
20070117
20081022
-
276.4
644
2
20070117
20081207
-
120.8
690
3
20070117
20090122
356.6
736
4
20070117
20090309
798.1
782
5
20070304
20070720
727.5
138
6
20070304
20080722
-
341.0
506
7
20070304
20090424
-
453.4
782
8
20070304
20090909
-
201.5
920
9
20070720
20071020
440.5
92
10
20070720
20100125
98.8
920
11
20070720
20100312
698.5
966
12
20071020
20080120
596.5
92
13
20071020
20100125
-
341.7
828
14
20071020
20100312
258.0
874
15
20071020
20100427
227.3
920
16
20071020
20100612
227.9
966
17
20080120
20100312
-
338.4
782
18
20080120
20100427
-
369.2
828
19
20080120
20100612
-
368.6
874
20
20080120
20100728
13.9
920
21
20080421
20080606
-
41.1
46
22
20080421
20101028
-
545.4
920
23
20080722
20090309
-
662.0
230
24
20080722
20090424
-
112.5
276
25
20080722
20090909
139.5
414
26
20080906
20081022
582.6
46
27
20080906
20081207
738.2
92
28
20081022
20081207
155.5
46
29
20081022
20090122
633.0
92
30
20081207
20090122
477.5
46
31
20090122
20090309
441.5
46
32
20090309
20090424
549.4
46
33
20090424
20090909
252.0
138
34
20090909
20091025
278.0
46
35
20100125
20100312
599.6
46
36
20100125
20100427
568.9
92
37
20100125
20100612
569.4
138
38
20100312
20100427
-
30.7
46
39
20100312
20100612
-
30.1
92
40
20100312
20100728
352.3
138
41
20100427
20100612
0.5
46
42
20100427
20100728
383.0
92
43
20100612
20100728
382.5
46
44
20100728
20101028
521.7
92Table 2
Basic information of the InSAR pairs for the case study in Cangzhou, China, corresponding to ENVISAT ASAR descending dataset.
Number
Master
Slave
Perpendicular baseline (m)
Temporal baseline (days)
1
20070223
20070608
-
163.2
105
2
20070223
20070817
-
150.1
175
3
20070608
20070817
13.0
70
4
20070608
20071026
-
163.2
140
5
20070608
20080314
-
82.6
280
6
20070817
20071026
-
176.2
70
7
20070817
20080314
-
95.7
210
8
20070817
20080523
4.3
280
9
20071026
20080104
-
241.5
70
10
20071026
20080314
80.5
140
11
20071026
20080523
180.5
210
12
20080314
20080523
100.0
70
13
20080314
20081010
79.4
210
14
20080314
20081219
-
105.3
280
15
20080523
20081010
-
20.6
140
16
20080523
20081219
-
205.3
210
17
20080523
20090123
134.7
245
18
20081010
20081219
-
184.7
70
19
20081010
20090123
155.3
105
20
20081010
20090717
126.1
280
21
20090123
20090717
-
29.2
175
22
20090123
20091030
-
189.5
280
23
20090717
20091030
-
160.3
105
24
20090717
20091204
32.4
140
25
20090717
20100319
21.2
245
26
20091030
20091204
192.7
35
27
20091030
20100319
181.5
140
28
20091030
20100702
231.7
245
29
20091204
20100319
-
11.1
105
30
20091204
20100702
39.0
210
31
20100319
20100702
50.1
105
32
20100702
20201015
273.7
105Figure 7
TCP-InSAR-derived relative LOS deformation velocity fields of Cangzhou. (a) ALOS PALSAR ascending track; (b) ENVISAT ASAR descending track. The triangles represent the arbitrarily selected reference points.
(a)
(b)Before integrating the deformation velocity fields acquired from cross-heading datasets, the results are resampled into the same grid. The black dot in Figure8 represents the identified location of the center of the largest deformation funnel in Cangzhou. In total, 38816 points are then selected based on an east deformation threshold of 1 mm/year, yielding the offsets for the calibrated ascending and descending LOS measurements, respectively.Figure 8
Actual 2D deformation velocity fields of Cangzhou estimated by the proposed approach. Black dots represent the identified location of the center of the largest deformation funnel. (a) The vertical component. (b) The east component. CZC: Cangzhou City; BTC: Botou City; CC: Cang County; QC: Qing County; MC: Mengcun County.
(a)
(b)Figures8(a) and 8(b) show the estimated actual vertical and east deformation velocity fields from 2007 to 2010 in Cangzhou, respectively, which are superimposed on the optical map of the Zhouqu area. Clearly, Cangzhou experienced obvious and large-scale vertical ground deformations, ranging from −23.4 to 20.9 mm/year. Ground subsidences concentrate in two areas: one is located in the southeastern of Qing County; the other is in the northwestern of Cang County. This subsidence can be ascribed to extensive groundwater pumping for irrigation and domestic usage [35]. However, the Cangzhou City experienced ground uplift instead of ground subsidence during 2007–2010, because most production wells in the urban area had been shut down, resulting in an increasing rise of the underground water level since 2005 [33]. Another ground uplift occurs between the Baotou City and Mengcun County, with smaller amplitude but larger coverage. Figure 8(b) shows nonnegligible east deformations in and around the area experiencing obvious ground subsidence or uplift, ranging from −22.0 to 12.0 mm/year. The most severe lateral deformation occurs in the Baotou City and the Qing County. It is also found that most of the lateral deformations move toward the center of subsidence and away from the center of uplift, especially those in the urban area of Cangzhou. This is also reasonable since the lateral deformations associated with underground exploitation are generally proportional to the gradient of vertical deformations, which can partly be responsible for the ground fissures in the Zhouqu area [38].Quantitative assessment, however, cannot be conducted in the Cangzhou experiment since the field observations are unavailable in the study. Nevertheless, the InSAR-derived relative LOS deformations in this study are comparable to other InSAR measurements reported in the previous studies [35, 39, 40]. The actual deformation results are also compatible with the geologic and hydrologic surveys [33, 41]. Compared to the relative LOS deformations retrieved from the single-track SAR dataset, the actual vertical and east deformations retrieved from the cross-heading SAR datasets are much more useful in the prevention of the potential geohazards as well as the future plan of groundwater exploitation in the Cangzhou area.


## 5. Conclusions

Estimation of actual deformation is of great importance for InSAR-derived deformation results that are inborn with relative values. In this study, we propose a novel approach that can obtain actual vertical and horizontal ground deformations associated with underground exploitations by integrating the ascending and descending InSAR-derived relative measurements. We conduct a series of simulated experiments to assess the feasibility and accuracy of the proposed method, which is then employed in the investigation of the actual ground deformations in Cangzhou, China.When GCPs and stable area are not available in the investigated area, the proposed method can be employed in the estimation of actual ground deformation velocities induced by underground exploitations. Although some model misfit errors are presented in the vertical and east estimations due to the neglecting of north deformation component, the results reveal that comparable accuracies with respect to the added InSAR noises can be achieved for both components. It is also found that the vertical estimation is superior in resisting the InSAR noises, especially when noises (STDs > 2 mm) are large.In the Cangzhou case, the actual vertical and east ground deformation velocities between 2007 and 2010 are retrieved by integrating the ALOS PALSAR ascending and ENVISAT ASAR descending measurements with the proposed approach. Ground subsidence of up to 23.4 mm/year has been detected in the southeast of Qing County and the northwest of Cang County, which can be ascribed to the groundwater pumping. The Cangzhou City experienced ground uplift of up to 20.9 mm/year during the study period, which results from a series of control measures such as shutting down wells, pumping-limit, and seawater irrigation. Meanwhile, considerable lateral deformations have been found in Cangzhou, highly related to the vertical deformations. These results are useful in the investigation of geohazards and groundwater in the Cangzhou area.More efforts can be made to improve the proposed method in the future. Firstly, the east deformation threshold (i.e., 1 mm) is empirically determined in the study; a variable one should be designed based on the noise levels of InSAR measurements. Secondly, resampling error is currently induced by the integration of cross-heading SAR datasets but can be avoided by identifying the homonymy points between ascending and descending datasets. Thirdly, the precision levels of different satellite observations are variant due to the differences of their wavelengths. The variances (i.e., weights) of the satellite observations should be introduced in the least squares adjustment by exploiting such as the coherences of the interferograms.



---

*Source: 3170506-2017-08-27.xml*