ULTRASHORT-BASELINE PERSISTENT SCATTERER RADAR INTERFEROMETRY FOR SUBSIDENCE DETECTION

This paper presents an improved approach of multi-temporal interferometric synthetic aperture radar (InSAR) for detecting land subsidence phenomena by using time series of high resolution SAR images. Our algorithm extends the capability of the temporarily coherent point (TCP) InSAR technique proposed previously to detect subsidence even in the case of a small number of SAR images available for a study area. The approach proposed in this paper is implemented by using the interferograms with ultrashort spatial baselines (USB) through several procedures, including selection of USB interferometric pairs, TCP identification, TCP networking and modeling, as well as TCP solution. As the topographic effects are negligible in the USB interferograms, an external digital elevation model is no longer necessary for differential processing, thus simplifying both TCP modeling and parameter estimating. The algorithm has been tested with the high resolution TerraSAR-X (TSX) images acquired over Tianjin (China), and validated by using the ground-based leveling measurements. The testing results indicate that the density and coverage extent of TCPs can be increased dramatically, and the quality of subsidence measurements derived by the USB-based TCPInSAR can be raised.


INTRODUCTION
Land subsidence is a common hazard that is usually caused by anthropic activities, which can bring severe consequences.Monitoring the progressive subsidence with InSAR can provide valuable data for assessing the impacts of land subsidence and for guiding compensation strategies.However, the subsidence measurements derived by InSAR are often degraded by two factors, i.e., spatial and temporal decorrelation (Massonnet & Feigl, 1998;Rosen et al., 2000;Liu, 2006;Zebker & Villaseno, 1992) and atmospheric artifacts, (e.g., Zebker & Villaseno, 1992;Buckley, 2000;Liu, 2003;Ding et al., 2004).
The persistent scatterer InSAR (PSI) technique is one of the advanced approaches of multi-temporal InSAR, which have been proposed in recent years to mitigate the technical limitations of InSAR, (Ferretti et al., 2000;Ferretti et al., 2001).The temporarily coherent point InSAR (TCPInSAR) is an alternative implementation of PSI, which was proposed by Zhang et al. (2011aZhang et al. ( , 2011b)).It has been successfully applied to extract deformation using ERS-1/2 SAR imagery with a resolution of decametres (Zhang et al., 2012).Without the need of phase unwrapping, the TCPInSAR approach is characterized by several features in TCP identification, TCP networking and TCP solution by a least squares (LS) estimator (Zhang et al., 2011b).As the selection of TCPs can be performed on a minimum of two SAR images, TCPInSAR is more advantageous when there are few SAR images available over an area of interest for deformation detection (Zhang et al., 2011b).
It is desired to extend the capability of TCPInSAR to detect surface deformation using the high resolution SAR images.This allows for extracting surface deformation in more detail (Wegmüller et al., 2010).When applying the TCPInSAR method to such high resolution SAR data, an external digital elevation model (DEM) with comparable resolution is usually not available for removing topography-related phase component for the purpose of deformation extraction.Utilization of a DEM with coarse resolution and low accuracy will result in underestimate or overestimate of deformation parameters (Massonnet & Feigl, 1998;Rosen et al., 2000;Liu, 2006).This paper presents an improved TCPInSAR approach to detect subsidence using time series of high resolution SAR images.It is implemented by using the interferograms with ultrashort spatial baselines (USB).As the USB-based interferograms are not sensitive to topographic effects, an external DEM is no longer required for differential processing, thus simplifying the TCP modelling and decreasing the number of parameters to be estimated by the TCP solution.The USB-based TCPInSAR algorithm will be tested with 23 high resolution TerraSAR-X (TSX) images acquired in a stripmap mode between 29 April 2009 and 11 November 2010 over Tianjin (China).The subsidence measurements derived by the USB-based TCP solution will be compared with those derived by the TCP solution using another set of interferometric pairs with long spatial baselines, and also validated by using the ground truth data obtained by precise leveling at seven leveling benchmarks and five man-made corner reflectors.

Selection of Interferometric Pairs with Ultrashort Baselines
The topographic phase contribution is modulated in the interferograms that are used for subsidence detection, which is a major concern for selecting the reasonable interferometric pairs (Massonnet & Feigl, 1998;Rosen et al., 2000;Liu, 2006).The conventional differential processing generally utilizes an accurate DEM of the imaged surface to remove the topographic fringes from the interferograms (Gabriel et al., 1989;Massonnet & Feigl, 1998;Rosen et al., 2000;Liu, 2006).For the high resolution SAR images acquired by a short wavelength radar system such as the TSX sensor, it is however a challenging task to obtain a DEM with both satisfactory accuracy and comparable resolution for removing topographic effects.It is a reasonable way to choose the interferometric pairs with ultrashort spatial baselines for minimizing topographic effects without needing an external DEM.
The criterion of ultrashort spatial baseline can be derived by analyzing the impact of topographic error onto the ground motion along radar line of sight (LOS).Assuming that the height error of the DEM used for differential processing is h  , the error l  in the resultant LOS motion can be expressed by (Massonnet & Feigl, 1998;Rosen et al., 2000) h R where  B denotes the perpendicular baseline length of an interferometric pair,  the radar look angle, and R the sensor-to- surface distance (i.e., slant range).It is clear that the contribution of the height error onto the LOS motion is proportional to both  B and h  .As an example, Figure 1

Temporarily Coherent Point Identification and Networking
The subsidence analysis concentrates on the temporarily coherent points (TCPs).The TCPs are defined as the surface targets that can retain satisfactory radar coherence in one or several interferograms with different temporal baselines (Zhang et al., 2011a).Different from the PS detection method proposed by Ferretti et al. (2001), we identify the TCPs by statistical calculation of the offsets in range and azimuth directions between two SAR images of each USB-based interferometric pair.Such offsets can be estimated through image matching on a pixel-by-pixel basis.The previous investigation indicated that the offsets estimated for strong scatterers (coherent) should be far less sensitive to the window size and the oversampling factor used for image matching than those for the distributed scatterers (incoherent) (Zhang et al., 2011a).
The TCPs can be selected through image matching between two SAR images (i.e., master and salve images) of each interferometric pair.For each pixel, we change the window size (e.g., varying from 5×5 to 125×125) for searching the corresponding points by correlation analysis from two SAR images to derive offsets at a sub-pixel level (Massonnet & Feigl, 1998;Balmer & Hartl, 1998;Rosen et al., 2000;Buckley, 2000).A set of offsets (OT j ) at a pixel (j) can be derived in this way for the N window sizes.Any pixel with a standard deviation of less than 0.1 pixels in offsets is determined as a TCP (Liu et al., 2009), i.e.,

 
1 2 ( ) 0.1 If the statistical result derived for a pixel cannot meet the requirement as given in (2), it will be identified as a false TCP.
After the selection of TCPs is completed for all the interferometric pairs, the two types of TCPs can be obtained, i.e. persistently coherent (PC) and temporarily coherent (TC) types.
In this work we concentrate on the PC type for extracting time series of deformations.
Once all the TCPs are identified, a network is formed by connecting the adjacent TCPs (each connection is referred to as an arc) for differential modeling.An approach of Delaunay triangulation has been widely used for this purpose (Liu et al., 2008;Liu et al., 2009;Liu et al., 2011).However, such a networking method cannot provide a robust framework for the multi-temporal InSAR analysis (Liu et al., 2009;Liu et al., 2011).To balance between the computation burden and the reliability for the TCP solution, we connect the TCPs using the local Delaunay triangulation method.For more details, see Zhang et al. (2011b).A grid with a given spacing (e.g., 100 m) is first imposed to the entire study area, and the local Delaunay triangulation is then performed around each grid node.All the TCPs within the circle of a given radius (e.g., 750 m) centred on the grid node are used for the local Delaunay triangulation.

Temporarily Coherent Point Modeling and Solution by Least Squares Estimator
Suppose that N USB-based interferograms with temporal baselines of t i (i=1, 2, … , N) have been generated from a set of M SAR images acquired over the study area.For a TCP at the pixel x in the i th interferogram with temporal baseline t i , the interferometric phase ) , ( i t x Φ can be written as If only vertical motion occurs in the study area, the deformation phase is a function of the vertical motion rates, the time intervals, and the radar look angle.As the M+1 SAR images can be ordered by their acquisition dates, there are M successive time intervals (i.e., ) corresponding to M vertical motion rates at the pixel x, which can be denoted by Suppose that the i th ( ) interferometric pair is formed by combining the k th and l th ( ) SAR images, the deformation phase term in (3) can be expressed by 1 1 ( ) ( ) where  stands for the radar wavelength (3.1 cm for the TXS case), and  for radar look angle (41° for the TXS case).
The equation ( 3) can be rewritten as where is the sum of two phase components related to atmospheric delay and noise.As the number of parameters to be estimated can be reduced in this way, it is different from the treatment as done by Zhang et al. (2011b).
In practical, our TCP modeling for vertical motion analysis is based on the concept of neighbourhood differencing along each arc of the TCP network formed by the local Delaunay triangulation method.For the i th interferogram, the phase increment (i.e., differential phase) between two TCPs (at the pixels x and y) of an arc can be expressed by is the sum of differential atmospheric phase and noise, and ΔV is the vector of differential vertical motion rates, i.e.,

 
where As a large number of arcs have very short distances, it is reasonable to assume that there is no phase ambiguity in the phase increment ) ( i t y, x, ΔΦ at these arcs (Zhang et al., 2011a(Zhang et al., , 2011b)).It means that for these arcs the phase unwrapping procedure is not required for the subsequent TCP solution.In addition, the differential modeling on each arc can cancel out or decrease the impacts of atmospheric delay and other systematic errors on deformation analysis, as the atmospheric effects generally exhibit high correlation in space.For an arc formed by connecting TCP x and TCP y, N equations can be written with the N interferometric pairs as follow The M differential vertical motion rates as given in (10) are the unknown parameters to be estimated.The equations ( 12) can be rewritten in the matrix form as follows where A stands for the observation matrix and each of its elements is ( ) j   or zero.The observation vector  and the constant vector W in ( 13) are expressed by where  denotes the estimated quantities, and P dd is the a priori weight matrix which can be obtained by taking the inverse of the variance matrix Q dd of the observations.The Q dd can be derived from the standard deviations in offsets estimated during the TCP identification.It should be noted that T dd A P A may be singular in the case of N<M.For solution safety, a minimum norm solution is obtained using singular value decomposition (SVD), i.e.,     The residuals of the measurements can be derived by After the LS solution for an arc of the TCP network, the outlier detector proposed by Zhang et al. (2011b) is applied to check if the phase measurements for the arc have phase ambiguities.
Both the LS residuals and variance components for the arc are used for such a quality check.For subsequent analysis, the unacceptable arcs with outliers should be discarded.After deriving the parameters (i.e., , valid arcs, we can determine the parameters (i.e., vertical motion rates j v  , j=1, 2, …, M) in an absolute sense of all the valid TCPs by a spatial integration method with a given reference point (Zhang et al., 2011b).Finally, the vertical motion time series at each TCP can be derived by successively adding the multiplications of the estimated vertical motion rates and the corresponding time intervals.

STUDY AREA AND DATA SOURCE
With an area of almost 12000 km 2 , Tianjin is one of four municipalities in China.It suffers water shortage due to its natural geographic condition and semi-arid climate (Enviro-Library, 2008).To meet the industrial and agricultural needs, a large amount of groundwater has been overpumped since 1920s, thus leading to severe land subsidence in many areas of Tianjin (Enviro-Library, 2008).As shown in the inset of Figure 2, we select the western part of Tianjin as the study area.Being close to the Bohai Bay, the study area of 144 km 2 belongs to the alluvial plains and appears to be generally about 2-3 m above sea level (Enviro-Library, 2008).

Distribution of TCPs Identified for the Study Area
The TCPs were identified using the 23 USB interferometric pairs with spatial baselines of less than 15 m (see Figure 3) by following the procedures as presented in Section 2.2. Figure 4 shows the distribution of the resultant 642058 TCPs (marked by red dots) which are superimposed onto the TXS amplitude image of the study area.
The spatial resolution of TSX imagery primarily governs the density of TCPs.The about 2-m pixel spacing of TSX imagery benefits distinguish individual surface objects.The high density of TCPs identified in built-up areas using the TSX images can be explained by that fact that many TCPs are formed by dihedrals or trihedrals related to buildings in the urban areas.In addition, the special shape and material of roofs of the buildings in Tianjin also contribute to the high density of TCPs.Most of the buildings in Tianjin are constructed in spires shape with clay or metal tile roofs.If the roof surface faces toward the incident radar echoes, it can be viewed as the perfect backscattering body, thus guaranteeing the satisfactory radar coherence even over months to years.

Subsidence Signature Derived by TCP Solution
The subsidence signatures at all the 642058 TCPs in the study area were estimated using the 23 USB interferograms by following the procedures as described in Section 2.3.We calibrated the subsidence time series using as the reference data the levelling measurements obtained at CR5 (see Fig. 2).The resultant subsidence rates at all the valid TCPs are shown in Fig. 5.It can be seen from Figure 5 that a subsiding trough marked by a smaller oval appears around CR3, BM2 and BM4 (at lower-right corner in Figure 5) and has a peak subsidence rate of 29 mm/yr, while a wider subsiding trough marked by a larger oval appears around the point P (at upper-left corner in Figure 5) and has a peak subsidence rate of 52 mm/yr.According to our field investigation, the former subsiding trough corresponds to an industrial park where groundwater has been exploited recently for manufacturing purpose, while the latter one corresponds to the thermal power plant where groundwater has been excessive extracted for electricity production in recent years.
We selected 29 April 2009 as the temporal reference to generate the subsidence time series of each TCP.Thus, the subsidence derived by the TCP solution was set to 0.0 for this date.As examples, Figure 6 shows the temporal evolution of subsidence (indicated by blue dots) derived using the 23 USB interferometric pairs at two BMs (BM1 and BM7), three CRs (CR1, CR3 and CR4), and the maximum subsidence point P.
For the BMs and CRs, the subsidence values obtained by leveling are also marked by red squares, which agree well with the results derived by the USB TCP solution (for further numerical comparison, see Section 4.3).It is clear that the linear trend of subsidence dominates at all the sample points, and the subsidence magnitude accumulated between April and August of 2009 is relatively smaller than that over other time spans.
From the subsidence time series of CR3 and P, it can be inferred that the accumulated peak subsidence of the industrial-park subsiding trough reaches up 40 mm in about one year and seven months, while that for the subsiding trough around the heat power plant reaches up to 100 mm (i.e., the maximum subsidence magnitude in the study area) in the same period of time.
For comparison purpose, we have also carried out the TCP solution using the 111 interferograms with long spatial baselines (LSB) by following the method proposed by Zhang et al. (2011).For subsidence estimation, a SRTM DEM over this study area was used to remove topographic effects from the interferograms (Liu et al., 2008;Liu et al., 2009;Liu et al., 2011), and the height errors were also considered as the parameters to be estimated in the LSB-based TCP modelling and solution (Zhang et al., 2011b).Figure 7 shows the distribution of subsidence rates derived in this way for all the 318310 TCPs (see Figure 7).Although the overall subsidencerate pattern derived by the LSB-based TCP solution (Figure 7) is very consistent with that derived by the USB-based TCP solution (Figure 5), the subsidence information gaps in Figure 7 are apparently wider than those in Figure 5.In particular, some subsidence points are lost in rural areas for the LSB case, which are however crucial for depicting the subsidence evolution in time and space.The quality assessment for the two types of TCP solutions will be presented in Section 4.3.
Figure 5. Subsidence rates estimated for all the valid TCPs using the 23 USB interferometric pairs.Both seven BMs and five CRs are annotated for numerical analysis.The maximum subsidence point detected in the study area is indicated by "P".The two subsiding troughs are marked by ovals.

Validation with Levelling Measurements
To assess the quality of the subsidence measurements, we compared the subsidence results derived by the USB-and LSBbased TCP solutions with the ground truth data, i.e., subsidence data obtained at seven BMs (BM1-BM7) and five CRs (CR1-CR5) by precise leveling.For better analysis, the subsidence results derived by the TCP solutions were interpolated for the dates of leveling campaigns if necessary.Table 1 lists three types of subsidence rates (i.e., L v , USB v , and LSB v ) derived by leveling, the USB-and LSB-based TCP solutions, respectively, at the BMs and the CRs.The results are also given in Figure 8.
The discrepancies (i.e., ) are also provided in Table 1.The statistical calculation indicates that the mean and SD of the discrepancies between leveling and the USB-based results are 0.4 mm/yr and ±2.3 mm/yr, respectively, while those between leveling and the LSB-based results are -0.3 mm/yr and ±4.0 mm/yr, respectively.It is evident that the accuracies in subsidence rates derived by the USB-based TCP solution are almost two times higher than those derived by the LSB-based TCP solution.
For further validation, we check the quality of the two types of TCP solutions by using as the reference data the leveling subsidence values (in absolute sense) obtained at the BMs and the CRs.The results derived by the TCP solutions generally agree well with the leveling results, and the quality of the USBbased TCP solution is apparently higher than that of the LSBbased TCP solution.For examples, Figure 6 shows the good agreement between leveling and the USB-based TCP solution at BM1, BM7, CR1, CR3 and CR4.The statistical calculation indicates that the mean and SD of the subsidence discrepancies between leveling and the USB-based TCP solution at BM1, BM7, CR1, CR3 and CR4 are (0.1 mm, ±2.2 mm), (0.4 mm, ±2.1 mm) and (-1.7 mm, ±2.8 mm), respectively, over three time spans, while those between leveling and the LSB-based TCP solution are (0.9 mm, ±5.2 mm), (-2.0 mm, ±5.1 mm) and (-1.5 mm, ±5.2 mm), respectively.The averaged mean and SD for the former are -0.4 mm and ±2.4 mm, respectively, while those for the latter are -0.9 mm and ±5.2 mm, respectively.It should be noted that the three time spans are April to September, 2009, September, 2009to April, 2010, and April to November, 2010.It is demonstrated that the accuracies in the subsidence measurements derived by the USB-based TCP solution are about two times higher than those derived by the LSB-based TCP solution.
The preceding assessments on the distribution of TCPs and the quality of subsidence measurements indicate that the USB-based TCPInSAR method is more advantageous than the LSB-based one.For the former case, since the USB-based interferograms with spatial baselines of less than 15 m are generated without needing an external DEM for removing topographic effects, the subsidence extraction is independent of topographic errors, thus resulting in higher accuracies of the subsidence measurements.For the latter case, the interferograms with spatial baselines of less than 250 m contain the height errors which cannot be completely compensated by modeling and parameter estimation, thus propagating into the subsidence measurements and decreasing their accuracies.

CONCLUSIONS
This paper presents a USB-based TCPInSAR methodology for land subsidence detection based on the time series of high resolution satellite SAR images.An external DEM is no longer needed for differential processing of the USB-based interferograms, thus simplifying the TCP modeling and decreasing the number of parameters to be estimated.The algorithm was tested using 23 high resolution TSX images acquired over Tianjin (China) between 29 April 2009 and 11 November 2010 to detect land subsidence related to groundwater overpumping.The subsidence results derived by the USB-based TCP solution were compared with those derived by the LSB-based TCP solution, and also validated by using the leveling data obtained at seven BMs and five CRs.
The comparison analysis indicates that the subsidence information gaps for the LSB case are apparently wider than those for the USB case.Some crucial TCPs are lost in rural areas for the LSB case, and the margins of the two subsiding troughs become ambiguous due to the obvious reduction of TCPs.The validation with leveling data obtained at seven BMs and five CRs shows that the accuracies in subsidence rates and subsidence measurements derived by the USB-based TCP solution are 2.3 mm/yr and 2.4 mm, respectively, which are about two times higher than those derived by the LSB-based TCP solution.

Figure 1 .
Figure 1.The LOS motion error is represented as a function of the perpendicular baseline length and the height error for the case of TSX system.
surface motion along radar LOS, atmospheric delay along radar LOS, and decorrelation noise, respectively.The topographic phase component is not included in (3) due to the USB-based interferogram used.
CRs during the period of 23 TXS acquisitions.The leveling dates for the four epochs are around 20 April 2009, 5 September 2009, 15 April 2010, and 30 October 2010, respectively.

Figure 3
Figure3shows the spatial and temporal baselines of the 23 USB interferometric pairs selected using 15 m as the threshold from the 23 TSX images, in which the color bar indicates the variation of spatial baseline lengths of all pairs, and each connection indicates the temporal baseline (i.e., time span) for the corresponding interferometric pair.It can be seen that the shortest temporal baseline is 11 days, and the longest one is 198 days.The spatial baselines range between -13 and 14 m.As the USB pairs are not sensitive to topographic effects and high-rise buildings, we carried out the differential processing without using an external DEM.

Figure 2 .
Figure 2. The location map of the study area and the amplitude image averaged from 23 TSX images.

Figure 3 .
Figure 3.The spatial and temporal baselines of the 23 USB interferometric pairs used for TCP solution.

Figure 4 .
Figure 4.The distribution of the TCPs identified from 23 USB interferometric pairs for the study area.

Figure 6 .
Figure 6.Time series of subsidence derived at two BMs, three CRs and the maximum subsidence point P using the 23 USB interferometric pairs.Note that yy.mm stands for year.month,and one can read 09.04 as April of 2009, for example.

Figure 7 .
Figure 7. Subsidence rates estimated for all the valid TCPs using the 111 LSB interferometric pairs.

Figure 8 .
Figure 8. Subsidence rates derived by levelling and the USB-/LSB-based TCP solutions at seven BMs and five CRs.

Table 1 .
Comparison between subsidence rates derived by leveling and the USB-/LSB-based TCP solutions at seven BMs and five CRs.(unit: mm/yr)