A Hierarchical Multi-Temporal InSAR Method for Increasing the Spatial Density of Deformation Measurements

Abstract

Point-like targets are useful in providing surface deformation with the time series of synthetic aperture radar (SAR) images using the multi-temporal interferometric synthetic aperture radar (MTInSAR) methodology. However, the spatial density of point-like targets is low, especially in non-urban areas. In this paper, a hierarchical MTInSAR method is proposed to increase the spatial density of deformation measurements by tracking both the point-like targets and the distributed targets with the temporal steadiness of radar backscattering. To efficiently reduce error propagation, the deformation rates on point-like targets with lower amplitude dispersion index values are first estimated using a least squared estimator and a region growing method. Afterwards, the distributed targets are identified using the amplitude dispersion index and a Pearson correlation coefficient through a multi-level processing strategy. Meanwhile, the deformation rates on distributed targets are estimated during the multi-level processing. The proposed MTInSAR method has been tested for subsidence detection over a suburban area located in Tianjin, China using 40 high-resolution TerraSAR-X images acquired between 2009 and 2010, and validated using the ground-based leveling measurements. The experiment results indicate that the spatial density of deformation measurements can be increased by about 250% and that subsidence accuracy can reach to the millimeter level by using the hierarchical MTInSAR method.

1. Introduction

The multi-temporal interferometric synthetic aperture radar (MTInSAR) methodology detects point-like targets (PTs) from SAR image time series to provide surface deformation information [1–7]. In recent years, it has been widely applied to investigate various deformation events related to volcano, earthquake, glacier, landslide, mining, groundwater exploitation, and so forth [3,8–13]. For MTInSAR analysis, the time series of satellite SAR images collected over a study area are used to analyze the spatiotemporal behavior of surface deformation in order to mitigate the negative effects of both temporal decorrelation [14] and atmospheric delay [15,16].

Several MTInSAR methods have been proposed by various research groups to overcome the technical limitations of the conventional InSAR [1–4,6,7,10,11,17]. Some methods, termed persistent scatterer (PS) InSAR, detect deformation only for point-like targets (PTs), such as rocks and man-made structures with the temporal steadiness of radar backscattering [3,9,13,18]. Some methods, termed the small baseline subset, detect deformation by using interferometric pairs with short spatial baselines to maximize interferometric coherence in the case of more distributed targets (DTs) available in a study area [1,11]. In addition, other methods can be applied to detect deformation for the specific study areas, including phase stability analysis [4], PS networking analysis [5,6], coherent pixels technique [2] and temporarily coherent point InSAR [17].

Previous studies demonstrated that all PTs and some DTs can remain temporally steady in radar backscatter [1,3,4,19], which can be used for deformation time series analysis. Ferretti et al. [3] reported that a PS corresponding to a PT can be identified by using the thresholding of the amplitude dispersion index (ADI), which can be derived by the statistical calculation of the time series of amplitude values at a pixel. However, the ADI analysis may result in the failure to identify a useful DT with the temporal steadiness of radar reflectivity [4,9]. Hooper [4] proposed analyzing the spatial correlation of phase data for raising reliability in PS solutions. The combined analysis of both PTs and DTs for deformation extraction is rarely considered in the literature. The spatial resolution of deformation measurements derived from MTInSAR processing is therefore not satisfactory for some study areas.

To increase spatial resolution and the coverage of deformation information, this paper presents an improved MTInSAR method by tracking both PTs and DTs with the temporal steadiness of radar reflectivity. As the Pearson correlation coefficient (PCC) works well for measuring the data correlation in time and space [20], we use the thresholding of both ADI and PCC for selecting the useful pixels corresponding to the PTs or DTs. To control error propagation, a hierarchical analysis strategy is applied to extract deformation rates at the useful pixels. For the pixels with lower ADI values, the deformation rates are first derived on an optimized network by a least squared (LS) estimator and a region growing method. For the pixels with higher ADI values, they are classified into several groups by the ADI intervals, and the deformation rates are then estimated through the multi-levels of processing.

For validation purpose, we test the proposed MTInSAR algorithm for subsidence detection over Tianjin in China using the 40 high resolution TerraSAR-X (TSX) images acquired between March 2009 and December 2010. For accuracy analysis, the ground truth data collected by precise leveling will be used to compare with the subsidence measurements derived from the MTInSAR solution.

This paper is organized as follows. Section 2 will present the modeling and hierarchical processing for the extraction of the deformation information. The experimental results and discussion will be given in Section 3. Some conclusions will be addressed in Section 4.

2. Hierarchical Processing for Deformation Extraction

Suppose that N differential interferograms can be generated from M SAR images available for a study area. We extract deformation rates at all the useful pixels by using a hierarchical processing approach. Different from the existing methods [1–3,18], we select the valid pixels with the temporal steadiness of radar reflectivity by using both ADI and PCC, taking into account both the temporal and spatial correlation of phase data used for deformation analysis. The PCC used to identify DTs will be addressed in this Section. It should be noted that the hierarchical processing strategy means that the pixels with lower ADI values are first treated, and the pixels with higher ADI values are then processed. In addition, the deformation time series at each useful pixel are derived by adding the linear motion component into the nonlinear motion components. The nonlinear motion components are estimated by using the spatiotemporal filtering method [21]. The primary phase models and procedures for extracting deformation information will be presented in this section.

2.1. Estimating Differential Deformation Rate between Two Valid Pixels

As a differential operation can be used to cancel out some spatially correlated errors [6,10,11], the basic phase modeling is performed along a link between two adjacent pixels. The differential phase between two pixels, x and y, in the k-th interferogram can be represented by:

$$\mathrm{\Delta }{\phi }_{(x,y)}^k=-\frac{4\pi }{\lambda }\mathrm{\Delta }{v}_{(x,y)}{T}^k-\frac{B_{\perp }^k}{R\text{sin}\theta }\mathrm{\Delta }{\varepsilon }_{(x,y)}+\mathrm{\Delta }{\phi }_{RE(x,y)}^k,$$

(1)

where T^k and $B_{\perp }^k$ denote the temporal and geometrical baseline, respectively; λ, R and θ are the radar wavelength (3.1 cm for the TSX sensor), the sensor-to-target distance and the incident angle, respectively; Δv and Δɛ are the differential deformation rate along radar line of sight (LOS) and the differential elevation error, respectively, between the two pixels; $\mathrm{\Delta }{\phi }_{RE(x,y)}^k$ is the residual phase consisting of increments of nonlinear deformation, the atmospheric phase and decorrelation noise.

If there are no phase ambiguities between two pixels, the least squared estimation can be applied to estimate Δv and Δɛ using N interferograms generated from N + 1 images through [17]:

$$\mathrm{\Delta }\mathbf{\Phi }=\mathbf{A}\left(\begin{array}{c}\mathrm{\Delta }v\\ \mathrm{\Delta }\varepsilon \end{array}\right)+\mathbf{W}$$

(2)

in which:

$$\begin{array}{l}\underset{N\times 1}{\mathrm{\Delta }\mathbf{\Phi }}={\left(\begin{array}{cccc}\mathrm{\Delta }{\phi }_{(x,y)}^1& \mathrm{\Delta }{\phi }_{(x,y)}^2& \cdots & \mathrm{\Delta }{\phi }_{(x,y)}^N\end{array}\right)}^T\\ \underset{N\times 2}{\mathbf{A}}=-{\left(\begin{array}{cccc}\frac{4\pi }{\lambda }{T}^1& \frac{4\pi }{\lambda }{T}^2& \cdots & \frac{4\pi }{\lambda }{T}^N\\ \frac{B_{\perp }^1}{R\text{sin}\theta }& \frac{B_{\perp }^2}{R\text{sin}\theta }& \cdots & \frac{B_{\perp }^N}{R\text{sin}\theta }\end{array}\right)}^T\\ \underset{N\times 1}{\mathbf{W}}={\left(\begin{array}{cccc}\mathrm{\Delta }{\phi }_{RE(x,y)}^1& \mathrm{\Delta }{\phi }_{RE(x,y)}^2& \cdots & \mathrm{\Delta }{\phi }_{RE(x,y)}^N\end{array}\right)}^T\end{array}$$

(3)

The most probable values of Δv and Δɛ are,

$$\left(\begin{array}{c}\mathrm{\Delta }\widehat{v}\\ \mathrm{\Delta }\widehat{\varepsilon }\end{array}\right)={\left({\mathbf{A}}^T\mathbf{A}\right)}^{-1}{\mathbf{A}}^T\mathrm{\Delta }\mathbf{\Phi }$$

(4)

where $\widehat{\cdot }$ denotes the estimated value. The vector of residual phases estimated can be represented by:

$$\widehat{\mathbf{W}}={\left(\begin{array}{cccc}\mathrm{\Delta }{\widehat{\phi }}_{\mathit{RE}(x,y)}^1& \mathrm{\Delta }{\widehat{\phi }}_{\mathit{RE}(x,y)}^2& \cdots & \mathrm{\Delta }{\widehat{\phi }}_{\mathit{RE}(x,y)}^N\end{array}\right)}^T=\mathrm{\Delta }\mathbf{\Phi }-\mathbf{A}{\left({\mathbf{A}}^T\mathbf{A}\right)}^{-1}{\mathbf{A}}^T\mathrm{\Delta }\mathbf{\Phi }$$

(5)

According to Zhang et al. [17], an easy and efficient outlier detection algorithm is used to reject the links with the effect of phase ambiguities,

$$\mathit{Max}\left(\mathrm{\Delta }{\phi }_{\mathit{RE}(x,y)}\right)>E\left(\mathrm{\Delta }{\phi }_{\mathit{RE}(x,y)}\right)+2\sigma \left(\mathrm{\Delta }{\phi }_{\mathit{RE}(x,y)}\right)$$

(6)

where Max(·), E(·) and σ(·) are the maximum value, expected value and standard deviation (SD) of a vector or matrix, respectively. When the Equation (6) is satisfied, the arc connecting two valid pixels is considered containing an outlier at the 95% confidence level. This arc is not used for further analyses.

2.2. Pearson Correlation Coefficient (PCC)

PCC was reported as a good indicator for measuring the data correlation in temporal and spatial domains [20]. A high PCC value assures the high spatial phase correlation between two points of a link. Besides, in the MTInSAR analysis, a high PCC value indicates that the dominant information in the phases of two points are spatially correlated, while the spatially uncorrelated terms, such as noise and part of the elevation errors [22], are considered as insignificant. To validate this, the modeling and simulation experiments will be carried out as follows.

For a link between two adjacent pixels, the PCC can be empirically estimated using the phase values obtained from N interferograms by:

$${\rho }_{(x,y)}=\left|\frac{\sum _{k=1}^N[{\phi }_x^k-\overline{{\phi }_x})({\phi }_y^k-\overline{{\phi }_y}]}{\sqrt{\sum _{k=1}^N{\left[{\phi }_x^k-\overline{{\phi }_x}\right]}^2\sum _{k=1}^N{\left[{\phi }_y^k-\overline{{\phi }_y}\right]}^2}}\right|$$

(7)

where ${\phi }_x^k$ and ${\phi }_y^k$ denote the phase values at the two pixels, respectively, obtained from the k-th interferogram; φ̄x and φ̄y are the mean phase values at the two pixels, respectively, derived from the N interferograms. The PCC varies within [0, 1], and a greater PCC means a higher spatiotemporal correlation between the phase of two adjacent pixels.

To check the usefulness of PCC, a simulation experiment was performed by using the temporal and geometric baselines given in

Table 1. The 40 TerraSAR-X (TSX) images and the interferometric parameters used in this study.

No. of Images	Imaging Dates	Tk (days)	$B_{\perp }^k$ (m)	No. of Images	Imaging Dates	Tk (days)	$B_{\perp }^k$ (m)
1	20090327 *	−231	42	21	20091205	22	127
2	20090407	−220	69	22	20091227	44	134
3	20090418	−209	−23	23	20100107	55	−25
4	20090429	−198	13	24	20100118	66	−27
5	20090510	−187	31	25	20100129	77	−7
6	20090521	−176	65	26	20100209	88	−383
7	20090623	−143	−76	27	20100220	99	−156
8	20090704	−132	−17	28	20100303	110	−152
9	20090715	−121	−33	29	20100314	121	−105
10	20090726	−110	−113	30	20100325	132	9
11	20090806	−99	139	31	20100405	143	−93
12	20090828	−77	−102	32	20100416	154	−127
13	20090908	−66	37	33	20100427	165	−36
14	20090919	−55	−64	34	20100621	220	19
15	20090930	−44	−182	35	20100702	231	−78
16	20091011	−33	−39	36	20100804	264	82
17	20091022	−22	−66	37	20100906	297	1
18	20091102	−11	120	38	20101009	330	158
19	20091113	0	0	39	20101111	363	−23
20	20091124	11	47	40	20101214	396	−93

^*For example, one can read 20090327 as 27 March 2009.

. For simulation purpose, we set Δv and Δɛ between two adjacent pixels to be −3.5 mm/yr and 6.7 m, respectively. A pixel with an ADI value of 0.25 was chosen as a reference pixel located at x, from the real data used in Section 4, thus obtaining a sequence of phase values for all the interferograms. The sequence of phase values at the adjacent pixel, y, was obtained by adding the relevant phase increments calculated using the given Δv and Δɛ by Equation (1) into the sequence of phase values at the reference pixel.

We added the normally distributed noises into the phase data. Four noise levels with an SD (σ_x) of 0.3, 0.5, 0.7 and 0.9 radians were considered for the phase values at the reference pixel, and the noise levels with SDs of 0–1 radians were considered for the phase values at the adjacent pixel. Five-thousand simulations were made for each case of noise levels and, thus, calculating the mean and SD of PCC. The statistical results for PCCs derived from the simulation experiment are shown in Figure 1a–d for the four cases of σ_x. For each case, the PCCs are represented as a function of σ_y. The solid curves in Figure 1 depict the mean values of PCCs, while the error bars denote the SDs of PCCs.

It can be seen from Figure 1 that the PCC values nonlinearly decrease with the increasing of σ_y. A comparison between four cases of σ_x of 0.3, 0.5, 0.7 and 0.9 radians (depicted by Figure 1a–d, respectively) at the reference pixel, x, indicates that the PCCs drop more significantly if the SD in the phase data at x is greater. The simulation results indicate that the PCC can be used to reflect the level of spatially uncorrelated noises and to assess the quality of phase values in the time domain. Besides, it is also visible that the uncertainties of PCCs increase with the increasing of SDs in the phase data at the two adjacent pixels. The uncertainties are mainly due to the phase ambiguities at the two pixels.

2.3. Estimating Deformation Rates at the Pixels with Lower ADI Values

To provide a reliable basis for the subsequent analysis, we first concentrate on estimating deformation rates in the LOS direction at the pixels with a higher signal-to-noise ratio (SNR) in phase data. As suggested by Hooper [4], we select the pixels with lower ADI values of ≤0.4 as the candidates for estimating deformation rates, thus obtaining a subset of potential pixels (SPP). As done by Liu et al. [6], any two adjacent pixels in the SPP are connected if their geometric distance is smaller than a given threshold, and the PCC between two pixels calculated using (7) is greater than a given threshold, thus forming a freely connected network (FCN) of the SPP. If a pixel in the SPP cannot be connected with any neighboring pixels, it is an isolated pixel and discarded for further processing. The estimation of LOS deformation rates is performed on the basis of the FCN using the LS solution and the region growing.

For each link of the FCN, the differential deformation rate between two pixels related to the link is estimated using the LS solution method, as described in Section 2.1, thus obtaining the relative deformation rates along all the links of the FCN. The subsequent processing is to estimate the absolute LOS deformation rates at the valid pixels in the FCN. Although Liu et al. [6] applied a LS solution for the entire FCN to estimate the absolute deformation rates, the computation cost is usually not acceptable, particularly for the processing of the high resolution SAR images. For example, in our experiment, which will be shown in Section 3, 91,610 PTs are maintained and 2,265,989 links are included in the FCN network. A design matrix for the LS solution by the method proposed by Liu et al. [6] requires the support of 773 GB of computer memory. Such a computation burden makes it impractical to apply the LS solution in high resolution MTInSAR processing. Although the LS estimation on the divided image blocks is possible [23], new problem arises, such as block-result mosaicking. We therefore propose to estimate the absolute LOS deformation rates through a region growing method to decrease the memory cost.

Starting from a reference pixel with known deformation information, the region growing is actually an operation of spatially integrating the relative deformation rates along the optimized paths of the FCN. It should be noted that the key to the region growing method is to search for the integral paths, thus ensuring the quality of such an integration solution. We take three constraints to find the nodes for the integral paths. If and only if a pixel meets the three constraints simultaneously, it can be accepted as a node of the integral paths. The first constraint is that the PCC between two adjacent pixels’ phase should be greater than a given threshold. The remaining two constraints are chosen to measure the spatial autocorrelation of the deformation rates and the elevation errors between the adjacent pixels, which are denoted by their localized root-mean-squared error (RMSE, i.e., δv and δɛ),

$$\begin{array}{l}\delta v=\sqrt{\frac{1}{L}\sum _{j=1}^L{\left(\mathrm{\Delta }{v}_{\left(x,y_j\right)}-\mathrm{\Delta }{\widehat{v}}_{\left(x,y_j\right)}\right)}^2}\\ \delta \varepsilon =\sqrt{\frac{1}{L}\sum _{j=1}^L{\left(\mathrm{\Delta }{\varepsilon }_{\left(x,y_j\right)}-\mathrm{\Delta }{\widehat{\varepsilon }}_{\left(x,y_j\right)}\right)}^2}\end{array}$$

(8)

where x is the coordinate of the pixel to be assessed; y_j are the adjacent pixels that have already been assessed and accepted; and L denotes the number of valid adjacent pixels in a window centered on x; Δv_(x,yj) = v_x – v_y; Δɛ_(x,yj) = ɛ_x – ɛ_y; Δv̂_(x,yj) and Δɛ̂_(x,yj) are the same as in Equation (4). According to the previous studies [6], δv and δɛ should not be greater than 5 mm/yr and 10 m, respectively, for optimizing the integral paths. The three aforementioned assessments are essentially indicators of phase quality. No matter how different the deformation patterns are, the high phase quality always ensures high values for the PCC and low values for the two localized RMSE indicators. Upon completion of region growing for the entire FCN, the valid PT pixels, their connectivity, deformation rates, as well as elevation error are determined, thus providing a basis for the subsequent hierarchical analysis.

2.4. Estimating Deformation Rates at the Pixels with Higher ADI Values

Different from the existing methods that concentrate on processing the PT candidates [10,21], the pixels that do not consist of dominant scatterers are also taken into account in our MTInSAR method. Therefore, we further treat the pixels with higher ADI values (>0.4) through the multi-levels of processing. It was demonstrated that some pixels with higher ADI values possess satisfactory SNR in the time series of phase values, which most likely correspond to the DTs with the temporal steadiness of radar reflectivity [4,19]. Some good examples of DTs include asphalt roads, concrete roads, roofs, non-cultivated lands and dessert areas with sufficient roughness [19]. Such kinds of targets can still be used for deformation extraction. We identify the useful pixels by mainly applying the spatiotemporal correlation analysis with the PCC thresholding, as discussed in Section 2.2.

To reduce the possibility of error propagation and to reduce the memory cost, all the pixels with ADI values >0.4 are processed hierarchically for the extraction of deformation rates. The mean ADI (m_ADI) and the SD of ADI (σ_ADI) are first calculated. Therefore, the upper bound of the ADI values is (m_ADI + 3σ_ADI) at a confidence level of 99.7%. According to the upper bound, all the pixels are classified into Q groups by the ADI intervals. The LOS deformation rates are estimated group-by-group using the procedures described in Section 2.3. The PTs belong to Group 0. To estimate the LOS deformation rates for the pixels in Group i, the valid results derived from the previous groups are used as the reference data.

For any pixel of interest (POI) in Group i, we first select the best neighboring pixel (BNP) from the valid pixels of the previous groups using two criteria, thus forming a link for estimating the deformation rate at the POI. The first criterion is that the geometric distance between the BNP and the POI should be shortest and at least less than a given threshold. The second criterion is that the PCC between the BNP and the POI should be greater than a given threshold. If the criteria cannot be met, the POI is discarded for further analysis. Otherwise, we then estimate Δv and Δɛ between the BNP and the POI using the LS solution method, as described in Section 2.1. The LOS deformation rate and the elevation error at the POI are derived by adding Δv and Δɛ into those at the BNP, respectively. The quality check is finally performed using the three constraints, as described in Section 2.3. If the quality check cannot be passed, the POI is viewed as an invalid pixel. In the same way, all the pixels in Group i are analyzed on a pixel-by-pixel basis, thus selecting the valid pixels and determining the deformation rates and the elevation errors at these pixels. After treating Group i, the same procedure are applied to Group i + 1, thus completing the analysis of all the Q groups.

For better understanding, Figure 2 shows the flowchart of the hierarchical processing procedures.

2.5. Extracting Nonlinear Deformation Components at the Useful Pixels

The nonlinear deformation components are extracted by using spatiotemporal filtering [21] at all the useful pixels, including both PTs and DTs. Previous studies indicate that the atmospheric perturbation and the nonlinear motion exhibits correlation in the space domain, and the former generally possesses a longer correlation distance than the latter [1,21]. In the time domain, the former is usually uncorrelated, due to different meteorological conditions on different SAR acquisition dates, while the latter is generally correlated due to temporal motion evolution [21,24]. As the nonlinear motion and the atmospheric artifact have different frequencies in time and space domain, the two components can be separated by spatiotemporal filtering [21].

Finally, for the given pixel of the k-th image, the full resolution deformation value, $S_{\mathit{full}}^k$, is the sum of the linear deformation components, $S_l^k$, and the nonlinear deformation component, $S_{nl}^k$:

$$S_{\mathit{full}}^k=S_l^k+S_{\mathit{nl}}^k$$

(9)

3. Experimental Results and Discussion

3.1. Study Area and Data Source

As shown in Figure 3, we selected the northwestern part of Tianjin in China as the area of investigation, which possesses complex hydrological settings and poor engineering geological properties [25–28]. It is reported that the land subsidence around Tianjin is ongoing, due to the overuse of groundwater [29]. In this study, we utilize 40 TSX images acquired between 27 March 2009, and 14 December 2010, over Tianjin for the detection of land subsidence. All the images were collected with an incidence angle of 41 degrees in HH polarization mode. Table 1 lists the acquisition dates of all the TSX images and the interferometric parameters used in this study. The original datasets were provided by Infoterra as single look complex images with pixel a spacing of 1.36 m in the slant range (2.07 m in the ground range) and 1.90 m in the azimuth.

Figure 3a shows that Tianjin is located in the eastern coastal region of the North China Plain, bordering with Beijing, Hebei and Bohai Bay. The entire coverage of the TSX scenes used for subsidence detection is marked by a larger filled rectangle in Figure 3a. We will conduct subsidence analysis in the typical area as marked by a small filled rectangle in Figure 3a. The amplitude image (5000 × 6200 pixels, equivalent to 10.4 × 11.7 km²) averaged from 40 TSX images of the selected study area is shown in Figure 3b with the annotation of seven leveling benchmarks (BM1–BM7) and the five man-made corner reflectors (CR1–CR5). It can be seen from Figure 3b that residential areas, cultivated lands, fishponds, rivers, railways and highways appear in the study area. The four epochs of the leveling campaigns (with an accuracy of 2 mm/km in height difference) were carried out for the BMs and CRs during the period of the 40 TSX acquisitions. The leveling dates of the four epochs are around 20 April 2009, 5 September 2009, 15 April 2010, and 30 October 2010, respectively. The leveling data will be used to validate the subsidence results derived from the MTInSAR processing.

The TSX image acquired on 13 November 2009, is selected as the unique master image for the interferometric combinations. The remaining 39 TSX images are used as the slave images, thus forming 39 interferometric pairs for subsidence analysis. Table 1 lists all the TSX images and the interferometric parameters (i.e., temporal and geometric baselines) used in this study for the MTInSAR processing. The 39 full-resolution differential interferograms were generated using the GAMMA DIFF processor through the two-pass differential processing approach. The digital elevation model (DEM, with a spatial resolution of 90 meters) derived from the Shuttle Radar Topography Mission (SRTM) is used to remove the topographic phases from the interferograms.

3.2. Subsidence Results and Analysis

Using 40 coregistered TSX images, the ADI values at all the pixels in the study area, as shown in Figure 3b, were first calculated by following the method by Ferretti et al. [3]. The statistical analysis indicates that the mean and SD of all the ADI values are 0.65 and 0.19, respectively, and an ADI value belongs to [0.08, 1.22] at a confidence level of 99.7%. We then extracted the deformation rates at the valid pixels through the hierarchical analysis by using the 39 TSX differential interferograms. A total of nine groups of pixels with ADI values of (0, 0.4], (0.4, 0.5], …, and (1.1, 1.2], respectively, were processed on a group-by-group basis. For the pixels in Group 0 with ADI values <0.4, the procedures as presented in Section 2.3, including the formation of FCN, LS-based estimation at each link, region growing and quality assessment, were applied to estimate the LOS deformation rates at the valid pixels. As the rigorous quality control was performed during the solution, the 91,601 pixels were actually determined as the valid pixels from all the 137,698 pixels in Group 0. The deformation rates can be easily converted to the subsidence rates by following the method by Liu et al. [29]. Figure 4 shows the subsidence rates estimated by the MTInSAR method at the 91,601 valid pixels, which mainly correspond to PTs with the temporal steadiness of radar backscattering [22]. It should be noted that all the subsidence results are referenced to CR5, and the spatial density of the deformation measurements is 753 km⁻².

As shown in Figure 4, most of the PTs are located in the built-up areas and along the road network, due to the availability of many natural and artificial hard objects, such as rocks in parks, buildings, bridges, iron fences and concrete bodies associated with urban infrastructures. These objects have less temporal variability of the radiometric property compared to incident radar echoes. However, very sparse PTs appear in farming lands and vegetated areas due to the temporal variability of radar reflectivity (e.g., caused by tilling, crop or vegetation growth and leaf fluttering related to wind effects). No PTs are available in fishponds, lakes and rivers, as expected, due to its specular behavior. The subsidence rates in the study area range between zero and 72 mm/yr. A subsiding trough with the maximum subsidence rate of 65 mm/yr is marked by an oval labeled by S₁ in the upper-left corner of Figure 4. According to our field investigation, the subsiding trough corresponds to a thermal power plant where groundwater has been over-pumped for electricity production in recent years.

As mentioned earlier, our region growing method is implemented through spatial integration by following the optimized paths to obtain the absolute subsidence rates, thus balancing between the solution reliability and the computation cost. The integral paths optimized using the three constraints (see Section 2.3) are crucial for the determination of the absolute subsidence rates at the valid pixels. As examples, four integral paths labeled by A, B, C and D are shown in Figure 5, which are used to calculate the subsidence rates at four PTs. Being close to the four integral paths, the four ovals labeled by S₂–S₅ are coincident with the rural areas where few PTs are available. Some fishponds around S₃ and S₄ cannot provide any PTs for subsidence tracking. It can be seen that the four integral paths go around the areas with very low radar coherence, without passing through them. It demonstrates that the integral paths selected by our method are reasonable and efficient for the determination of the absolute subsidence at the valid pixels.

To expand the subsidence information in the study area, we continued to estimate the subsidence rates at the valid pixels from Group 1 to Group 8 by the hierarchical processing discussed in Section 2.4.

Table 2. The number of valid pixels in nine groups obtained by the hierarchical processing.

Group No.	Number of Valid Pixels	Group No.	Number of Valid Pixels	Group No.	Number of Valid Pixels
0	91,601	3	2707	6	3364
1	177,805	4	1149	7	2196
2	40,551	5	1253	8	563

lists the number of valid pixels obtained for each of Groups 0–8, at which the subsidence information have been extracted. The expansion situation of subsidence rates from Group 0 to 8 is shown in Figure 6a–i for the study area. It should be noted that each sub-figure depicts the combined subsidence rates obtained from the previous and current groups. For example, the results in Figure 6a are obtained from Group 0, while the results in Figure 6i are from Groups 0–8.

As shown in Table 2, the total number of the valid pixels obtained from Groups 0–8 reaches up to 321,189; the corresponding spatial density is 2640 km⁻². While the number of valid pixels obtained from Groups 1–8 is 229,588, being 251% more than that (91,601) obtained from Group 0. The greatest contribution is provided by Group 1 with ADI values of (0.4, 0.5], resulting in the largest number (i.e., 177,805) of valid pixels. Group 2 with ADI values of (0.5, 0.6] provides the second contribution, resulting in 40,551 valid pixels. For Groups 3–8, the numbers of valid pixels range between 563 and 3364. Groups 0–2 provide a contribution of 96.5% to the total valid pixels in this study area. As shown in Figure 6, the spatial resolution and coverage of subsidence data are significantly enhanced by a combination of results derived from the various groups. Although the overall subsidence-rate pattern derived from Group 0 (Figure 4) is very consistent with those derived from the multiple groups (Figure 6), the point density of Group 0 is less than that of the multiple groups, according to the aforementioned data. In particular, some valid pixels in rural areas cannot be identified from Group 0, but can be detected from other Groups, which are crucial for revealing the subsidence evolution in the areas with low radar coherence. In addition, a closer inspection with Figure 4 and 6 indicates that the integrity and margin of the subsidence trough marked by S₁ have been improved through the multi-level processing.

For a better visualization, Figure 7a shows the enlarged map of Figure 6i, and Figure 7b shows comparisons between Figure 6a and Figure 6i in the four selected areas labeled by S₆, S₇, S₈ and S₉, respectively. The number of the valid points in the four areas is increased from [545, 321, 1,386, 454] to [8,972, 2,460, 5,449, 1,426], respectively. Our field investigation for the areas indicates that some valid pixels increased from the multi-level processing are coincident with the DTs along the asphalt road (S₆), on the building tops (S₇, S₈), in bare lands (S₈) and around the overhead-bridge system.

As the typical DTs, like asphalt roads, roofs and non-cultivated lands, can be seen in the cyan rectangle area marked in Figure 7a, we also utilized the Stanford Method for PS (StaMPS) to calculate the subsidence rates in this area for comparison purposes. By following the standard PS processing chains, we have identified a total number of 4909 valid points for the selected area, which is less than that (5395) from our MTInSAR method. Figure 8 shows the detailed comparison between the two types of subsidence rates derived from our method and the StaMPS for the selected area. Figure 8a–c depicts the subsidence rates by our method and the StaMPS and their comparison at the common points, respectively. The colorbar in Figure 8c indicates the discrepancies between the two types of subsidence rates at the common points. From Figure 8c, it can be seen that the subsidence rates by our method are in good agreement with those by the StaMPS. The statistical analysis indicates that the correlation of the two types of results is 0.83, thus meaning a strong consistency between our method and the StaMPS [30].

The nonlinear subsidence values at all the valid pixels are finally estimated using the spatiotemporal filtering method presented in Section 2.5. As examples, Figure 9 shows the time series of subsidence values (by red dots) at BM6, CR2, CR3 and CR4 derived from the MTInSAR solution. The subsidence values (denoted by green squares) obtained by leveling agree well with the results derived from the MTInSAR solution.

Previous studies reported that the covering layer on the bed rock in the study area belongs to the Neogene and quaternary sedimentary formation [25,27]. This layer mainly consists of grits and loose or semi-loose argillaceous sediments. Such strata possess high a compression ratio and are very suitable for the evolution of land subsidence, due to autologous gravity, surface loading and groundwater exploitation. Our field investigation indicates that a huge amount of groundwater has been pumped from the deep wells in the study area. We believe that this might be the primary cause for land sinking. Some efficient measures should be taken to optimize the groundwater use planning and decrease the groundwater use in the study area.

To assess the quality of the subsidence measurements, we compared the subsidence results derived from the MTInSAR solution with the ground truth data, i.e., subsidence data obtained at BM1–BM7 and CR1–CR5 by precise leveling (see Figure 10). For better analysis, the subsidence results derived from the MTInSAR solution were interpolated for the dates of leveling campaigns if necessary. Figure 10a shows the two types of subsidence rates derived from leveling and the MTInSAR solution, respectively, at the BMs and the CRs, while Figure 10b–d shows the two types of subsidence values over three successive time spans derived from leveling and the MTInSAR solution, respectively, at the BMs and the CRs. The statistical calculation indicates that the RMSE between two types of subsidence rates is 2.5 mm/yr, and the RMSEs between two types of subsidence values are 1.8, 3.6 and 3.8 mm over the three time spans, respectively. This indicates that the accuracy in subsidence values derived from the MTInSAR solution can reach up to the millimeter level.

4. Conclusions

A hierarchical MTInSAR method is presented in this paper to provide the high spatial density of subsidence measurements for a suburban area in Tianjin, China. Different from our previous studies presented in [6] and [29], we have proposed the improved MTInSAR for processing both the PTs and DTs in the study area, thus resulting in more details of the subsidence measurements. In the improved MTInSAR method, two important parameters, including ADI and PCC, are used to identify all the useful pixels from SAR image time series. All the image pixels are classified into several groups by thresholding the ADI values and processed group-by-group using the hierarchical processing strategy. The method was tested using 40 TSX images collected between 2009 and 2010 over Tianjin. Seven BMs and five CRs were used to validate the subsidence measurements.

In our experiment, nine groups of image pixels are considered. Group zero is PTs, which are expected to have a dominant scatterer within the pixels. Its number is 91,601. The corresponding spatial density is 753 km⁻². After the hierarchical processing, the spatial density is dramatically increased to 2640 km⁻². The density of valid pixels obtained from Groups 1–8 is 251% higher than that obtained from Group 0. Groups 0–2 provide a contribution of 96.5% to the total valid pixels in this study area. The valid pixels increased from Groups 1–8 are coincident with the DTs along the asphalt road, on the building tops, in bare lands and around the overhead-bridge system. The subsidence rate map derived from all the valid pixels reveals a tiny subsidence center with a subsidence rate of about 72 mm/yr. The validation has been carried out using 12 leveling points. The RMSEs of the subsidence values and the subsidence rates derived from the proposed MTInSAR method are 2–4 mm and 2.5 mm/yr, respectively.

Our experiments show that the proposed MTInSAR method is applicable in the areas with the availability of PTs. The typical areas include the urban areas, suburban districts, rocky areas, and so forth. For some areas, like vegetated or cultivated lands with a scarcity of PTs, our method may fail in detecting surface deformation. In addition, the comparison analysis indicates a strong consistency between our method and the StaMPS. By providing the high spatial density of the deformation data, the proposed MTInSAR method is useful for detecting the minor deformation bowls and assessing land deformation.