Adaptive Covariance Matrix Estimation for Multi

IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING

1

Adaptive Covariance Matrix Estimation for Multi-Baseline InSAR Data Stacks Michael Schmitt, Student Member, IEEE, Johannes L. Schönberger, and Uwe Stilla, Senior Member, IEEE

Abstract—For many multidimensional applications of synthetic aperture radar (SAR) imaging, the estimation of the covariance matrix for each resolution cell is a critical processing step. The context of this work is the application of covariance matrix estimation for multi-baseline interferometric SAR data sets. In order to ensure local stationarity, which is needed for an unbiased estimation, adaptive techniques are necessary. In this paper, a new approach for adaptive covariance matrix estimation is proposed and evaluated based on measures known from the field of image processing. The procedure is centered around the idea of checking whether the neighboring pixels belong to the same statistical distribution as the currently investigated pixel by applying a threshold to the respective probability density function. All inlier pixels are then used to estimate the complex covariance matrix of the reference pixel. From this covariance matrix, both amplitude and interferometric phase values are extracted, which are then combined for all pixels in the stack in order to employ techniques for the evaluation of filtering efficiency that are typically used in image denoising research. It is found that the proposed algorithm provides high filtering efficiency and good detail preservation at the same time. Apart from that, it is found to be particularly suitable for small-sized stacks of coregistered SAR imagery. Index Terms—Adaptive filtering, covariance matrix estimation, multilooking, synthetic aperture radar (SAR).

I. I NTRODUCTION

F

OR many multidimensional applications of synthetic aperture radar (SAR) imaging, the estimation of the covariance matrix for each resolution cell is a critical processing step. Examples include SAR polarimetry (PolSAR), SAR interferometry (InSAR), polarimetric SAR interferometry (PolInSAR), or SAR tomography (TomoSAR) [1]–[7]. This paper focuses on stacks of interferometric SAR data as they are used as input to both multi-baseline or multitemporal InSAR applications such as deformation estimation [8] or change detection [9] or multilooking-based TomoSAR. Whereas the task of covariance matrix estimation is not very challenging for low- or mediumresolution data of natural environments and rural terrain, the

Manuscript received September 29, 2013; revised November 7, 2013 and December 17, 2013; accepted January 23, 2014. This work was supported by the German Research Foundation (DFG) under Project STI 545/4-1. M. Schmitt and U. Stilla are with the Department of Photogrammetry and Remote Sensing, Technische Universität München, 80333 Munich, Germany (e-mail: [email protected]; [email protected]). J. L. Schönberger is with the Department of Photogrammetry and Remote Sensing, Technische Universität München, 80333 Munich, Germany and also with the University of North Carolina at Chapel Hill, Chapel Hill, NC 27599 USA (e-mail: [email protected]). Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org. Digital Object Identifier 10.1109/TGRS.2014.2303516

advent of decimeter-resolution sensors has significantly risen the demand for sophisticated processing strategies: While formerly fixed-sized boxcar windows were usually sufficient, the hypothesis of statistical homogeneity of neighboring pixels cannot be considered valid for complex scenes such as urban areas anymore [10], [11]. Therefore, many papers have been published on the topic of adaptive filtering of different kinds of SAR data. Beginning with the problem of simple speckle filtering for intensity imagery [12]–[14], the research has gradually been extended toward interferometric and polarimetric SAR data with the goal of an unbiased estimation of phase, coherence, and polarimetric scattering information [15]–[17]. One of the most efficient filters in this context was proposed by Deledalle et al. [18]: This algorithm, which utilizes a nonlocal estimation framework [19], allows for the simultaneous extraction of all relevant information of a pair of coregistered InSAR images. Although the results are quite promising, however, NL-InSAR’s most severe drawback is that it cannot be used for data stacks containing more than two single-look complex (SLC) SAR images, as the derivation of the algorithm is strictly based on bivariate data sets. As an alternative, more recently, some first papers have been published about the adaptive filtering of multidimensional SAR stacks. One of the first approaches in this context was designed for multitemporal data, which were realized by detection of 3-D adaptive neighborhoods in order to properly consider changes in the scene [20]. The next generation, now aiming at the filtering of multi-baseline InSAR data, was presented by Ferretti et al. [21], who proposed DespecKS, an algorithm embedded in their SqueeSAR framework that uses a twosample Kolmogorov–Smirnov test in order to evaluate if two stack pixels within a predefined search window belong to the same statistical distribution. Parizzi and Brcic [22] further investigated this ansatz with respect to different goodness-of-fit tests such as Kullback–Leibler divergence, Anderson–Darling test, or generalized likelihood ratio test. Finally, most recently, Navarro-Sanchez and Lopez-Sanchez [23] extended the DespecKS concept to polarimetric InSAR data by proposing a likelihood ratio test for testing the equality of complex Wishart matrices. Although all these formulations show promising results, they all suffer from one certain disadvantage: They typically work only for stacks of at least eight images and up [24]. If, however, just standard InSAR pairs or stacks with a limited number of images, e.g., acquired by single-pass multi-baseline systems as they are frequently equipped on airborne platforms, are to be processed, alternative methods are necessary. A first approach for this problem field has been already presented

0196-2892 © 2014 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

2


in [25]. However, like all other mentioned filters, this method also relies on the assumption that the amplitude values of the pixels can be used as a hint for changes in their phase values as well; the determination of statistical similarity is solely based on amplitude information, whereas phase and coherence information is neglected. If an object in the scene causes a backscattering amplitude similar to the background, its phase values will be joined with the phase of this background during filtering, potentially causing blurring or even the loss of the object’s phase information. This paper therefore proposes a new approach for the adaptive estimation of covariance matrices for stacks of multibaseline InSAR data, fulfilling two main requirements: 1) The algorithm is supposed to work for stacks of arbitrary size, i.e., from two conventional InSAR images to large multi-baseline InSAR stacks. 2) The homogeneity of the pixels shall not be determined by analyzing the amplitudes only; instead, the whole complex information is to be considered. The general idea of our method is based on analyzing whether resolution cells belong to the same statistical distribution as the currently investigated center pixel, which is realized by a simple maximization of the probability density function of the center pixel. While we assume circularly symmetric complex Gaussian single-look observations in this context, the authors of a recently published related paper employ a Wishart test in order to evaluate the similarity of several covariance matrices derived from the original measurements [26]. For the assessment of the adaptivity and filtering efficiency of our algorithm, we use amplitude and interferometric phase data, which can be extracted from the complex covariance matrix of each pixel, as a proxy. In this way, evaluation methods known from image processing can be applied conveniently. The tests are carried out on both simulated and real test data. II. A DAPTIVE C OVARIANCE M ATRIX E STIMATION Here, we propose a new approach for an adaptive estimation of the covariance matrices of the pixels in a multi-baseline InSAR data stack. In analogy to many competing methods, the procedure is centered around a sliding window operation (see Fig. 1). This arbitrarily large window is called search window in Fig. 1 and defines the maximum size of the homogeneous pixel neighborhood. In addition to that, a smaller central window, which does not need not be rectangular necessarily, is defined. Using the pixels in this central window, a robust covariance estimator is applied in order to generate an initial guess for the covariance matrix of the respective center pixel from its direct enclosing neighbors. Afterward, the probability density function for each pixel in the search window is evaluated with respect to this initial covariance matrix, and a threshold is applied such that pixels belonging to the distribution of the center pixel are separated from pixels belonging to a different distribution. The adaptivity test is then finalized by a connectivity check in order to ensure that only pixels connected to the center pixel are considered for the neighborhood. Finally, the sample covariance matrix we are looking for is estimated from all connected inliers, i.e., from all pixels that were detected to be part of the homogeneous neighborhood.

Fig. 1. Sketch of the sliding-window-based adaptive filter. Within an arbitrarily large search window around the pixel of interest, all pixels are thresholded whether they belong to the probability distribution defined by the robust covariance matrix initially estimated from the pixels in the central window. A subsequent connectivity check then separates connected inliers from unconnected inliers.

A. Statistical Properties of SLC InSAR Observations The whole procedure is centered around the assumption that each SLC observation z agrees at least approximately with Goodman’s model [27], where both the real and imaginary parts follow a zero-mean Gaussian distribution (∼ N (0, σ)) and are statistically independent. From this, it follows that each pixel vector z ∈ CN in a stack of N coregistered SAR images follows the multivariate probability density function f (z|C) =

1 exp(−zH C−1 z) π N det(C)

(1)

which is fully characterized by its complex covariance matrix C = E{zzH } ∈ CN ×N .

(2)

Although De Zan [28] has investigated this assumption in comparison with a constant plus Gaussian model, which might be considered more appropriate for heterogeneous scenes such as urban areas, where often quasi-deterministic point scattering might occur, he came to the conclusion that (1) is a sufficient approximation. This view is also seconded by other authors working with X-band SAR [6], [29], whereas for smaller wavelengths (e.g., millimeter-wave SAR), the validity of the assumption should be even better [30]. Therefore, our idea is now to check every pixel in the search ˆ init ), where window whether they fit to the distribution N (0, C ˆ Cinit is an initial and robust guess of the covariance matrix of the center pixel’s distribution. B. Robust Estimation of the Initial Covariance Matrix The initial guess for the covariance matrix of the center pixel’s distribution is calculated from its enclosing 4- or 8-neighborhood. However, it is possible—and in heterogeneous areas such as urban scenes even probable—that not all resolution cells in this neighborhood really belong to a statistically homogeneous population of backscattering observations. ˆ init is employed in order to Hence, a robust estimation of C

SCHMITT et al.: ADAPTIVE COVARIANCE MATRIX ESTIMATION

3

minimize the potential bias. An additional advantage of this robust initialization is that the effect of single pixels affected by strong speckle is mitigated. The utilized M -estimator, which has been proposed in [31] and has been recently discussed in further detail in [32], is iteratively defined as ˆ init,k+1 = C

1 Linit

L init

H ˆ −1 w zH l Cinit,k zl zl zl

(3)

l=1

where Linit is the initial number of looks, i.e., the size of ˆ init,0 is the standard sample the initial neighborhood, and C covariance matrix estimated from the classic 3 × 3 boxcar window. Since Linit needs to be larger than N in order to ensure ˆ init,0 is nonsingular and can be inverted in (3), a larger that C initial neighborhood has to be employed for stacks with more than nine images. w(x) is a robust weighting function, which is descending to zero, such that a highly deviating observation zl ˆ −(1/2) zl 2 = zH C ˆ −1 zl receives smaller weights with large C l init init in the estimation. The general idea behind this robust estimation procedure is the assumption that the random vectors z in the local window are assumed to follow a complex multivariate t-distribution Ctk,ν (μ, C) with ν > 0 degrees of freedom, i.e., 2k+ν s − 2 (4) f (z|μ, C) = c|C|−1 1 + 2 ν where s = (z − μ)H C−1 (z − μ), and c is a normalizing constant. With ν = 1, this distribution is called the multivariate complex Cauchy distribution, which is a prominent robust heavy-tailed alternative for the Gaussian distribution, which is obtained for ν → ∞. Due to this reason, the complex multivariate t-distributions are very useful for analyzing the robustness of multivariate statistics, since a decrement of ν yields a distribution with an increased heaviness of the tails. Therefore, the weighting function is chosen by 2N + ν . w(x) = wν (x) = ν + 2x

(5)

Usually, the parameter is selected in the range 1 ≤ ν ≤ 5 if the goal is to ensure a certain robustness against outliers and noisy data. As a stopping criterion for the iterative estimation of (3) ˆ init,k+1 − C ˆ init,k < ε C

(6)

is chosen, where ε > 0 is a small number. The benefit achieved by relying on a robust estimation of the initial covariance matrix instead of just the standard sample estimate is illustrated in Fig. 2. C. Similarity and Connectivity Testing After a robust initial guess of the distribution’s covariance matrix has been attained, the statistical similarity to the corresponding central population is determined for each pixel within the larger search window. A stack pixel z is assumed to be part of the distribution of the central pixel if a certain threshold on the probability density is exceeded, i.e., ˆ init > εPDF . (7) f z|C

Fig. 2. Comparison of conventional and robust covariance matrix estimation. The Frobenius norm of the difference between the estimated covariance matrix and the theoretical covariance matrix is plotted with respect to a growing percentage of outlier observations. The benefit of robust covariance matrix estimation is clearly visible.

The choice of this threshold is dependent on the noise level and the dimensionality of the observed data and is further investigated in Section III-A1. After applying the thresholding, a final 4- or 8-connectivity test is carried out in order to ensure that all inliers are spatially connected to the central pixel. This step is carried out in order to enhance the probability that all pixels of the adaptively chosen neighborhood really belong to an area of homogeneous backscattering with the same statistical characteristics. The main advantage of this inspection is that the search window can, in principle, be chosen arbitrarily large because spatial similarity, i.e., all pixels really belonging to the same backscattering phenomenon, is ensured by the pixel connectivity. D. Final Covariance Matrix Estimation After a set of homogeneous pixels has been determined by the procedure described in the previous sections, the final complex covariance matrix is ultimately estimated from all inlier samples by ˆ = 1 zl zH C l L L

(8)

l=1

where L now is the final number of looks. Due to the predefined fixed size of the search window, L cannot exceed the number of pixels inside this search window. The most interesting characteristic of the complex covariance matrix of InSAR stacks is that the covariance information is directly related to the denoised interferometric measurements, i.e., √ √ ⎡ γ12 I1 I2 . . . γ1N √I1 IN ⎤ I1 √ ∗ I1 I2 I2 . . . γ2N I2 IN ⎥ ⎢ γ12 ⎥ . (9) ⎢ C=⎣ .. .. .. .. ⎦ . . . . √ √ ∗ ∗ IN γ1N I1 IN γ2N I2 IN . . . Ii denote the (despeckled) intensities of all N acquisitions in the stack, whereas E{zi zj∗ } = |γij | exp (jφij ) γij = Ii Ij

(10)

is the complex coherence, which itself is composed of the magnitude of coherence (or correlation) |γij | and the interferometric phase φij between acquisitions i and j.

4


This means that all relevant information contained in an InSAR data stack can be deduced just by estimating the covariance matrix for every pixel. This information contains despeckled intensities, all possible interferogram combinations, as well as all related coherence magnitudes. Therefore, the feasibility of the adaptivity checking procedure proposed in this paper will be proven by evaluating the resulting intensity, phase, and coherence images in an image processing manner.

III. E XPERIMENTAL R ESULTS For the evaluation of the adaptive covariance matrix estimation procedure proposed in this paper, experiments on both simulated and real test data are carried out. In order to get an impression about both filtering efficiency and adaptivity, the amplitude map of the master image and the interferogram corresponding to the longest available baseline combination are extracted from the covariance matrices of all pixels in the stack, as described in Section II-D. Therefore, evaluation methods known from the field of conventional image denoising can be employed for evaluation.

Fig. 3. Optimal choice of degrees of freedom ν for the robust initial covariance matrix estimation and optimal choice of probability density threshold εPDF for similarity checking. (a) Plot for fixed stack size (two images) with respect to growing noise level. (b) Plot for fixed noise level (STD = 0.5) with respect to growing number of images in the stack.

A. Tests on Simulated InSAR Data The tests on simulated InSAR data aim at an empirical determination of optimal parameter settings for the algorithm, as well as a comparison with state-of-the-art InSAR filters under laboratory conditions in order to evaluate the theoretical filtering efficiency with respect to the noise level and the number of images in the stack. The latter evaluation is carried out in terms of both the full covariance matrices and image information extracted from them, enabling the utilization of image processing methods for the assessment. The simulated data are synthesized from a true set of amplitude and phase images according to the multiplicative speckle noise model discussed, e.g., in [33]. The primary noise level is assumed to be equal for the real and imaginary parts, resulting in Rayleigh distributed amplitude and Gaussian distributed interferometric phase. Finally, in addition to this SAR inherent speckle, the signal is augmented with Gaussian distributed thermal noise of the SAR sensor. Both speckle and thermal noise are added in order to describe the overall noise level used for the following evaluations. 1) Determination of Optimal Parameter Settings: From Section II-B and C, it is obvious that the proposed algorithm depends on two main parameters, which have to be tuned manually: 1) the degrees of freedom ν for robust estimation of the initial covariance matrix and 2) the probability density threshold εPDF for determination whether a pixel belongs to the same distribution as the center pixel. Fig. 3 shows two plots analyzing the optimal parameter choice with respect to the noise standard deviation and the number of images in the stack. Obviously, ν does not depend on the stack size but only on the noise level and is optimally chosen between 2.5 and 3.9, which corresponds to [31], which suggested 1 ≤ ν ≤ 5. Due to the small variations, the choice of ν is not very critical; if the

Fig. 4. Filtering efficiency comparison with respect to noise standard deviation. The normalized bias between estimated and true covariance matrices averaged over all pixels is plotted. The number of images in the stack was set to N = 2.

Fig. 5. Filtering efficiency comparison at a constant noise level and with respect to growing stack size. Again, the normalized bias between estimated and true covariance matrices averaged over all pixels is plotted.

noise level is not known a priori, we suggest ν = 3 as a good compromise. The decadic logarithm of εPDF stays at approximately −10 independently of the noise level, whereas it slowly decreases for larger stacks. However, also the choice of the threshold is not very critical. Therefore, in the remainder of this paper, we employ a constant threshold of εPDF = 10−10 . If large


Fig. 6.

5

Filtering results for a simulated InSAR pair. The boxes around numbers 1, 2, and 3 are enlarged for better interpretability.

stacks, such as typically encountered in the context of repeatpass spaceborne SAR, are to be processed by the presented approach, the threshold should be lowered accordingly. Considering these empirical findings, the method proposed in this paper does not necessarily require manual settings or parameter tunings and can be applied to different kinds of InSAR data with predefined default parameters. Only if the highest precision is required and, e.g., very large stacks are to be processed or very large noise levels are present, a refinement of the parameters can be advisable. 2) Comparison With Other Filtering Techniques: In the presented experiments, the proposed estimation procedure is compared with four specifically chosen state-of-the-art multilooking procedures: • Boxcar multilooking: Boxcar filtering is the common reference procedure that has been efficiently used for lowand medium-resolution SAR and InSAR data for decades [34]. Its application to multi-baseline stacks is straightforward. For the experiments in this paper, a 15 × 15 window was employed. • DespecKS: Being the first multilooking approach specifically designed for stacks of SAR imagery, this method can be considered a direct benchmark [21]. However, its filter-

Fig. 7. Quantitative comparison of filtering results on a simulated InSAR pair. (a) Amplitude and (b) interferometric phase data. For the evaluation of the phase, shadow regions have been excluded.

ing capacities break down for stacks smaller than about eight acquisitions. It is based on a goodness-of-fit test evaluating the statistical similarity between the amplitude

6


Fig. 8. Filtering results for a simulated multi-baseline stack with four images. The boxes around numbers 1, 2, and 3 are enlarged for better interpretability.

values of two neighboring stack pixels. In addition, for DespecKS, a 15 × 15 search window was used. • NL-InSAR: Although NL-InSAR is designed for just single pairs of InSAR data,1 it can be considered the most efficient multilooking algorithm at the moment [18]. The MATLAB implementation provided at http://www.math.u-bordeaux1.fr/~cdeledal/nlinsar is therefore used for comparison. Please note that, although NL-InSAR is not meant to be a window-based algorithm (instead, all pixels are supposed to be considered in a nonlocal manner), in the provided implementation, a 21 × 21 window is used due to computational reasons. • PCA-TV: Recently, we have already proposed a different approach for the filtering of multi-baseline InSAR data with the special motivation to provide a tool for small stacks consisting of only three to five acquisitions, as they are frequently encountered with airborne single-pass missions [25]. With this method, a hitherto existing gap in the literature has been tackled for the first time. Just as DespecKS, the PCA-TV filter employs amplitude values, which are preprocessed by principal component analysis and total variation denoising, as hint for pixel similarity. Again, a 15 × 15 search window was employed during the experiments for this paper. 1 Only recently, a new preprint extending the NL-InSAR principle to multimodal data has come to our knowledge [35], which could not be considered for the investigations in this paper.

Filtering Efficiency: First, the theoretical estimation efficiency was investigated by Monte Carlo experiments on the covariance matrices of simulated InSAR data pairs. For the proposed probabilistic similarity determination, the optimal parameter settings as determined in Section III-A1 (εPDF = 10−10 , ν = 3) were used. The mean normalized bias ˆ − C0 F /N C

(11)

ˆ to the true covariance of the estimated covariance matrix C matrix C0 with respect to a growing noise level is shown in Fig. 4. In this context, · F denotes the Frobenius norm, i.e.,

n

m AF := |aij |2 (12) i=1 j=1

where m and n are the number of rows and columns of A = [aij ], and N is the number of images in the stack. While the bias reduction of the proposed method is close the bias reduction offered by NL-InSAR and PCA-TV, Boxcar filtering and DespecKS show comparably lower efficiency. This is caused by the lack of adaptivity in the Boxcar case and the fact that DespecKS requires large sample numbers for the goodness-of-fit testing. In order to prove the efficiency also for larger stacks with more than two images, additional tests have been carried out on simulated stacks of different sizes using a constant noise


7

level. The corresponding plot is shown in Fig. 5. Note that in this case, NL-InSAR is not included due to its restriction to the two-image case. Again, PCA-TV and probabilistic similarity determination show comparably high efficiency (i.e., strong bias reduction), whereas Boxcar filtering and DespecKS are again limited by nonadaptivity and reliance on large sample numbers, respectively. 3) Image-Based Evaluation: Qualitative results achieved by extracting amplitude and phase values (corresponding to the longest available baseline) for every pixel in the stack are visualized in Fig. 6. The corresponding quantitative results can be found in Fig. 7. The results of Boxcar filtering, DespecKS, and NL-InSAR are as expected: The Boxcar approach provides high filtering efficiency but completely blurs all image details. DespecKS is not able to filter the data since two samples per pixel are too few for statistical goodness-of-fit testing. Eventually, NL-InSAR provides the benchmark result. Both PCA-TV and the newly proposed probabilistic filter provide strong filtering efficiency and detail preservation. However, the PCA-TV result is disturbed by single dark spots and an unsatisfactory result for the phase of letter 3, which is probably caused by the fact that two images do not allow for a sufficient separation between signal and noise components. In contrast, the probabilistic filter leads to a result comparable with NL-InSAR in quality. Since the actual goal of this paper is to propose a filtering technique for stacks of multi-baseline InSAR data rather than just standard InSAR pairs, the evaluation has been additionally carried out on simulated stacks containing four images. Note that NL-InSAR is not part of the experiment anymore since it is restricted to conventional InSAR data. Apart from that, the results in Figs. 8 and 9 lead to an impression similar to the previously investigated two-image case: While Boxcar filtering blurs the results, DespecKS still suffers from the low number of images. The PCA-TV result shows fewer dark spot outliers, but the proposed probabilistic approach still creates the best results considering the desired tradeoff between filtering efficiency and detail preservation.

Fig. 9. Quantitative comparison of filtering results on a simulated stack consisting of four images. (a) Amplitude and (b) interferometric phase data. For the evaluation of the phase, shadow regions have been excluded. TABLE I T EST DATA S YSTEM PARAMETERS

B. Tests on Real InSAR Data In spite of all simulated test results, in the end, the quality of a filter can only be evaluated based on real data experiments. For the investigations in this section, test data acquired by the German MEMPHIS sensor in 2011, courtesy of the Fraunhofer Institute for High Frequency Physics and Radar Techniques (FHR), have been investigated [36]. For detailed sensor specifications, see Table I. 1) Filtering Efficiency: The crucial capability of any filter is high filtering efficiency, i.e., a combination of variance reduction and mean preservation. Since these measures, which are commonly estimated for homogeneous image patches, tend to provide a biased view in some cases, we refer to the speckle suppression index (SSI) and the speckle suppression and mean preservation index (SMPI), as proposed in [37], in order to get a more reliable evaluation. The SSI is calculated by SSI =

σ f μ0 · μf σ 0

(13)

Fig. 10. Quantitative comparison of filtering results for a four-image MEMPHIS stack. (a) Amplitude and (b) interferometric phase data.

where μ0 , μf , σ0 , and σf denote the means and the standard deviations of the original and the filtered images, respectively. If the filter is efficient in speckle reduction, the SSI is usually less than 1.

8


Fig. 11. Filtering results for a multi-baseline stack with four images, acquired by the airborne MEMPHIS sensor. The interferometric phase data of the longest baseline are shown.

Unfortunately, the SSI tends to be unreliable if the filter overestimates the mean value. Therefore, the more sophisticated SMPI, which is calculated by SMPI = (1 + |μ0 − μf |) ·

σf σ0

(14)

is used additionally. The lower the SMPI values are, the more efficient the filter performs. Fig. 10 shows the standard deviation (STD), SSI, and SMPI for amplitude and phase data extracted from the complex covariance matrices that were estimated for each pixel. Again, NL-InSAR was not considered. It can be easily observed that the real data results agree with the results achieved with simulated data: Both PCA-TV and the newly proposed probabilistic filter achieve filtering efficiency

comparable with the optimal Boxcar method. The relatively bad result for DespecKS is again caused by the low number of images in the stack, providing only four samples per pixel for the goodness-of-fit test. 2) Filter Adaptivity: In order to ensure local stationarity for the estimation of the covariance matrix of a pixel in an InSAR stack, all sample pixels must belong to the same statistical distribution. Therefore, adaptive pixel selection is needed. Unfortunately, this adaptivity can hardly be evaluated quantitatively. Therefore, we choose to provide a visual comparison of both amplitude and phase images for the assessment. In this way, the adaptivity of the proposed filtering method can be analyzed by its capability to preserve fine image details such as edges or strong point scatterers. Fig. 11 shows the corresponding


9

Fig. 13. Mean and standard deviation of the coherence magnitudes within a small homogeneous window on the lawn area in the real test data set.

IV. D ISCUSSION

Fig. 12. Coherence maps estimated for the longest baseline of the real test data set.

comparison for the different filtering approaches. As expected, it can be seen that the Boxcar approach completely blurs the image and destroys fine details, which proves the lack of adaptivity. DespecKS, PCA-TV, and the probabilistic filter all preserve the exemplary chosen facade scatterers in a satisfying manner. Looking at the phase images, however, it becomes obvious that the PCA-TV approach tends to stronger filtering and less detail preservation in comparison with the probabilistic filter. 3) Coherence Estimation: A very important task for interferometric SAR applications or change detection tasks is the estimation of the coherence between the individual acquisitions. As shown in Section II-D, also this measure can be easily extracted from the complex covariance matrix of each pixel. In order to prove the effect of the proposed adaptive estimation procedure on coherence estimation, the results for the longest baseline available in the real test data set are shown in Fig. 12, whereas the mean and standard deviations of the coherence magnitudes within the same homogeneous window already used for amplitude and phase evaluation in Section III-B1 are displayed in Fig. 13. It can be seen that PCA-TV and the probabilistic filter achieve almost the same coherence quality as the Boxcar filter while still preserving finestructured details in the data.

From both the results of the simulated and real data experiments, several insights can be acquired. They are summarized in Table II. It is obvious that most of the filtering methods provide advantages and disadvantages at the same time: While Boxcar filtering generally shows very high filtering efficiency, it naturally destroys fine details and is not adaptive in any way. NL-InSAR is probably the most adaptive and efficient filter known at the moment, but it is not able to process multi-baseline InSAR stacks with more than two images. In analogy to that, DespecKS is a very promising approach for stack filtering, but it requires a rather high number of acquisitions in order to make the goodness-of-fit testing work. PCA-TV filtering closes the gap between NL-InSAR and DespecKS, providing an efficient and adaptive filtering technique also for small-sized stacks with about three and more images. However, its results are weaker in case covariance matrices for only a standard interferometric image pair; for only a standard interferometric image pair have to be estimated; this is caused by the need to compress signal components in the first principal component as a preprocessing step. In contrast, the newly proposed probabilistic filter provides high filtering efficiency and strong detail preservation independent of the stack size while still being optimally suitable for small-sized stack cases. This is caused by the filter principle, which is centered around the multidimensional probability density function of each pixel. Therefore, the general design is independent of the number of samples per pixel. Apart from that, the fact that the two necessary parameters (ν for the robust initial covariance matrix estimation and εPDF for the probability density thresholding) can be set globally for a certain data configuration makes it a powerful and still easy to use tool. Nevertheless, a more sophisticated fully automatic choice of these parameters completely independent of framework conditions such as stack size or dynamic range of the images should be investigated in the future. In addition, also another drawback of the method has to be mentioned: Since it is based on a sliding window operation and a three-stage covariance matrix estimation process, it is rather computationally expensive. V. C ONCLUSION In this paper, a new procedure for the adaptive estimation of complex covariance matrices of the resolution cells in multibaseline InSAR data sets has been proposed. By application of a thresholding operation on the probability density function

10


TABLE II F ILTER E VALUATION

defined by an initial robust guess of the reference pixel’s covariance matrix for each neighboring pixel, a set of statistically homogeneous inlier pixels is determined. From these, the final covariance matrix can be estimated. Experiments on both simulated and real test data were able to show both the high filtering efficiency and the excellent detail preservation capabilities of the proposed algorithm. The main advantages of this probabilistic filter is that both amplitude and phase values are used for the homogeneity check and that it is applicable to InSAR stacks of arbitrary size, starting from conventional two-image InSAR pairs up to large repeat-pass data collections.

ACKNOWLEDGMENT The authors would like to thank Dr. S. Stanko and T. Brehm of the Fraunhofer Institute for High Frequency Physics and Radar Techniques, as well as Dr. E. Meier and C. Magnard of Remote Sensing Laboratories, University of Zurich, for providing the MEMPHIS test data. R EFERENCES [1] D. Just and R. Bamler, “Phase statistics of interferograms with applications to synthetic aperture radar,” Appl. Opt., vol. 33, no. 20, pp. 4361– 4368, Jul. 1994. [2] R. J. A. Tough, D. Blacknell, and S. Quegan, “A statistical description of polarimetric and interferometric synthetic aperture radar data,” Proc. Math. Phys. Sci., vol. 449, no. 1937, pp. 567–589, Jun. 1995. [3] S. R. Cloude and K. P. Papathanassiou, “Polarimetric SAR interferometry,” IEEE Trans. Geosci. Remote Sens., vol. 36, no. 5, pp. 1551–1565, Sep. 1998. [4] C. Lopez-Martinez and X. Fabregas, “Polarimetric SAR speckle noise model,” IEEE Trans. Geosci. Remote Sens., vol. 41, no. 10, pp. 2232– 2242, Oct. 2003. [5] F. Gini and F. Lombardini, “Multibaseline cross-track SAR interferometry: A signal processing perspective,” IEEE Aerosp. Electron. Syst. Mag., vol. 20, no. 8, pp. 71–93, Aug. 2005. [6] F. Baselice, A. Budillon, G. Ferraioli, and V. Pascazio, “Layover solution in SAR imaging: A statistical approach,” IEEE Geosci. Remote Sens. Lett., vol. 6, no. 3, pp. 577–581, Jul. 2009. [7] J.-S. Lee and E. Pottier, Polarimetric Radar Imaging—From Basics to Applications. Boca Raton, FL, USA: CRC Press, 2009. [8] B. Pinel-Puyssegur, R. Michel, and J.-P. Avouac, “Multi-link InSAR time series: Enhancement of a wrapped interferometric database,” IEEE J. Sel. Topics Appl. Earth Observ. Remote Sens., vol. 5, no. 3, pp. 784–794, Jun. 2012. [9] R. Zandona Schneider and D. Fernandes, “Entropy concept for change detection in multitemporal SAR images,” in Proc. 4th Eur. Conf. Synthetic Aperture Radar, 2002, pp. 221–224. [10] R. Touzi, A. Lopes, J. Bruniquel, and P. Vachon, “Unbiased estimation of the coherence from multi-look SAR data,” in Proc. IEEE Int. Geosci. Remote Sens. Symp., 1996, vol. 1, pp. 662–664. [11] G. Gao, “Statistical modeling of SAR images: A survey,” Sensors, vol. 10, no. 1, pp. 775–795, Jan. 2010.

[12] J.-S. Lee, “Digital image enhancement and noise filtering by use of local statistics,” IEEE Trans. Pattern Anal. Mach. Intell., vol. PAMI-2, no. 2, pp. 165–168, Mar. 1980. [13] V. S. Frost, J. A. Stiles, K. Shanmugan, and J. C. Holtzman, “A model for radar images and its application to adaptive digital filtering of multiplicative noise,” IEEE Trans. Pattern Anal. Mach. Intell., vol. PAMI-4, no. 2, pp. 157–166, Mar. 1982. [14] D. T. Kuan, A. A. Sawchuk, T. C. Strand, and P. Chavel, “Adaptive noise smoothing filter for images with signal-dependent noise,” IEEE Trans. Pattern Anal. Mach. Intell., vol. PAMI-7, no. 2, pp. 165–177, Mar. 1985. [15] J.-S. Lee, K. Papathanassiou, T. Ainsworth, M. Grunes, and A. Reigber, “A new technique for noise filtering of SAR interferometric phase images,” IEEE Trans. Geosci. Remote Sens., vol. 36, no. 5, pp. 1456–1465, Sep. 1998. [16] G. Vasile, E. Trouve, M. Ciuc, and V. Buzuloiu, “General adaptiveneighborhood technique for improving synthetic aperture radar interferometric coherence estimation,” J. Opt. Soc. Amer. A, Opt. Image Sci., vol. 21, no. 8, pp. 1455–1464, Aug. 2004. [17] G. Vasile, E. Trouve, J.-S. Lee, and V. Buzuloiu, “Intensity-driven adaptive-neighborhood technique for polarimetric and interferometric SAR parameters estimation,” IEEE Trans. Geosci. Remote Sens., vol. 44, no. 6, pp. 1609–1621, Jun. 2006. [18] C.-A. Deledalle, L. Denis, and F. Tupin, “NL-InSAR: Nonlocal interferogram estimation,” IEEE Trans. Geosci. Remote Sens., vol. 49, no. 4, pp. 1441–1452, Apr. 2011. [19] X. Yang and D. Clausi, “Structure-preserving speckle reduction of SAR images using nonlocal means filters,” in Proc. 16th IEEE Int. Conf. Image Process., 2009, pp. 2985–2988. [20] M. Ciuc, P. Bolon, E. Trouve, V. Buzuloiu, and J. Rudant, “Adaptiveneighborhood speckle removal in multitemporal synthetic aperture radar images,” Appl. Opt., vol. 40, no. 32, pp. 5954–5966, Nov. 2001. [21] A. Ferretti, A. Fumagalli, F. Novali, C. Prati, F. Rocca, and A. Rucci, “A new algorithm for processing interferometric data-stacks: SqueeSAR,” IEEE Trans. Geosci. Remote Sens., vol. 49, no. 9, pp. 3460–3470, Sep. 2011. [22] A. Parizzi and R. Brcic, “Adaptive InSAR stack multilooking exploiting amplitude statistics: A comparison between different techniques and practical results,” IEEE Geosci. Remote Sens. Lett., vol. 8, no. 3, pp. 441–445, May 2011. [23] V. D. Navarro-Sanchez and J. M. Lopez-Sanchez, “Spatial adaptive speckle filtering driven by temporal polarimetric statistics and its application to PSI,” IEEE Trans. Geosci. Remote Sens., vol. 52, no. 8, pp. 4580– 4589, Aug. 2014. [24] M. Stephens, “Use of the Kolmogorov–Smirnov, Cramer–Von Mises and related statistics without extensive tables,” J. R. Stat. Soc. B, Methodol., vol. 32, no. 1, pp. 115–122, 1970. [25] M. Schmitt and U. Stilla, “Adaptive multilooking of airborne single-pass multi-baseline InSAR stacks,” IEEE Trans. Geosci. Remote Sens., vol. 51, no. 1, pp. 305–312, Jan. 2014. [26] S.-W. Chen, X.-S. Wang, and M. Sato, “PolInSAR complex coherence estimation based on covariance matrix similarity test,” IEEE Trans. Geosci. Remote Sens., vol. 50, no. 11, pp. 4699–4710, Nov. 2012. [27] J. Goodman, “Statistical properties of laser speckle patterns,” in Laser Speckle and Related Phenomena, ser. Topics in Applied Physics, J. Dainty, A. Ennos, M. Francon, J. Goodman, T. McKechnie, G. Parry, and J. Goodman, Eds. Berlin, Germany: Springer-Verlag, 1975, pp. 9–75. [28] F. De Zan, “Optimizing SAR interferometry for decorrelating scatterers,” Ph.D. dissertation, Politecnico di Milano, Milan, Italy, 2008. [29] Y. Wang, X. X. Zhu, and R. Bamler, “Retrieval of phase history parameters from distributed scatterers in urban areas using very high resolution SAR data,” ISPRS J. Photogramm. Remote Sens., vol. 73, pp. 89–99, Sep. 2012.


[30] M. Schmitt, C. Magnard, S. Stanko, C. Ackermann, and U. Stilla, “Advanced high resolution SAR interferometry of urban areas with airborne millimeterwave radar,” Photogramm. Fernerkundung Geoinf., vol. 2013, no. 6, pp. 603–617. [31] E. Ollila and V. Koivunen, “Robust antenna array processing using M-estimators of pseudo-covariance,” in Proc. 14th IEEE Pers., Indoor Mobile Radio Commun., 2003, vol. 3, pp. 2659–2663. [32] A. M. Zoubir, V. Koivunen, Y. Chakhchoukh, and M. Muma, “Robust estimation in signal processing. A tutorial-style treatment of fundamental concepts,” IEEE Signal Process. Mag., vol. 29, no. 4, pp. 61–80, Jul. 2012. [33] J. A. Richards, Remote Sensing With Imaging Radar. Berlin, Germany: Springer-Verlag, 2009. [34] R. F. Hanssen, Radar interferometry—Data Interpretation and Error Analysis. New York, NY, USA: Academic, 2001. [35] C.-A. Deledalle, L. Denis, F. Tupin, A. Reigber, and M. Jäger, NL-SAR: A unified non-local framework for resolution-preserving (Pol)(In)SAR denoising, Tech. Rep. HAL, hal-00844118. [Online]. Available: http:// www.math.u-bordeaux1.fr/~cdeledal/publis_all.php [36] H. Schimpf, H. Essen, S. Boehmsdorff, and T. Brehm, “MEMPHIS—A fully polarimetric experimental radar,” in Proc. IEEE Int. Geosci. Remote Sens. Symp., 2002, pp. 1714–1716. [37] A. Shamsoddini and J. Trinder, “Image texture preservation in speckle noise suppression,” in Proc. Int. Archives Photogramm., Remote Sens. Spatial Inf. Sci., 2010, vol. 38, no. 7A, pp. 239–244.

Michael Schmitt (S’08) was born in Munich, Germany, in 1984. He received the Dipl.-Ing. (Univ.) degree in geodesy and geoinformation from the Technische Universität München (TUM), Munich, Germany, in 2009. He is currently working toward the Ph.D. degree in airborne SAR interferometry using multi-baseline and multi-aspect millimeterwave data in the Department of Photogrammetry and Remote Sensing, TUM, where he has been a full-time research assistant since May 2009. Besides SAR interferometry, his main interests are (statistical) signal processing and parameter estimation theory applied to sensor data fusion and the extraction of information from imagery.

Johannes L. Schönberger was born in Munich, Germany, in 1991. He received the B.Sc. degree in geodesy and geoinformation from the Technische Universität München (TUM), Munich, Germany, in 2013. He is currently working toward the M.Sc. degree in geodesy and geoinformation at TUM and the University of North Carolina, Chapel Hill, NC, USA. In addition to SAR, his main research interests are image processing, parameter estimation theory, and computer vision.

11

Uwe Stilla (M’04–SM’09) was born in Cologne, Germany, in 1957. He received the Dipl.-Ing. degree in electrical engineering from Gesamthochschule Paderborn, Paderborn, Germany, in 1980 and the Dipl.-Ing. degree in biomedical engineering and the Ph.D. degree in the field of pattern recognition from the University of Karlsruhe, Karlsruhe, Germany, in 1987 and 1993, respectively. From 1990 until 2004, he was with the Institute of Optronics and Pattern Recognition (FGAN-FOM), a German research establishment for defense-related studies. Since 2004, he has been a Professor with the Technische Universität München, Munich, Germany, where is currently the Head of the Department of Photogrammetry and Remote Sensing and the Director of the Institute of Photogrammetry and Cartography. He is the Vice Dean of the Faculty of Civil, Geo and Environmental Engineering and the Dean of Student Affairs of the Bachelor’s and Master Program “Geodesy and Geoinformation” and the international Master Programs “Earth Oriented Space Science and Technology (ESPACE)” and “Cartography.” He has been the Chair of the ISPRS Working Group III/VII “Pattern Analysis in Remote Sensing.” He is also the Principal Investigator of the International Graduate School of Science and Engineering (IGSSE), Vice President of the German Society of Photogrammetry, Remote Sensing and Geoinformation (DGPF), member of the Scientific Board of German Commission of Geodesy (DGK), and member of the Commission for Geodesy and Glaciology (KEG) of the Bavarian Academy of Science and Humanities. He has been the Organizer and Chair of the conferences “Photogrammetric Image Analysis (PIA),” “City Models, Roads and Traffic (CMRT),” “GRSS/ISPRS Joint Urban Remote Sensing Event (JURSE 2011),” “Earth Observation and Global Changes (EOGC 2011),” and the “IEEE-GRSS Remote Sensing Summer School (RSSS12).” He has published more than 300 contributions. His research focuses on image analysis in the field of photogrammetry and remote sensing.