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.

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 [59], 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) [1722]. In addition, 3D measurements provided by sparse GPS sites can also be employed to aid the LOS measurements for mapping spatial continuous 3D deformations [2325].

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 [1720]. 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 point with respect to a reference point:where is the InSAR-derived relative LOS deformation measurement at point ; is the actual LOS deformation at point ; 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]:where , , and are the relative deformations at point in the vertical, east, and north directions, respectively; , , and are the projection coefficients of the InSAR-derived LOS measurement at point in the vertical, east, and north directions, respectively:where and are the incidence and azimuth angles (clockwise from the north) of point , respectively.

By providing the InSAR-derived LOS relative measurements from ascending and descending tracks, (2) can be rewritten aswhere the subscripts and indicate ascending and descending tracks, respectively; and are the errors of the ascending and descending InSAR-derived LOS relative measurements at the point , 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 asThe 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., ) to construct a joint model:where is the east deformation threshold, usually very small (e.g., 1 mm/year); represents the number of selected points; and are the offsets for the calibrated ascending and descending measurements, respectively. It can be observed that the number of observations (i.e., ) is always larger than the number of unknowns (i.e., ). 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:

The flowchart of the approach is shown in Figure 1. 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.

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 Figures 2(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.

The identified center of the largest deformation funnel is shown by black dots in Figure 4. 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.

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 Figure 5, 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.

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 Figure 6, 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].

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 Figure 6. 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.

Before integrating the deformation velocity fields acquired from cross-heading datasets, the results are resampled into the same grid. The black dot in Figure 8 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.

Figures 8(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.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

The ALOS PALSAR and ENVISAT ASAR data are provided by the Japan Aerospace Exploration Agency (JAXA) and European Space Agency (ESA), respectively. This work was supported by the National Basic Research Program of China (no. 2013CB733303), the Nature Science Foundation of China (nos. 41404011, 41674010, and 51178063), the Science and Technology Project of Hunan Province (nos. 2014FJ3068, 2016SK2002), the Project of Education Department of Hunan Province (nos. 13C1037, 16K053), Major Projects of High Resolution Earth Observation (Civil Part) (Grant no. 03-Y20A11-9001-15/16), the Demonstration System of High-Resolution Remote Sensing Application in Surveying and Mapping of China (no. AH1601-8), and the Advanced Project of Civil Aerospace Research of China: Earth Application and Key Technology Research of 20 m-Resolution Geosynchronous SAR Satellite.