Interferometric SAR Phase Denoising Using Proximity-Based K-SVD Technique

This paper addresses the problem of interferometric noise reduction in Synthetic Aperture Radar (SAR) interferometry based on sparse and redundant representations over a trained dictionary. The idea is to use a Proximity-based K-SVD (ProK-SVD) algorithm on interferometric data for obtaining a suitable dictionary, in order to extract the phase image content effectively. We implemented this strategy on both simulated as well as real interferometric data for the validation of our approach. For synthetic data, three different training dictionaries have been compared, namely, a dictionary extracted from the data, a dictionary obtained by a uniform random distribution in [−π,π], and a dictionary built from discrete cosine transform. Further, a similar strategy plan has been applied to real interferograms. We used interferometric data of various SAR sensors, including low resolution C-band ERS/ENVISAT, medium L-band ALOS, and high resolution X-band COSMO-SkyMed, all over an area of Mt. Etna, Italy. Both on simulated and real interferometric phase images, the proposed approach shows significant noise reduction within the fringe pattern, without any considerable loss of useful information.


Introduction
Interferometric Synthetic Aperture Radar (InSAR) [1][2][3][4][5][6][7][8][9] is a consolidated remote sensing technique with broad applications in the field of Earth and environmental sciences. In the last two decades, InSAR has been playing a significant role for measuring several geophysical quantities, including land topography, surface deformations, land changes, water levels, ocean currents, soil moisture, glacier dynamics and vegetation properties [10]. Basically, InSAR uses two coherent SAR images to form an interferogram that can be acquired either from two different antennas on the same platform or from different passes of the same antenna at different times. In general, SAR interferograms are affected by several decorrelation effects, depending on different noise sources, which collectively produce interferogram phase noise. Decorrelation stems from system noise, processing errors and other internal and external factors (e.g., atmospheric fluctuations) [6,9,11]. In the last decade, several techniques have been proposed in the literature for getting rid of the decorrelation effects in the interferometric phase noise. Non-adaptive filtering methods, including the mean filtering technique proposed by Rosen [9], are not so effective for InSAR interferograms. In fact the adoption of some fixed windows for filtering can induce phase distortions due to the periodic character of interferograms not being considered. Another well known filtering method is the adaptive noise filtering proposed by Lee [12]. This method uses suitably selected windows, whose orientations better fit the fringes. Although it has an advantage

Fundamental of Denoising Problem in Interferometry
To set the stage, let us start by considering one single SAR interferogram I, obtained from two co-registered SAR acquisitions S 1 and S 2 in InSAR data processing. Let φ 1 and φ 2 represent the phase information of the acquired images S 1 and S 2 , such that: where φ m = φ 1 − φ 2 represents the interferometric phase [8,51]. In general, measurements of the interferometric phase are both noisy (due to various decorrelation effects such as sensor noise, temporal and geometrical fluctuations), and affected by undeterminate additive multiples of 2π (periodic functions, wrapped phase). Consequently retrieving the noise free φ m from noisy data is a very difficult inverse problem [40]. Given a phase φ ∈ , the corresponding wrapped interferometric phase φ 2π can be expressed as and (a mod b) is the remainder of a/b. The main goal of the interferometric phase estimation problem is to figure out the 2D phase map φ 2π from the observed 2D map φ z given by where ν represents noise [1][2][3][4][5][6][7][8][9]52]. Our aim is evaluating an estimateφ 2π of φ 2π from φ z , capitalizing on the fact that the original noise-free phase image should be smooth within the fringes. We propose here a modified K-SVD algorithm for interferometric phase denoising, named Proximity-based K-SVD (ProK-SVD).
The conventional K-SVD algorithm has been originally formulated in connection with (digitized) optical image denoising [34]. The semantic and structure of differential interferograms is quite different and hence the performance of K-SVD denoising is not obvious, as discussed below.

The Proximity-Based K-SVD Method
For the sake of the casual reader we summarize here the basics of sparse representations, and the rationale of the conventional K-SVD algorithm. Signal processing techniques for denoising problems require that the chosen representation should efficiently separate signal and noise. Representing a signal translates into the choice of a dictionary, a set of elementary signals or atoms [21]. Orthogonal dictionaries (bases) have been widely used due to their mathematical simplicity and general applicability. However, orthogonal dictionaries trade generality for compactness of representation. This led to the development of new classes of (overcomplete) dictionaries, which allow to represent specific classes of signal in more compact way. Let us consider the dictionary D = [d 1 d 2 · · · d L ] ∈ N×L , where the columns constitute the dictionary atoms, and L ≥ N. Representing a signal x ∈ N , using this dictionary, can follow of two alternative paths: either the analysis path, where the signal is represented via its inner products with the atoms, γ a = D T x (5) or the synthesis path, where it is represented as a linear combination of the atoms, In the general case (where the dictionary is not a basis), analysis and synthesis of a signal may differ very much. If D is overcomplete, the family of representations satisfying (6) is actually infinitely large and we can seek the most informative representation of the signal with respect to some cost function C(γ): In the present context, C(γ) will promote the sparsity of the representation. Then, the dictionary is updated assuming known and fixed coefficients. Given a set of samples X = [x 1 x 2 , · · · x n ], the goal of sparse representation is to find a dictionary D and a sparse matrix Γ which minimize the representation error, where {γ i } represents the columns of Γ , and the l 0 sparsity measure · 0 counts the number of non-zeros in the representation, and · F is the Frobenius distance [32,46,47]. The resulting optimization problem is combinatorial and highly non-convex. K-SVD [34] finds a numerical solution to this optimization problem rather than using matrix inversion for dictionary update, changing atom-by-atom via a simple and efficient process. Applying the original K-SVD algorithm proposed by [34] to differential radar interferogram denoising requires further modification of the algorithm.
Following a local approach in [34] and by using the efficient code developed by Rubinstein in [53] modified for our purposes as we shall explain better later, the original two dimensional N × N interferogram is reduced to a one dimensional array whose elements are patches of size n × n, with n N to be scanned sequentially, with partial overlap (see Figure 1 ), and local sparsification on each patch is applied. Mathematically this can be described bŷ argmin φ 2π (i) (9) whereφ 2π (i) for i = 1, . . . , h represents the sparsification of the ith patch viz.: with a dictionary (matrix) of size D(i) ∈ n 2 ·k (with k n 2 ), andα(i) defined such that with φ z (i) being the ith (noisy) patch. To set the patch dimension Q we estimate the number of principal components of the original image. This latter is partitioned into a variable number of patches, and in each patch the number of Principal Component Analysis (henceforth PCA) is computed. Remarkably, irrespective of the partitioning, each patch can be considered as a realization of the stochastic process characterizing the interferogram structure. To estimate the optimal PC-patch dimension we perform an unsupervised training Principal Component Analysis on a simulated interferogram dataset shown in Figure 2.  The dataset consists of 64 blocks of 4096 elements each. The normalized mean square error between the simulated noise-free data and its patched Q-components approximant for varying Q are shown in Figure 3. The typical PCA behaviour is observed namely the presence of a knee in the curve for Q ≈ 60. Accordingly, we will use 64 principal components to describe each block. The next step is to reduce each patch to a one dimensional array, as in [53]. However, instead of scanning the patch in columns, appending columns one after another, in order to preserve spatial correlation of the data, we introduce a proximity concept, assuming that proximity implies similarity. Each patch is accordingly scanned as exemplified in Figure 4.
Log-lin plot of normalized mean square error when approximating the no ipal components.  We will show that this choice will give better performance than the K-SVD column-by-column scanning. K-SVD is applied to each patch in two steps [54]: (1) a block-coordinate minimization algorithm and (2) a search of optimalα i . Orthonormal matching pursuit [32,54] is used, selecting one atom at a time, and stopping when the error Dα − φ z (i) 2 2 goes below a fixed threshold. Given allα(i), α is updated. The block diagram of the whole processing chain of the proposed ProK-SVD algorithm for SAR interferogram denoising is shown in Figure 5. An initial dictionary is selected to start the process of phase denoising on a single wrapped interferogram from a stack of data. Next the ProK-SVD method is applied in two steps, i.e., sparse representation and dictionary update. When the algorithm meets the preset threshold T, the process stops and the denoised phase map is obtained.
In the next section numerical experiments on simulated as well as real data are described and commented. Figure 5. Block diagram of ProK-SVD SAR interferogram denoising algorithm. In the first step an unsupervised PCA analysis is implemented aimed at decomposing the original image into P patches of size n × n. In the second step each patch X k is re-arranged in a 1D array using proximity, an initial dictionary D (0) k is chosen, and K-SVD algorithm is run (up to a suitable maximum number of iteration) to obtain a denoised version of the patchX k , and the related (optimal) dictionary D k . Finally, all denoised patches are re-assembled to form the denoised (full) interferogramX.

Figure 5.
Block diagram of ProK-SVD Synthetic Aperture Radar (SAR) interferogram denoising algorithm. In the first step an unsupervised PCA analysis is implemented aimed at decomposing the original image into P patches of size n × n. In the second step each patch X k is re-arranged into a 1D array using proximity, an initial dictionary D (0) k is chosen, and the K-SVD algorithm is run (up to a suitable maximum number of iteration) to obtain a denoised version of the patchX k , and the related (optimal) dictionary D k . Finally, all denoised patches are re-assembled to form the denoised (full) interferogramX.

Results and Discussion
We discuss here the performance of our ProK-SVD algorithm, using different simulated as well as real interferometric data (provided by CNR-IREA: the ALOS and COSMO-SkyMed data in the frame of the MED-SUV project (http://med-suv.eu/); the ENVISAT and ERS data in the frame of the ASI, DCP and MIUR project "A multidisciplinary study on the preparatory phases of an earthquake" (http://www.irea.cnr.it/en/index.php?option=com_k2&view=item&id=545: a-multidisciplinary-study-on-the-preparatory-phases-of-an-earthquake&Itemid=166).).

Simulated Data
For interferogram simulation, we follow the procedure described in [55], using two SAR acquisitions of the same area, and a Digital Elevation Model (DEM). Specifically, we use a Shuttle Radar Topography Mission (SRTM) DEM of three arcseconds (i.e., 90 m spatial resolution), and two ERS-sensor data of the city of Rome, acquired on December 1st 1996 and on 9 June 1996. In Figure 6a we show a 100 × 100 close-up of the whole simulated 512 × 512 interferogram.
In Figure 6b we added zero mean white Gaussian noise with a standard deviation of 0.5 to the simulated fringe pattern. For simplicity we confine our investigation here to additive Gaussian noise, but the proposed method does not rely on any specific assumption about the noise statistics. For producing Figure 6, the initial dictionaries D k were chosen as the columns of the X k patch. Figure 7 shows a close-up of Figure 6 comparing visually the denoising performance of K-SVD and ProK-SVD.
The final dictionaries {D k } are shown in Figure 8, together with some close-ups whereby the effect of proximity can be qualitatively grasped.  For each (n × n) patch we may estimate the local Peak Signal-to-Noise Ratio (PSNR k ) [56] PSNR k = 20 · log 10 max(X k ) X k beeing the true andX k the estimated phase. The map of the (local) PSNR k for K-SVD and ProK-SVD is displayed in Figure 9 for the simulated interferogram shown in Figure 6. The related (empirical) distribution functions of the PSNR k for the K-SVD and ProK-SVD denoised interferograms are shown in Figure 10. It is seen that proximity boosts the performance of K-SVD in the range [5][6][7][8][9][10] dB.  The average of the quantities (12) over all patches making up our test-image is given in Table 1. Overall, the proximity strategy helps extracting the structure of the original phase map in a more faithful way. Accordingly, we will adopt the ProK-SVD in our further experiments on real data. Different choices of the initial dictionaries D The ODCT dictionary has a strong energy compaction property, tending to concentrate the signal features in a few low-frequency components. The RD taken from a uniform distribution in the interval [−π, π] provides a structure-free playground. The ODCT and RD do not depend on the content and noise level of the image. The DD dictionary on the other hand is built from the noise corrupted interferogram.
The performance of the above dictionaries for different noise levels is compared in Figure 11 for the considered simulated interferogram in terms of the average (over patches) PSNR. It is seen that the ODCT and RD behave almost equivalently better than the DD dictionary. A visual comparison of the denoised interferograms obtained using the DD, ODCT, and RD initial dictionaries, for various added noise levels, is shown in Figure 12. It is visually evident that for low to moderate noise levels the random dictionary provides the best (smoothest, more faithful) reconstruction. For any fixed dictionary the reconstruction becomes worse as the noise level is increased, in particular loosing contrast.

212
On the basis of the above, we finally illustrate the performance of ProK-SVD on real SAR 213 interferometric data, using the RD initial dictionary.

214
On purpose we use interferometric data pairs from various SAR platforms over the area of the 215 Etna Volcano (Italy), namely archived data from the ERS, ENVISAT, ALOS and COSMO-SkyMed SAR 216 sensors, with varying spatial and temporal baselines summarized in Table 2. Data having smaller 217 Figure 12. On top the original noisy-free map is shown. The second row displays the original phase map corrupted by additive noise with given standard deviation. The subsequent rows describe the reconstruction obtained using ProK-SVD for different initial dictionary D.
The choice of the dictionary has little impact on the number of iterations required in the K-SVD step of the algorithm, as illustrated in Figure 13, in terms of the l 2 reconstruction error in Equation (12)

Real Data
On the basis of the above, we finally illustrate the performance of ProK-SVD on real SAR interferometric data, using the RD initial dictionary.
On purpose we use interferometric data pairs from various SAR platforms over the area of the Etna Volcano (Italy), namely archived data from the ERS, ENVISAT, ALOS, and COSMO-SkyMed SAR sensors, with varying spatial and temporal baselines summarized in Table 2. Data having smaller spatio-temporal baseline, e.g., in our case COSMO-SkyMed and ENVISAT have lower interferometric noise. Conversely data having large spatial temporal baselines, i.e., ALOS and ERS have higher noise. As a first benchmark we consider Goldstein filtering [7], widely used in the conventional SAR denoising techniques, with the power factor α = 0.5 (As shown in Figure 4 that represents the ratio of the energies of the noisy phase map and the energy of the (fiducial) noise 228 removed by the de-noising process.

229
In Fig. 15 we extend the comparison to other denoising techniques 3 . It is noted that the In order to preserve phase circularity, the mean, median, and wavelet filters have been applied to the complex interferogram and then the wrapped filtered phase has been extracted. Assessing the quality of image denoising algorithms in the case where no noise-free version is available is an extensively studied (and still open) issue. Several quality metrics have been proposed in the literature, as discussed e.g., in [57]. We adopt the signal-to-distorsion ratio (SDR), discussed, e.g., in [58], defined as that represents the ratio of the energies of the noisy phase map and the energy of the (fiducial) noise removed by the de-noising process.
In Figure 15 we extend the comparison to other denoising techniques (in order to preserve phase circularity, the mean, median, and wavelet filters have been applied to the complex interferogram and then the wrapped filtered phase has been extracted). It is noted that the ProK-SVD entails minimum smoothing, and yields a higher SDR while producing effective denoising.
In Table 3 the various denoising techniques are accordingly compared in terms of the average SDR. Table 3. Signal-to-Distortion Ratio (SDR) values relevant to ProK-SVD, K-SVD, Goldstein filter with parameter α = 0.5, non-local filter, wavelet filter, median filter, mean filter (ENVISAT SAR sensor data-set).    All zoomed views display significant reduction of the noise level compared to the original data. The ProK-SVD technique is seen to preserve effectively the local features within the fringe pattern without any sensible loss of valuable information, and without introducing any artifacts. This is further illustrated by the values of SDR collected in Table 4. Summing up, the ProK-SVD approach performs well for all SAR sensors considered, in terms of noise suppression and lack of artifacts.

Conclusions
In this paper we addressed the interferometric phase image denoising problem by solving the sparse and redundant representations problem over a trained dictionary. The key new idea is to apply a proximity-based modified K-SVD algorithm to the noisy interferograms so as to obtain both a sparse representation and an updated dictionary together. This approach, remarkably, does not require any fine tuning of the relevant parameters, nor any a priori information to work reliably. We tested the proposed algorithm on both simulated and real SAR interferometric data, from low resolution (ERS and ENVISAT), medium resolution (ALOS), and high resolution (COSMO-SkyMed) data. We discussed the choice of the initial dictionaries, referring to three particular cases, namely, random, ODCT and data-driven, and evaluated their performances on simulated data with additive (Gaussian) noise with varying sigma values. The random dictionary was found to yield the best performance. The proposed technique was capable of effectively retrieving the fringe pattern from the noise interferograms, without introducing significant artifacts in a wide range of signal to noise ratios.
As far as computational complexity and burden are concerned, assuming a dictionary of dimension N × L and T iterations (see Equation (8)), the ProK-SVD algorithm is dominated by the K-SVD stage, requiring T 2 L + 2NL floating point operations, as demonstrated in [53]. We plan to implement the algorithm using parallel (GPU) architectures for further optimization.