HIGH RESOLUTION DEFORMATION TIME SERIES ESTIMATION FOR DISTRIBUTED SCATTERERS USING TERRASAR-X DATA

: In recent years, several SAR satellites such as TerraSAR-X, COSMO-SkyMed and Radarsat-2 have been launched. These satellites provide high resolution data suitable for sophisticated interferometric applications. With shorter repeat cycles, smaller orbital tubes and higher bandwidth of the satellites; deformation time series analysis of distributed scatterers (DSs) is now supported by a practical data basis. Techniques for exploiting DSs in non-urban (rural) areas include the Small Baseline Subset Algorithm (SBAS). However, it involves spatial phase unwrapping, and phase unwrapping errors are typically encountered in rural areas and are difficult to detect. In addition, the SBAS technique involves a rectangular multilooking of the differential interferograms to reduce phase noise, resulting in a loss of resolution and superposition of different objects on ground. In this paper, we introduce a new approach for deformation monitoring with a focus on DSs, wherein, there is no need to unwrap the differential interferograms and the deformation is mapped at object resolution. It is based on a robust object adaptive parameter estimation using single look differential interferograms, where, the local tilts of deformation velocity and local slopes of residual DEM in range and azimuth directions are estimated. We present here the technical details and a processing example of this newly developed algorithm.


INTRODUCTION
Various differential SAR interferometric stacking techniques have been developed that exploit phase stable i.e. persistent scatterers (PSs). PS interferometry provides a parametric estimation of the displacement and 3D location based on the assumption of one or two dominant scatterers in the resolution cell (Ferretti et al., 2000, Ferretti et al., 2001, Kampes, 2006, Adam et al., 2008. However, it has some limitations in nonurban (rural) areas due to low density of PSs, their inhomogenous spatial distribution, phase ambiguities and atmospheric effects. There is an increasing focus on utilizing distributed scatterers (DSs) to extract geophysical parameters of interest (i.e. LOS deformation and DEM error) for surfaces characterized by fields, forests, soil and rock surfaces. Application examples include mines, volcanoes, oil/gas/water reservoirs and seismic zones. Distributing scattering mechanism involves a coherent sum of many independent small scatterers (no dominant scatterer) within a resolution cell (Goodman, 1976).
Techniques such as the Small Baseline Subset Algorithm (SBAS) and SqueeSAR have been proposed to process DSs. SBAS makes use of spatially unwrapped small baseline differential interferograms (Berardino et al., 2002, Mora et al., 2003. These are linked using the Singular Value Decomposition (SVD) method and a minimum-norm least squares (LS) solution is obtained (Golub and Loan, 1996). Essentially, the phase is averaged in an estimation window to reduce the phase noise (Zebker and Villasenor, 1992). However, the drawback of the straightforward rectangular estimation window is a loss of resolution and a superposition of different objects on ground. In addition, phase unwrapping (Constantini, 1998, Eineder et al., 1998 is an important step in SBAS and phase unwrapping errors are often encountered in natural terrains. There might be several decorrelated areas (e.g. trees, soil, water etc.) separating the coherent patches and the relative values in the different coherent patches can have unknown integer multiples of π 2 phase offsets.
SqueeSAR, on the other hand, makes use of all possible interferograms and employs an adaptive spatial multilooking to reduce the phase noise and estimate the covariance matrix for each DS (Ferretti et al., 2011, Zan et al., 2008. Afterwards, a phase triangulation algorithm is applied to each DS to retrieve the optimized phase values for the SAR images. The DSs are then processed jointly with the PSs using the traditional PS interferometric chain. However, this technique can be computationally expensive. With respect to the above-mentioned techniques, we have developed and implemented an alternative method for high resolution deformation monitoring of DSs. The proposed method performs an object adaptive parameter estimation which is based on two principles. First, the high resolution of satellites such as TerraSAR-X leads to many resolution cells covering a homogenous object in non-urban areas. This object can be described by a single deformation parameter. The object's area (pixels) thus needs to be identified for an optimal estimation of the model parameter, namely, the LOS deformation velocity. Second, we just concentrate on small baseline differential interferograms to reduce the effects of topography on the DSs, and mainly to reduce the computational complexity.
Practically, the algorithm involves the identification of the object's area (pixels) by a similarity test algorithm using a stack of SAR amplitude images (Parizzi and Brcic, 2011). The differential interferometric phase values of the object's pixels (calibrated with respect to a reference pixel which lies within the object boundary, so that, the atmospheric effects and orbital errors are cancelled out) are then exploited for parameter estimation using a search algorithm in the solution space. For each object patch, a periodogram approach is applied in the spatial domain (object area) and time domain (small baseline differential interferometric stack) for a robust estimation of the local tilts of displacement velocity and local slopes of residual DEM in range and azimuth directions with respect to the reference pixel. Finally, since the independent estimated neighboring patches are close and deformation is assumed to be smooth, a 2D-model-deformation integration can be performed to get the absolute deformation. The new concept with respect to the existing algorithms is that there is no need for spatial phase unwrapping and the deformation time series is estimated at full resolution. Even in the presence of high phase noise, the algorithm compensates DEM errors and additionally, atmospheric artifacts are removed automatically for each patch as the phase values are referred to one pixel within the patch. Objective of this paper is to present the developed algorithm and demonstrate it using TerraSAR-X high resolution spotlight data of Lueneburg in Germany.

METHODOLOGY
Assuming that we have N SAR images and M single look small baseline differential interferograms available, the implemented methodology for deformation estimation of DSs involves the following steps:

Identification of Homogenous Patches
We start with the identification of independent homogenous patches using a stack of coregistered and calibrated SAR amplitude images. These statistically homogenous patches are dependent on the reflectivity, have typically a constant local slope and clear boundaries (e.g. fields, roads etc.).
Various statistical tests have been proposed in recent years to identify homogenous pixels based on the amplitude of coregistered and calibrated stack of SAR images. These include the Kullback-Leibler Divergence test (Bishop, 2006), the Kolmogorov-Smirnov (KS) test (Papoulis and Pillai, 1984) (used in the SqueeSAR approach (Ferretti et al., 2011)) and the Anderson-Darling (AD) test (Pettitt, 1976). The AD test has been proven to be the most effective statistical test to identify if two pixels arise from the same distribution (Parizzi and Brcic, 2011). It is a non-parametric test i.e. we do not assume that the samples belong to a defined probability distribution. Instead, using the amplitude of the stack of SAR images, we obtain the empirical cumulative distribution functions of amplitudes for the two pixels (points) under consideration. The distance between the distributions, with weighting given to tails (higher order moments of the distribution are taken into consideration), tells us if the two points statistically arise from the same distribution. For a set of points a and b , the AD statistic 2 A is: where, N is the number of SAR amplitude images, are empirical cumulative distribution functions of amplitudes x for points a and b , If the AD statistic is less than a threshold value, the two pixels are assumed to belong to the same homogenous area.
In this step, we divide the area into rectangular blocks. The block size is chosen in such a way that the atmosphere is mitigated if we subtract the phase of the reference pixel from the phases of the other pixels in the block, and also considering that the patch size should be large enough to provide us with a reliable estimation. Then, within each rectangular block, homogenous patch pixels are identified based on a certain criterion which gives us the best possible estimate of the DS inside the block. In practice, the estimation takes advantage of large homogenous areas. It is based on the fact that the large number of samples improves the precision and the large spatial extension increases the sensitivity of the estimation. This is the reason that we consider a minimum patch size. Also, the average coherence of the patch should be larger than a certain threshold. The coherence of the pixels can be calculated using an adaptive multilooking algorithm Brcic, 2011, Goel andAdam, 2011a). For each detected patch, a reference pixel is selected. Fig. 1 shows an example of the identification of homogenous patches for a small region in the town of Lueneburg in Germany. 17 TerraSAR-X images (high resolution spotlight mode) of the test site from 2010-2011 were used. In this demonstration, the region was divided into blocks of 40 × 40 pixels. A minimum patch size of 20 pixels and a coherence threshold of 0.3 was applied.

Tilt and Slope Estimation
The single look small baseline differential interferograms can now be exploited for estimation of tilts of deformation velocity and slopes of residual DEM for each patch with respect to the reference pixel of that patch. First, the phase values of each patch are corrected with respect to the reference pixel so that the atmospheric effects and orbital errors are negligible. Second, these are then used for parameter estimation using a 4D periodogram.
For each object patch, to estimate the local tilts of deformation velocity where, L is the number of homogenous pixels for the patch, (3) where, λ is the transmitted wavelength, k t B is the temporal baseline for interferogram k . The height conversion factor for an interferogram k is given by: where, k B ⊥ is the perpendicular baseline for interferogram k , R is the sensor-target distance, θ is the local incidence angle (for flat terrain).
Further on, we make the parameter estimation robust by averaging all the periodograms for a patch to reduce side lobes.
We finally get an averaged periodogram ξ as follows: The local tilts of the deformation velocity

Network Inversion
Since the independent estimated neighboring patches are close and deformation is assumed to be smooth, a 2D-modeldeformation integration can finally be performed to get the absolute deformation. A simple 1D case is shown in Fig. 2.

APPLICATION TEST CASE AND RESULTS
The town of Lueneburg in Germany has been used to demonstrate our technique. It is situated in the German state of Lower Saxony. The old part of this town lies on a salt dome. As a result of constant salt mining dating back to the 19th century, various areas of the town have experienced a gradual or high subsidence, became unstable and had to be demolished. The sinking still continues even today. Many ground stations have been established since 1946 to monitor the deformation, but due to the changing subsidence patterns and locations, space-borne differential SAR interferometric technique is better suited for deformation mapping of Lueneburg (Goel et al., 2011).
We used 17 high resolution spotlight TerraSAR-X images of Lueneburg from October, 2010 to September, 2011 with a look angle of 29.6 degrees and 'HH' polarization. We generated 89 small baseline differential interferograms. Fig. 3 shows the identification of homogenous patches for a part of Lueneburg which has deformed highly in the considered time period. In Fig. 4, we present deformation estimation results for this region using a different approach so that we can compare the new technique with them. The approach which we used for comparison is object adaptive phase filtering, followed by an L1 norm based SBAS technique (Goel and Adam, 2011b). Fig. 5 shows the preliminary tilt estimation results for the deformation velocity in range and azimuth directions in mm/year/pixel using the newly developed technique. We can see that the results compare well with the SBAS results. Fig. 6 shows the preliminary slope estimation results for the residual DEM in range and azimuth directions in m/pixel. The comparison of the results with ground truth data is foreseen in the future. We are in contact with the local government of Lueneburg for the levelling data.

CONCLUSION
A new concept has been developed to estimate deformation of DSs at full resolution without any phase unwrapping. Even in the presence of high phase noise, the algorithm compensates DEM errors and additionally, atmospheric artifacts are removed automatically for each patch as the phase values are referred to one pixel within the patch. A demonstration has been provided using TerraSAR-X data of Lueneburg in Germany. Future work will concentrate on network inversion algorithms to retrieve the absolute LOS deformation velocities and residual DEM.