Bilateral Distance Based Filtering for Polarimetric SAR Data

Abstract

This paper introduces a non-linear Polarimetric SAR data filtering approach able to preserve the edges and small details of the data. It is based on exploiting the data locality in both, the spatial and the polarimetric domains, in order to avoid mixing heterogeneous samples of the data. A weighted average is performed over a given window favoring pixel values that are close on both domains. The filtering technique is based on a modified bilateral filtering, which is defined in terms of spatial and polarimetric distances. These distances encapsulate all the knowledge in both domains for an adaptation to the data structure. Finally, the proposed technique is employed to process a real RADARSAT-2 dataset.

1. Introduction

Polarimetric Synthetic Aperture Radar (PolSAR) is a multidimensional SAR type based on exploiting the vectorial nature of the electromagnetic waves. Different wave polarization states are employed for the transmitted and the received radar echoes in order to achieve a more comprehensive characterization of the target. Additionally, it allows the exploration of the complete space of polarization states allowing, for instance, optimization procedures [1]. In the last decade, PolSAR has demonstrated its usefulness for Earth surface characterization and study, being able to retrieve some biophysical and geophysical information from the area under observation [2–6].

Usually in SAR systems, the resolution cell dimensions are larger than the wavelength and, consequently, the measured echo is the coherent sum of all the individual target contributions. Since this contribution may be constructive or destructive, for distributed scatterers, the measured echo depends on the spatial distribution and importance of the individual targets within each resolution cell. As a result, homogeneous areas of the scene appear with a grainy appearance on SAR data, which is referred to as speckle. Note that the speckle is a true electromagnetic measure but, due to its complexity, it can only be characterized statistically.

In order to mitigate the effects of the speckle, some filtering is applied to the data. This speckle filtering process is needed for a proper target response estimation and for a reliable information extraction from PolSAR data. A classical Boxcar linear filter may be applied to the data to reduce the effect of the speckle but it results into spatial resolution loss and, additionally, it will not be valid near contours, for instance, due to the mixture of heterogeneous samples. Consequently, a proper filtering technique should employ only homogeneous samples of the data for speckle filtering. However, SAR images are strongly inhomogeneous, as they reflect the complexity of the scene and, therefore, some sort of spatial adaptation is needed. Recently, some state-of-the-art speckle filtering techniques have been developed in order to adapt to the spatial contours of the scene. In [7], this adaptation is performing by selecting the filtering window from a set of predefined oriented windows. On a more recent contribution [8], the region growing technique is employed to find an arbitrary homogeneous neighborhood for each pixel. Similarly, some edge-preserving filtering techniques have been defined for optical image denoising, as the bilateral filter [9] or the non-local filtering [10]. Additionally, some of these techniques have been already adapted to process SAR [11] and PolSAR data [12], respectively. They are based on a weighted average favoring similar pixels. The weights are calculated on a pixel-based approach for the bilateral filtering and in a patch-based approach for the non-local means. Non-local means filtering seems to be more robust to noise but, since it requires an homogeneous neighborhood, it may be suboptimal under some circumstances [13].

These techniques are based on identifying homogeneous samples for each pixel over the image. Note that this information is related with the scene itself, and, usually, we have no access to it, so this information should be inferred from the data. In this case, it is assumed that homogeneous pixels will have a similar radar response. Consequently, a distance or similarity measure has to be defined over the polarimetric space. In fact, this is an aspect in common with a wide number of techniques for PolSAR data processing as, for instance, classification [14] or segmentation [15]. In this paper, a speckle filtering technique is proposed based on a modification of the bilateral filter defined in terms of similarity measures.

This paper is organized as follows: Section 2 introduces the bilateral filtering modification based on distances that is proposed. Section 3 describes the iterative weight refinement process that will be employed to achieve a more reliable estimation of the filter weights. On Section 4 the proposed filtering technique is applied to RADARSAT-2 real data and it is evaluated and compared with other state-of-the-art speckle filtering techniques. Finally, Section 5 presents the conclusions.

2. Bilateral Filter

PolSAR systems measure the complex reflectivity of the scene, collecting the scattering matrix S for each resolution cell [2]. When assuming a monostatic radar configuration, this information may be vectorized onto the target vector k [16]:

$$\mathbf{k}=\left[\begin{array}{c}{S}_{\mathit{hh}}\\ \frac{1}{\sqrt{2}}\left(S_{\mathit{hv}}+S_{\mathit{vh}}\right)\\ S_{\mathit{vv}}\end{array}\right]$$

(1)

where the subscripts denote the transmitted and received polarization, respectively, and h and v stands for horizontal and vertical wave polarization states, respectively.

As mentioned before, the speckle is a true electromagnetic measure but it has to be characterized statistically. If the resolution cell is much larger than the wavelength, and there is no dominant scatterer within, then, by applying the central limit theorem, it may be assumed that the measured reflectivity is following a zero-mean complex Gaussian distribution [17]:

$$p_{\mathbf{k}}(\mathbf{k})=\frac{1}{{\pi }^3\left|\mathbf{C}\right|}\text{exp}\left(-{\mathbf{k}}^{†}{\mathbf{C}}^{-1}\mathbf{k}\right)$$

(2)

where ^† is the complex conjugate transpose of a vector and C represents the covariance matrix.

In this case, the distribution Equation (2) may be completely characterized by the second order moment, which corresponds to the covariance matrix C for multidimensional data. This information may be estimated from data as the average sample covariance matrix, which corresponds to the Maximum Likelihood Estimator (MLE):

$$\mathbf{Z}={\langle {\mathbf{kk}}^{†}\rangle }_n=\frac{1}{n}\sum _{i=1}^n{\mathbf{k}}_i{\mathbf{k}}_i^{†}$$

(3)

where k_i represents the scattering vector of the i-th pixel, n represents the number of pixels averaged and ^† represents the complex conjugate transpose. This estimated covariance matrix Z may be statistically determined by the Wishart distribution [18–20]:

$$p_{\mathbf{z}}(\mathbf{Z})=\frac{n^{3n}{\left|\mathbf{Z}\right|}^{n-3}}{{\left|\mathbf{C}\right|}^n{\tilde{\Gamma }}_3(n)}\text{etr}\left(-n{\mathbf{C}}^{-1}\mathbf{Z}\right)$$

(4)

where etr(X) is the exponential of the matrix trace, being only valid for n ≥ 3, as it will be discussed later, and

$${\tilde{\Gamma }}_3(n)={\pi }^3\prod _{i=1}^3\Gamma (n-i+1)$$

(5)

However, the sample covariance matrix as defined in Equation (3) only makes sense over homogeneous areas. The Boxcar or multilook filter, for instance, assumes spatial locality and defines a window around each pixel assuming that all the samples within the window are homogeneous. The bilateral filter [9] introduces also the domain locality. For PolSAR data one would assume that homogeneous samples are close in both the spatial and the polarimetric domains. Accordingly, the averaging in Equation (3) is replaced by a weighted average depending on the samples closeness in both domains [21,22]. Then, the estimated covariance matrix Ẑ^ij at position (i, j) may be expressed as:

$${\widehat{\mathbf{Z}}}^{\mathit{ij}}=\frac{1}{k(i,j)}\sum _{m,n\in V(i,j)}{\mathbf{Z}}^{\mathit{mn}}{w}_s(i,j,m,n)w_p\left({\mathbf{Z}}^{\mathit{mn}},{\mathbf{Z}}^{\mathit{ij}}\right)$$

(6)

where V (i, j) represents the local window around the pixel located at (i, j), Z^ij represents the estimated covariance matrix of the input image at (i, j), w_s and w_p are the weighting functions on the spatial and polarimetric domains, respectively, and k(i, j) is the normalization factor:

$$k(i,j)=\sum _{m,n\in V(i,j)}{w}_s(i,j,m,n)w_p\left({\mathbf{Z}}^{\mathit{mn}},{\mathbf{Z}}^{\mathit{ij}}\right)$$

(7)

The keystone of the filter, then, is the definition of the w_s and w_p weighting functions. They should be defined in order to exploit the spatial and polarimetric locality of the data, that is, they should fulfill the following properties:

• w_s and w_p are within the interval [0, 1], presenting higher values for closer values in the spatial and polarimetric domains, respectively.
• Consequently, it is assumed that

  $$w_s(i,j,i,j)=1,\forall i,j$$

  (8)

  $$w_p(\mathbf{Z},\mathbf{Z})=1,\forall \mathbf{Z}$$

  (9)
• As symmetry is assumed in data closeness, w_s and w_p should be symmetric, that is,

  $$w_s(i,j,m,n)=w_s(m,n,i,j)\forall i,j,m,n$$

  (10)

  $$w_p\left({\mathbf{Z}}^{\mathit{mn}},{\mathbf{Z}}^{\mathit{ij}}\right)=w_p\left({\mathbf{Z}}^{\mathit{ij}},{\mathbf{Z}}^{\mathit{mn}}\right)\forall i,j,m,n$$

  (11)

Consequently, w_s and w_p should depend on the samples closeness and similarity on the spatial and polarimetric domains, respectively. As mentioned in Section 1, some PolSAR data processing techniques are defined in terms of a polarimetric distance or dissimilarity measure. In this paper, we propose to define the weighting functions w_s and w_p in terms of a spatial and a polarimetric distance d_s and d_p:

$$w_s(i,j,m,n)=\frac{1}{1+\frac{d_s^2(i,j,m,n)}{{\sigma }_s^2}}$$

(12)

$$w_p\left({\mathbf{Z}}^{\mathit{mn}},{\mathbf{Z}}^{\mathit{ij}}\right)=\frac{1}{1+\frac{d_p^2\left({\mathbf{Z}}^{\mathit{mn}},{\mathbf{Z}}^{\mathit{ij}}\right)}{{\sigma }_p^2}}$$

(13)

where σ_s and σ_p control the weights sensitivity in the spatial and polarimetric domains, respectively. Mathematically, d_s and d_p are metric functions that quantify the sample closeness on their domains, following the properties:

• d(x, y) ≥ 0. Non-negativity.
• d(x, y) = 0 if and only if x = y. Identity of indiscernibles.
• d(x, y) = d(y, x). Symmetry.
• d(x, z) ≤ d(x, y) + d(y, z). Triangle inequality.

Different similarity measures on the polarimetric space are defined and analyzed in [15,23]. These measures may be classified on full matrix and diagonal measures, depending on the use of all the elements of the estimated covariance matrix, or just the diagonal ones, respectively. Only full matrix measures employ all the polarimetric information under the Gaussian hypothesis. However, the matrix estimated from original pixels Z = kk^† is singular, having rank equal to 1, and these measures need full-rank matrices. This problem also affects the Wishart distribution Equation (4) since it is not defined for singular matrices as |Z| = 0. Some techniques circumvent this problem through a small initial multilook filtering as a regularization process but it reduces the spatial resolution and blurs small details and contours. The option considered in this paper is to apply a diagonal measure, since they do not have this limitation. For diagonal measures, the off-diagonal elements of the covariance matrix are considered to be 0. This automatically leads to a full-rank matrix avoiding the regularization process and its associated spatial resolution loss. Note that this idea may be applied to obtain a diagonal measure from any full matrix measure. Nevertheless, this approach has the limitation of being insensitive to the off-diagonal information.

In the proposed technique, diagonal measures are considered since it is focused on speckle filtering while also preserving the spatial resolution and small details of the image as much as possible. On the spatial domain the euclidean distance is proposed:

$$d_s^2(i,j,m,n)={\left(i-m\right)}^2+{\left(j-n\right)}^2$$

(14)

On the polarimetric domain two different measures are proposed, one based on the Wishart distribution and other based on the hermitian positive definite matrix cone geometry. The revised Wishart dissimilarity measure [14] is based on a test statistic for the hypothesis that the two covariance matrices belong to the same distribution, with some modifications to obtain a symmetric measure. In this paper, as mentioned before, the diagonal Wishart is employed [15] in order to avoid the need for an initial filtering:

$$d_{\mathit{pw}}^2\left({\mathbf{Z}}^{\mathit{mn}},{\mathbf{Z}}^{\mathit{ij}}\right)=\left(\sum _{k=1}^3\left(\frac{{\left(Z_{\mathit{kk}}^{\mathit{mn}}\right)}^2+{\left(Z_{\mathit{kk}}^{\mathit{ij}}\right)}^2}{Z_{\mathit{kk}}^{\mathit{mn}}{Z}_{\mathit{kk}}^{\mathit{ij}}}\right)-6\right)$$

(15)

where Z_ij is the index notation for the (i, j)th element of the 3 by 3 matrix Z.

The geodesic similarity measure is based on the Riemannian geometry of the hermitian positive definite cone of matrices [24]. For the proposed method, as in the previous case, a modification of the diagonal version is employed [23]. This measure is defined in terms of the natural logarithm of the diagonal elements quotient. This fact reduces considerably the rate of growth of the similarity measure having not enough separability between different matrices to prevent the heterogeneous region mixture. Thus, a modification is proposed to revert this rate of growth reduction caused by the logarithm:

$$d_{\mathit{pg}}^2\left({\mathbf{Z}}^{\mathit{mn}},{\mathbf{Z}}^{\mathit{ij}}\right)=\text{exp}\left(\sqrt{\sum _{k=1}^3{\text{ln}}^2\left(\frac{Z_{\mathit{kk}}^{\mathit{ij}}}{Z_{\mathit{kk}}^{\mathit{mn}}}\right)}\right)-1$$

(16)

where ln represents the natural logarithm.

When dealing with real PolSAR data, consideration must be given to the system noise, caused by thermal noise, discretization errors, etc. [25]. This noise is usually modeled as a zero-mean uncorrelated Gaussian noise, as opposed to the speckle noise, that is following a Wishart distribution. To address this issue, a modification d_pt is proposed to Equations (15) and (16), taking into account the power of the system noise ${\sigma }_t^2$:

$$d_{\mathit{pt}}^2\left({\tilde{\mathbf{Z}}}^{\mathit{mn}},{\tilde{\mathbf{Z}}}^{\mathit{ij}}\right)=d_p^2\left({\mathbf{Z}}^{\mathit{mn}}+{\sigma }_t^2\mathbf{I},{\mathbf{Z}}^{\mathit{ij}}+{\sigma }_t^2\mathbf{I}\right)$$

(17)

where I represents the 3 by 3 identity matrix and $d_p^2$ may be $d_{\mathit{pw}}^2$ or $d_{\mathit{pg}}^2$.

Consequently, the d_pt distance will have a lower sensitivity to the changes on the retrieved power as it gets close to the system noise power ${\sigma }_t^2$. The value of ${\sigma }_t^2$ may be obtained from the sensor parameters or, alternatively, it may be estimated from the data. In [25] it is proposed to estimate it as the mean power of a dark area on the image, where no backscatter is received from the scene and, then, the measured power corresponds to the system noise. Following this strategy, an automatic ${\sigma }_t^2$ parameter is proposed, dividing the scene into reasonably sized blocks to avoid noise effects, for instance 9 × 9 pixel blocks, and taking as ${\sigma }_t^2$ the minimum power of all the channels and blocks of the entire image. However, it is preferred to employ the sensor information when available, as this method may induce wrong values when no dark areas are present on the image.

Note that an important property of the proposed filtering in Equation (6) is that, as opposite to the multilook filtering, the number of averaged pixels is different for every pixel of the image. This fact may cause some difficulties for further data analysis and processing that assume that all pixels are equally filtered. However, this information is still available as the k parameter Equation (7) that may be obtained and taken into account since it contains the amount of filtering for each pixel. Additionally, due to the adaptive nature of the bilateral filter, it holds also information related with the image structure.

3. Iterative Weight Refinement

As stated in Section 1, usually SAR images are strongly contaminated by speckle noise. Then, a given pixel of the image will have a large component induced by the speckle that will also contaminate all the pixel-based comparisons of the filter defined in Equation (6). As a consequence, the computed w_p weights will be also affected by noise resulting in a noisy image. To reduce this effect, an iterative weight refinement procedure is proposed, employing the filtered images to compute a more reliable estimation of the filter weights. The filter weights at each iteration are computed over the filtered image of the previous iteration, as represented on Figure 1.

In this diagram, “Img 0” represents the original image and “Img i” represents the output result obtained after iteration i. CBL represents a cross-bilateral filter [26] based on the filter defined in Equation (6); the weighted average is performed over the input image in, whereas the weight computation is performed over the reference image ref:

$${\widehat{\mathbf{Z}}}^{\mathit{ij}}=\frac{1}{k(i,j)}\sum _{m,n\in V(i,j)}{\mathbf{Z}}_{\text{in}}^{\mathit{mn}}{w}_s(i,j,m,n)w_p\left({\mathbf{Z}}_{\text{ref}}^{\mathit{mn}},{\mathbf{Z}}_{\text{ref}}^{\mathit{ij}}\right)$$

(18)

where ${\mathbf{Z}}_{\text{in}}^{\mathit{ij}}$ and ${\mathbf{Z}}_{\text{ref}}^{\mathit{ij}}$ represent the estimated covariance matrix of the input and reference images at position (i, j), respectively.

It is worth noting that the proposed scheme shown in Figure 1 defines iterative weight modification, not iterative filtering. Iterative filtering may produce stronger filtering but usually tends to produce overfiltering by mixing similar non-homogeneous areas and blurring some details due to the propagation of errors. In contrast to iterative filtering, the proposed technique processes the original image at each iteration. This scheme avoids the propagation of errors and allows a more reliable estimation of the polarimetric weights w_p based on filtered images. Additionally, iterative filtering has the disadvantage of merging samples with different amounts of filtering through the iterations. This fact will produce a k parameter that does not correspond to the real number of samples averaged from the original image, making the analysis and interpretation of the obtained results harder.

4. Results

The proposed filtering technique has been employed to process the real RADARSAT-2 image presented in Pauli RGB composition in Figure 2a. This image is represented in slant range radar geometry to avoid the geocoding distortion before filtering and, thus, the axes of the image represent range and azimuth directions, respectively, with a pixel spacing of approximately 4.7 by 4.7 m. It corresponds to a Fine Quad-Pol image from a test-site in Flevoland, the Netherlands, that was acquired on 4 April 2009 employing the beam FQ13, during the ESA AgriSAR 2009 campaign. An optical image of this region with the acquisition area marked in red is represented in Figure 3 [27]. The AgriSAR 2009 campaign was devoted to analyze agricultural fields evolution employing polarimetric SAR data. Consequently, the image is mainly composed by agricultural fields, but it also contains sea surface and an urban area at the bottom part of the image. This image has been filtered with the proposed technique, employing the d_pw similarity measure, a 11 × 11 pixel local window V, σ_s = 3 and σ_p = 0.6. These values have been chosen experimentally, as they provide a good level of speckle filtering while also preserving the details of the image. The effect of these parameters in the filtering process will be studied later. Additionally, the automatic method for ${\sigma }_t^2$ calculation has been applied, as stated before. In this case, the original image contains dark areas, as seen on Figure 2a, corresponding to closed and calm water, where this parameter may be estimated correctly. The filtered image is presented on Figure 2b after employing the weight refinement scheme presented on Figure 1 with 5 iterations.

To see clearly the edge and small detail preservation of the proposed technique, Figure 4e shows an area from this result, comparing it with the classical 7 × 7 multilook, represented on Figure 4b, the refined Lee filter [7], on Figure 4c, and the IDAN [8] filter, on Figure 4d. For the computation of the refined Lee and IDAN results the PolSARPro [28] implementation has been employed. This area contains some forest and agricultural areas in the top part and urban and sea areas in the bottom part. Additionally, a zoom is shown from the central part of the image, in order to see these results at pixel level. As it may be seen, the multilook filter blurs the edges on the image, mixing inhomogeneous samples of the image, and enlarges small point scatterers, which may also be seen as a spatial resolution loss. The refined Lee filter may achieve a good adaptation to the contours of the image, but it produces some artifacts around it, distorting the resulting image. The IDAN filter achieves a good level of filtering and preserves some contours, but it blurs small details, as point scatterers, of the image. On the contrary, the proposed filtering technique has a good edge preservation without introducing artifacts, that may be clearly seen over the coastline, and also preserves the small details within the urban area, which is particularly visible at the zoom of the images. Figure 4f shows the k parameter for the last iteration. Although this parameter may not be regarded as the Equivalent Number of Looks (ENL) [29] it is related to it, in a similar manner to the multilook filtering window size. As explained before, the proposed method employs also the polarimetric locality to avoid mixing inhomogeneous samples of the image. This effect is also reflected on Figure 4f, where small values may be observed near contours. On the other hand, larger values of the k parameter are observed over homogeneous areas, in the order of 40, similarly to the 7 × 7 multilook filter, that would have an equivalent k parameter constant over the image and equal to 7 × 7 = 49 averaged pixels.

As mentioned before, the iterative weight refinement scheme, presented on Figure 1, has been employed, performing 5 iterations, to obtain a more reliable estimation of the filter weights. The effect of this refinement scheme over the weights is reflected indirectly in Figure 5, where the k parameter is represented for each iteration of the iterative scheme. In each image a different color scale has been employed in order to increase the image contrast, as it may be seen on the corresponding color bar. Note that the output at the last iteration corresponds to Figure 4f and, therefore, it is not represented on Figure 5. It may be seen that, due to the large amount of speckle noise present on the image, small values of k are observed at the first iteration, indicating that the weights w_p present low values. Additionally, the noise of the original image is also affecting the computed weights which are also noisy as may be seen on the zoom images. In successive iterations, larger values of k are obtained in homogeneous areas, due to the increased level of speckle filtering, resulting in a more clear detection of the contours and homogeneous areas. This evolution may also be seen on Figure 6, where the histograms of the k parameter over the image are represented for all the iterations. This iterative process may be continued, as a convergence trend is observed. However, small differences were reported in the following iterations and, for computational simplicity, we have decided to stop it after the 5th iteration.

In the previous examples only the Pauli RGB representation and the k parameter have been compared. To analyze the preservation of the polarimetric information, the Entropy (H) and averaged alpha angle (ᾱ) polarimetric decomposition parameters [16] are represented on Figure 7. These parameters are obtained from the eigen-decomposition of the estimated covariance matrix Z, depending, thus, on the complete covariance matrix. A comparison is presented between the proposed method and the 7 × 7 multilook, refined Lee and IDAN filters. As it may be seen on Figure 7, qualitatively similar values are obtained in all cases for large areas, meaning that the estimated covariance matrix by these methods have similar polarimetric information to the one obtained with the multilook over homogeneous areas. However, the improved spatial resolution preservation is clearly stated for small details of the image, specially over the urban area, as seen on the zoomed area, where the multilook filter enlarges those details according to its window. Refined Lee, shown on Figure 7b,f, preserves the small details of the image whereas the IDAN filter, Figure 7c,g, almost destroys some small details that appear as small blue and red dots in H and ᾱ, respectively, when employing the refined Lee or the proposed method. The proposed filter is able to preserve the shape and contours of the small details, similarly to the refined Lee, while also may attain a higher amount of filtering over homogeneous areas while not introducing distortion to the image.

Additionally, to make a quantitative evaluation,

Table 1. Mean estimated values over homogeneous areas for different filtering techniques. For the proposed technique, an 11 × 11 local window has been employed with 5 weight refinement iterations and σ_s = 3, σ_p = 0.6.

Region	Filtering	C11	C22	C33	|ρ13|	arg(ρ13)(°)	H	A	ᾱ(°)
	Original	2.302 × 10−1	1.125 × 10−1	1.906 × 10−1	-	-	-	-	-
Z1	7 × 7 Multilook	2.305 × 10−1	1.130 × 10−1	1.933 × 10−1	0.4194	1.670	0.8504	0.2599	42.87
Forest	Refined Lee	1.879 × 10−1	9.965 × 10−2	1.623 × 10−1	0.3971	2.873	0.8442	0.3063	44.51
5,000	IDAN	1.471 × 10−1	7.353 × 10−2	1.252 × 10−1	0.2993	2.339	0.8878	0.2513	46.03
pixels	Proposed dpw	2.227 × 10−1	1.090 × 10−1	1.869 × 10−1	0.3865	1.569	0.8698	0.2244	43.46
	Proposed dpg	2.210 × 10−1	1.082 × 10−1	1.853 × 10−1	0.3869	1.550	0.8661	0.2290	43.52
	Original	1.787 × 10−2	1.280 × 10−3	2.552 × 10−2	-	-	-	-	-
Z2	7 × 7 Multilook	1.801 × 10−2	1.271 × 10−3	2.552 × 10−2	0.7443	3.088	0.4351	0.6407	19.30
Water	Refined Lee	1.492 × 10−2	1.255 × 10−3	2.099 × 10−2	0.7062	2.535	0.4701	0.6438	21.89
3,900	IDAN	1.158 × 10−2	8.901 × 10−4	1.667 × 10−2	0.6263	2.473	0.5248	0.6967	24.32
pixels	Proposed dpw	1.749 × 10−2	1.263 × 10−3	2.479 × 10−2	0.7094	2.872	0.4689	0.6381	20.59
	Proposed dpg	1.732 × 10−2	1.254 × 10−3	2.458 × 10−2	0.7079	2.871	0.4693	0.6392	20.75
	Original	1.862 × 10−1	3.195 × 10−2	1.771 × 10−1	-	-	-	-	-
Z3	7 × 7 Multilook	1.855 × 10−1	3.170 × 10−2	1.758 × 10−1	0.7489	8.423	0.5444	0.2998	22.87
Crop	Refined Lee	1.434 × 10−1	2.862 × 10−2	1.393 × 10−1	0.7059	7.478	0.5781	0.3653	25.86
2,200	IDAN	1.020 × 10−1	1.948 × 10−2	1.006 × 10−1	0.6187	8.641	0.6514	0.3783	27.84
pixels	Proposed dpw	1.765 × 10−1	3.087 × 10−2	1.679 × 10−1	0.7062	8.474	0.5880	0.2933	24.50
	Proposed dpg	1.757 × 10−1	3.058 × 10−2	1.671 × 10−1	0.7072	8.470	0.5808	0.2987	24.53

shows a comparison of some elements of the covariance matrix estimated Z for the 7 × 7 multilook filter and the previously analyzed methods. Three homogeneous areas have been manually selected, that are marked on Figure 4a, corresponding to a forested area, a sea area and an agricultural crop. Table 1 shows the mean estimated values for the diagonal elements of the Z matrix averaged over those areas for original and filtered images employing the 7 × 7 multilook, refined Lee and IDAN filters and the proposed method employing both dissimilarity measures defined in Equations (15) and (16). Additionally, for the filtered images the mean value of the magnitude and phase of the correlation coefficient ρ₁₃ over HH and VV polarization states and Entropy (H), Anisotropy (A), and averaged alpha angle (ᾱ) are represented. Note that these values may not be estimated over the original image since some filtering has to be applied to the data in order to have full-rank matrices [16]. As it may be seen from Table 1, all the adaptive filters analyzed seem to underestimate the covariance matrix diagonal elements to some extent. However, it is worth to notice that the proposed method has the lower figures of bias in all cases. Nonetheless, all the methods obtain similar values to the multilook in terms of the H/A/ᾱ decomposition parameters.

In order to compare the amount of speckle filtering among different techniques, the equivalent (or effective) number of looks (ENL) is a well-known parameter. It describes the amount of averaging performed to SAR data under the complex Gaussian assumption [30]:

$$\mathit{ENL}=\frac{E{\left\{I\right\}}^2}{\mathit{Var}\left\{I\right\}}$$

(19)

where E{·} denotes the expectation operator and Var{·} denotes the variance of the intensity image I.

Usually, the ENL is estimated over an homogeneous area of the image, assuming that the speckle is fully developed and that no texture is present, that is, assuming a constant radar cross section [30]. In this scenario, the estimated ENL L̂ may be expressed as [29]:

$$\widehat{L}=\frac{{\langle I\rangle }^2}{\langle I^2\rangle -{\langle I\rangle }^2}$$

(20)

where 〈·〉 denotes the averaging over the homogeneous area selected.

However, this measure was defined for single polarization SAR data and, when employing PolSAR data, different ENL values are obtained for each polarimetric channel. Some extensions to PolSAR data were proposed and analyzed in [29] based on the matrix trace moments (TM) L̂_TM and on the Maximum Likelihood (ML) based on the Wishart distribution Equation (4) L̂_ML:

$${\widehat{L}}_{\mathit{TM}}=\frac{\text{tr}{\left(\langle \mathbf{Z}\rangle \right)}^2}{\langle \text{tr}(\mathbf{ZZ})\rangle -\text{tr}(\langle \mathbf{Z}\rangle \langle \mathbf{Z}\rangle )}$$

(21)

$$\langle \text{ln}\left|\mathbf{Z}\right|\rangle -\text{ln}|\langle \mathbf{Z}\rangle |-\sum _{i=0}^2{\Psi }^{(0)}\left({\widehat{L}}_{\mathit{ML}}-i\right)+3\text{ln}\left({\widehat{L}}_{\mathit{ML}}\right)=0$$

(22)

where tr(·) denotes the matrix trace, | · | represents the matrix determinant and Ψ⁽⁰⁾(·) refers to the digamma function [29]. No analytical solution has been derived for L̂_ML in Equation (22) and, then, this equation has to be solved numerically.

To asses quantitatively the effectiveness of the different filtering techniques in terms of speckle noise reduction,

Table 2. Estimated equivalent number of looks (ENL) for different filtering techniques. For the proposed technique, an 11 × 11 local window has been employed with 5 weight refinement iterations and σ_s = 3.

Region	Filtering	EN L(C11)	EN L(C22)	EN L(C33)	EN LTM(Z)	EN LML(Z)
	Original	0.9070	0.8309	0.9111	0.9241	-
	7 × 7 Multilook	15.61	11.40	16.45	16.54	17.24
Z1	Refined Lee	12.04	7.982	12.23	12.25	12.45
Forest	IDAN	8.417	6.565	8.975	15.92	-
5,000	Proposed dpw, σp = 0.6	8.610	7.764	10.20	13.24	18.52
pixels	Proposed dpw, σp = 0.9	20.20	15.04	23.24	24.31	26.57
	Proposed dpg, σp = 0.6	7.860	6.993	9.090	11.97	17.19
	Proposed dpg, σp = 0.9	15.92	12.68	18.50	20.73	24.12
	Original	1.018	0.9865	1.075	1.053	-
	7 × 7 Multilook	23.25	15.67	21.63	22.73	20.74
Z2	Refined Lee	15.54	11.06	15.09	15.23	13.66
Water	IDAN	8.407	8.180	9.027	10.30	-
3,900	Proposed dpw, σp = 0.6	7.297	17.32	7.516	8.085	21.48
pixels	Proposed dpw, σp = 0.9	28.96	22.89	26.44	29.09	33.81
	Proposed dpg, σp = 0.6	6.631	14.29	6.845	7.375	19.73
	Proposed dpg, σp = 0.9	18.52	20.55	17.73	19.40	29.79
	Original	0.7275	0.7319	0.8457	0.8090	-
	7 × 7 Multilook	8.372	7.522	15.27	11.97	15.00
Z3	Refined Lee	6.798	6.364	9.467	8.526	10.62
Crop	IDAN	5.348	5.707	7.364	8.285	-
2,200	Proposed dpw, σp = 0.6	2.947	4.889	4.464	4.254	12.97
pixels	Proposed dpw, σp = 0.9	6.810	8.061	12.93	10.59	19.77
	Proposed dpg, σp = 0.6	2.911	4.561	4.308	4.134	12.55
	Proposed dpg, σp = 0.9	5.540	7.164	9.755	8.425	17.84

shows the ENL estimated over the homogeneous areas shown on Figure 4a employing the estimators presented before. The columns 3 to 5 show the ENL L̂ estimated according to Equation (20) for each of the polarimetric channels, representing the C_ii components of the covariance matrix. The 6th column represent the ENL L̂_TM estimated employing the full covariance matrix according to Equation (21), whereas the 7th column represent the L̂_ML according to Equation (22). In the case of the original image and the IDAN filtered image it was not possible to compute the L̂_ML estimator due to the presence of singular matrices.

As it may be seen, there are significant differences according to the distinct ENL estimators. The performance of these estimators has been studied in detail in [29], where the advantages of taking into account the full covariance matrix in L̂_TM and L̂_ML are discussed. Moreover, the L̂_ML estimator is labeled as the best one in terms of variance and robustness to texture. Analyzing the values obtained according to this estimator, it may be seen that similar speckle reduction to the 7 × 7 multilook filter may be obtained with the proposed technique employing σ_p = 0.6 and even larger filtering is obtained when increasing this parameter. Results obtained for different similarity measures are very close in terms of ENL, obtaining a slightly higher level of filtering for d_pw than for d_pg. It is worth to mention that the ENL measure only refers to the amount of speckle filtering over homogeneous areas. The main advantage of the proposed technique lies in its ability to achieve a filtering level comparable to the 7 × 7 multilook while also retaining the original resolution and preserving the small details of the image.

As mentioned in Section 2, the parameters σ_s and σ_p control the weights sensitivity in the spatial and polarimetric domains, respectively. To analyze the effect of these parameters on the filtering process, the same real image, presented in Figure 2a, has been processed employing different values of sigma σ_s and σ_p. Some results are represented on Figure 8 for the same detail area of Figure 4, showing on Figure 8a,b, and Figure 8d,e changes on the σ_s and σ_p parameters, respectively. As it may be seen when comparing these results with the one presented on Figure 4e, reducing the value of σ_s or σ_p increases the weights sensitivity, resulting in a decrease of the speckle filtering level. Similarly, increasing the values for σ_s or σ_p has the opposite effect. Additionally, when comparing Figure 8b,c,e,f the effect of changing the dissimilarity measure may be observed. In this case, the dissimilarity measure has a smaller impact, since the obtained results are very similar in both cases.

However, these parameters have a different impact on the filtering process, as they are referring to different domains. In order to analyze this impact, Figure 9 represents the 50-bin histograms of the k parameter, for different values of σ_s and σ_p. As it may be seen, when changing the value of σ_s, on Figure 9a–c, the spatial decay is changed but the sensitivity on the polarimetric domain remains unaltered. In terms of the k parameter, this implies a larger amount of filtering per pixel, increasing the values of k as σ_s increases. Similarly, the σ_p parameter also affects the k parameter in the same way, as represented on Figure 9d–i. However, note that changing σ_p changes significantly the shape of the k histogram which is not the case for the σ_s parameter. Moreover, the σ_p parameter has an important role to play in terms of the contour preservation. In fact, the selection of this parameter is a compromise between the contour preservation and the amount of speckle filtering, as increasing σ_p may result in a heterogeneous sample mixture over areas with subtle contrasts, as it may be observed on Figure 8e,f. It is worth noting that the amount of filtering can not be increased indefinitely, as it is firstly limited by σ_s, which defines a spatial decay, and ultimately it is limited by the local window V, defining an upper bound for k. Furthermore, note that the proposed filter tends to the mulitlook filter over the V window when σ_s → ∞ and σ_p → ∞. In this sense, it may be seen as a generalization of the multilook filter.

In the previous examples, the proposed method has been analyzed from the point of view of the speckle filtering application. However, it may also be useful for other applications as, for instance, matrix regularization. As mentioned before, the sample covariance matrix for original pixels Z = kk^† has rank equal to one, whereas for distributed scatterers the covariance matrix defining their statistical distribution has full-rank. Moreover, the Wishart distribution, that describes the statistics of the sample covariance matrix Z under the Gaussian hypothesis, is not defined for singular matrices. This may be a problem for some techniques that assume full-rank matrices as, for instance, PolSAR data filtering and segmentation [15,21] or classification [14]. In order to circumvent this issue, these techniques apply an initial filtering stage to the data, usually a small multilook. The ultimate objective of this stage is not to reduce the effect of the speckle, but to avoid having singular matrices. Furthermore, the mentioned applications are based on similarity measures or distances over the covariance matrix space, as in the case of the proposed technique, making it a natural preprocessing stage as it is based on parallel concepts. For Binary Partition Tree (BPT) [31] based PolSAR speckle filtering and segmentation [15,32] an initial 3 × 3 multilook filter is applied for this purpose. However, this initial step implies a small spatial resolution loss, as seen previously. The proposed method may be applied instead of the multilook filter in order to preserve, as much as possible, the original spatial resolution of the data. Note that, in this case, the amount of filtering needed is much lower than in previous examples, just to regularize the obtained sample covariance matrices.

In order to see the benefits of employing the proposed method as a covariance matrix regularization instead of the 3 × 3 multilook, Figure 10 shows a crop area of the results of applying the same BPT-based filtering to the image presented in Figure 2a with these initial regularization steps. Figure 10a shows an area of the image obtained after applying a 3 × 3 multilook to the original data, whereas Figure 10b shows the same result after applying the proposed technique employing 3 weight refinement iterations, σ_s = 2 and σ_p = 0.4 with a 5 × 5 local window V. As in the previous examples, the automatic method for ${\sigma }_t^2$ calculation has been applied. A BPT homogeneity based pruning has been applied to segment the data into homogeneous regions, as described in [15], employing the geodesic d_sg similarity function defined and analyzed in [23]. Results are shown on Figure 10c,d after processing the previously commented full images employing a BPT pruning factor of δ_p = −2 dB for both cases. As it may be seen, better resolution preservation is achieved when employing the proposed method as a matrix regularization process, on Figure 10d, as some point scatters of the image are maintained with its original size whereas if the 3 × 3 multilook is applied, as shown on Figure 10c, they are enlarged according to the filtering window. This is particularly evident over urban areas, having high density of point scatters and small structures.

5. Conclusions

In this paper, a novel Polarimetric Synthetic Aperture Radar (PolSAR) speckle filtering technique is proposed by exploiting the polarimetric locality, as well as its spatial locality in terms of polarimetric and geometric distances, respectively. Then, the spatial averaging is replaced by a weighted averaging favoring similar pixels in both spatial and polarimetric domains. This approach is based on a modified bilateral filter where a new iterative scheme is introduced to mitigate the noise over the estimated weights, leading to a large filtering over homogeneous areas, comparable to the classical multilook filter, while also preserving contours and small spatial details of the scene. In addition, the use of a polarimetric distance depending only on the diagonal elements of the covariance matrix, allows to consider the proposed technique as a regularization pre-processing step for other applications that require full-rank matrices, producing better results than a small multilook filtering for this purpose.

The technique has been applied to process real RADARSAT-2 data, leading to a better preservation of the polarimetric and the spatial information than other polarimetric approaches. In terms of polarimetric information, the retrieved values of Entropy, Anisotropy and Mean Alpha Angle present, in general, a lower bias. Additionally, the proposed approach allows a better preservation of spatial details associated to point targets presenting low Entropy.

A possible drawback of the proposed technique is that, as opposed to the multilook filter, it is employing a different number of averaged samples for each pixel of the image. However, this information is related to the k parameter that is obtained by the technique. Furthermore, this parameter may be also useful for information extraction as, for instance, contour detection, since it is strongly related with the scene structure.