Regional Land Subsidence Analysis in Eastern Beijing Plain by InSAR Time Series and Wavelet Transforms

Land subsidence is the disaster phenomenon of environmental geology with regionally surface altitude lowering caused by the natural or man-made factors. Beijing, the capital city of China, has suffered from land subsidence since the 1950s, and extreme groundwater extraction has led to subsidence rates of more than 100 mm/year. In this study, we employ two SAR datasets acquired by Envisat and TerraSAR-X satellites to investigate the surface deformation in Beijing Plain from 2003 to 2013 based on the multi-temporal InSAR technique. Furthermore, we also use observation wells to provide in situ hydraulic head levels to perform the evolution of land subsidence and spatial-temporal changes of groundwater level. Then, we analyze the accumulated displacement and hydraulic head level time series using continuous wavelet transform to separate periodic signal components. Finally, cross wavelet transform (XWT) and wavelet transform coherence (WTC) are implemented to analyze the relationship between the accumulated displacement and hydraulic head level time series. The results show that the subsidence centers in the northern Beijing Plain is spatially consistent with the groundwater drop funnels. According to the analysis of well based results located in different areas, the long-term groundwater exploitation in the northern subsidence area has led to the continuous decline of the water level, resulting in the inelastic and permanent compaction, while for the monitoring wells located outside the subsidence area, the subsidence time series show obvious elastic deformation characteristics (seasonal characteristics) as the groundwater level changes. Moreover, according to the wavelet transformation, the land subsidence time series at monitoring well site lags several months behind the groundwater level change.


Introduction
Land subsidence in metropolitan area has been observed as series of disastrous phenomena over decades worldwide [1].For example, more than 60 countries [2] as well as more than 45 cities in China [3] have suffered from land subsidence.The spatial variety of land subsidence in urban areas often leads to damage to urban infrastructures including buildings, airports, railways, highways, subways, and other underground facilities.Land subsidence in urban areas is mainly caused by human activities, including over-exploitation of groundwater [4] or oil/gas [5], underground works, and massive construction of high buildings.As reported, Beijing is plagued by a serious shortage of water resources, and 2/3 of water demand is supplied by groundwater [6].Nevertheless, 87.6% of the whole plain area is suffering from excessive groundwater withdrawal.The long-term and high-intensity groundwater consumption has resulted in water level falling and local city sinking [7].
Land subsidence, excessive groundwater withdrawal and thick haze are the new unbearable pains of Beijing.Beijing has suffered from land subsidence for decades [8].The first record of subsidence in Beijing was discovered by the survey and mapping departments in 1935, near Xidan.The measured maximum accumulative displacement was only 58 mm until 1952.Nevertheless, since the 1950s, with the rapid development of the light industry in eastern Beijing, many wells were drilled and lots of water was unboundedly pumped from underground for industrial consumption.The local land subsidence area has gradually formed since.As reported by the China Geological Survey (CGS), the groundwater level in Beijing has dropped rapidly and has exceeded natural recharge [9] since the 1970s, due to industrial development and the demands of a growing population.Especially since 1999, groundwater overdraft due to drought has resulted in rapid development of the land subsidence in Beijing.
A necessary step to perform a proper analysis of the land subsidence is to obtain accurate measurements of the actual amount of subsidence at certain intervals.Compared with the traditional leveling and global positioning system (GPS) measurement methods, the interferometric synthetic aperture radar (InSAR) technique can obtain a wide spatial range of surface deformation information with millimeter-scale precision [10].In this case, InSAR has been developed as a new geodetic technique over the past few decades.InSAR has been used for the measurement of topography and surface deformation, such as seismic deformation [11], volcanic-tectonic deformation [12], glacial kinematics [13], landslide [14], delta sinking [15], and land subsidence [16][17][18][19].However, conventional InSAR can only extract the phase shifts between SAR image acquisition dates, rather than the time-series deformation rates, because it relies only on two SAR images with a short time span [20], and not on a stack.In addition, the application of conventional InSAR is mainly limited by temporal decorrelation, spatial decorrelation, and atmospheric delay, which are caused by different satellite observation positions, long image-acquisition intervals, and atmospheric fluctuations [21], respectively.All these issues can be addressed using time-series InSAR analysis methods.Temporal InSAR techniques such as persistent scatterer interferometry (PSI) [22,23] and Small BAseline Subset (SBAS) [24,25] have been successfully used to detect, map and analyze long-term and slow land subsidence [26][27][28][29][30]. Previous study showed that the land subsidence in Beijing was comprehensively affected by over-exploitation of groundwater [31][32][33][34], urban construction [35], and geological fault structures, and that the mechanisms of land subsidence are complex in different areas.
In urban areas, land subsidence is a serious threat to urban infrastructure, high-speed railway and the utilization of underground space, thus restricting the sustainable development of society.Therefore, the study of the regional subsidence evolution in Beijing Plain is of great significance.As mentioned above, previous studies focused on time series InSAR technology to obtain deformation measurement combined with groundwater level changes to analyze the statistical correlation between them.In this study, a newly proposed InSAR time series analysis approach is implemented on 47 Envisat ASAR images and 33 TerraSAR-X images to estimate the time-series deformation of the east Beijing Plain (Figure 1a).Extra leveling data are used to calibrate the InSAR based measurements to ensure the reliability of the observation.Then, wavelet multi-scale analysis is applied to seize the responses of the time-varying surface deformation to hydraulic head levels changes in Beijing Plain.This paper is organized as follows.We introduce the history of land subsidence in Beijing, and summarize the study area (Section 2).Then, the InSAR time series technique and the wavelet transform used in this study are presented in Section 3. We detail the subsidence rate and time series analysis from InSAR results, compare these with hydraulic head levels, and then present the results of the wavelet transform for both the land subsidence and hydraulic head level (HHL) time series (Section 4).Finally, discussions from this study are given in Section 5, and our concluding remarks are presented in Section 6.
summarize the study area (Section 2).Then, the InSAR time series technique and the wavelet transform used in this study are presented in Section 3. We detail the subsidence rate and time series analysis from InSAR results, compare these with hydraulic head levels, and then present the results of the wavelet transform for both the land subsidence and hydraulic head level (HHL) time series (Section 4).Finally, discussions from this study are given in Section 5, and our concluding remarks are presented in Section 6.

Study Area
Behind Yanshan Mountains to the north and Taihang Mountains to the southwest and facing the sea to the east, Beijing Plain lies on the northern edge of the North China Plain and covers an area of over 6300 km 2 (Figure 1a).The climate is temperate semi-humid continental monsoon type: the annual average temperature is approximately 10-12 • C, and the mean annual precipitation is 510 mm.The Beijing Plain is a typical piedmont alluvial diluvial plain formed from sediment carried by the Yongding, Chaobai, and Wenyu Rivers.The Quaternary sediments from the piedmont to the plain areas are generally divided into three parts: alluvial-pluvial fan in the top, alluvial-pluvial fan in the middle, and fringe part of alluvial-luvial fan and alluvial plain area [36].In the middle part are 2-3 layers of sand gravel strata, and, in the fringe and alluvial plain area, it transfers gradually into the multi-layer structure of the coarse sand, medium sand, and fine sand [37].The soft sediment in the pluvial-alluvial plain contains several aquifers that are compressible in natural.It is confirmed that extraction of groundwater results in the compaction of alluvial aquifer systems and generates surface deformation [10,38].Previous study showed that the spatial distribution of land subsidence in Beijing Plain was consistent with that of the quaternary compressible layers [39].The thickness of the quaternary alluvial fan of the Yongding River is more than 340 m.One of the subsidence centers, Houshayu (marked by red circle in Figure 1a), is located in the Wenyu River alluvial fan with a thickness of compressible sediments over 150 m.
Since 1999, as shown in Figure 1a, six sinking bowls have gradually formed in the eastern Beijing Plain.Several active faults (Liangxiang-qianmen-shuyi fault, Nankou-Sunhe fault, Nanyuan-tongxian fault and Huangzhuang-Gaoliying fault) cross throughout the whole sinking area.Nevertheless, the Changping bowl, the Chaoyang bowl, and the Tongzhou bowl are connected to each other, and the sinking rate shows an accelerating tendency.By late 2011, an area of more than 3900 km 2 , which is about 60% of the Beijing Plain, has been affected by land subsidence, with a cumulative subsidence over 100 mm, while the maximum accumulative displacement in this area reached 1302 mm [40].

InSAR Datasets
Envisat is an Earth observation satellite launched by ESA (European Space Agency) on 1 March 2002, with a revisit frequency of 35 days and orbit altitude of 800 km.ASAR is one of the sensors carried on Envisat and measures the radar backscatter of the Earth's surface at C-band.In this study, a stack of 47 Envisat ASAR images acquired between June 2003 and August 2010 from descending track was provided by the ESA.The maximum temporal baseline is 1960 days, and the perpendicular baseline ranges from 39.8 m to 1011.8 m (Figure 1b).The SAR image stack covers an area of approximately 10,000 km 2 (green box in Figure 1a), including most of the Beijing Plain except the Daxing County.
Another stack of 33 TSX/TDX SAR images recorded at 100-200 MHz was acquired between April 2010 and November 2013 from ascending track and was provided by the German Aerospace Center (German: Deutsches Zentrum für Luft-und Raumfahrt e.V., abbreviated as DLR).The maximum temporal baseline is 682 days, and the perpendicular baseline ranges from 3.1 m to 509.0 m (Figure 1b).The stack covers an area of approximately 1500 km 2 (red box in Figure 1a), part of the eastern Beijing Plain.The stacks were single-polarization HH mode images recorded in stripmap mode with a resolution of up to 3 m.
TerraSAR-X, launched on 15 June 2007, is the first in the world to provide services using X-band microwave sensors with a revisit frequency of 11 days and orbit altitude of 514 km.As a 1-m resolution class SAR satellite, TerraSAR-X delivers Earth observation data for scientific, institutional and commercial users worldwide.TanDEM-X (TerraSAR-X add-on for Digital Elevation Measurement) is the name of the twin satellite of TerraSAR-X.It is a second, almost identical spacecraft to TerraSAR-X.The new mission is to generate a global digital elevation model (DEM) with 12 m resolution and a vertical accuracy better than 2 m.

Hydraulic Head Level Data
Hydraulic head levels from observation wells provide direct measurements of the fluid pressure at depth [17].The groundwater database maintained by CGS provides observation well locations, measurements, and other information on well status.Annual Beijing Water Resources Bulletin published by Beijing Water Authority (BWA) also provides auxiliary information about phreatic water.In this study, more than 30 wells selected for time series analysis are scattered around the plain and not necessarily located near zones of maximum subsidence or uplift.In Section 4.3, we examine the head level data in conjunction with vertical surface deformation data obtained from InSAR time series.

InSAR Time Series Analysis
To measure the time series surface deformation in Beijing Plain, we implement the Wavelet Based InSAR (WabInSAR) approach [41].WabInSAR uses multiple-master interferograms within a certain perpendicular baseline, similar to SBAS approach [24].We firstly coregister the SAR images to the same master image within each stack.Then two sets of interferograms are generated with respect to the pre-defined perpendicular and temporal baseline thresholds (Figure 2c,d).The 90 m Shuttle Radar Topography Mission (SRTM) DEM is used to remove the flat earth effect and the phase due to topography.WabInSAR implements a statistical approach to identify elite (i.e., less noisy) pixels by examining complex interferometric phase noise of the dataset.Then, the algorithm performs a variety of wavelet-based filters to reduce the effects of topography and the orbit errors and correlated atmospheric delay [42].An iterative approach, incorporating a minimum cost flow approach used for phase unwrapping followed by Legendre polynomial wavelets filter used for estimating the high-frequency nuisances, is used for phase unwrapping.Following the inversion of the unwrapped interferograms and the generation of the time series of the surface deformation, the effect of temporally uncorrelated atmospheric delay is then removed by using a continuous wavelet transform.Through a reweighted least square approach, WabInSAR inverts the interferometric dataset and generates a uniform time series of the line of sight (LOS) surface deformation and uses these values to fit a linear velocity.
In this study, we apply WabInSAR technique to process the two stacks of the Envisat and TSX/TDX satellites for InSAR time series generation.The vertical displacement is considered to be more appropriate for comparison with ground leveling measurements.Thus, the measurements in line of sight (LOS) are directly projected into the vertical direction assuming that the detected movements are mostly in vertical direction during the period spanned by SAR image acquisition in this study.The assumption is also supported by previous works, in which GPS [43] and leveling [44] investigations show that this area presents a low relative horizontal movement (1.57-1.93mm/year).Then, the InSAR time series estimates are validated by independent ground leveling data.

Wavelet Transform for Time-Series Analysis
Signal decomposition algorithms, such as principal component analysis (PCA), the Fourier transform (FT), have been used to separate temporal component from geodetic time series [11,[45][46][47].Rudolf et al. [48] used PCA to extract the dominant temporal behavior from InSAR time series, while Chaussard et al. [49] performed PCA to extract the spatiotemporal variability of deformation and to investigate the seasonal deformation and water level changes.Despite its advantages for signal decomposition, as previous study shows, PCA is not suitable for dealing with signals with low amplitude [50].The FT can identify different frequencies through decomposing a sequence of deformation data in the time domain, but it gives no information regarding where those components appear in timeline [51].
To overcome these limitations, we apply complex continuous wavelet transform (CWT) which provides the time-frequency representation in different scale.Wavelet analysis has been widely used in geophysical time series analysis, in particular for SAR interferometry and further processing [52].This includes SAR filters, orbital error removal, topography-correlated atmospheric delay reduction, interferogram inversion and time series generation, a new approach for InSAR time series analysis [41], and analyzing the spatiotemporal characterization of land subsidence and the correlation between vertical displacement and hydraulic head level time series in Phoenix [17].
In this study, the complex Morlet wavelet is chosen to serve as mother wavelet in the procedure.The wavelet transform analysis is performed by using the free Matlab package created by Grinsted et al. [49] (refer to: [53]).In mathematics, the Morlet wavelet is a wavelet composed of a complex exponential multiplied by a Gaussian window, which is closely related to human perception.In this case, The Morlet wavelet transform analysis is widely used in signal analysis, heart disease diagnosis, and so on.For a more detailed description of the Morlet wavelet, refer to: [54].Before performing CWT, preparatory work should be done to ensure the InSAR time series has the same interval.Besides, the hydraulic head levels are also interpolated in this study to consistent with the dense InSAR time series.Then, we generate the wavelet power spectrum of the interpolated InSAR and the synchronous hydraulic head level time series.Finally, we perform cross wavelet transform (XWT) and wavelet transform coherence (WTC) to investigate the relationship between InSAR and hydraulic head level time series.Nevertheless, XWT and WTC will expose their common power and relative phase in time-frequency space [47].
Remote Sens. 2018, 10, x FOR PEER REVIEW 6 of 18 This includes SAR filters, orbital error removal, topography-correlated atmospheric delay reduction, interferogram inversion and time series generation, a new approach for InSAR time series analysis [41], and analyzing the spatiotemporal characterization of land subsidence and the correlation between vertical displacement and hydraulic head level time series in Phoenix [17].In this study, the complex Morlet wavelet is chosen to serve as mother wavelet in the procedure.The wavelet transform analysis is performed by using the free Matlab package created by Grinsted et al. [49] (refer to: [53]).In mathematics, the Morlet wavelet is a wavelet composed of a complex exponential multiplied by a Gaussian window, which is closely related to human perception.In this case, The Morlet wavelet transform analysis is widely used in signal analysis, heart disease diagnosis, and so on.For a more detailed description of the Morlet wavelet, refer to: [54].Before performing CWT, preparatory work should be done to ensure the InSAR time series has the same interval.Besides, the hydraulic head levels are also interpolated in this study to consistent with the dense InSAR time series.Then, we generate the wavelet power spectrum of the interpolated InSAR and the synchronous hydraulic head level time series.Finally, we perform cross wavelet transform (XWT) and wavelet transform coherence (WTC) to investigate the relationship between InSAR and hydraulic head level time series.Nevertheless, XWT and WTC will expose their common power and relative phase in time-frequency space [47].

InSAR Measurements in Two Time Periods
We generate time series for the Envisat ASAR (June 2003-August 2010) stack and TerraSAR-X (April 2010-November 2013) stack using the method described in Section 3.1.Interferogram pairs with perpendicular, temporal, and Doppler separations less than the threshold values were selected for processing (Figure 2c,d).For Envisat ASAR dataset (Figure 2c), the maximum perpendicular absolute baseline difference (500 m), maximum absolute temporal difference (400 days), and coherence threshold (0.5) were applied according to the data availability and the environmental settings in the study site.For TSX/TDM dataset (Figure 2d), the threshold values were 500 m, 300 days, and 0.5, respectively.
Interferograms of TerraSAR-X SAR stack cannot cover the entire study area but provide good coverage of the sinking area in Beijing Plain. Figure 2a,b shows the linear land subsidence velocities (cm/year) from two different time periods.Negative values indicate downward movement (subsidence, while positive values indicate uplift.The black circle with label indicates the location for each ground leveling benchmark used in this study.The black dashed circles in Figure 2b highlights areas with significant unwrapping errors due to sparse elite points.Overall, displacement rate ranges from −12.86 cm/year subsidence to 0.92 cm/year uplift with a standard variance of 69.8 mm/year relative to the reference point (Figure 1a), between June 2003 and November 2013.
During June 2003-August 2010, as can be seen from Figure 2a, significant land subsidence occurs in the eastern Beijing Plain, especially in Changping, Chaoyang, Tongzhou, Shunyi and Pinggu County.Regarding to the displacement rate during April 2010-November 2013, both of sinking rate and range increase are observed near Tongzhou and Chaoyang County.As a result, the maximum subsidence rate of 12.86 cm/year is observed in Chaoyang funnel.
We plot the deformation profiles through time for every red dashed line feature in two different periods (Figures 3 and 4).As can be seen in the figures, profiles across Pinggu funnel, Chaoyang funnel and Changping funnel, labeled PG, CY-1, CY-2 and CP, respectively, are characterized by steady subsidence.Nevertheless, the CY-1 profile indicates a small zone of uplift to the north and east, while TZ-1/2 profiles crossing Tongzhou funnel feature more cumulative subsidence, and indicate the Tongzhou funnel is made up of independent subsidence bowls.The SY profile crossing Shunyi funnel has significant difference in two different time periods.As mentioned in Section 3.1, for individual interferograms, we identify pixels by examining phase noise to find elite pixels, and then use an iterative MCF approach to estimate phase value of the network of every pixel.This strategy would miscalculate cycles of phase unwrapping, such as, if an isolated elite point is too far away from its neighbor.Unfortunately, we note that, in Shunyi funnel, there are several areas with significant unwrapping errors due to sparse elite points during the SAR acquisitions period (2010-2013).

Comparison with Leveling Measurements
Interpolation in the time domain was used for quantitative comparison of the cumulative displacement time-series measured from both datasets.The SBAS-measured displacements were interpolated to the dates of the ground leveling surveys.To make the results from two measurements are comparable, we choose benchmark No. 139 (marked by filled orange circle with black border) close to the InSAR reference point as the leveling reference point.Then, the average displacement of pixels lying <300 m from each benchmark (marked by black circles in Figure 2b) was obtained as the corresponding InSAR measurement, assuming that the surface does not change significantly over this distance.
Then, a scatter plot generated from interpolated InSAR-measured time-series and ground leveling surveys at the 40 benchmarks is shown in Figure 5a.The differences between the two measurements vary from −10.8 mm/year to 7.6 mm/year, and the standard deviation (SD) is less than 6 mm/year, with a correlation coefficient (R 2 ) of 0.82. Figure 5b shows that the time-series measurements from both techniques match reasonably well; however, the temporal sampling rate of the leveling surveys (annual) is much lower than that of the SAR acquisitions (almost monthly).Nevertheless, there are strong fluctuations detected by InSAR between 2004 and 2009 at benchmark BJ046 and BJ163.Similar subsidence trends are observed at the benchmark BJ014, BJ064, BJ111, and BJ167, all of which are located at the edges of the sinking funnels (Figure 2).While the points located in the center of the funnels, such as the benchmark BJ094, shows continually sinking trend.These features are not obvious in ground leveling surveying due to the poor temporal sampling.
Remote Sens. 2018, 10, x FOR PEER REVIEW 8 of 18 SAR acquisitions (almost monthly).Nevertheless, there are strong fluctuations detected by InSAR between 2004 and 2009 at benchmark BJ046 and BJ163.Similar subsidence trends are observed at the benchmark BJ014, BJ064, BJ111, and BJ167, all of which are located at the edges of the sinking funnels (Figure 2).While the points located in the center of the funnels, such as the benchmark BJ094, shows continually sinking trend.These features are not obvious in ground leveling surveying due to the poor temporal sampling.

Temporal Evolution of Land Subsidence and Groundwater
As mentioned in Section 2.1, a large sinking bowl (orange-red areas in Figure 6a) with accumulated subsidence of more than 50 cm (from 2003 to 2010) were formed since the Changping funnel, the Chaoyang funnel, and the Tongzhou funnel have connected each other.A case study in the northern Beijing Plain showed that the subsidence in Beijing is mainly caused by intense groundwater extraction, apart from the hydrogeological conditions; however, the uneven distribution of the land subsidence has been strongly affected by the presence of compressible

Temporal Evolution of Land Subsidence and Groundwater
As mentioned in Section 2.1, a large sinking bowl (orange-red areas in Figure 6a) with accumulated subsidence of more than 50 cm (from 2003 to 2010) were formed since the Changping funnel, the Chaoyang funnel, and the Tongzhou funnel have connected each other.A case study in the northern Beijing Plain showed that the subsidence in Beijing is mainly caused by intense groundwater extraction, apart from the hydrogeological conditions; however, the uneven distribution of the land subsidence has been strongly affected by the presence of compressible

Temporal Evolution of Land Subsidence and Groundwater
As mentioned in Section 2.1, a large sinking bowl (orange-red areas in Figure 6a) with accumulated subsidence of more than 50 cm (from 2003 to 2010) were formed since the Changping funnel, the Chaoyang funnel, and the Tongzhou funnel have connected each other.A case study in the northern Beijing Plain showed that the subsidence in Beijing is mainly caused by intense groundwater extraction, apart from the hydrogeological conditions; however, the uneven distribution of the land subsidence has been strongly affected by the presence of compressible deposits [31].As illustrated in Figure 6a, the distribution of subsidence matches that of the sinking groundwater table; areas with high subsidence rates are always located in areas with a low hydraulic head level (such as groundwater depression cone).Nevertheless, the spatial distributions of the accumulated subsidence and the groundwater depression cones are not fully aligned with each other, regarding to the profile A-A' (Figure 6b).
In Beijing, deep water has been extracted at high rates for industrial and agricultural uses.However, the reported groundwater extraction volumes cannot be used because of the large number of unregistered wells.Thus, we plot the monitoring results from two selected wells (marked by red diamonds in Figure 6a).Well 2 is located at the edge of Tianzhu-tongzhou funnel while Well 1 is located near the center of the funnel.As can be seen in Figure 6b, the time-series hydraulic head level presents certain seasonal characteristics, especially for Well 2. For the location of Well 1, the measured subsidence increases along with the hydraulic head level.
In north China, agricultural irrigation is the main water consumption in winter and spring.Since the Wheat-Maize double cropping system plays a very important role in agriculture production in North China Plain, irrigation of winter wheat in Beijing Plain is implemented according to the solar term.In this case, the water level at Well 2 located near farmland presents strong seasonal characteristics.The groundwater table declines year by year, simply because the water increasing in summer cannot afford the decrease in winter.Nevertheless, InSAR time-series measurements at different locations are not strongly affected by the seasonal change of the groundwater.The inner periodicity in time domain and the relationship between the subsidence and groundwater level are analyzed in detail based on wavelet transform in the next section.
Remote Sens. 2018, 10, x FOR PEER REVIEW 10 of 18 deposits [31].As illustrated in Figure 6a, the distribution of subsidence matches that of the sinking groundwater table; areas with high subsidence rates are always located in areas with a low hydraulic head level (such as groundwater depression cone).Nevertheless, the spatial distributions of the accumulated subsidence and the groundwater depression cones are not fully aligned with each other, regarding to the profile A-A' (Figure 6b).In Beijing, deep water has been extracted at high rates for industrial and agricultural uses.However, the reported groundwater extraction volumes cannot be used because of the large number of unregistered wells.Thus, we plot the monitoring results from two selected wells (marked by red diamonds in Figure 6a).Well 2 is located at the edge of Tianzhu-tongzhou funnel while Well 1 is located near the center of the funnel.As can be seen in Figure 6b, the time-series hydraulic head level presents certain seasonal characteristics, especially for Well 2. For the location of Well 1, the measured subsidence increases along with the hydraulic head level.
In north China, agricultural irrigation is the main water consumption in winter and spring.Since the Wheat-Maize double cropping system plays a very important role in agriculture production in North China Plain, irrigation of winter wheat in Beijing Plain is implemented according to the solar term.In this case, the water level at Well 2 located near farmland presents strong seasonal characteristics.The groundwater table declines year by year, simply because the water increasing in summer cannot afford the decrease in winter.Nevertheless, InSAR time-series measurements at different locations are not strongly affected by the seasonal change of the groundwater.The inner periodicity in time domain and the relationship between the subsidence and groundwater level are analyzed in detail based on wavelet transform in the next section.

Wavelet Transform Results
We implement CWT analysis on hydraulic head level and InSAR time series for two wells to identify their temporal characteristics, respectively.The CWT analysis of the hydraulic head level and InSAR time series at two wells are shown in Figures 7 and 8, respectively.The peaks and troughs are denoted in blue and red in the contour map. Figure 7a-d shows that there are two main periods

Wavelet Transform Results
We implement CWT analysis on hydraulic head level and InSAR time series for two wells to identify their temporal characteristics, respectively.The CWT analysis of the hydraulic head level and InSAR time series at two wells are shown in Figures 7 and 8, respectively.The peaks and troughs are denoted in blue and red in the contour map. Figure 7a-d shows that there are two main periods under different scales from January 2005 to December 2011 for HHL time series at Well 1.As can be seen in Figure 7a, there are about 2.5 circles around scale 50 and >3 circles withn scale 30-40 from January 2005 to December 2011.Nevertheless, the high coefficient modulus of the time series around scale 50 (Figure 7b) indicates a more significant periodicity.To determine the scales, we plot the wavelet variance in Figure 7c.It can be clearly seen that there are two peaks appearing at scale 37 and 51.Then, we plot the whole signal of CWT coefficient in Figure 7d.Consequently, there are 2.5 cycles indicating a period of about 33.6 months at scale 51 and 3.5 cycles with a period of about 24 months at scale 37. Figure 7e-h shows there is no constant period for InSAR time series at Well 1.The time interval between adjacent peak and trough increases along with scale, as can be seen from Figure 7e.Besides, taking coefficient curve in scale 60 as example, the amplitude shows significant increase over the same time interval (Figure 7h).
Figure 8a-d shows that there is one main period at scale 18 from January 2005 to December 2011 for HHL time series at Well 2. It can be seen in Figure 8d that there are seven cycles with a period of about 12 months at scale 18 and are 14 cycles with a period of about six months.The higher coefficient at scale 18 (Figure 8b) indicates that the HHL time series at Well 2 have stronger annual variation characteristics apart from the feature of seasonal variation.Figure 8e-h shows that there are no permanent periods for InSAR time series at Well 2 from June 2003 to November 2013.As can be seen, there are 5.5 cycles with a not strong period (because of the low coefficient in Figure 8f) of about 15.6 months at scale 33.Nevertheless, the amplitudes show significant increase after 2011 at scales 33 and 59 in Figure 8h.As can be seen in Figure 7a, there are about 2.5 circles around scale 50 and >3 circles withn scale 30-40 from January 2005 to December 2011.Nevertheless, the high coefficient modulus of the time series around scale 50 (Figure 7b) indicates a more significant periodicity.To determine the scales, we plot the wavelet variance in Figure 7c.It can be clearly seen that there are two peaks appearing at scale 37 and 51.Then, we plot the whole signal of CWT coefficient in Figure 7d.Consequently, there are 2.5 cycles indicating a period of about 33.6 months at scale 51 and 3.5 cycles with a period of about 24 months at scale 37. Figure 7e-h shows there is no constant period for InSAR time series at Well 1.The time interval between adjacent peak and trough increases along with scale, as can be seen from Figure 7e.Besides, taking coefficient curve in scale 60 as example, the amplitude shows significant increase over the same time interval (Figure 7h). Figure 8a-d shows that there is one main period at scale 18 from January 2005 to December 2011 for HHL time series at Well 2. It can be seen in Figure 8d that there are seven cycles with a period of about 12 months at scale 18 and are 14 cycles with a period of about six months.The higher coefficient at scale 18 (Figure 8b) indicates that the HHL time series at Well 2 have stronger annual variation characteristics apart from the feature of seasonal variation.Figure 8e-h shows that there are no permanent periods for InSAR time series at Well 2 from June 2003 to November 2013.As can be seen, there are 5.5 cycles with a not strong period (because of the low coefficient in Figure 8f) of about 15.6 months at scale 33.Nevertheless, the amplitudes show significant increase after 2011 at scales 33 and 59 in Figure 8h.

XWT and WTC on Hydraulic Head Level and InSAR Time Series
The connection between the portrayed patterns is quite difficult to see, therefore we use cross wavelet transform and wavelet transform coherence.Cross wavelet and wavelet coherence analysis are powerful techniques for testing proposed linkages between two time-series.To untangle the relation between the hydraulic head level and land subsidence in Beijing Plain, we apply XWT and CWT on hydraulic head level and accumulated subsidence (AS) time series from Wells 1 and 2. The wavelet power spectrum (WPS) shows the distribution of a frequency component with respect to the time, as illustrated in Figure 9.The thick contour is the 5% significance level against red noise, while the thinner black line indicates the cone of influence (COI).Areas outside the COI are shown as lighter shadows and are not reliable in signal analysis.
The continuous wavelet power spectrums of the HHL and AS from Well 1 are shown in Figure 9a,b, respectively.The HHL time series has high power in the 16-20-month band in the period from 2005 to 2008, although the power is under the 5% significance level.The CWT power of the HHL time series from Well 2 shown in Figure 9e has 8-14-month band in the period from 2006 to 2011.However, for the AS time series from two wells, the powers are quite low and therefore there are no high power above the 5% significance level in the whole period (Figure 9b,f).
The XWT power of the HHL and SA for the selected wells are shown in Figure 9c,g, respectively.The relative phase relationship vector is shown as arrows (with in-phase pointing right, anti-phase

XWT and WTC on Hydraulic Head Level and InSAR Time Series
The connection between the portrayed patterns is quite difficult to see, therefore we use cross wavelet transform and wavelet transform coherence.Cross wavelet and wavelet coherence analysis are powerful techniques for testing proposed linkages between two time-series.To untangle the relation between the hydraulic head level and land subsidence in Beijing Plain, we apply XWT and CWT on hydraulic head level and accumulated subsidence (AS) time series from Wells 1 and 2. The wavelet power spectrum (WPS) shows the distribution of a frequency component with respect to the time, as illustrated in Figure 9.The thick contour is the 5% significance level against red noise, while the thinner black line indicates the cone of influence (COI).Areas outside the COI are shown as lighter shadows and are not reliable in signal analysis.
The continuous wavelet power spectrums of the HHL and AS from Well 1 are shown in Figure 9a,b, respectively.The HHL time series has high power in the 16-20-month band in the period from 2005 to 2008, although the power is under the 5% significance level.The CWT power of the HHL time series from Well 2 shown in Figure 9e has 8-14-month band in the period from 2006 to 2011.However, for the AS time series from two wells, the powers are quite low and therefore there are no high power above the 5% significance level in the whole period (Figure 9b,f).
The XWT power of the HHL and SA for the selected wells are shown in Figure 9c,g, respectively.The relative phase relationship vector is shown as arrows (with in-phase pointing right, anti-phase pointing left, and HHL leading AS by 90 • pointing straight down), indicating the phase difference between the HHL and AS.Here, we notice that there is a relatively high common power in the 8-16-month band from 2006 to 2011 standing out above the 5% significance level (Figure 9g).We therefore speculate that there is another link between HHL and AS than that implied by the XWT power.
WTC power can reveal areas with high common power and can describe the coherence of the time series in time frequency space.The statistical significance level of the WTC is estimated by Monte Carlo methods.The WTC of the HHL and SA for the selected wells are shown in Figure 9d,h, respectively.Compared with the XWT power, there are more areas with higher power standing out above the 5% significance level.Nevertheless, these areas show different phase relationship between HHL and AS in different periods.Figure 9d shows the HHL and AS are in-phase with significant common power in the ~2 band from 2007 to 2009, but are anti-phase between June 2010 and November 2010.A similar situation is found in Figure 9h from the time series of Well 2. There are areas with significant coherence between the HLL and AS time series in ~8-month band in several periods, although the cross-phase angle ranges from −84 pointing left, and HHL leading AS by 90° pointing straight down), indicating the phase difference between the HHL and AS.Here, we notice that there is a relatively high common power in the 8-16month band from 2006 to 2011 standing out above the 5% significance level (Figure 9g).We therefore speculate that there is another link between HHL and AS than that implied by the XWT power.WTC power can reveal areas with high common power and can describe the coherence of the time series in time frequency space.The statistical significance level of the WTC is estimated by Monte Carlo methods.The WTC of the HHL and SA for the selected wells are shown in Figure 9d,h, respectively.Compared with the XWT power, there are more areas with higher power standing out above the 5% significance level.Nevertheless, these areas show different phase relationship between HHL and AS in different periods.Figure 9d shows the HHL and AS are in-phase with significant common power in the ~2 band from 2007 to 2009, but are anti-phase between June 2010 and November 2010.A similar situation is found in Figure 9h from the time series of Well 2. There are areas with significant coherence between the HLL and AS time series in ~8-month band in several periods, although the cross-phase angle ranges from −84° ± 5° to 164° ± 17°.

Other Factors Concerned with Land Subsidence in Beijing
As mentioned in Section 2.1, Beijing Plain is a typical piedmont pluvial-alluvial plain composed of several alluvial-pluvial fans and an alluvial plain area.The aquifer on the top of the alluvial fans consists of cobbles and gravels; the middle part is full of sand gravel strata, while in the fringe part and alluvial plain area it gradually converted to the multi-layer structure of the coarse sand, medium sand, and fine sand.In this case, land subsidence in Beijing Plain is increasingly serious and its causes are complicated, due to its complex geological conditions and growing anthropogenic intervention.
In Beijing Plain, the scope, range and rate of land subsidence have the direct relation with overexploitation of groundwater.The poroelastic response to groundwater withdrawal is detected as surface deformation due to elastic recoverable subsidence or inelastic, permanent compaction [17].Then, inelastic, permanent subsidence occurs where the aquifer system presents a continuous compaction associated with the groundwater level decrease.In this case, the land subsidence in these areas present an increasing trend in the rate.Apart from the Groundwater exploitation, previous works from our group indicate that the distribution and development trend of the local land subsidence in eastern plain area are controlled by geological structures [19].As can be seen in Figure 1a, there are several active high angle normal faults in the Beijing plain.Tectonic movements cause tectonic units to decline slowly under the regional stress field, while the relative displacement of the hanging walls and the footwalls results in differences in the land subsidence patterns on both sides of the fault.However, in this case, obvious horizontal deformation may be formed under the brute force of active faults, and the issue should be further analyzed through offset-tracking analysis [55].Besides, land subsidence also have connection the quaternary compressible layers [56] and LU/LC [57].Nevertheless, by creating an index-based built-up index (IBI) from landsat TM image, the time-varying static load field has been generated according to the period of the acquisitions used for InSAR measurements in [35].The result showed that there is a positive correlation between building density and land subsidence, especially for the areas with high density buildings or extensive transportation networks.

Conclusions
Land subsidence in Beijing Plain has connection with hydraulic head level falling caused by over-exploitation of groundwater.We implement the Wavelet Based InSAR approach on two SAR image stacks to investigate the long-term displacement in eastern Beijing Plain.The InSAR-measured sinking rates and accumulated displacements are in good agreement with results estimated from ground leveling surveys.Then, through the distribution of subsidence areas and groundwater funnel, the InSAR based time series and the monitoring well based groundwater level changes, the spatial correlations and responses between land subsidence field and groundwater flow field are analyzed.The results show that the distribution of the subsidence centers in the northern Beijing Plain is consistent with that of the groundwater drop funnels, with a similar downward trend over the whole observation time.Continuous wavelet transform provides a powerful tool for signal decomposition and extracting periodic components of a time series.According to the analysis of well based results located in different areas, the long-term groundwater exploitation in the northern subsidence area has led to the continuous decline of the water level, resulting in the inelastic and permanent compaction, while, for the monitoring wells located outside the subsidence area, the subsidence time series show obvious elastic deformation characteristics (seasonal characteristics) as the groundwater level changes.According to the wavelet transformation, the land subsidence time series at monitoring well site lags several months behind the groundwater level change.

Figure 1 .
Figure 1.(a) The location of Beijing Plain and the coverage of the SAR sensors.The light gray dashed lines are the main faults published by CGS, and the labels indicate the name of the active faults mentioned in this study.The red circle with cross is the reference point.The red circle is the location of Houshayu, while the black bold plane represents the Beijing Capital International Airport.(b) Temporal-spatial baseline distributions of Envisat ASAR and TerraSAR-X image stacks in this study.TSX/TDM is the abbreviation of TanDEM-X/TerraSAR-X.

Figure 1 .
Figure 1.(a) The location of Beijing Plain and the coverage of the SAR sensors.The light gray dashed lines are the main faults published by CGS, and the labels indicate the name of the active faults mentioned in this study.The red circle with cross is the reference point.The red circle is the location of Houshayu, while the black bold plane represents the Beijing Capital International Airport.(b) Temporal-spatial baseline distributions of Envisat ASAR and TerraSAR-X image stacks in this study.TSX/TDM is the abbreviation of TanDEM-X/TerraSAR-X.

Figure 2 .
Figure 2. Line of sight velocity maps and connections from two satellite tracks: (a) Envisat ASAR descending track (June 2003-August 2010) with maximum subsidence of 10.06 cm/year and maximum uplift of 0.64 cm/year.(b) TSX/TDM ascending track (April 2010-November 2013) with maximum subsidence of 12.86 cm/year and maximum uplift of 0.92 cm/year.InSAR time series profiles for two periods are marked by red dashed line and are detailed in Figures 3 and 4, respectively.(c) Connection graph generation from Envisat ASAR interferogram pairs.(d) Connection graph generation from TSX/TDM interferogram pairs.

Figure 2 .
Figure 2. Line of sight velocity maps and connections from two satellite tracks: (a) Envisat ASAR descending track (June 2003-August 2010) with maximum subsidence of 10.06 cm/year and maximum uplift of 0.64 cm/year.(b) TSX/TDM ascending track (April 2010-November 2013) with maximum subsidence of 12.86 cm/year and maximum uplift of 0.92 cm/year.InSAR time series profiles for two periods are marked by red dashed line and are detailed in Figures 3 and 4, respectively.(c) Connection graph generation from Envisat ASAR interferogram pairs.(d) Connection graph generation from TSX/TDM interferogram pairs.

Figure 3 .
Figure 3. Profiles of InSAR vertical displacement time series (June 2003-August 2010) marked by red dashed line in Figure 2a.The gaps occur where noise covers up the pixels throughout the interferograms, and the pixels are excluded as described in Section 3.1.

Figure 3 .
Figure 3. Profiles of InSAR vertical displacement time series (June 2003-August 2010) marked by red dashed line in Figure 2a.The gaps occur where noise covers up the pixels throughout the interferograms, and the pixels are excluded as described in Section 3.1.

Figure 4 .
Figure 4. Profiles of InSAR vertical displacement time series (April 2010-November 2013) marked by red dashed line in Figure 2b.The gaps occur where noise is statistically significant throughout the interferograms, as described in Section 3.1.

Figure 5 .
Figure 5.Comparison of InSAR and ground leveling-derived velocities during 2003-2013: (a) the black line represents the one-to-one function; and (b) comparison between the time series subsidence measured by InSAR and by ground leveling surveys at the benchmarks "BJ046", "BJ163", and "BJ094".

Figure 4 .
Figure 4. Profiles of InSAR vertical displacement time series (April 2010-November 2013) marked by red dashed line in Figure 2b.The gaps occur where noise is statistically significant throughout the interferograms, as described in Section 3.1.

Figure 4 .
Figure 4. Profiles of InSAR vertical displacement time series (April 2010-November 2013) marked by red dashed line in Figure 2b.The gaps occur where noise is statistically significant throughout the interferograms, as described in Section 3.1.

Figure 5 .
Figure 5.Comparison of InSAR and ground leveling-derived velocities during 2003-2013: (a) the black line represents the one-to-one function; and (b) comparison between the time series subsidence measured by InSAR and by ground leveling surveys at the benchmarks "BJ046", "BJ163", and "BJ094".

Figure 5 .
Figure 5.Comparison of InSAR and ground leveling-derived velocities during 2003-2013: (a) the black line represents the one-to-one function; and (b) comparison between the time series subsidence measured by InSAR and by ground leveling surveys at the benchmarks "BJ046", "BJ163", and "BJ094".

Figure 6 .
Figure 6.Spatial and temporal characteristics of land subsidence and hydraulic head levels: (a) InSARmeasured accumulated subsidence between June 2003 and August 2010 with the average hydraulic head level during the same period; and (b) monitoring results for confined hydraulic head level and InSAR-measured subsidence for selected wells.The illustration at bottom details the profile A-A' in Figure 6a.

Figure 6 .
Figure 6.Spatial and temporal characteristics of land subsidence and hydraulic head levels: (a) InSAR-measured accumulated subsidence between June 2003 and August 2010 with the average hydraulic head level during the same period; and (b) monitoring results for confined hydraulic head level and InSAR-measured subsidence for selected wells.The illustration at bottom details the profile A-A' in Figure 6a.
Remote Sens. 2018, 10, x FOR PEER REVIEW 11 of 18 under different scales from January 2005 to December 2011 for HHL time series at Well 1.

Figure 7 .
Figure 7. CWT on HHL and InSAR time series for Well 1: (a) the contour lines of the real part of the CWT analysis for HHL time series; (b) the coefficient modulus of the HHL time series in different scales; (c) wavelet variance curve of HHL time series; and (d) the CWT coefficients curve of HHL time series in different scales.(e-h) The same as (a-d) but for InSAR time series.

Figure 7 .
Figure 7. CWT on HHL and InSAR time series for Well 1: (a) the contour lines of the real part of the CWT analysis for HHL time series; (b) the coefficient modulus of the HHL time series in different scales; (c) wavelet variance curve of HHL time series; and (d) the CWT coefficients curve of HHL time series in different scales.(e-h) The same as (a-d) but for InSAR time series.

Figure 8 .
Figure 8. CWT on HHL and InSAR time series for Well 2: (a) the contour of the real part of the CWT analysis for HHL time series; (b) the coefficient modulus of the HHL time series in different scales.(c) wavelet variance curve of HHL time series; and (d) the CWT coefficients curve of HHL time series in different scales.(e-h) The same as (a-d) but for InSAR time series.

Figure 8 .
Figure 8. CWT on HHL and InSAR time series for Well 2: (a) the contour of the real part of the CWT analysis for HHL time series; (b) the coefficient modulus of the HHL time series in different scales.(c) wavelet variance curve of HHL time series; and (d) the CWT coefficients curve of HHL time series in different scales.(e-h) The same as (a-d) but for InSAR time series.

Figure 9 .
Figure 9. (a) Wavelet power spectrum (Morlet wavelet) of the HHL at Well 1; (b) wavelet power spectrum (Morlet wavelet) of the AS at Well 1; (c) the cross-wavelet power between the HHL and AS at Well 1; and (d) the wavelet coherency and phase between HHL and AS at Well 1. (e-h) The same as (a-d) but for Well 2. Arrows indicate the phase difference between the HHL and AS.The thick black contour is the 5% significance level using the red noise model, while the thin black line indicates the cone of influence (COI).

Figure 9 .
Figure 9. (a) Wavelet power spectrum (Morlet wavelet) of the HHL at Well 1; (b) wavelet power spectrum (Morlet wavelet) of the AS at Well 1; (c) the cross-wavelet power between the HHL and AS at Well 1; and (d) the wavelet coherency and phase between HHL and AS at Well 1. (e-h) The same as (a-d) but for Well 2. Arrows indicate the phase difference between the HHL and AS.The thick black contour is the 5% significance level using the red noise model, while the thin black line indicates the cone of influence (COI).