The performance analysis based on SAR sample covariance matrix.

Multi-channel systems appear in several fields of application in science. In the Synthetic Aperture Radar (SAR) context, multi-channel systems may refer to different domains, as multi-polarization, multi-interferometric or multi-temporal data, or even a combination of them. Due to the inherent speckle phenomenon present in SAR images, the statistical description of the data is almost mandatory for its utilization. The complex images acquired over natural media present in general zero-mean circular Gaussian characteristics. In this case, second order statistics as the multi-channel covariance matrix fully describe the data. For practical situations however, the covariance matrix has to be estimated using a limited number of samples, and this sample covariance matrix follow the complex Wishart distribution. In this context, the eigendecomposition of the multi-channel covariance matrix has been shown in different areas of high relevance regarding the physical properties of the imaged scene. Specifically, the maximum eigenvalue of the covariance matrix has been frequently used in different applications as target or change detection, estimation of the dominant scattering mechanism in polarimetric data, moving target indication, etc. In this paper, the statistical behavior of the maximum eigenvalue derived from the eigendecomposition of the sample multi-channel covariance matrix in terms of multi-channel SAR images is simplified for SAR community. Validation is performed against simulated data and examples of estimation and detection problems using the analytical expressions are as well given.

Introduction

Multi-channel systems with random nature appear in a wide range of fields in the literature. For several of them, the central limit theorem applies and their random behavior can be modeled by Gaussian statistics, being thus statistically fully described by the first and second order moments. In the case of multi-channel SAR systems, the assumption of zero-mean multivariate complex Gaussian distribution is frequently valid for geophysical media, being thus fully described by the complex Hermitian covariance matrix. This is the statistical case treated in this work.

Since a SAR signal corresponds to the superposition of the scattered fields from all scatters inside a resolution cell, physical parameter estimation may be performed by exploiting the multi-channel covariance matrix. For that, eigendecomposition theorems have been widely used in the literature, decomposing the whole covariance matrix into elementary quantities, in order to provide a better physical interpretation of the data.

The eigendecomposition of the multi-channel SAR covariance matrix leads to useful information in a large range of SAR applications, including Ground Moving Target Indication (GMTI) [1], polarimetric SAR (PolSAR) [2], interferometric SAR (InSAR) [3], polarimetric interferometric SAR (PolInSAR) [4], change detection [5], target detection [6] and filtering [7,8]. Particularly, the covariance matrix maximum eigenvalue has been proved to be a key parameter in various areas. In [5] and [9], the change detection problem was elaborated using the maximum eigenvalue of the interferometric covariance matrix. The eigendecomposition of the SAR covariance matrix with application to SAR GMTI is used in and to calculate the probability of moving target detection in homogeneous and heterogeneous terrain. Regarding polarimetry, the eigendecomposition of the coherency (which is also a covariance) matrix is the most common tool among incoherent target decomposition theorems for the interpretation of the earth’s surface. The maximum eigenvalue of the coherency matrix is also often used in polarimetry, due to the close relationship to the dominant scattering mechanism of the scattering process [2,10,11].

Although there is a wide area of applications regarding the maximum eigenvalue of the multi-channel covariance matrix, a certain lack of information in the literature concerning its statistical characterization in terms of multi-channel SAR may be verified. The statistical description can be useful to understand bias effects and to analyze estimation as well as detection problems. The estimated (sample) covariance matrix follows a Wishart distribution, which has been of interest along the years since its first derivation. In areas as communication systems, several works using such results have been published [12-16]. The application of the Wishart distribution in remote sensing was however considerably late introduced [17].

The cornerstone study considering the statistical description of the covariance matrix eigendecomposition in polarimetry has been carried out in [18]. The majority of the analysis in [18] was performed on the basis of numerical methods. In this paper the results of [18] is supported by addressing analytical solutions. Additionally, it is derived an exact closed form expressions for the Moment Generating Function (MGF) of the sample covariance matrix maximum eigenvalue.

The validation of the derived expressions is performed using simulated data, demonstrating its agreement with the theory. The effect of the number of samples and underlying correlation scenario of the sample covariance matrix on the bias in the estimation of the maximum eigenvalue is also investigated. Finally, examples of applications are as well given, in the fields of estimation and detection theory.

The next section reminds the basics of the statistical description of multi-channel SAR systems and of the covariance matrix eigendecomposition. It also includes the derived theorems presenting the statistical founds. Section 3 includes numerical examples validating the theorems as well as an analysis of their behavior. Sections 4 and 5 show their usage for the analysis of the estimation bias and for the elaboration of detection problems, respectively. Section 6 concludes the paper with discussions and directions for future work. The detailed statistical derivations can be found in the Appendix.

Statistical Characteristics of the Maximum Eigenvalue of the Multi-Channel Sample Covariance Matrix

Preliminaries

The individual elements of a m-dimensional multi-channel system may be organized in a m-dimensional vector k. As the elements are assumed to follow zero mean complex Gaussian distributions, the vector k is said to follow a m-dimensional multivariate normal distribution, with zero mean and true covariance matrix Σ among the vector elements, and is represented by [19,20]. For zero mean Gaussian statistics the covariance matrix fully describes the data, playing a key role in several application fields. In practical situations however, the true covariance matrix Σ is unknown and has to be estimated by its maximum likelihood estimator (MLE), the sample covariance matrix , where n is the number of estimation samples and † is the transpose conjugate operator. In the SAR context, the number of independent samples is also called looks. The elements of Z follow a m dimensional complex Wishart distribution with n degrees of freedom and true covariance matrix Σ, represented by and defined as [20]: where Γ(·) is the gamma function and etr(·) is the exponential trace of a matrix.

The spectral theorem from linear algebra allows the decomposition of the full rank m-dimensional covariance matrix in a set of m one-rank covariance matrices using its eigenvalues and eigenvectors. Accordingly, the decompositions of the true covariance matrix and its estimator are given by with their real non-negative eigenvalues l and λ, and respective complex eigenvectors e, e′, for i = {1, ...m}.

Sample Maximum Eigenvalue Statistical Description

The following Theorem III concerns the derived Moment Generating Function (MGF) of the maximum eigenvalue of the sample covariance matrix Z of a multi-channel statistical system, which is a critical step in removing the bias of the largest eigenvalue. The Theorem III, which is derived due to the previous works (Theorems I and II), is the main result of the paper. It is detailed in the appendix and will be used for the elaboration of illustrative estimation and detection problems in Sections 4 and 5.

Throughout the paper, | · | is the matrix determinant, denotes the estimator of the random matrix X formed from a sample of size n.

Theorem I: Let be a m-dimensional complex vector whose elements follow zero mean Gaussian distributions with associated m × m covariance matrix Σ. Let Σ have l ≤ .... ≤ l1 eigenvalues. Then the Cumulative Density Function (CDF) of the maximum eigenvalue λ of the sample covariance matrix 〈kk†〉, with the assumption m ≤ n, is given by where Ψ(x) is a m × m matrix with its (i, j)th element , and γ is the incomplete Gamma function (Equation (2.42) in [21]).

Proof: [14] gave the closed expression of CDF of the largest eigenvalue of the complex Wishart matrices. Here, it is written in the form of SAR covariance matrix after small continuation, see also Appendix A.

Theorem II: Let be a m-dimensional complex vector whose elements follow zero mean Gaussian distributions with associated m × m covariance matrix Σ. Let Σ have l ≤ .... ≤ l1 eigenvalues. Then the Probability Density Function (PDF) of the maximum eigenvalue λ of the sample covariance matrix 〈kk†〈, with the assumption of m ≤ n, is given by where Ω(x) is an m × m matrix with its i, jth elements , and Ψ(x) and 𝒮 are defined in Equation (3).

Proof: Equation (4) is obtained by differentiating Equation (3) with respect to x using (Equation (9) in [22]) Here it can be noted that when the true covariance matrix Σ is diagonal, the PDF of the largest eigenvalue reduces to Ermolaev and Rodyushkin result [23].

Theorem III: Let be a m-dimensional complex vector whose elements follow zero mean Gaussian distributions with associated m × m covariance matrix Σ. Let Σ have l ≤ .... ≤ l1 eigenvalues. Then for any positive integers s, the sth moment of the maximum eigenvalue λ of the sample covariance matrix 〈kk†〉, with the assumption of m ≤ n, is given by where 𝒮 indicates the constant term as in Theorem I and F(· · ·) is the hypergeometric function of several variables with , u = (2 + n − π(k))/2, , , , , . Here, the sum is computed over (m − 1)! permutations of π. S denotes the set of all m! permutations of the set S = {1, 2, ...,m}, and sgn(π) denotes the sign of the permutation π : +1 if π is an even permutation and −1 if it is odd.

Proof: See Appendix section B. In the appendix the moment of sample maximum eigenvalue is also given in more friendlier form for programming obtained by splitting the hypergeometric function of several variables into a sum of its variables.

Dependence of the Covariance Matrix Eigenvalues

Before starting the following analysis, it would be nice to illustrate the dependencies of the eigenvalues on the covariance matrix parameters, which is necessary for the further understanding of the addressed topics. For that, a two-dimensional covariance matrix is used, originated from the expected values of the outer product of the vector k = [k1 k2], i.e., where , are the variance or power of k1 and k2, respectively, ρ is the absolute and ϕ the phase of their complex correlation coefficient, and T means transpose.

The eigenvalues of the true covariance matrix can be shown to be given by [5]: The behavior of l1 and l2 is presented in Figure 1, as a function of a normalized power ratio b and a correlation ρ. Note their mutual dependence, making clear that b and ρ directly determine the behavior of both eigenvalues.

Figure 1.

Two-dimensional system configuration. (a) True eigenvalue l1 and (b) true eigenvalue l2 as a function of correlation and normalized power ratio b. (c) Counter plots of l1.

Validation and Analysis of the Theoretical Expressions

This section aims to validate the theorems mentioned in the previous section using simulated data. The simulated data have been generated using different multi-dimensional configurations, where the correlation between channels have been generated using the well known Mahalanobis transformation [24].

Figure 2(a) shows the comparison of the Equation (4) with simulations. The theoretical PDF curves clearly agree with the histograms obtained from simulated data. As expected, the PDFs become narrower with increasing n, indicating less variance around the true value of the maximum eigenvalue l1 = 1.8. This behavior can be better seen in Figure 2(b) where the distribution of λ as a function of the number of estimation samples n is presented. Note also that the expected value of the distributions seems to change for different n.

Figure 2.

(a) Comparison between the theoretical distributions and the histograms obtained from simulated data of the maximum sample eigenvalue. (b) The distribution of the maximum sample eigenvalue as a function of the number of samples, for a three-dimensional system. In both case, the powers are given by σ = σ = σ = 1 and correlations by ρ = 0, ρ = 0.8 and ρ = 0. When n → ∞, λ = l1 = 1.8.

Figure 3(a) shows the variation of the histogram mean as a function of n for a two-dimensional case with fix correlation ρ = 0.2. In Figure 3(a), the theoretical expected value of λ has been also over plotted for comparison. The curves match well, which validates Theorem III. The same has been carried out for the second order moment, i.e., the variance of λ, and is presented in Figure 3(b). The agreement between the theoretical variance of λ and the variance of the simulated data can again be confirmed. Also, as already observed in Figure 2(a), the lower the number of samples n, the higher the variance of λ.

Figure 3.

Theoretical versus simulation results (a) for the first order statistics (mean), and (b) for the second order statistics (variance), of the sample maximum eigenvalue of a two-dimensional system with powers σ = σ = 1 and correlation ρ = 0.2.

The impact of the correlation and the number of samples on the behavior of the third (skewness) and fourth (kurtosis) order moment of the maximum sample eigenvalue λ is presented in Figure 4. Skewness is a measure of how symmetrical the distribution is with respect to its mean. Figure 4(a,b) indicates that the skewness of λ converges to zero for increasing n and increasing ρ, expressing the tendency to a symmetrical distribution in that cases.

Figure 4.

Theoretical results for the third (skewness) and the fourth (kurtosis) order statistics of the sample maximum eigenvalue of a 2D system having powers σ = σ = 1 and correlations ρ = {0.2, 0.3, . . ., 0.9} versus the number of samples n = {2, 3, . . ., 62}.

Kurtosis is a measure of the peakedness of the distribution. Figure 4(c,d) shows that for increasing n and ρ, the kurtosis of the λ distribution tends to three, which corresponds to the kurtosis of the normal Gaussian distribution.

Estimation Bias

Figure 5 shows the expected value (first order moment) of the maximum sample eigenvalue λ as a function of the true eigenvalues l for a fixed number of samples n = 3 and k = 1. The variation of the underlying correlation of the true covariance matrix ρ, which changes the values of l, has been as well indicated in a color code. Observe that the higher the correlation is, the higher l becomes.

Figure 5.

Effects of the correlation between channels on the expected value of the maximum sample eigenvalue keeping n = 3 and k = 1.

The bias of the estimator of a certain parameter is defined as the difference between the expected value of the estimator and the true value of the parameter. In the present case, the bias is hence given by E[λ] − l. Notice thus from Figure 5 that the bias of λ becomes very strong for low values of the correlation l (or equivalently, low values of ρ).

The effect is emphasized in Figure 6(a,b), where the bias of λ is presented as a function of l and l2, for n = 3 and n = 16, respectively. The values of the underlying correlation ρ have been also indicated. The bias is small either when the correlation is high or the number of samples is sufficiently large. When the correlation between channels is low, the bias becomes very significant and a large number of samples is necessary to decrease the estimation bias.

Figure 6.

The bias, E[λ] − l, of the sample maximum eigenvalue with the number of samples 3 and 16 in various correlated channels having a standard deviation of σ = σ = 1.

Analysis of Detection Problems

Detection theory is a means to quantify the ability of a procedure to detect a parameter (or signal) immersed in a noise environment. A decision has to be taken, in order to say yes, there is a signal, or no, there is no signal. Such decision is usually taken under the application of a decision threshold. Noise contributes off course negatively inducing wrong decisions. For the quantification of detection accuracy, some quantities are usually defined as the probability of detection (PD) and probability of false alarm (PFA), allowing a detection problem analysis.

The maximum eigenvalue of covariance matrices has been frequently used in the literature for the elaboration of certain detection problems. In the Ground Moving Target Indication (GMTI) area, for instance, the signal to clutter plus noise ratio has been studied under the distribution of the sample eigenvalues, which allowed the implementation of a constant false alarm rate detector [9]. The eigenvalues of the covariance matrix of a SAR image pair have been also used in order to formulate a problem in the change detection area [5]. Another example is the determination of the existence of just a single dominant scattering mechanism in a SAR polarimetric acquisition, which also has been evaluated using a threshold in the maximum sample eigenvalue of the polarimetric covariance matrix [11]. Target detection and polarimetric filtering represent other fields of application in which the maximum eigenvalue can be used.

Having the closed form expressions of the PDF (Equation (4)) and/or CDF (Equation (3)) of the maximum sample eigenvalue of the covariance matrix, the PD and PFA when applying a threshold in λ can be analytically computed, allowing a complete detection problem analysis.

Detection of a Dominant Maximum Eigenvalue

Since several detection problems rely in fact on the choice of a dominant or non-dominant l by applying a threshold in λ, the following problem with two hypotheses is elaborated for i = {2, 3, . . ., m}.

In one application area H0 may mean that just a single scattering mechanism is present inside a resolution cell of polarimetric SAR data, in other application H0 can mean that a change happened or a moving target is present in the scene.

Since l is not achievable, a threshold 𝒯 has to be used in λ, originating the following PD and PFA as a function of the decision threshold 𝒯.

Hence, for a two-dimensional system configuration, the distribution p (x; H0) is the distribution of λ when l1 = 2 and l2 = 0. For a three-dimensional system, p (x; H0) is determined evaluating the distribution of λ when l1 = 3 and l2 = l3 = 0. Higher dimensional system configurations follow the same rule.

On the other hand, the distribution p (x; H1) is given by l1 = l2 = 1 and l1 = l2 = l3 = 1 for a two- and three-dimensional system, respectively.

For each value of 𝒯, there exists a pair (PFA, PD). The curves of PD versus PFA are called Receiver Operating Characteristic (ROC) curves and express the detection performance. The more a ROC curve bends toward the upper left, the better is the detection performance since a higher PD and lower PFA is achieved.

Figure 7(a) shows the detection performance as a function of the dimension of the multichannel system, i.e., m = 2, e.g., interferometric, m = 3, e.g., polarimetric and m = 6, e.g., polarimetric-interferometric system. For all cases, with fixed number of samples n = 6, one can realize that the performance of detection is significantly improved as the number of SAR images increases. Figure 7(b) shows the ROC curves for m = {2, 3} and for different number of samples n. It can be seen that for both multidimensional system configurations the number of samples increases the detection performance.

Figure 7.

Maximum sample eigenvalue detection performance as a function of the number or channels m, samples n and correlation ρ. (a) ROC curves for different multidimensional systems, m = {2, 3, 6} and n = 6. (b) ROC curves for the two- and three-dimensional cases, and for different number of samples n. (c) ROC curves for the two-dimensional case, n = 3, and for different correlation.

In order to state how correlation effects the detection performance, the detection problem is reformulated as follows In this way, every time that correlation is greater than zero the eigenvalues have different values, and one is larger than the others. For the two-dimensional case, for instance, p (x; H0) is the distribution of λ when l1 > l2, which changes for different values of the correlation ρ. Figure 7(c) presents the ROC curves for this case. For small correlation between channels, the detector suffers from a significant false-alarm rate. As expected, when correlation increases the ROC curve bends toward the upper left, indicating better detection performance. The limit is reached when ρ = 1, meaning that l1 = 2 and l2 = 0, which corresponds to the cases of the previous detection problem.

Target Detection Using Polarimetric Matched Filter

In this section the application of the expressions for the detection of specific targets using the Polarimetric Matched Filter (PMF) concept is illustrated. For that, a short review on PMF is required [6,7].

In a polarimetric acquisition, the dimension of the system corresponds to the number of channels, which are in general four, but three for the backscattering reciprocal case. The measured vector k is hence three- or four-dimensional. The elements of the vector k may be weighted in a given way, in order to satisfy a certain condition. The weighting can be accounted for by making y = h†k, where h is a complex vector with same dimension as k. A condition usually aimed to be satisfied in the literature is the maximization of the quadratic detector |y|2 (Equation (59) in [6]). The optimal weighting vector h in order to detect a distributed target with the target vector k is dependent on the target characteristics. Accordingly, the expected value of |y|2 is given by where Σ = E{kk†} is the target covariance matrix. Regarding Rayleigh quotient (Theorem 15.91 in [25]), for any m-dimensional complex vector x and a given m × m Hermitian matrix A, x†Ax ≤ ||x||2a, where a is the maximum eigenvalue of A. The equality is valid if x is along the direction of the a eigenvector U (||U|| = 1). Hence, for a given distributed target with covariance matrix Σ, the optimum weighting vector h is given by the eigenvector of the maximum eigenvalue of Σ.

The target vector k is a random sample and has alone no physical meaning. Therefore, usually multi-look processing is performed making, for n samples (or looks) The target is assumed to be immersed in polarization independent clutter. In this way, the detection procedure is thus evaluated by choosing h as the eigenvector corresponding to the maximum eigenvalue of the covariance matrix Σ of the target to be detected, and applying a detection threshold y̅2 > 𝒯.

In the presence of target y̅2 is maximum, given by y̅2 = ||h̅||2λ, where h̅ is the eigenvector of λ. Hence, assuming without loss of generality that h̅ has unitary length (||h̅||2 = 1), detection performances can be made using the derivations given in this work, as done in the previous section, where a threshold has been used in order to detect a dominant maximum sample eigenvalue (λ = y̅2 > 𝒯).

Note that the probabilities of detection and false alarm defined in last section can be more simply evaluated by

A three-dimensional polarimetric target detection problem was formulated making k = [S S S], where S is the polarimetric scattering matrix element in channel i, and using the target covariance matrix structure where , and ρ is the complex correlation coefficient between S and S [6,26]. Figure 8 shows the ROC curves for different number of samples for three different kinds of targets. Figure 8(a) corresponds to an azimuthal symmetric target having covariance matrix with parameters ɛ = 1, γ = 0.5 and ρ = 0. Figure 8(c) corresponds to a reflection symmetric target having covariance matrix with parameters ɛ = 1, γ = 0.8 and ρ = 0.8, while Figure 8(b) to a reflection symmetric target with parameters ɛ = 0, γ = 0.8 and ρ = 0.8 in its covariance matrix. For all three cases, S = 1 and the PFA was evaluated using the polarization independent clutter having ɛ = 1, γ = 1 and ρ = 0. Note that the curves vary not just with the number of used samples n but are also different for different targets having different Σ. This means that some types of targets are easier to detect than others, when using the quadratic detector described here and when the clutter is polarization independent. When the target covariance matrix is similar to the one of the clutter, the detection performance weakens (Figure 8(a)). On the other hand, when the covariance matrix of the target is significantly different from the clutter one, the detection performance improves (Figure 8(c)).

Figure 8.

Detection performance of three different distributed scatterers using the PMF concept for different number of samples. (a) Azimuthal symmetric scatterer with ɛ = 1, γ = 0.5 and ρ = 0. (b) Reflection symmetric scatterer with ɛ = 1, γ = 0.8 and ρ = 0.8. (c) Reflection symmetric scatterer with ɛ = 0, γ = 0.8 and ρ = 0.8. The PFA is evaluated using the polarization independent clutter with ɛ = 1, γ = 1 and ρ = 0.

Conclusions and Discussion

In this paper a depth statistical analysis of the maximum eigenvalue of the eigendecomposition of the sample covariance matrix in terms of SAR applications is presented. The proposed analysis are supported by simulation results via several examples. The results are based on a exact closed-form expressions of PDF, CDF and MGF. In this study, existing density functions of the sample maximum eigenvalue were extended and/or implemented into multi-channel SAR system in order to obtain a simple expression of the sample eigenvalues giving a way to fruitful applications. From these closed-form expressions, it has been possible to develop new algorithms for unbiased calculations of parameters extracted from multi-channel SAR covariance matrix. In addition to these implementations, closed-form expressions were developed for the MGF of the sample maximum eigenvalue, which can be critical in the area of bias removal and detection performance analysis. This new closed-form expressions of the MGF can be also interesting for other application areas like MIMO systems (Multiple-Input Multiple-Output). Apart from estimation theory analysis including the MGF, the detection problem of the sample maximum eigenvalue has been also discussed.