Recent Surface Deformation in the Tianjin Area Revealed by Sentinel-1 A Data

In this study, we employed multitemporal InSAR (Interferometric Synthetic Aperture Radar) (MT-InSAR) to detect spatial and temporal ground deformations over the whole Tianjin region in the North China Plain area. Twenty-five ascending Sentinel-1A terrain observation by progressive scans (TOPS) synthetic aperture radar (SAR) scenes covering this area, acquired from 9 January 2016 to 8 June 2017, were processed using InSAR time series analysis. The deformation results derived from Sentinel-1A MT-InSAR were validated with continuously operating reference stations (CORS) at four sites and four stations of the Crustal Movement Observation Network of China (CMONOC). The overall results show good agreement, demonstrating the suitability of applying Doris with Sentinel-1A data to high-resolution monitoring of surface deformation. Significant deformation variations have been observed in different parts of Tianjin. These gradually increased from the central part of the metropolitan area to the nearby suburbs. The deformation rate of the main urban area is well-balanced and it is also relatively linear, with uplifting rates ranging from 0 to 20 mm/yr. However, due to the diversity of the geological conditions and anthropogenic activities, remarkable signs of subsidence were found in several parts of Tianjin. In particular, the south-western part of Wuqing District and western part of Beichen District showed subsidence rates of up to −136 mm/yr. We also found that, in addition to groundwater over-exploitation and lithological characteristics, additional factors also influence ground subsidence, including dynamic loads (e.g., railways), static loads (e.g., urban construction), and groundwater recharging.


Introduction
Ground subsidence has always been one of the most severe and widespread geological hazards in the cities of China [1].Over 70,000 km 2 in 17 provinces in China have been directly affected by subsidence [2,3].The potential consequences of land subsidence mainly include degradation of the aquifer systems and damage to the utility infrastructure, buildings, railroads, highways, and bridges [4].In 2010, Herrera et al. [5] selected the TerraSAR-X data acquired from July 2008 to September 2009 to study the deformation of the city of Murcia.Detailed discussions were given with respect to the main factors that control the subsidence mechanisms: the thickness of the compressible layer, the presence of pumping wells, and the water table variation.Local analysis of several buildings and infrastructures was presented.The results demonstrated that surface deformation variations were correlated with seasonal variations of groundwater table level.The rapid response of ground deformation to pore pressure changes is explained by the drainage capacity of the aquifer system, and particularly by the presence of sandy layers in shallow acquifers, that increase its permeability.Differential subsidence is caused when different types of foundations of buildings and infrastructures are jointed.With the rapid development of the economy, the amount of groundwater exploitation in cities has also risen sharply.As a result, a series of problems related to urban land surface subsidence have arisen.Surface subsidence is a slowly changing geological hazard that can cause a series of secondary chain of processes that are closely related to important factors such as human production, economic development, and urban development [6].The North China Plain (NCP), Yangtze River Delta region and the Fenwei basin in China are the three major areas in which serious surface subsidence has occurred over the past few decades.The NCP includes three subsidence centres, in Tianjin, Cangzhou and the north-eastern corner of Beijing.The land subsidence in Tianjin began in the 1920s, accompanying large-scale urbanization and industrialization in Tianjin.After 1949, land subsidence in Tianjin accelerated and the maximum value exceeded 3.1 m by 2005.This made it the most serious area of subsidence in China [1].As one of the most serious affected areas in Tianjin, the special location of Wangqingtuo Town cannot be uniformly managed in administrative divisions [7].
To control the rate of land subsidence and reduce its harmfulness effectively, subsidence monitoring is extremely important.This is also basic work that is indispensable for studying the mechanism and trend of development of land subsidence.Traditional methods of land subsidence measurement include levelling, bedrock measurement and stratified measurement.A bedrock standard is a special monitoring point buried in intact underground bedrock by drilling.It is the marker point from the underground bedrock to the surface and acts as a starting point or elevation control point for leveling (i.e., the reference point for stratified target monitoring).Stratified marks are a special monitoring points buried at different depths of the ground layer by drilling methods according to the nature of the soil layer.With the development of modern surveying and mapping technology, GPS surveying has also become an important method of subsidence monitoring.Although traditional levelling and GPS data are able to provide reliable observations of ground subsidence, they are limited by low spatial sampling density and high cost.It is hard to obtain subsidence information with consistent spatial distribution at large scale [6].
In order to overcome these problems, an advanced remote sensing technique, referred to as differential interferometric synthetic aperture radar (D-InSAR), was introduced as a way to monitor deformation over wide areas with centimetre to millimetre accuracy [8].Repeat-pass satellite InSAR is potentially a unique tool for low cost, precise generation of digital elevation models (DEM) and broad-coverage surface deformation monitoring [9].The D-InSAR technique uses two SAR scenes of the same area taken at different times and calculates an interferogram that shows the line-of-sight (LOS) ground-surface displacement (range change) between two time periods [8,10,11].It has been proven to be an effective remote sensing tool for mapping urban land subsidence and is widely used to monitor surface subsidence, landslides, volcano dynamics, earthquakes, snow changes, etc. [12][13][14][15].However, D-InSAR is susceptible to atmospheric disturbance, and to temporal and spatial decorrelation.These make it difficult to use for long time-interval surface monitoring with high precision [16][17][18].Temporal decorrelation makes D-InSAR measurements unfeasible over vegetated areas and where the electromagnetic profiles and/or positions of the scatterers change with time within the resolution cell.Geometric decorrelation limits the number of image pairs suitable for interferometric applications and prevents full exploitation of the data sets available.Atmospheric inhomogeneities create an atmospheric phase screen (APS) superimposed on each SAR image that can seriously compromise accurate deformation monitoring [18].
In summary, under the influence of decorrelation factors and atmospheric delay effects, the accuracy in deformation monitoring achieved with the traditional D-InSAR method is generally not high enough, ranging from the decimetre to centimetre level.However, the accuracy of deformation monitoring often needs to be higher, so D-InSAR has substantial limitations in urban infrastructure monitoring applications.To overcome temporal and spatial decorrelation and the atmospheric delay of D-InSAR and to improve the applicability of the D-InSAR technique, multi-temporal InSAR (MT-InSAR) was developed to extract and analyse high precision and high density measurement points using multi-view SAR images.In 2001, Ferretti et al. [9] proposed a kind of MT-InSAR, the 'permanent scatterers' InSAR (PS-InSAR) method-known as the first persistent scatterers interferometry (PSI) technique [19,20]-a powerful technique for wide-area deformation monitoring with millimetre accuracy [21].Reliable deformation and height correction information can be obtained by adopting stable scattering targets (buildings, structures, etc.) from a long temporal series of interferometric stacks.Later, Berardino et al. [22] proposed a small baseline subset InSAR (SBAS-InSAR) method to study low-resolution, large-scale deformation.SBAS-InSAR is another TS-InSAR method that effectively mitigates decorrelation phenomena by analysing distributed scatterers (DSs) with high coherence based on an appropriate combination of interferograms.These can reduce the geometrical and temporal decorrelation more than with PS-InSAR [22][23][24].
There has been some debate about the relative merits of the PS and SBAS methods.However, because they are optimized for different models of ground scattering, the two approaches are in fact complementary, at least in the usual case where a data set contains pixels with a range of scattering characteristics; the scattering characteristics of the terrain are often between these two models.Recently, both of these two advanced InSAR techniques have been widely used in monitoring urban ground deformation around the world.For instance, Perissin and Wang [25] used PS-InSAR to monitor subsidence and landslides in urban areas of China (e.g., Shanghai, Tianjin, and Badong) and near the Three Gorges Dam.Heleno et al. [26] used PS-InSAR to detect and measure ground subsidence in Lisbon, Portugal.Castellazzi et al. [27] applied the SBAS-InSAR method to monitor the spatiotemporal variations of displacements and fracturing from 2012 to 2014 in Mexico.Bai et al. [4] presented results based on SBAS-InSAR observations to monitor the ground deformation in the Wuhan region of China.
In 2008, Hooper et al. [28] presented a new algorithm (released software package 'StaMPS/MTI' (Stanford Method for Persistent Scatterer/Multi temporal InSAR)) that combined both the PS and SBAS approaches to maximize the spatial sampling of a useble signal and get a higher signal-to-noise ratio than with either approach alone.Improvement of the spatial sampling is important, not only because the resolution of any deformation signal is increased, but it also enables more reliable estimation of integer phase-cycle ambiguities present in the data (phase unwrapping).
In this paper, we present the first investigation of land deformation in the Tianjin area based on the Sentinel-1A terrain observation by progressive scans (TOPS) data with Delft object-oriented radar interferometric software (Doris ) at large scale.The StaMPS approach was applied to 24 ascending Sentinel-1A scenes acquired from 9 January 2016 to 8 June 2017 to retrieve and characterize spatiotemporal variations of ground deformation in the Tianjin region.Moreover, the corresponding validation, detailed illustration, and analysis provided were done using the results of the MT-InSAR and GPS time series.Finally, the potential factors influencing subsidence and analysis of the mechanisms were investigated from the viewpoint of geological tectonic movement, urban development and various external factors.

Study Area
Tianjin is located at the junction of three tectonic units, with the east-west trending Yanshan mountain structural high to the north, the north-easterly Taihang mountain structural high to the southwest and the North China rift-depression basin to the south and southeast.The city is also on the western shore of the Bohai Gulf [45].Tianjin is separated by the Baodi-Ninghe Fault into two parts, with the northern part on the Yanshan Folded Belt.The southern part is in the Bohai Bay Basin, which contains three tectonic divisions: the Jizhong Depression, Cangxian Uplift, and Huanghua Depression.The Jizhong Depression and Cangxian Uplift are separated by a Paleogene lacuna.The Cangxian Uplift and Huanghua Depression are separated by the Cangdong Fault [46].To control further deterioration of the land subsidence condition in Tianjin, during 1995, a special group named the 'Tianjin Office of Land Subsidence' was set up.To date, 2100 benchmarks (BMs), 2 bedrock marks, 14 layerwise marks, 503 groundwater monitoring wells, 11 continuously operating reference stations (CORS), and 35 GPS stations have been built to monitor land subsidence [47,48].
The study area of our research covers almost all districts, except for a few small parts of Tianjin.The overall study area was approximately 170 by 167 km (Figure 1).It covered Ninehe and Jinghai Counties, Jixian, Baodi, Wuqing, Beichen, Dongli, Jinnan, and Binhai new districts (i.e., Hangu, Tanggu, and Dagang).Figure 2 illustrates the tectonic framework of the study area.

Sentinel1-A Datasets
In order to investigate the latest land deformation over the whole of Tianjin, 25 ascending Sentinel-1A TOPS SAR scenes with VV polarisation from the European Space Agency (ESA) were processed.As we considered only the ascending orbit dataset, negative values could only show subsidence in flat plain sectors; this was not valid for hilly and mountainous sectors, where the geometry (slope and aspect) of the mountain and hill sides influences the LOS velocity values.Except for a small part in the northeastern corner which is mountainous area, the study area is a flat plain.Therefore, the deformation results are reliable in the area we discuss below.Data were acquired from 9 January 2016 to 8 June 2017.The Sentinel-1A TOPS mode data consist of series of bursts with mutual overlaps, which provides a wide swath with resolution of 5 × 20 m [49].To improve the solution efficiency, swath1 and swath2 were selected for each of the scenes.Detailed parameters of the Sentinel-1A data are given in Table 1.
The Shuttle Radar Topography Mission (SRTM) data product released in 2015 by the U.S. Geological Survey (USGS) with a spatial resolution of 1 arc-second (30 m) was used to remove the topographic phase [4,27] (https://lta.cr.usgs.gov/SRTM1Arc).Finally, time series for 2015 to 2018 from four continuously operating reference stations (CORS) and four stations of the Crustal Movement Observation Network of China (CMONOC) were utilised to evaluate the results derived from Sentinel-1A MT-InSAR data.

MT-InSAR Data Processing with the Application of StaMPS
In the traditional PS-InSAR technique [9], PS candidates (PSC) are selected based on thresholding pixel amplitude dispersion with time, and a temporal linear deformation model is assumed.For high signal-to-noise ratio (SNR) targets, this PS method and its variant [50,51] have been successfully used for picking bright PS, e.g., urban areas with much man-made construction.Without prior assumptions about its temporal nature, StaMPS takes spatial correlation of phase into consideration and succeeds in picking PSC with low SNR scatterers in the natural terrain [20,49].The StaMPS software package, which contains the PS, SBAS, and combined MT-InSAR algorithms, can be download from http://www.hi.is/~ahooper/stamps [52].The general processing procedures used for StaMPS-PS are briefly introduced in the following steps; a more detailed description can be found in [20,28,52].

Interferogram Formation
Firstly, Doris software was used to generate the interferograms needed [53,54].Generally, at least 25 interferograms are required based on the functional temporal model [28,55].However, according to Hooper et al. [28], if StaMPS is applied, the number of the required interferograms can be reduced to approximately 12.In this study, 25 scenes were applied to generate the interferograms (Figure 3).In order to minimize the temporal and spatial decorrelation, the scene acquired on 22 December 2016 was selected as the master image, and all of the remaining scenes were co-registered and resampled in relation to the master image.Next, the topographic phase was removed from the interferograms using the SRTM data.

PS Pixel Selection
Next, PS methods of the StaMPS software package have been used to process the images [20,28,46].Initially, StaMPS selected a subset of pixels that included almost all of the PS pixels based on an analysis of their amplitude dispersion.The amplitude dispersion index, D A , was defined by [9] as: where σ A is the standard deviation and µ A is the mean of a series of amplitude values.Generally, a threshold value of D A is used to select PSC.For different study areas, [9,28] adopted threshold values of 0.25 and 0.40, respectively.Because most of the area under investigation was covered by natural vegetation, the threshold value we adopted was 0.40.After selection of a subset of pixels as the initial PSC, StaMPS estimated the phase stability for each of them using phase analysis.The wrapped phase of the x-th point in the i-th interferogram, can be illustrated as follows [9,28,56]: where φ de f ,x,i is the phase change due to surface deformation of pixels in the direction of LOS, φ α,x,i is the phase of atmospheric delay between different passes, φ orb,x,i is the residual phase caused by satellite orbit error, φ ε,x,i is the phase due to look angle error (DEM error), n x,i is noise (e.g., thermal effect, coregistration error, azimuth phase center uncertainty, etc.).Based on the error characteristics of the variables in Equation ( 2), [28] defined a phase stability indicator γ x (3) to evaluate whether the pixel is a PS: where N is the number of interferograms, φ x,i is a wrapped estimate of the spatially correlated parts on the right side of Equation ( 2).φ u ε,x,i is the estimate of the uncorrelated part of φ ε,x,i .Once this solution has converged to the estimated stability of each pixel and stops iterating, the resulting PS pixels was selected.Typically, PS pixel densities of 5-10 px per km 2 is required to properly estimate and remove atmospheric heterogeneity.With the different threshold values of the amplitude dispersion we selected, each of the PS densities could be raised to 145 px per km 2 by StaMPS.

Three Dimensional Phase Unwrapping and LOS Deformation Rates Deriving
Phase unwrapping is the process of recovering data from only known modulo 2π radians phase to unambiguous phase values.For large displacement of this area, a full three-dimensional (3-D; two spatial dimensions, as with conventional InSAR, and one temporal) unwrapping algorithm is required to derive the deformation rates [57].Before unwrapping, the spatially uncorrelated term DEM error was subtracted.After unwrapping, there are still four error terms in Equation ( 2).The spatially correlated portion of these errors are assumed to be temporally uncorrelated.Thus, high-pass filtering the unwrapped data in the temporal dimension and low-pass filtering in the spatial dimensions was applied to estimate the spatially correlated error [57].Then least-squares inversion was applied to derive the LOS displacement time series for PS points [52], and all the PS points were calibrated to a stable CORS reference station located in the north-east corner of Tianjin (Figure 12).

Precision Evaluation
At present, there are two methods to assess the accuracy of the PSInSAR technology.One is represented by Politecnico di Milano, which uses the SNR to calculate the noise of the signal [58][59][60].Another method comes from the Delft University of Technology and establishes the statistical relationship of various phase components.The accuracy of the PS points is finally determined based on the error propagation principle [53,61].As shown in Equation (3), the PS points were identified based on coherence in the temporal dimension.Therefore, the first method was adopted for accuracy assessment.
For temporally correlated radar targets, noise is mainly caused by the uncorrelated noise of the PS points.Generally, the residual phases are considered to conform to a normal distribution.The relationship between the temporal coherence and the phase variance of PS points can be expressed as: Equation (4): where the standard deviation of the phase σ ω can be obtained by the temporal coherence value γ: where γ is the estimated value of γ x , which can be obtained from Equation (3).The velocity accuracy of the PS points can be illustrated as where λ is the C-band microwave wavelength used by the Sentinel satellite.T i is the i-th acquisition date of the images and T is the mean value of T i .The values mentioned in this section are all 'precision of inner coincidence', which are obtained based on the phase stability.This mainly reflects the dispersion of the phase for PS points in the time dimension.The distribution of the calculated temporal coherence γ for all of the PS pixels are illustrated in Figure 4 and the values of the statistics are summarized in Table 2.The relationship between the temporal coherence γ and the theoretical value of the velocity accuracy for PS points is also given as in Figure 5.In theory, we can find that the accuracy of the velocity σ v can reach up to about 2 mm/yr when the coherence reaches 0.7 for each PS point.Although this velocity accuracy is the theoretical limit, it can also explain the reliability of the derived velocity for most of the PS points to some extent.In order to access the accuracy of the derived velocity, the standard deviation of the velocity for PS points in the direction of LOS have been shown in Figure 6.It shows that the accuracy rate is very high in most areas, except for the southern area and a small part of Tangshan.This was mainly caused by phase unwrapping errors, especially in the edge region where 3-D phase unwrapping could not be well constrained.Results indicate that the minimum, maximum, and mean standard deviation of the derived velocity for PS points are 1.841 mm/yr, 12.124 mm/yr, and 4.761 mm/yr, respectively.Root-mean-square error reached 5.015 mm/yr.For 60.720 % of the PS points, the standard deviation was less than 5 mm/yr, and over 99.31 % of all PS points retain a relatively high precision with the error lower than 10 mm/yr.Furthermore, we also compared the results between the theoretical and calculated statistics of the standard deviation, as shown in Figure 7.For the theoretical value of the standard deviation, most of the PS points range from 0.5 mm/yr to 1.5 mm/yr, while the calculated standard deviations are mostly between 2.5 mm/yr and 4 mm/yr.These inner coincidence results suggest that the deformations derived by InSAR are reliable.Theoretical statistics of the PS points

Validation with GPS Time Series
To further verify the accuracy of the InSAR derived results, a comparative analysis was performed by comparing the results with the GPS time series.Before comparison, all of the GPS time series were projected along the LOS direction and all the PS pixels located within a circle with 50 m radius around the corresponding GPS station were selected.Time series of the two CORS stations (i.e., XQYY and CH01) and the corresponding PS points were all referenced to the JIXN CORS station, as can be seen from Figure 10.This time series is relatively stable and the velocity in the direction of LOS was only 1.9 mm/yr.The other GPS time series were also referenced to a stable IGS (International GPS Service) station.In addition, all of GPS time series and the InSAR results were referenced to the data of the master image (22 December 2016).Due to various random noises in the GPS time series, a Db3 wavelet was introduced to filter the original data [48], as shown in Figures 8 and 9.
The Daubechies extremal phase wavelets (dbN) refers to a particular family of wavelets.N is the order of the wavelet and refers to the number of vanishing moments.Basically, the higher the number of N, the smoother the wavelet (and longer the wavelet filter).A wavelet with N vanishing moments is orthogonal to polynomials of degree at most N.In this paper, the Db3 wavelet was used to filter the original GPS time series.Figure 9 indicates that the SNR of the time series improved after the filtering was implemented.After the above data processing, we compared the deformation time series of GPS and InSAR, as shown in Figure 10.Since time series of XQYY and CH01 are too short, it is not reliable for us to compare the results with a limited time series.Therefore, these two stations were removed for validation.A curvilinear regression analyzing method was also applied for quantitative assessment the accuracy of PS time series.The results were that root mean square errors for TJWQ, TJBD, JIXN and HETS were 4.35 mm, 4.69 mm, 7.27 mm and 5.03 mm, respectively.Generally, the GPS and InSAR-derived deformation time series may not be consistent with each other.Since there was no Sentinel-1B data in Tianjin area during the investigated time period, there is a long time interval for other scenes, and only Sentinel-1A data was used to derive the deformation velocity, which may have led to three-dimensional phase unwrapping errors.Due to the long time interval of the Sentinel data and the severe subsidence situation in Tianjin, it could lead to poor coherence between scenes and phase unwrapping errors.Furthermore, a linear model was adopted to fit the Sentinel-1A image, which might generate unwrapping errors and lead to misleading results [48].Despite these issues, the long-term trends for both time series show good consistency.These results suggest that it is reliable enough for further research on the mechanism of deformation in Tianjin.

Time Series of Deformation
According to the processing procedure mentioned above, the cumulative time series of vertical deformation for Tianjin and its surrounding area were derived, as shown in Figure 11.To obtain a better surface deformation trend in this area, the oldest image, acquired on 9 January 2016, was selected as the time reference.From the spatial distribution of land subsidence, the center of the subsidence area exhibited two large shifts.In the early 1980s, with the completion of the Luanhe-Tianjin water diversion project, subsidence funnels were transferred from the central city to four districts in the vicinity (i.e., Beichen, Dongli, Jinnan, and Xiqing) and to Jinghai District.In recent years, in additional to Jinghai district, with the increasing control of ground subsidence in the four districts, the subsidence funnels have been transferred to Hangu district, and the junction of Wuqing, Beichen and Xiqing districts.
As illustrated in Figure 12, heterogeneous land deformation patterns are widely found in different districts of Tianjin and its surrounding cities.The deformation rates are relatively linear in the north part of Tianjin (i.e., Jixian, north Wuqing, Baodi, and north Ninghe).Actually, this is related to the geological conditions of these areas, where the groundwater is mainly stored in an aquifer with small storage, so that it has not been withdrawn.The deformation rates in the main urban areas of Tianjin are well-balanced and stable, and exhibit a uplifting trend.The average deformation rate along the direction of LOS reached 7.98 mm/yr, which indicates that the subsidence of Tianjin has been well controlled in recent years.Additionally, we found that the subsidence rate gradually increased from central Tianjin toward the surrounding areas.These subsidence areas include mainly three parts (Region 1 to 3 in Figure 12 ).Among these three subsidence funnels, Region 1 has the widest coverage and the most serious deformation rates.Additionally, two subsidence funnels have been detected in this continuous deformation region, one is located in the town of Wangqingtuo in the southern part of Wuqing District.The other is in Shengfang Town, east of Bazhou City in Hebei Province.This shows that the most serious area of subsidence is in Wangqingtuo Town, where the maximum subsidence rate reached −136 mm/yr.It has been difficult to determine the complex causes and to take effective measures to control this subsidence [7].Detailed discussion and explanation of possible causes for the various deformation patterns of Wangqingtuo and the town surrounding it will be presented in Section 5.
The maximum subsidence rates in the LOS direction in Shengfang Town reached −110 mm/yr.Due to the influence of these two major subsidence funnels, the areas of subsidence continue to expand outwards.Most of the northern and southern parts of Beichen District, central Wuqing District, and the north-eastern corner of Jinghai District are subject to different degrees of subsidence.In particular, the subsidence rate in the direction of LOS of Donggugang Town in south Langfang City is more serious (−35 mm/yr to −95 mm/yr), and located between the two subsidence funnels.Yangfengang Town has also been seriously affected.
As shown in Region 2, there are still large areas of subsidence in the northern and eastern parts of Jinghai District.Actually, as a region in Tianjin with a high degree of agricultural development, the amount of groundwater exploitation in Jinghai District is still increasing.Likewise, due to the rapid development of industrial production in the southern part of Xiqing District, a large area of subsidence is found in Wangwenzhuang Town (i.e., Region 3).The cities surrounding Tianjin, such as Langfang (i.e., Region 4) and Tangshan (i.e., Region 5, industrial area) exhibit the highest subsidence rates, −75 mm/yr and −84 mm/yr in LOS direction, respectively.In particular, the support of aquaculture policies has driven further increase in the amount of groundwater exploitation and the average subsidence velocity has reached −21 mm/yr in the east of Hangu District in LOS direction.

Subsidence Investigation and Their Causes
In order to better study the latest deformation velocities and trends of the serious subsidence funnels, and analyze its subsidence mechanism, the subsidence funnel of Part 1 in Figure 12 was enlarged as shown in Figure 13.It can be observed that most parts of the subsidence areas are located in and around Shengfang Town, Donggugang Town, and Wangqingtuo Town, and are distributed on both sides of the Rongcheng-Wuhai (i.e., Rongwu) highway and of the Tianjin-Baoding Intercity Railway.We can also see that the coverage area of Shengfang Town is not only the largest but also the one with the greatest area of subsidence.The most severe subsidence area was located in the western and north-eastern parts of Shengfang Town and in small portions of the eastern parts, showing an obvious regional subsidence pattern.The most serious subsidence area detected was located in Wangqingtuo Town, where the maximum subsidence velocity reached −136 mm/yr.The intercity railway from Tianjin to Baoding and the highway from Rongcheng to Wuhai are also affected by serious subsidence, and show obvious subsidence velocities in different sections.In order to figure out the subsidence mechanism, areas A, B, and C were enlarged, a well station has also been chosen to illustrate the variation of groundwater level, as shown in Figures 14-16.The subsidence in Tianjin might be caused by urban exploitation and geological influence.For Shengfang Town, urban development may play the most important role.The groundwater storage in Shengfang Town is abundant, but the shallow water is salt water, not suitable for irrigation.The deep groundwater has low mineralization and is widely exploited.There are 26 different types of industrial districts with more than 5200 industrial enterprises in this town.These are mainly furniture manufacturing and steel factories.As shown in Figure 14, the factories are mainly distributed in the north part of the North Outer Ring Road and east of Shengda Road.Most of these are furniture factories and there is a dense residential area located west of Shengda Road.A distinct subsidence funnel was detected in Zones a1 to a6 and in the Blue Line Zone.It can be seen that the subsidence rates in the a1, a2, and a3 commercial areas and near the new high-rise buildings in the south part of the North Outer Ring Road are much more serious than in the surrounding area (−110 mm/yr to −83 mm/yr).There has been little research on Shengfang Town, especially geological strata in this area.We have obtained a small amount of soil data from the Government development document.The document indicates that the soil in the north part of Zhongting River is lacustrine sediment, mainly fine silt and fine medium sand.Typically, in the centre of Shengfang Town, the soil is miscellaneous fill and the bearing capacity of such foundations is small, merely about 6 to 8 ton/m 2 .This weakness could easily lead to subsidence in major loading areas with intensive construction of buildings.
The north of the North Outer Ring Road is the old industrial district, most of the factories in this area are operating poorly, and some have already faced production suspension or closed down.However, the factories east of Shengda Road are operating well and most of these factories are new enterprises.Figure 14 illustrates that a high density of PS pixels were finally identified and less serious subsidence rates are detected in the new industrial district (−66 mm/yr to −10 mm/yr) than the old area (−110 mm/yr to −60 mm/yr).Furthermore, serious subsidence funnels were detected in Zones a4, a5, and a6.After investigating, we found that these factories are mainly for metal processing and glass manufacturing, so excessive groundwater exploitation may be the main cause of the differences in subsidence.In order to monitor the variation of groundwater level caused by the over exploitation of groundwater, a number of well stations have been built in different areas.A well station has been chosen to illustrate the variation of groundwater level, as shown in Figure 15.  Figure 15 shows that, apart from the visually prominent seasonal oscillations, a long-term linear trend was also detected in the time series.It indicates that the groundwater level has decreased year by year, mostly caused by the consumption of groundwater.This is even though the well station is located in Xiqing District, about 45 kilometers away from Shengfang Town.However, both of these areas are industrial cities, and Figure 12 shows that the subsidence in Shengfang Town and Wangqingtuo Town are the most severe in the study area.Therefore, the groundwater level variation in Xiqing District could also represent the groundwater level variation trend in Shengfang Town and Wangqingtuo Town to some extent.
To better study the surface deformation for the whole town, a profile line from point P1 to point P2 along the direction of the town was selected.The mean LOS deformation velocity of the profile is shown in Figure 14.It indicates that ground subsidence rates change with the different types of land use, as indicated by good correlation.In residential areas mainly located east and west of Shengfang Town, the subsidence rates range from −10 mm/yr to −60 mm/yr.The subsidence rates in the commercial areas are more serious than in the residential areas, with subsidence rates ranging from −30 mm/yr to −90 mm/yr.The most serious subsidence area is the industrial area, where the maximum subsidence rate reached −110 mm/yr.In 2017, Guo et al. [48] used a total of 27 Sentinel-1A images over Tianjin acquired between 31 May 2015 and 13 May 2016 to derive the subsidence of Tianjin.They found that the maximum subsidence rates have reached about −80mm/yr.The subsidence rate is further accelerating during 2016 to 2017.This result shows that the ground subsidence there is mainly caused by surface loads and groundwater extraction.In Shengfang Town, due to the large number of factories, the subsidence caused by over exploitation of groundwater is more serious.

Influence of Geological Factors on Subsidence
Region B and C have also been enlarged in Figure 16.We can see different subsidence rates on both sides of Line L1 and Line L2.In Donggugang Town, some research has indicated that a rock layer was found when drilling at a depth of 100 m.Thus, different geological conditions may cause the differential subsidence rates on two sides of Line L1 [62].Although much research has been carried out on subsidence in Wangqingtuo Town, little attention has focused on these differences.Here, hydrological and geological factors were analysed to discover the reason for the differences.The geological and hydrological map of Wangqingtuo Town is shown in Figure 17.
Figure 17 shows that the aquifer can be divided into four layers [62], with the first to fourth layers under Wangqingtuo Town at depths of 30 to 70 m, 180 to 220 m, 280 to 320 m and 370 to 430 m, respectively.A single well has also been used to assess the water abundance, and the experimental evidence suggests that the water yield property of the first layer is less than 100 m 3 /d and that the other three layers are in a relatively strong class (1000 to 3000 m 3 /d).This means that more space is needed between wells when drilling for equal amounts of groundwater from the different aquifers.This may lead to more serious subsidence while extracting shallow groundwater from the first aquifer than while extracting the deep phreatic water in Wangqingtuo Town.The areas in Wangqingtuo Town have the same geological conditions, so the differences of subsidence rates between the residential area and the industrial area may be due to the depth of the wells.For factories which require a large amount of water, it may be that only the deep, plentiful phreatic water can meet its production needs.In contrast, the shallow groundwater could meet the demand for water consumption in residential areas.This could be a good explanation for the differences in subsidence.Of course, with further improvement of the water supply system, it is also possible that the industrial water is derived from wells located in the main urban area.These are also likely to lead to the present subsidence differences.In order to get better and rational interpretations, more field trips need to be conducted and more relevant data need to be obtained.Excessive extraction of groundwater is the main contributor to the subsidence in Tianjin.However, different geological conditions can also contribute to different degrees of subsidence.As Figure 16 illustrates, there are distinct differences in the subsidence rates within Donggugang Town and Wangqingtuo Town.Subsidence is more serious in Wangqingtuo Town than in Donggugang Town, and the maximum difference could reach 46 mm/yr.Generally, sediments underground to depths of about 20 to 30 m (e.g., Holocene) have great influence on the degree of subsidence, especially where shallow groundwater is withdrawn.The surface sediment of Wangqingtuo Town is alluvial-diluvial, which is composed of fine sand and silty sand with high porosity.The surrounding areas are formed by alluvium and lacustrine deposits, which are mainly composed of clay with low porosity.Therefore, subsidence is more likely to happen in Wangqingtuo Town when the same amount of groundwater is withdrawn in the two areas [62].
Additionally, since 2006, the groundwater exploitation in Wangqingtuo Town has extended from the third aquifer to the first aquifer and the amount of groundwater mining has also been increasing.Particularly in 2008, the amount of groundwater mining of the second aquifer was increased by about 26% relative to 2007, while use of the shallow aquifer was increased to 900,000, accounting for 36% of the total volume.Therefore, it is suggested that the continuing and increasing consumption of the shallow groundwater could accelerate the subsidence rates in Wangqingtuo Town, and that different geological conditions also have an important effect on surface subsidence.

Conclusions
In this study, we first applied StaMPS MT-InSAR methodology to Sentinel-1A TOPS SAR scenes collected from 9 January 2016 to 8 June 2017, to detect spatial and temporal ground deformation across the city of Tianjin.The evaluation of the results of the deformation in the Tianjin area suggests that the accuracy is very high in most of the area.The root mean square (RMS) reached 5.015 mm/yr and the standard deviation reached 5 mm/yr for 60.720 % of the PS points, except for the edge region where three-dimensional phase unwrapping cannot be well constrained.Furthermore, the InSAR-derived deformation time series show good agreement with the GPS time series, especially for the long-term trends, and present good consistency with each other.
This study shows that significant deformation has occurred in different parts of the Tianjin area, gradually increasing from the central part of the metropolitan areas toward the suburbs.The deformation rate of the main urban area is well-balanced and it is also relatively stable, with deformation rates ranging from 0 to 20 mm/yr in the LOS direction.This indicates that the Tianjin ground subsidence has been well controlled due to groundwater recharge by the government and factory migration in recent years.Remarkable subsidence funnels have also been detected in the eastern part of Tianjin, especially in Shengfang Town, Donggugang Town, and Wangqingtuo Town, with maximum subsidence rates reaching −120 mm/yr.Subsidence funnels in Wangqingtuo Town are mainly caused by intensive construction of buildings with heavy surface loading as well as industrial production with over exploitation of groundwater.
From the difference of subsidence rates in both Donggugang Town and Wangqingtuo Town, we draw the conclusion that different geological conditions also contribute to different degrees of subsidence, and that the maximum difference could reach 46 mm/yr.Shallow deposits formed by alluvial-diluvial sediment are more likely to exhibit subsidence than that by alluvium and lacustrine deposits in the case where the same amount of groundwater is withdrawn.Furthermore, obvious subsidence differences between the residential area and industrial area may demonstrate that differences in the depth of the wells used could also cause different subsidence rates, even in the same geological region.
Author Contributions: All authors contributed to the manuscript and discussed the results.T.Z.performed all the data processing and analysis, and contributed to the manuscript of this paper.W.-B.S. designed and performed the scheme of this paper, and provided critical comments, re-analysis and interpretation, and contributed to final version of this paper.W.W. provided the technical guidance on Sentinel-1A data processing, B.Z. and Y.P. contribution to the discussion.All authors discussed the results and implication and comments of the manuscript at all stage.

Figure 1 .
Figure 1.This satellite image shows the location of the study area in our research and the red line is the provincial boundary.The pink rectangular frame denotes the coverage of the Sentinel-1A image and the red triangles represent the distribution of GPS stations.

Figure 2 .
Figure 2. The tectonic framework of the study area.Piedmont plain, alluvial fan and flood plain, and coastal plain in North China Plain are bounded by thick red and blue lines.Major faults and subsidiary faults are shown in thick yellow and black solid lines, respectively.(Revised from [46])

Figure 3 .
Figure 3. Temporospatial baseline distribution of ascending Sentinel-1A interferograms.The vertical axis shows the space perpendicular baseline in metres, the horizontal axis shows the image acquisition dates.All of the baselines represent interferometric pairs used for the time-series analysis.The green triangle represents the master image and the blue asterisks represent the slave images.

Figure 4 .
Figure 4. Distribution of the calculated temporal coherence γ for all of the PS pixels: The X-axis and Y-axis denote the coordinates of PS pixels in the direction of Azimuth and Range, respectively.The color scale shows the temporal correlation ranges from '0' to '1'.

Figure 5 .
Figure5.The blue line shows the relationship between the temporal coherence γ and the theoretical value of the velocity accuracy σ v for PS points (i.e., left-hand vertical axis).The histogram represents the statistic for the calculated coherence of PS points (i.e., right-hand vertical axis).

Figure 6 .
Figure 6.Distribution of the standard deviation of the velocity for PS points in the direction of LOS.

Figure 7 .
Figure 7. Theoretical and calculated statistics of the standard deviation of the mean velocity for PS points in the direction of LOS.Blue histogram denotes the theoretical statistics (left axis) and the red histogram denotes the calculated statistics (right axis).

Figure 8 .Figure 9 .Figure 10 .
Figure 8.Comparison of the original and filtered GPS vertical deformation time series.The grey line denotes the original GPS time series and the red line denotes the time series filtered with a Db3 wavelet.

Figure 11 .
Figure 11.Time series of the vertical deformation for Tianjin from 26 February 2016 to 8 June 2017 referenced temporally to the earliest image (9 January 2016) and referenced spatially to a stable CORS station (JIXN) located in the north-east corner of Tianjin.

4. 4 .
Figure 12 shows the annual surface deformation rates derived by the PS technique along the LOS direction of ascending pairs.The blue triangle denotes the reference points (TJA2), located in the north-east of Tianjin, north of Baodi District.The positive and negative signs indicate uplift and subsidence, respectively, in the direction of LOS.A total of approximately four million PS pixels were ultimately identified.Due to the difference in the distribution of man-made objects, significant gaps exist between different parts of Tianjin.The distribution of PS points with higher density are approximately located in the southern part of Wuqing District, the southwest part of Beichen District, in central Tianjin, in central south Dongli District, in north and central Binhai New District, in Jinnan District, and in Jinghai District.Other areas with more vegetation have fewer PS points.

Figure 12 .
Figure 12.Linear deformation trend derived by PS technique along the LOS direction of ascending pairs, blue triangle denotes the reference point, red triangles denote the GPS stations and the double circles denote the well station ZDK.Each administrative region and severely deformation areas are marked in black font.1-5 represent five major subsidence areas which enclosed by red lines.Two surface uplifting areas are also represent by blue ovals.

2 T
ia n B a o Ra il way Rongwu Hi g h w a y

Figure 13 .NorthFigure 14 .
Figure 13.The mean LOS deformation velocity of the subsidence funnel over Part 1 as shown in Figure 12: The white boxes A, B, and C represent Shengfang, Donggugang and Wangqingtuo Town, respectively.The red lines L1 and L2 represent the borders of different subsidence velocities.

Figure 15 .
Figure 15.Groundwater level variation of the well station ZDK.The location of this station has been given in Figure 12.Red dotted line denotes the time span for the Sentinel data acquisition.Vertical axis denotes the measured values, which is contrary to the actual groundwater level.

Figure 16 .
Figure 16.The mean LOS deformation velocity map of Shengfang Town and Wangqingtuo Town: The black dotted line from point P3 to P4 represent the profile of the Jintong Road.

Figure 17 .
Figure 17.(a) Geological map of Wangqingtuo Town, (b) Geological profile map (above) and Hydrological profile map (below).The black box in panel b denotes the span in the direction from north-west to south-east, W1 denotes Wangqingtu Town, W2 and W3 denote Xiqing District, the red dotted line L3 in the panel a denotes the profile from W1 to W3.This figure is revised from [62,63].

Table 1 .
Parameters of the Sentinel-1A data.

Table 2 .
Statistics table for the calculated temporal coherence.