NEW SAR INTERFEROGRAM DENOISING METHOD VIA SPARSE RECOVERY BASED ON L 0 NORM

The goal of this paper is to estimate a denoised phase image from the observed noisy SAR interferogram. We proposed a linear model to obtain a sparse representation of the interferomteric phase image. The main idea is based on the smoothness property of the phases inside interferometric fringes which leads to get a sparse image when applying the gradient operator twice, along x or y direction, on the interferogram. The new sparse representation of the interferometric phase image allows to transform the denoising problem to an optimization one. So the estimated interferogram is achieved using the approximate message passing algorithm. The proposed approach is validated on different cases of simulated and real interferograms.


INTRODUCTION
The SAR interferogram (InSAR) filtering is a fundamental step before phase unwrapping process.Several filters have been proposed last decades, many of them are spatial filters (Lee et al., 1998), (Baran et al., 2003), (Abdallah and Abdelfattah, 2013) or filters opertaed in wavelet domain (Lopez-Martinez and Fabregas, 2002), (Bian and Mercer, 2011) and (Chang et al., 2011).Since the compressed sensing (CS) technique was invented by Donoho in 2006(Donoho, 2006), it has been used in the SAR images context.Most of the proposed appraoches were intended to the polarimetric SAR data such as (Si et al., 2009), (Lin et al., 2010) and (Li et al., 2012).Other approaches are dedictaed to the object detection in the SAR images (Anitori et al., 2013) and (Ahmad et al., 2013).In the case of the InSAR denoising, few filtering algorithms are proposed and they are mostly based on the sparsity context rather than the CS technique.For example, mary et al. proposed a sparse denoising approach of the interferometric data by extracting geometrical features of the images (Mary et al., 2010).In (Hao et al., 2013), the authors used a sparse representation and a generated dictionary to estimate the noise free phase from a noisy InSAR.In this paper, we propose a new approach for interferometric phase image denoising using sparse coding.Our idea is based on the fact that the original phase image without noise should be smooth within the fringes.This property allows us to obtain a sparse representation of the image to estimate if we apply the gradient operator.The elements of the result image will be close to 0 except the pixels in the edges of fringes.Next, the filtered interferometric image can be estimated by solving the l0 minimization problem of the obtained sparse representation.This minimization problem is achieved using the approximate message passing algorithm (Vila and Schniter, 2012).The rest of this paper is organised as follows.Section 2 is dedicated to present the proposed sparse model and the formulation of the optimization problem, section 3 presents the experimental results and discussions and we finish by concluding in section 4.

THE PROPOSED SPARSE MODEL FOR INTERFEROGRAM DENOISING
In this section we present a model of the interferometric phase image based on a sparse representation.The observed noisy interferogram ϕ can be considered as a M × N image as following (Lee et al., 1998): where φ denotes the interferometric phase image without noise and n is the noise.The denoising problem aims to estimate the phase φ from the observed one ϕ.Our goal is to rewrite the expression of φ with a sparse representation.For this reason, we assume that the phase image φ should be smooth inside the fringes i.e the difference between two neighbours phases in the same fringe should be small.In fact, the phase stationarity in homogeneous areas of the interferogram is ensured by the removal of the orbital component from the phase images before the generation of the interferometric phase image (Deledalle et al., 2011).But, in case of almost reliefs, the variation of the phase within a given interferometric fringe could be important and our hypothesis is not necessary valid.For this reason, and in order to insure the applicability of the sparse representation we apply the gradient operator twice.So we use the gradient operators twice along x or y direction to compute the phase difference between neighbours pixels as follows: We obtain the second order derivative (∂ 2 ) of the interferometric phase image φ using a simple linear formulation: This contribution has been peer-reviewed.doi:10.5194/isprsarchives-XL-7-37-2014where G is a M × M matrix having the following structure: An alternative of this sparse representation could be the Laplacien of the interferometric phase.However, the linear expression of that is very complicated.So, as mentioned above Sx and Sy can be considered as sparse matrices whose almost of their elements are close to 0 except the pixels located at the edges of the fringes.Since the treatments along two directions are symmetric to each other, in the rest of this paper we will interested only on y direction.So now we can express the unknown phase image φ with a sparse representation using Γ = G −1 which is expressed as a lower triangular matrix: The next step is to construct the vectors ϕ col and φ col containing respectively the elements of ϕ and φ arranged column by column.By the same way, we generate the vectors S col and n col from Sy and n respectively.So the size of these vectors is M N × 1.Also, we define the diagonal matrix Ψ ∈ M N ×M N as following: According to (3) we can write: And the final sparse representation of (1) becomes: According to (7), the estimate of the noise free phase φ is given by φcol = Ψ Ŝcol when Ŝcol is the solution of the following l0 minimization problem with constrain: where S col 0 denotes the norm l0 that gives the number of nonzero elements of s and is a given parameter to control the error rate.
By minimizing the number of nonzeor of S col , the homogeneous regions in the interferogram are filtered and only the phase jumps between different fringes are taken into account.The minimization problem ( 9) is NP-hard meaning that it is very difficult to reach the exact solution.But the problem can be solved using l1 approximation or if we use a greedy algorithm such as the Matching Pursuit (Mallat and Zhang, 1993) or the Basis Pursuit (Chen et al., 2001).In this paper we used a recent greedy algorithm called Approximate Message Passing (AMP) algorithm proposed by Vila et al. (Vila and Schniter, 2012) to solve the l0 minimization problem in (9).The choice of AMP is based on the fact that it is a fast algorithm with high probability to find a sparsest solution (Mohimani et al., 2010).

Results from simulated interferograms
We begin our validation tests with five different simulated interferograms with 512 × 512 pixels size representing respectively a cone, a pyramid, a Gaussian variation and two phase images generated using meshgrid function in MATLAB as shown in Fig. 2. First, we filtered the simulated interferograms affected by a white Gaussian noise with σ 2 = 0.8.Fpr , we assume that the error rate between the estimated phase φ and the observed one ϕ for one pixel should not exceed π 1000 .So for the whole phase image with size N × N the global error rate will be compute as the sum of local error rate note that more details about the relationship between the noise variance and the estimaton error threshold is given in (Hongxing et al., 2013).
For the filtering results, we notice that the proposed approach reduce noise inside fringes very well.It denoise perfectly the homogeneous areas of the interferograms and therefore it gives lower Mean Square Error (MSE) values with respect to noise level as shown in Table 1.Thus, the great advantage of the proposed algorithm is that it preserves the edges of the fringes forms.
As quantitative comparison, we compute the MSE between the original phase image without noise and the filtered interferograms using FAMM (Abdelfattah and Bouzid, 2008), EWMF (Abdallah andAbdelfattah, 2013), WInPF (Lopez-Martinez andFabregas, 2002), NL-InSAR (Deledalle et al., 2011) and the proposed algorithm.Table 1 shows the different MSE values when varying the noise level.We notice that the proposed filtering approach gives, globally, the best MSE values except in few cases.

Results from real interferograms
For the real SAR interferograms, we used four different phase images produced from SLC pairs of Radarsat-2, ERS-2, Envisat and Cosmo-SkyMed satellites respectively.All informations about  1 shows that the proposed approach gives lower error rate for the four real InSAR except in the case of ERS satellite data.

CONCLUSION
The main idea of the interferometric phase filter proposed in this paper is to smooth the homogeneous regions within the interfrogram's fringes.The low phase difference between pixels located in the same fringe of the interferogram leads to obtain a sparse signal when we applying the gradient operator.This transform the denoising problem to a l0 minimization problem which can be solved using greedy algorithm.We used in this paper the Approximate Message Passing (AMP) algorithm.The proposed approach is tested and validated on simulated and real interferograms and compared with two recent filters FAMM and EWMF.
As future work, we will apply this approach to estimate the true unwrapped phase value from the observed noisy interferogram.This will be done by twice applying the gradient operator and the Laplacien simultaneously.

Figure 1 :
Figure 1: Flowchart of the proposed approach for SAR interferometric denoising via sparse recovery.

Figure 2 :
Figure 2: Simulated interferograms filtering results.From top to bottom: cone, pyramid, Gaussian variation, first meshgrid function and second meshgrid function.From left to right: the original interferograms without noise, the corresponding noisy interferograms with σ 2 = 0.8 and the filtered results with the proposed approach.

Figure 3 :
Figure 3: Real interferograms filtering results.The first column using Radarsat-2 interferogram and the second column using ERS satellite and from top to bottom: the original interferogram produced from ASTER DEM, the observed noisy interferogram, filtered with the proposed approach and the filtering error with respect to column 1. these four InSAR are available in (Ben Abdallah and Abdelfattah, 2013).All tested real images are 512 × 512 pixels size.Similar to the simulated data, we consider the interferogram generated by wrapping the Digital Elevation Model (DEM) provided by ASTER Global DEM available in the Land Processes Distributed Active Archive Center website (Land Processes Distributed Active Archive Center, 2013) as the original interferogram without noise.Note that since the four satellites, mentioned above, having different spatial resolutions, we should make under or over sampling to their InSARs to obtain the same resolution of ASTER's interferogram (30 × 30 meters).Next, we compare the filtering results of the four InSARs with the reference one generated from ASTER's DEM.In Fig.3, we illustrate our estimation result and the corresponding error map.For the statistical comparison, the second part of the Table1shows that the proposed approach gives lower error rate for the four real InSAR except in the case of ERS satellite data.
The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, VolumeXL-7, 2014   ISPRS Technical Commission VII Symposium, 29 September -2 October 2014, Istanbul, Turkey

Table 1 :
Mean Square Error (MSE) computed using different filtering algorithms.The first part of table: MSE of simulated interferograms with different noise variances.The second part: MSE between interferograms obtained using ASTER's DEM and the filtered ones.