Outlier Rejection by Means of Median Matrices for Polarimetric SAR Covariance Symmetry Classification

This letter exploits the intrinsic selectivity properties of the median to enhance the covariance symmetry classification in polarimetric synthetic aperture radar (PolSAR) images. More in detail, the median matrices are utilized to properly detect and remove outliers in the data belonging to a reference window, in turn used to estimate the covariance structure of the pixel under test. Hence, the scene is classified in terms of the structures assumed by the covariance under specific symmetric scattering mechanisms. To do this, for each pixel under test, the data in a reference window are filtered through the application of a generalized inner product (GIP)-based procedure involving the median matrix in its computation. The filtered data are then used as input to a model order selection (MOS)-based procedure for the final scene classification. Tests conducted on L-band real-recorded SAR data show the effectiveness of the devised framework.


I. INTRODUCTION
R ECENT years have seen an ever wide growth in the implementation of polarimetric synthetic aperture radar (PolSAR) platforms, followed by a strong increment in the availability of fully polarimetric data. These achievements have led to the developments of many research topics based on the exploitation of PolSAR, especially in terms of scene segmentation and/or land-cover classification, e.g., [1], [2], [3], [4], [5], [6], [7], [8], [9], [10], [11], [12] and references therein. In this context, the main utilized tools are the polarimetric covariance, coherence, and Muller matrix [3], [4], [5], [13] that allow to fully exploit the polarimetric information. These matrices are the basis in many processes aimed at describing the scattering mechanisms characterizing the observed scene. Moreover, as shown in [2], the covariance assumes particular forms when some symmetric structures in the target characterize the polarimetric returns. This structure, if properly revealed, can be exploited to obtain more accurate covariance estimates to be used in place of the sample one to improve the final scene classification.
To properly detect the specific structure assumed by the covariance (as well as the coherence), a methodology exploit- Manuscript  ing the model order selection (MOS) has been designed in [6]. In fact, the MOS rules [14], [15] were introduced since the resulting problem was designed as a multiple hypothesis test with nested instances. Hence, the method in [6], for each pixel under test, assuming homogeneity, extracts a set of pixels through a sliding window to be used in the next inferences. However, homogeneity assumption often fails due to the presence of outliers in the set of data. Sources for outliers are point-like targets in the reference data, strong speckle, power variations, or discontinuities in the scene. This drawback has been dealt in [16], where pixels in the reference window are filtered before the computation of the MOS detectors exploiting a generalized inner product (GIP)-based scheme involving geometric barycenters. However, thanks to the intrinsic selectivity properties of the median [17], [18], in this letter, following the lead of [16], a more robust algorithm for symmetry classification in the presence of outliers is proposed. As a matter of fact, it exploits the median matrices of elementary covariances related to the pixels in the reference window to perform the outlier cancellation in the extracted data. The performances of the proposed method are assessed on the fully polarimetric L-band ElectroMagnetic Institute SAR (EMISAR) data in comparison with its homogeneous counterpart not performing any data screening.

II. PROBLEM FORMULATION AND ALGORITHM DESCRIPTION
The proposed method exploits the intrinsic robustness of the median to outliers in the data spatially close to the pixel under test to improve the covariance symmetry classification capabilities. The proposed scheme depicted in Fig. 1 starts with the acquisition of the fully PolSAR image. Then, for each pixel under test, a neighbor window is extracted to perform inference. Hence, data in the selected window are filtered using a GIP-based outlier cancellation procedure that makes use of the estimated median matrices [18]. The filtered data are finally used to classify the covariance symmetry as described in [6] for the homogeneous environment.

A. PolSAR Datacube
The pipeline proposed in this letter uses three polarimetric channels, say S HH , S HV , 1 and S VV , with channel S VH involved in the noise estimation process. Hence, for each pixel, N = 3 complex polarimetric returns are stored in x l,m , l = 1, . . . , L and m = 1, . . . , M to obtain an L × M × N datacube X (with L and M the azimuth/range sizes). Then, for each pixel under test, a window B of size K = W 1 × W 2 ≥ N is extracted and its pixels are organized in the random matrix, R, whose columns, r 1 , . . . , r K , are modeled as independent and identically distributed (i.i.d.) zero-mean circular complex Gaussian random vectors with covariance M.

B. Median Matrices
In this section, the covariance is estimated as the median matrix related to a distance defined in the space of positive definite matrices. The median is derived from a set of elementary covariance estimates M k , k = 1, . . . , K , that is is a specific distance defined in the space of positive definite matrices, with A ≻ 0 (resp. ⪰) denoting that A is positive definite (resp. semidefinite). The elementary covariance estimates M k , k = 1, . . . , K , are computed from the data r k in the reference window [18], [21], [22], that is , U k the unitary matrix of the eigenvectors of r k r H k with the first eigenvector corresponding to the eigenvalue ∥r k ∥ 2 , and σ 2 0 denoting the thermal noise power. Moreover, diag (·) is the diagonal matrix whose diagonal is the entry vector, the symbol (·) H denotes the conjugate transpose operator, and ∥·∥ is the Frobenius norm. The procedure for noise power estimation has been derived in [16] and exploits the residual mismatches between the cross-polarization channels, S HV and S VH , due to thermal noise variations only, that iŝ with |·| representing the modulus of its complex argument.
As proved in [18], the optimal solution to the optimization problem (1) can be found solving the following equivalent convex optimization semidefinite programming (SDP) problem with vec(·) the column vector obtained lining-up the columns of its matrix argument, and C the set of Complex vectors of size N 2 .
In this letter, we refer to the Log-Euclidean distance between two positive definite Hermitian matrices A and B, that have shown the best selective capabilities, 2 i.e., with tr {·} the trace of its matrix argument. 3 Therefore, with distance (5), substituting M k with log M k and M with log M, and denoting by M the corresponding optimal solution to (4), the Log-Euclidean median estimator is Before concluding this Section, it is worth to recall that the computational complexity for solving the SDP is O(N 3.5 log(1/η)), with η a desired accuracy [23].

C. Filtering Data in the Reference Window for Outlier Removal
The filtering operation for outlier removal consists in evaluating the GIPs for each pixel in the reference window, so as to remove those sharing the highest κ 0 values [24], [25]. Following the lead of [18], the GIP is computed with the above-described median matrices in place of the classic sample covariance matrix (SCM). More in detail, the method consists in evaluating the following quadratic form for each single data, r 1 , . . . , r K , contained in the window under analysis, that is After the computation of ρ (d) k , k = 1, . . . , K , the κ 0 data sharing the highest GIP values are removed. Obviously, κ 0 is a tuning parameter that allows to rule the trade-off between the amount of data to excise to ensure homogeneity and that for the next covariance estimations. Setting a high κ 0 value leads to a strong screening of the data, but this in turn reduces the number of available data for the sample estimate (recall that it should be K − κ 0 ≥ N to guarantee that the sample matrix is a good covariance estimate [26]). In this letter, κ 0 is set with the adaptive selection procedure devised in [16]. It consists in choosing it as the value to which the corresponding ordered GIP values contain a preassigned percentage of the whole energy (that is set to the 20% in the tests of Section III to be conservative with respect to the amount of data cancellation).

D. Classification of Dominant Covariance Symmetry Structures
The classification of the dominant covariance symmetry structure (derived from a symmetry property in the polarimetric returns [2]) for each pixel of the considered synthetic aperture radar (SAR) image is performed through the adaptive approach designed in [6]. In particular, due to the fact that the methodology in [6] assumes that data in the reference window are homogeneous, herein they are substituted by the data filtered through the median-matrix-based approach as described above. As to the classification process, in what follows four different covariance symmetries [2] are considered, and hence the resulting problem is formulated as a multiple hypothesis test comprising both nested and not nested hypotheses, that is H 1 : no symmetry H 2 : reflection symmetry H 3 : rotation symmetry H 4 : azimuth symmetry. (8) It worth recalling that the above-mentioned covariance structures can be found in several practical situations. As an example, reflection symmetry with respect to a vertical plane can be observed on forest, snow, etc. [2]. Rotation symmetry can be seen in medium like the Earth's ionosphere [2]. Conversely, azimuth symmetry can be experienced in the presence of vertical tree trunks when signals are scattered from a randomly oriented distribution of dipoles, e.g., randomly oriented branches [2], [27].
Due to the presence of nested hypotheses 4 in (8), in [6], a modified version of generalized maximum likelihood (GML) based on the use of the MOS rules [14], [15] is introduced to efficiently solve the detection problem. Hence, the considered with M (n(h)) the ML estimate of M comprising n(h) parameters, viz., n(1) = 9, n(2) = 5, n(3) = 3, and n(4) = 2 [6]. Moreover, fR(R|M) is the complex multivariate probability density function (pdf) of the filtered observable matrixR, whose columns are modeled as zero-mean complex circularly symmetric Gaussian vectors [6]. More in detail, (9) is evaluated under each hypothesis, h = 1, . . . , 4, and the order corresponding to the minimum between the four statistics is selected, that iŝ In other words, for each pixel under test, the selected structure is the one associated with Hĥ. Note that, the term n(h) η(n(h), K − κ 0 ) is the penalty coefficient that is aimed at penalizing overfitting [14], that depends on the specific MOS. In the tests conducted in this letter, we analyze and compare the Bayesian information criterion (BIC) [14], [15] and Hannan-Quinn information criterion (HQC) [28], for which the term η(n(h), K − κ 0 ) becomes log(K − κ 0 ) and 2 log(log(K − κ 0 )), respectively.

III. RESULTS AND DISCUSSION
In this section, the results of the application of the proposed framework on measured data are analyzed and discussed. To do this, the L-band (1.25 GHz) EMISAR coherent polarimetric sample data 5 are used. The image refers to a scene  of Foulum Area (DK), containing buildings, vegetation, and water. The conducted analyses concentrate on a portion of the entire SAR image of size 401 × 401 pixels, whose span (i.e., the total power of the polarimetric image |S HH | 2 + |S HV | 2 + |S VH | 2 + |S VV | 2 ) is shown in Fig. 2 as ground truth. Results of the tests are shown in terms of classified symmetries in Fig. 3. More in detail, the proposed detectors based on Bayesian information criterion (BIC) and HQC exploiting the Log-Euclidean median in the outlier cancellation process (shortly denoted as median BIC and median HQC) are compared with their homogeneous counterparts [6] (indicated as BIC, and HQC) as well as the classifiers based on the exploitation of the competing barycenters [16]. Results are obtained extracting for each pixel K = 25 data with a 5 × 5 window. Moreover, for each pixel a specific color is assigned to each symmetry as detailed in the legend of the figure.
Observing the colored maps, as expected, the vegetated areas (essentially forests) show a dominant azimuth symmetry due to the reflection of the tree trunks. Moreover, reflection symmetry is observed over crops and bare fields, with some rotations shown by man-made objects in urban areas. Finally, the absence of symmetry is mainly present over the lake in the lower left corner of the image, and on a cultivated field in the bottom center, that probably contains a cultivation with different structural characteristics than the others. Comparing the subplots, it is also interesting to see that the proposed median-based estimator tends to provide a sharper separation of the different areas in the classified image that also share more homogeneity with respect to the standard counterpart. Moreover, the median and barycenter approaches share similar classification capabilities as shown in Table I, even tough, the median has a general strong selective capabilities. In fact, pixels classified as having azimuth symmetry (that is the most challenging case) are the 26.50% with the median and 26.02% with the barycenter. Finally, for this specific scenario, the BIC detectors works better than the HQC counterparts.
To demonstrate further the benefits of the proposed approach, Fig. 4 shows a quantitative comparison between the H/α classification [3], [5] and the one exploiting our algorithm as pre-processing stage.
More in detail, the coherence matrix used to perform the H/α decomposition is that provided by the proposed pipeline exploiting the median BIC selector in place of the classic sample one. Comparing the two figures, it is evident that the proposed pre-processing allows to better emphasize areas representing the same scene with the boundaries that appears sharpened. As an example, with the proposed pipeline forests are mostly classified as class 2, representing random anisotropic scatterer, whereas with the standard approach forest pixels share a mixed classification among classes 2 and 5 (anisotropic particles). Additional, trees are classified with class 4 (double reflection propagation effects), which is more evident in Fig. 4(b).

IV. CONCLUSION
A polarimetric covariance symmetries classification pipeline has been devised to properly account for the presence of outliers in the data used during the adaptation process. The method exploits the selectivity properties of median matrices to reject outliers in the reference data. After that, the remaining data are exploited in the MOS-based detector to classify the covariance symmetry for each pixel under test. Results have shown the effectiveness of the proposed framework on measured L-band data, also in comparison with its counterpart assuming homogeneity. Possible future studies could consider the optimal selection of the reference window size, the estimation of the pdf of κ 0 , and its adaptive selection.