AN IMPROVED MULTI-TEMPORAL INSAR METHOD FOR INCREASING SPATIAL RESOLUTION OF SURFACE DEFORMATION MEASUREMENTS

The multi-temporal interferometric synthetic aperture radar (InSAR) technology has proven very useful in extracting surface deformation with time series of SAR images over a study area. To increase spatial resolution of deformation information, this paper presents an improved multi-temporal InSAR (MTI) method by tracking both the point-like targets (PTs) and the distributed targets (DTs) with temporal steadiness of radar backscattering. The valid pixels corresponding to PTs and DTs are identified mainly by thresholding of the amplitude dispersion index (ADI) and the Pearson correlation coefficient (PCC). To efficiently reduce error propagation, a hierarchical analysis strategy is applied to extract deformation rates at the valid pixels. For the pixels with lower ADI values, the deformation rates are estimated on an optimized network by a least squared 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 estimated through the multi-levels of processing. The nonlinear deformation values at all the valid pixels are estimated by spatiotemporally filtering and spatially integrating. The proposed MTI algorithm has been tested for subsidence detection over Tianjin in China using the 40 high resolution TerraSAR-X images acquired between 2009 and 2010, and validated by using the ground-based leveling measurements. The testing results indicate that the spatial resolution and coverage of subsidence data can be increased dramatically by the hierarchical analysis, and the accuracy in subsidence values derived from the MTI solution can reach up to a millimeter level. Corresponding Author: Dr. Guoxiang Liu, E-mail: rsgxliu@swjtu.edu.cn


INTRODUCTION
The multi-temporal interferometric synthetic aperture radar (InSAR) technology has proven very useful in extracting surface deformation due to tectonic movements and anthropogenic activities.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 (Casu et al., 2008;Colesanti et al., 2003;Ferretti et al., 2001;Kampes, 2006;Lanari et al., 2004;Tizzani et al., 2007;Wegmüller et al., 2010).Although all the methods are proposed to concentrate on the PTs, some previous studies demonstrate that all PTs and some DTs can remain temporal steadiness of radar backscattering (Berardino et al., 2002;Ferretti et al., 2000;Hooper, 2008).However, 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 MTI processing is therefore not satisfactory for some studies areas.To increase spatial resolution and coverage of deformation information, this paper presents an improved MTI method by tracking both PTs and DTs with temporal steadiness of radar reflectivity.The two kinds of targets are identified using both ADI and PCC.The proposed MTI algorithm were tested for subsidence detection over Tianjin in China using the 40 high resolution TerraSAR-X (TSX) images acquired between March of 2009 and December of 2010.For accuracy analysis, the ground truth data collected by precise leveling were used to compare with the subsidence measurements derived from the MTI solution.

HIERARCHICAL PROCESSING FOR DEFORMATION RATES
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 (Berardino et al., 2002;Blanco-Sá nchez et al., 2008;Colesanti et al., 2003;Ferretti et al., 2001;Hooper, 2008;Kampes, 2006;Mora et al., 2003), we select the valid pixels with temporal steadiness of radar reflectivity by using two criteria (i.e., ADI and PCC), taking into account both temporal and spatial correlation of phase data used for deformation analysis.To efficiently reduce error propagation, the pixels with lower ADI values are first treated, and those with higher ADI values are treated through the multi-levels of processing.The primary models and procedures for extraction of deformation rates will be presented in this Section.

Pearson Correlation Coefficient
Rodgers and Nicewander reported that the PCC is a good indicator for measuring the data correlation in time and space (Rodgers and Nicewander, 1988).We introduce it into our MTI analysis.For a link between two adjacent pixels, the PCC ( ( , ; , ) i i j j x y x y


) can be empirically estimated using the phase values obtained from N interferograms by respectively, derived from the N interferograms.The PCC varies within [0, 1], and a greater PCC means a higher spatiotemporal correlation between two adjacent pixels.
To check the usefulness of PCC, a simulation experiment was performed.We set v  and  between two adjacent pixels to be -3.5 mm/yr and 6.7 m, respectively.A pixel with ADI value of 0.25 was chosen as a reference pixel (x i , y i ) from the real data used in Section 3, thus obtaining a sequence of phase values for all the interferograms.The sequence of phase values at the adjacent pixel (x j , y j ) was obtained by adding the relevant phase increments calculated using the given v  and  into the sequence of phase values at the reference pixel.We added the normally distributed noises into the phase data.Four noise levels with standard deviations (SDs, ) , ( i i y 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 ( ) , (  It can be seen from Fig. 1 that the PCCs nonlinearly decrease with increasing of SDs in the phase data at the adjacent pixel (x j , y j ).Comparison between four cases of noise levels of 0.3, 0.5, 0.7 and 0.9 radians (depicted by Fig. 1(a), (b), (c) and (d), respectively) at the reference pixel (x i , y i ) indicates that the PCCs drop more significantly if the SD in the phase data at (x i , y i ) is greater.It is also visible that the uncertainties (i.e., SDs) of PCCs increase with increasing of SDs in the phase data at the two adjacent pixels.These demonstrate that the PCC can be used to efficiently measure the spatiotemporal correlation of phase data between two adjacent pixels.In addition, the further simulation shows that the PCC directly derived from phase data is very sensitive to the ambiguity of phase values, which is obviously advantageous if compared with other measures such as MCF.

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 LOS direction at the pixels (primarily corresponding to PTs) with higher signal to noise ratio (SNR) in phase data.As suggested by Hooper et al. (2008), 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. (2009), we create a freely connected network (FCN) of the SPP.If any pixels in the SPP cannot be connected with any neighboring pixels, they are the isolated pixels and discarded for further processing.The estimation of LOS deformation rates is performed on the basis of the FCN using the LS solution (Zhang et al., 2011) and the region growing (Mora et al., 2003).

Estimating Deformation Rates at the Pixels with Higher ADI Values
To extend the coverage of deformation measurements, 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 time series of phase values, which most likely correspond to the DTs with temporal steadiness of radar reflectivity (Colesanti et al., 2003;Ferretti et al., 2011;Hooper, 2008).Some good examples of DTs include asphalt roads, concrete roads, roofs, noncultivated lands and dessert areas with sufficient roughness (Ferretti et al., 2011).Such kind 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.1.
To reduce possibility of error propagation, all the pixels with ADI values of  0.4 are processed hierarchically for extraction of deformation rates.These pixels are first classified into Q groups by the ADI intervals, i.e., (0.4, 0.5], …, (0.3+0.1×i, 0.4+0.1×i],(0.3+0.1×Q, …, 0.4+0.1×Q], in which i=1, 2, …, Q.The LOS deformation rates are separately estimated for the order of Group 1, 2, …, and Q.It is pointed out here that the subset of valid pixels determined using the procedures as described in Section 2.2 belongs to Group 0. Our testing shows that the hierarchical processing of 8 to 10 groups is suitable for the TSX imagery.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 (e.g., 25 pixels for the TSX case).The second criterion is that the PCC between the BNP and the POI should be highest and at least greater than a given threshold of 0.75.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 by Zhang et al. (2011).The LOS deformation rate and the elevation error at the POI are derived by adding v  and  into those at the BNP, respectively.

Spatiotemporal Filtering for Nonlinear Deformation
Once the LOS deformation rates and the elevation errors at all the valid PTs and DTs are derived from the hierarchical processing, further processing concentrates on separating nonlinear motion from atmospheric artifact.Suppose that the numbers of all the valid pixels and links determined from the Groups 0-Q are S and P. Then P phase gradients can be generated in both range and azimuth directions.We use the spatiotemporal filtering (Ferretti et al., 2000) on the phase gradients instead of the unwrapped phase.Then the phase gradients are integrated to extract both the nonlinear deformation and the atmospheric phase screen.Finally, the whole deformation evolution is calculated by using

Study Area and Data Source
As shown in Fig. 2, we select the northwestern part of Tianjin in China as the investigating area, which possesses the complex hydrological settings and the poor engineering geological properties (Hu et al., 2009;Hu et al., 2004;Lixin et al., 2011).
It is reported that the land subsidence around Tianjin is ongoing due to overuse of groundwater (Liu et al., 2011).In this study, we will utilize 40 TSX images collected between 27 March 2009 and 14 December 2010 over Tianjin for detection of land subsidence.The TSX image collected on 13 November 2009 is selected as the unique master image for the interferometric combinations.
The 39 full resolution differential interferograms were generated using the GAMMA DIFF processor through the twopass differential processing approach (Wegmüller et al., 2010).

Subsidence Results and Analysis
The statistical analysis indicates that the mean and SD of all the ADI values in this experiment 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 40 TSX differential interferograms.A total of 10 groups of pixels with ADI values of [0, 0.4], (0.4, 0.5], …, and (1.2, 1.3], respectively, were processed on a group-by-group basis.In Group 0, as the rigorous quality control was performed during the solution, the 91601 pixels were actually determined as the valid pixels from all the 137698 pixels with ADI values of (0.4, 0.5].The deformation rates can be easily converted to the subsidence rates by following the method by Liu et al. (2011).Fig. 3 shows the subsidence rates estimated by the MTI method at the 91601 valid pixels which mainly correspond to PTs with temporal steadiness of radar backscattering (Hooper, 2008).
As shown in Fig. 3, most of PTs are located in the built-up areas and along the road network due to 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 the less temporal variability of radiometric property to incident radar echoes.However, very sparse PTs appear in farming lands and vegetated areas due to 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.
The subsidence rates in the study area range between 0 and 72 mm/yr.A subsiding trough with the maximum subsidence rate of 65 mm/yr is marked by an oval labeled by S 1 in the upperleft corner of Fig. 3.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.
Figure 3. Distribution of subsidence rates estimated at 91601 valid pixels identified from Group 0 with ADI values of (0.4, 0.5].An oval labeled by S 1 is marked for a subsidence trough that corresponds to a thermal power plant. 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 9 by the hierarchical processing.For a pixel in Group 1, its BNP is first selected from the valid pixels of the previous groups by using two criteria related to geometric distance and PCC.The geometric distance between the BNP and the pixel should be shortest and at least less than 25 pixels, while the PCC between the BNP and the pixel should be highest and at least greater than 0.75.The pixel is discarded for further analysis if the criteria cannot be met.Otherwise, the relative deformation rate and relative DEM error between the BNP and the pixel are estimated by the LS solution.The subsidence rate and the elevation error at the pixel are derived by adding the relative values at the BNP, respectively.If the quality check (see Section II-C) cannot be passed, the pixel is viewed as the invalid pixel.In the same way, all the pixels in Group1 were processed on a pixel-by-pixel basis, thus selecting the valid pixels and determining the subsidence rates and the elevation errors at these pixels.After treating Group 1, the same procedures were carried out for the order of Group 2, 3, …, and 9, on a group-by-group basis.Table 1 lists the number of valid pixels obtained for each of Groups 0-9, at which the subsidence information have been extracted.The expansion situation of subsidence rates from Group 1 to 9 is shown in Fig. 5(a)-(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 Fig. 5(a) are obtained from Groups 0-1, while the results in Fig. 5(i) are from Groups 0-9.
As shown in Table 2, the total number of the valid pixels obtained from Groups 0-9 reaches up to 322301.The number of valid pixels obtained from Groups 1-9 is 230700, being 252% more than that (91601) 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., 177805) of valid pixels.Group 2 with ADI values of (0.5, 0.6] provides the second contribution, resulting in 40551 valid pixels.For Groups 3-9, the numbers of valid pixels range between 563 and 3364.Groups 0-2 provide a contribution of 96.2% to the total valid pixels in this study area.As shown in Fig. 4, 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 (Fig. 3) is very consistent with those derived from the multiple groups (Fig. 4), the subsidence information gaps in Fig. 3 are apparently wider than those in Fig. 4. 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 Fig. 3 and 4 indicates that the integrity and margin of the subsidence trough marked by S 1 have been improved through the multi-levels of processing.

Group
No.
Number   ) shows comparisons between subsidence rates derived from Group 0 and Groups 0-9 in the four selected areas labelled by S 6 , S 7 , S 8 and S 9 , respectively.It is clear that the density and spatial coverage of the subsidence rates in the areas are apparently increased by the hierarchical processing.
Our field investigation for the areas indicates that some valid pixels increased from the multi-levels of processing are coincident with the DTs along the asphalt road (S 6 ), on the building tops (S 7 , S 8 ), in bare lands (S 8 ) and around the overhead-bridge system.Previous studies reported that the covering layer on the bed rock in the study area belongs to the neogene and quaternary sedimentary formation (Enviro-Library, 2008;Hu et al., 2004).This layer mainly consists of grits and loose or semi-loose argillaceous sediments.Such strata possess high 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, which is 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 MTI solution with the ground truth data, i.e., subsidence data obtained at BM1-BM7 and CR1-CR5 by precise levelling (see Fig. 7).For better analysis, the subsidence results derived from the MTI solution were interpolated for the dates of leveling campaigns if necessary.Fig. 7(a) shows the two types of subsidence rates derived from levelling and the MTI solution, respectively, at the BMs and the CRs, while Fig. 7(b)-(d) show the two types of subsidence values over three successive time spans derived from Group 0. Groups 0-2 provide a contribution of 96.2% to the total valid pixels in this study area.It is found that most of valid pixels increased from the multi-levels of processing are coincident with the DTs along the asphalt road, on the building tops, in bare lands and around the overhead-bridge system.The subsidence rates in the study area range between 0 and 72 mm/yr, and a subsiding trough with the peak subsidence rate of 65 mm/yr is revealed around the thermal power plant.The validation with levelling data at seven BMs and five RCs shows that the accuracies in subsidence rates and in subsidence values derived from the MTI solution can reach up to 2.5 mm/yr and 2-4 mm, respectively.
phase values at the two pixels, for the phase values at the adjacent pixel.5000 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 Fig.1(a)-(d) for the four cases of Fig.1depict the mean values of PCCs, while the error bars denote the SDs of PCCs.

Figure 1 .
Figure 1.The variation of PCCs represented as SDs in the simulated phase data at the two adjacent pixels.The solid curves depict the mean values of PCCs, while the error bars denote the SDs of PCCs.

Fig. 2
Fig. 2(a) shows our study area of Tianjin, China.In Fig. 2(b) seven levelling benchmarks (BM1-BM7) and the five manmade corner reflectors (CR1-CR5) are also plotted for validation purpose.The four epochs of levelling campaigns (with 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 levelling dates of the four epochs are around 20 April 2009, 5 September 2009, 15 April 2010, and 30 October 2010, respectively.

Figure 2 .
Figure 2. The location map and the study area around Tianjin in China.The coverage of the TSX scenes and the study area selected are marked in (a) by a larger and smaller filled rectangle, respectively.The amplitude image averaged from 40 TSX images of the study area is shown in (b) with annotation of seven BMs by red triangles and five CRs by green squares.

Figure 4 .
Figure 4.The subsidence rates expanded for the study area.(a)-(i) depicts the combined subsidence rates obtained from the previous and current groups, respectively.For example, (a) is obtained from Groups 0-1, and (i) from Groups 0-9.

Figure 5 .
Figure 5. (a) distribution of subsidence rates at all the valid pixels, and (b) comparisons between subsidence rates derived from Group 0 and Groups 0-9 in four areas labeled by S6, S7, S8 and S9, respectively.The nonlinear subsidence values at all the valid pixels were finally estimated through the spatiotemporal filtering and the spatial integration.As examples, Fig.6shows the time series of subsidence values (by red dots) at BM6, CR2, CR3 and CR4 derived from the MTI solution.The subsidence values (denoted by green squares) obtained by levelling agree well with the results derived from the MTI solution.

Figure 6 .
Figure 6.Time series of subsidence values (by red dots) derived from the MTI solution at BM6, CR2, CR3 and CR4.The subsidence values obtained by leveling are indicated by green squares.

Table 1 .
The number of valid pixels in 10 groups obtained by the hierarchical processing.