ADAPTIVE SPATIAL FILTERING OF INTERFEROMETRIC DATA STACK ORIENTED TO DISTRIBUTED SCATTERERS

Standard interferometry poses a challenge in non-urban areas due to temporal and spatial decorrelation of the radar signal, where there is high signal noise. Techniques such as Small Baseline Subset Algorithm (SBAS) have been proposed to make use of multiple interferometric combinations to alleviate the problem. However, the interferograms used in SBAS are multilooked with a boxcar (rectangle) filter to reduce phase noise, resulting in a loss of resolution and signal superstition from different objects. In this paper, we proposed a modified adaptive spatial filtering algorithm for accurate estimation of interferogram and coherence without resolution loss even in rural areas, to better support the deformation monitoring with time series interferometric synthetic aperture radar (InSAR) technique. The implemented method identifies the statistically homogenous pixels in a neighbourhood based on the goodness-of-fit test, and then applies an adaptive spatial filtering of interferograms. Three statistical tests for the identification of distributed targets will be presented, applied to real data. PALSAR data of the yellow river delta in China is used for demonstrating the effectiveness of this algorithm in rural areas. * Corresponding author. Tel.: +86 010 64838047 E-mail address: zhangyj@irsa.ac.cn (Y. Zhang), xiechou@irsa.ac.cn (C. Xie).


INTRODUCTION
InSAR is a microwave remote sensing technique that uses satellite images to measure surface deformation over large areas with millimetre precision (Rosen, Hensley et al., 2000).With the availability of various synthetic aperture radar (SAR) images and continual observations over the same area, long time deformation can be extracted by using stacking techniques such as Persistent Scatterer Interferometry (PSI) (Ferretti, Prati et al., 2000;2001) and Small Baseline Subset Algorithm SBAS (Berardino, Fornaro et al., 2002;Mora, Lanari et al., 2002).PSI processes differential interferograms with respect to a common master image, aims to identify coherent targets exhibiting high phase stability during the whole time period of observation.These phase stable points, slightly affected by temporal and geometrical decorrelation, are called persistent scatterers (PSs).They often correspond to one or two dominant scatterers in a resolution cell and are typically characterized by high reflectivity values created by dihedral reflection or simple single-bounce scattering (Perissin and Ferretti, 2007).PSI is powerful in urban areas where there are a lot of man-made structures which behave as good PS.While in non-urban areas, characterized by vegetated or low reflectivity homogenous regions, such as cropland, volcanoes, mines, reservoirs, the spatial density of PSs extracted by PSI is low (< 10 PS/sqkm), which is a key limitation for its application in rural areas (Ferretti, Fumagalli et al., 2011).Therefore, the deformation exploiting of distributed scatterers (DSs, or Gaussian scatterers) is gathering more and more attention.DS belongs to areas of moderate coherence in some interferometric pairs of stack, where there is sufficiently high number of random small scatterers within a resolution cell with no dominant scatterer, and follows the complex circular Gaussian distribution (Bamler and Hartl, 1998).Various techniques, such as SBAS and SqueeSAR (Ferretti, Fumagalli et al., 2011) have been proposed to process DS, which are widespread in rural areas.
SBAS implements a combination of the multilooked differential interferograms, which are properly selected in order to minimize the temporal and spatial baselines, hence extenuating the decorrelation phenomena (Lanari, Casu et al., 2007).The multilooking filters the interferogram with a rectangle window to increase the Signal-to-Noise Ratio (SNR), and also brings out several drawbacks.First is the lower resolution.Secondly, superposition of different objects mixes various radar signals from scatterers with different properties, which results in severer temporal decorrelation.In addition, undulating terrain and abrupt changes in deformation are smoothed out, which cause spatial decorrelation.All of these break the stationary random condition of Gaussian distribution and lead to inaccurate estimation of interferogram and coherence (Bamler and Hartl, 1998;Jiang, Ding et al., 2013).While, new technique has been developed to select scatterers with stable signals and to compensate the topographic phase (Goel and Adam, 2012).
In contrast, SqueeSAR detects DSs based on their amplitude statistics (adaptive spatial filtering), exploits all possible interferograms and estimates the optimal wrapped phase of each DS with phase triangulation algorithm (temporal filtering), which is based on the error theory.Then the DSs are processed with conventional PSI chain jointly with the PSs.However, the Kolmogorov-Smirnov test used in SqueeSAR is too simple and not suitable for the long tail distribution which applies in radar signal.The linear model assumption is not always the case, and improper for the highly non-linear deformation.Besides, the phase triangulation algorithm can be computationally expensive (Ferretti, Fumagalli et al., 2011), which should be taken into consideration in practice.
International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XL-7/W1, 3rd ISPRS IWIDF 2013, 20 -22 August 2013, Antu, Jilin Province, PR China The main objective of this work is to find a proper method to identify and make use of DSs, which behave as statistically homogenous pixels (SHP) under the assumption that their radar signal and geophysical parameters of interest are linear correlated.The paper is organized as follows.Section 2 describes the statistics theory for SHP identification, and the algorithm for adaptive spatial phase filtering.Experimental results for a test scene located at Yellow River Delta and the comparison of different statistics are drawn in Section 3. Finally, we wrap-up our findings and propose further steps in the conclusion.

METHODOLOGY
Suppose we have a stack of N complex SAR images coregistered to subpixel accuracy and M single look differential interferograms available.The adaptive spatial filtering involves the following steps:

Identification of statistically homogenous pixels
For each pixel, we wish to identify the surrounding pixels which present similar statistical behaviour.Assuming the process is stationary and continuous over time, we can get N observations by sampling the process temporally, to test the degree of similarity between pixels.Since the radar signal of distributed scatterers follows the complex circular Gaussian distribution with unknown parameters according to the central limit theorem, the identification algorithm detects SHP by testing if the two corresponding random processes belong to the same distribution in statistics, where it is referred to a nonparametric goodness-of-fit (GOF) test, and many methods have been developed (Conover, 1980).Then the problem can be defined as a null hypothesis test, 0 : p q H F F  , that the two distribution F p and F q are the same, versus the alternative, 1 : p q H F F  , that they are not (Parizzi and Brcic, 2011).
Considering the uncertainty of the radar signal's distribution, we choose three statistical tests based on the amplitude of coregistered and calibrated stack of SAR images, namely: the Kolmogorov-Smirnov (KS) test, the Cramer-von Mises (CVM) test and the Anderson-Darling (AD) test.These tests are nonparametric, i.e. the samples are not assumed to follow any defined probability distribution.

Kolmogorov-Smirnov Test
The Kolmogorov-Smirnov test is the first goodness-of-fit test for general distribution (Kvam and Vidakovic, 2007;Stephens, 1970).For two-sample version with the same sample size, the test statistics D N is defined as the maximum vertical distance between the cumulative distribution functions (cdfs): max ( )-( ) where () p Fx and () q Fx are the empirical cdfs of the amplitude at pixel p and q, respectively.The null hypotheses H 0 will be rejected at significant level  if the test statistics D N exceeds K  , where K  is the 1- quantile of the two-sided KS test for two samples of equal size N (Conover, 1980).

Cramer-von Mises Test
The Cramer-von Mises test is another nonparametric GOF test.It measures the weighted distance between the empirical cdfs based on a squared-error function (Conover, 1980;Kvam and Vidakovic, 2007).In our two samples with equal size case, the test statistics is defined as where () p Ri and () q Ri are the ranks, in the combined ordered sample, of the ith smallest of the amplitudes at point p and q, respectively.Reject H 0 at the significant level if 2 N  exceeds the 1- quantile 1 w   , the quantiles based on asymptotic distribution are given by Anderson and Darling (Anderson and Darling, 1952), the exact quantiles for 8 N  are given by Burr (Burr, 1963).

Anderson-Darling Test
Anderson-Darling test is a specification of CVM test with the weighted functional . Compared to KS test, AD test puts more weight on the tails of the distributions, which is very important in radar application where the distribution's tail plays an important role, leading to a lower rate of the second kind of error in the hypotheses (Kvam and Vidakovic, 2007).The test statistics is defined as A  of the two-sample AD test.The asymptotic distribution was given by Anderson (Anderson and Darling, 1954) for N , while the quantiles' approximations for finite sample size are given by Pettitt (Pettitt, 1976).
Based on the GOF tests, SHP was selected to mitigate the temporal decorrelation phenomenon.However, the radar returns of distributed scatterers are still affected by atmospheric phase delay, which is spatially smooth with a correlation length of 1~3km (Hanssen, 2001).It is necessary to limit the scale of the identified SHP, according to the geometric resolution of the image.
The identification algorithm of SHP was designed as follows: 1) For each pixel P, define an estimation window centred on P.
2) Apply the two-sample GOF test between every pixel within the estimation window and centre pixel P at a given level of significance; select all the pixels that can be considered as statistically homogenous.3) Abandon the statistically homogenous image pixels that are not connected to P directly or through other SHP.4) Connect SHP to pixel P for the pixelwise post processing, such as interferometric phase filtering and coherence estimation.
International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XL-7/W1, 3rd ISPRS IWIDF 2013, 20 -22 August 2013, Antu, Jilin Province, PR China  of earth between rice-fields, which shows higher coherence and tendency of coherence loss affected by temporal decorrelation.

Improved interferogram and coherence estimation
Based on the SHP identified in the step above, we estimate the M filtered interferograms and their coherence.This step was performed pixelwise with an adaptive filter by using the SHP (5) For identified SHP of pixel P, we estimated its covariance matrix C(P) given by (Ferretti, Fumagalli et al., 2011): where H indicates Hermitian conjugation and T indicates transposition, () dP is the N temporal complex SAR data vector, () i dP is the complex reflectivity value of the ith image on pixel P.
The covariance matrix C transferred to a coherence matrix  after the normalization with amplitude data.The off-diagonal elements of  are an estimation of the complex coherence value for all possible interferograms of data stack, while the phase values of the matrix correspond to the adaptive filtered interferometric phase , jk  (Ferretti, Fumagalli et al., 2011):

RESULTS AND DISCUSSION
To assess the improvements related to the new algorithm, we selected Yellow River Delta located in northeast Shandong Province of China as test site.Yellow River Delta is a typical fan-shaped delta, with natural slope 1/8000~1/12000.The widespread tidal flats and seasonal inundated wetlands with 73.6% coverage of reeds made it a unique estuary ecosystem (Xie, Shao et al., 2013).The test site is visualized in Landsat TM image in Fig. 4a.13 ALOS PALSAR images in fine mode acquired from June, 2007 to October, 2009 with an incident angle of 38.73 °and HH polarization were used.The wavelength of PALSAR is 23.6 cm, longer than X or C band microwave, resulting in a higher penetration depth, which have advantage in rural areas.An area of approximately 20 km×25 km has been processed.The image resolutions in range and azimuth directions are 9.39 m and 3.14 m respectively.
The data stack was pre-processed with Gamma software.All algorithms have been implemented in Matlab for the experiment.Fig. 1 shows the identification results of SHP (in red) for the point (in green) based on three goodness-of-fit tests according to Eq. ( 1), ( 2), (3).For a detect window of 23×23 pixels (in black), 106 pixels were selected as SHP based on KS test and CVM test, as shown in Fig. 1b and Fig. 1c, 105 pixels were selected as SHP based on AD test, as shown in Fig. 1d.We can see that the identified SHPs lay on the same paddy field, and there is no significant difference among these GOF tests' identification results.The coherence matrix on two identified DSs are illustrated in Fig. 2, calculated according to Eq. ( 6), ( 7), (8).It depicts all the possible combinations of the acquisitions we have, with color coded according to its average coherence.Compared to paddy field, earth bank shows higher average coherence, which are mainly affected by temporal decorrelation and present a seasonal behaviour due to annual planting and irrigation.The matrix can be a good reference for interferometric data pairs' combination.
In order to test the effectiveness of algorithm, an area distributed with small paddy fields and citizen buildings was chosen and adaptive spatial filtering algorithm was performed according to Eq. ( 4), ( 5), in comparison with the conventional boxcar multilooking.Fig. 3a illustrates the Landsat TM image of this area and the 1×5 multilooked interferogram formed with the data acquisitions of August 13 th , 2013 and September 28 th , 2013.Since there was rice field irrigation with water drawn from Yellow River between the acquisitions dates, water level changes along the azimuth direction and interferometric fringe was detected within the single rice field.Fig. 3b is the coherence and interferogram after conventional 4×20 pixels boxcar multilooking.Fig. 3c is the coherence and interferogram after the 4×20 pixels adaptive spatial filtering based on AD test.We can clearly distinguish the rice fields from the bank of earth between them (on the upper part of image); building's outline was retained with sharp edge, along with the road connected to it (on the lower part of image).And the object resolution was preserved.
In Fig. 4 and Fig. 5, the full area results are presented.3114121 pixels have been selected as SHP, account for 92.68% of all pixels in the site.Fig. 4a is the corresponding Landsat TM image; Fig. 4b shows the coherence map, with an average value of 0.2805, higher than the average coherence after boxcar multilooking 0.2676.Because of the pixelwise filtering, the computational complexity is increased with the estimation window size.There is a slight reduction of resolution due to the adaptive filtering, but the object resolution was preserved.

CONCLUSION
In this paper, we developed an adaptive spatial filtering algorithm oriented to the distributed scatterers, which are widely distributed in rural areas.The identification of statistically homogenous pixels by using three different goodness-of-fit tests on real SAR data proved its effectiveness.
Compared with the typical rectangle multilooking, this new technique improves the interferogram and coherence estimation with an adaptive spatial filter, preserves the detailed information and increased the average coherence.
Future work will be devoted to the optimization of interferometric pairs' combination and more robust and computationally cheap DS identification technique.

Fx
is the empirical cdf of the pooled distribution of the two samples.Like the CVM test, the null hypotheses will be rejected at level  if test statistics

Fig. 1 .
Fig. 1.Identification results of statistically homogenous pixels (SHP) (in red) for the green pixel with a detection window of 23×23 pixels by goodness-of-fit test.(a) is 1×5 multilooked interferogram enclosed in black rectangle, (b) is a zoom-in of the black rectangle with SHP identified by Kolmogorov-Smirnov test, (c) Cramer-von Mises test, (c) Anderson-Darling test.

Fig. 2 .
Fig. 2. Two examples of Coherence matrix of the PALSAR data stack on the tested distributed scatterers in Yellow River Delta, China.It shows the average coherence of all interferometric combination.(a) is a rice field, shows low coherence, (b) is the low bank of earth between rice-fields, which shows higher coherence and tendency of coherence loss affected by temporal decorrelation.
identified.The adaptive filtered interferogram value , () jk IP from the jth SAR image S j and the kth SAR image S k on pixel point P was estimated as(Goel and Adam, 2012): the conjugation,  is the set of SHP identified and ref  is the reference phase, i.e. the topographic phase, the orbit phase, etc.
Fig. 3. Comparison of boxcar multilooking and adaptive spatial filtering algorithm, applied on the interferometric data pair acquired on August 13 th , 2013 and September 28 th , 2013.(a) is Landsat TM image and its 1×5 multilooked interferogram, (b) is coherence and interferogram estimate after 4×40 pixels boxcar multilooking, (c) is coherence and interferogram estimate after 4×20 pixels adaptive spatial filtering based on Anderson-Darling test.The vertical interferometric phase change within the paddy field was caused by irrigation.

Fig. 4 .
Fig. 4. (a) is Landsat TM image of Yellow River Delta, (b) is the coherence after the 4×20 pixels adaptive spatial filtering based on Anderson-Darling test