Polarization Optimization for the Detection of Multiple Persistent Scatterers Using SAR Tomography

Abstract

The detection of multiple interfering persistent scatterers (PSs) using Synthetic Aperture Radar (SAR) tomography is an efficient tool for generating point clouds of urban areas. In this context, detection methods based upon the polarization information of SAR data are effective at increasing the number of PSs and producing high-density point clouds. This paper presents a comparative study on the effects of the polarization design of a radar antenna on further improving the probability of detecting persistent scatterers. For this purpose, we introduce an extension of the existing scattering property-based generalized likelihood ratio test (GLRT) with realistic dependence on the transmitted/received polarizations. The test is based upon polarization basis optimization by synthesizing all possible polarimetric responses of a given scatterer from its measurements on a linear orthonormal basis. Experiments on both simulated and real data show, by means of objective metrics (probability of detection, false alarm rate, and signal-to-noise ratio), that polarization waveform optimization can provide a significant performance gain in the detection of multiple scatterers compared to the existing full-polarization-based detection method. In particular, the increased density of detected PSs at the studied test sites demonstrates the main contribution of the proposed method.

1. Introduction

Synthetic Aperture Radar (SAR) images have been used and studied in interferometric analyses for decades to map the elevations and displacements of the earth’s surface [1]. The temporal decorrelation of SAR image stacks, acquired over repeated orbits, limits interferometric analyses to specific pixels, i.e., persistent scatterers (PSs), since the latter exhibit scattering behaviour that is stable over time. Typically, the selection of valid pixels and the detection of reliable scatterers in each resolution cell is an intense research topic, where the goal is to find as many valid scatterers as possible at the best spatial resolution. With respect to interferometric techniques, SAR Tomography (TomoSAR) has introduced a different strategy for processing a SAR data stack and has supported the evolution of interferometry in three-dimensional (3D) [2,3], four-dimensional [4,5], and five-dimensional [6,7]) applications. In contrast to interferometry, TomoSAR offers the possibility of identifying multiple permanent scatterers superimposed in a single resolution cell, thus significantly increasing the density of detected scatterers.

In analogy to time series interferometry, e.g., PS interferometry [8], TomoSAR relies on the persistence of the scatterers during the observation period. Although the tomographic model of distributed sources [9] can mitigate the short-term effects of decorrelating sources, distinguishing between reliable scatterers and false alarms is a daunting task, especially when dealing with data obtained over a long time interval.

The key point in TomoSAR PS detection is to identify and increase the number of detectable stable scatterers in the presence of noise. At the same time, the possibility of signal misinterpretation, i.e., false alarms [5], needs to be minimized. Currently, efforts are still underway to develop modern and advanced methods to further increase the detection probability of persistent scatterers.

In the framework of SAR tomography, the detection of reliable scatterers can be generally addressed with detection methods that control the false alarm rate (FAR), such as the generalized likelihood ratio test (GLRT). GLRT has been used extensively, not only in tomographic PS recognition but also in various target recognition applications [10,11]. Within the GLRT scheme, the detection of individual or single scatterers present in a resolution cell is straightforward [12,13]. The critical aspect is to distinguish between multiple scatterers and especially to distinguish between single and double scatterers. Among the different strategies, there are two successful methods commonly used in the literature for characterizing multiple scatterers lying in the same range azimuth pixel: (1) sequential GLRT with cancellation (SGLRTC) [12] and (2) support GLRT [13], which is based on estimating the support of the unknown signal that best fits the data. SGLRTC adopts a cancellation strategy that first cancels the dominant scatterer and then evaluates the presence of the secondary scatterer in the residual signal. Support GLRT is a sequential test in which the presence of scatterers is detected first, followed by the discrimination of multiple scatterers. In other words, testing the signal strength in sub-spaces with higher dimensionality is achieved by sequentially adding a single direction, starting with the first one. When comparing the two approaches, SGLRTC is more efficient in terms of computational requirements, while support GLRT has a super-resolution capability that enables the detection of multiple layovered scatterers at distances below the Rayleigh resolution. To improve the computational efficiency of support GLRT, ref. [14] proposed a fast implementation that is able to maintain the super-resolution capability.

Spatial averaging or multi-looking operations, through the estimation of a second-order data covariance matrix, represent a solid and significant strategy for improving the detection capability of typical GLRT-based methods over the scene and, in particular, over regions characterized by a low signal-to-noise ratio (SNR) [15,16]. The use of the covariance matrix of the data can improve the detection probability at the cost of a small reduction in spatial resolution. Indeed, by trading off spatial resolution, undeniable improvements have been reported in the context of SAR tomography for the detection of multiple permanent scatterers. In addition, the potential of polarimetric SAR (PolSAR) data is expected to further improve detection performance. The authors of [17,18] have shown that increasing the length of the observation vector as input to support GLRT [13] using dual polarization data can mitigate the problem of a reduced number of available baselines. To take full advantage of PolSAR data, however, a method has been developed in [19] (denoted as scattering property-based support GLRT) that incorporates polarimetric information into the detection test. The aim is not only to improve detection performance but also to allow the test to identify PS backscattering mechanisms. Indeed, this method can be considered as a first attempt at the classification of PSs besides the detection and the extraction of their parameters related to the elevation, displacement, or thermal dilation in 3D, 4D, and 5D applications.

This paper presents a comparative study to answer the question of which polarization design for a radar antenna provides the best detection probability. To this end, we introduce an extension of the scattering support GLRT [19] that takes the dependence on the polarization basis into account. The analysis is based upon the evaluation of the detection performance in all synthesized polarizations obtained from polarimetric measurements on a linear orthonormal basis. The main consequence of the analysis is the complete characterization of the studied scatterers. This allows the selection of the optimal polarization, i.e., the polarization that provides the highest detection probability. The identification of the best polarization basis is, in fact, the main contribution of the proposed method, which consequently leads to an increase in the density of the detected permanent scatterers and provides a better possibility of modeling the three-dimensional structures of the sensed objects and their possible slow deformation and thermal dilation. Additionally, the proposed method paves the way for comprehending the detection performance in TomoSAR data in terms of the transmitted and received polarizations used, thus complementing polarimetric-based PS detection studies.

The paper is organized as follows. Section 2 begins with the basic principles of SAR polarimetry, while an overview of existing scattering support GLRT is given in Section 3. After these introductory descriptions, the proposed method is presented in Section 4. Experimental results using simulated and real data sets from the Experimental SAR System (ESAR) of the German Aerospace Center and Uninhabited Aerial Vehicle SAR (UAVSAR) are presented in Section 5. The discussion and conclusions regarding the reported results are provided in Section 6 and Section 7.

2. SAR Polarimetry

The polarization information in SAR images is related to both the geometrical structure and the orientation of the observed objects [20]. The characteristic information of any target sensed by a full polarimetric radar sensor can be represented in a $2\times 2$ scattering matrix $\mathbf{S}$.

$$\mathbf{S}=\left[\begin{array}{cc}{x}_{hh}& x_{hv}\\ x_{vh}& x_{vv}\end{array}\right],$$

(1)

where the complex scattering coefficient $x_{hv}$ represents the electromagnetic wave emitted through the polarization channel v and received on channel h. In the case of interchanging the role of the transmitting and receiving antennas, i.e., with the reciprocity theorem, the scattering matrix is symmetric, i.e., $x_{hv}=x_{vh}$. Once the scattering matrix in (1) is known, a key point is that radar cross-sections can be synthesized for any arbitrary combination of transmit and receive polarizations [21]. In this study, we will consider this point to determine the optimal transmit and receive polarizations that offer the best PS detection performance.

Let us now consider the geometry of the multi-baseline, multi-temporal polarimetric SAR imaging system that operates at wavelength $\lambda$ and range distance r. In this geometry, referring to a resolution cell, the stack data vector $\mathbf{x}$, which collects the information from N polarimetric sensors, can be represented as Equation (2) [22,23]. The data stack is acquired under interferometric conditions and co-registered with respect to a given master image and compensated for atmospheric effects and possible small-scale nonlinear deformations. Each polarimetric sensor is assumed to acquire images in three linear polarization channels (hh, hv, vv) with reciprocal condition. The multi-baseline imaging geometry thus yields $3N$ image channels or $\mathbf{x}\in {\mathbb{C}}^{3N\times 1}$, where ${\mathbb{C}}^{3N\times 1}$ represents the set of $3N\times 1$ complex vectors.

$$\mathbf{x}\left(l\right)=\mathbf{A}({\Theta }_M,{\mathbf{K}}_M){\mathbf{s}}_M\left(l\right)+\mathbf{n}\left(l\right)$$

(2)

In (2), ${\mathbf{s}}_M\left(l\right)\in {\mathbb{C}}^{M\times 1}$ is the lth realization, i.e., look, of the complex reflectivity vector of M superimposed scatterers into the resolution cell, which is invariant with respect to the different polarizations. Moreover, the matrix $\mathbf{A}({\Theta }_M,{\mathbf{K}}_M)$ equals to $[\mathbf{a}({\mathbf{\theta }}_1,{\mathit{k}}_1)\dots \mathbf{a}({\mathbf{\theta }}_M,{\mathit{k}}_M)]\in {\mathbb{C}}^{3N\times M}$ with $\mathbf{a}({\mathbf{\theta }}_m,{\mathbf{k}}_m)={\mathbf{k}}_m\otimes \mathbf{a}\left({\mathbf{\theta }}_m\right)$, where $\mathbf{a}\left({\mathbf{\theta }}_m\right)$ represents the steering vector of the mth superimposed scatterer, where $m\in \{1,2,\dots ,M\}$, and ${\mathbf{\theta }}_m$ is the parameter vector representing the mth scatterer position in the elevation direction in the 3D case, or the position and mean deformation velocity in the 4D case, or the position, mean deformation velocity, and thermal dilation coefficient in the 5D case:

$$\begin{array}{c}\mathbf{a}\left({\mathbf{\theta }}_m\right)=exp\left(i\frac{4\pi }{\lambda }\frac{\mathbf{b}}{rsin\left(\beta \right)}{z}_m\right)\left(\mathit{in}\mathit{3}\mathit{D}\mathit{case}\right)\hfill \\ \mathbf{a}\left({\mathbf{\theta }}_m\right)=exp(i\frac{4\pi }{\lambda }\frac{\mathbf{b}}{rsin\left(\beta \right)}{z}_m+i\frac{4\pi }{\lambda }\mathbf{t}{v}_m)\left(\mathit{in}\mathit{4}\mathit{D}\mathit{case}\right)\hfill \\ \mathbf{a}\left({\mathbf{\theta }}_m\right)=exp(i\frac{4\pi }{\lambda }\frac{\mathbf{b}}{rsin\left(\beta \right)}{z}_m+i\frac{4\pi }{\lambda }\mathbf{t}{v}_m+i\frac{4\pi \mathbf{T}}{\lambda }{\gamma }_m).\left(\mathit{in}\mathit{5}\mathit{D}\mathit{case}\right)\hfill \end{array}$$

(3)

In (3), $z_m$, $v_m$, and ${\gamma }_m$ denote the elevation, velocity, and thermal dilation, respectively, of the mth scatterer, and the vectors $\mathbf{b}$, $\mathbf{t}$, and $\mathbf{T}$ indicate the spatial, temporal, and thermal baselines in multi-baseline multi-temporal acquisition mode. Moreover, $\beta$ is the incident angle, $exp$ is the exponential operator, and i is the standard imaginary unit. Additionally, the operator ⊗ indicates the Kronecker product, and ${\mathit{k}}_m\in {\mathbb{C}}^{3\times 1}$ is the unitary polarimetric target vector representing the polarimetric backscattering pattern of the mth scatterer.

From a statistical point of view, the data vector $\mathbf{x}$ in (2) is typically modeled as a zero-mean complex circular Gaussian random vector, with covariance matrix:

$$\mathbf{R}=E\left\{\mathbf{x}{\mathbf{x}}^{†}\right\}=\mathbf{A}({\Theta }_M,{\mathbf{K}}_M)diag\left\{\right|{\mathbf{s}}_M{|}^2\}\mathbf{A}{({\Theta }_M,{\mathbf{K}}_M)}^{†}+{\sigma }_n^2{\mathbf{I}}_{3N},$$

(4)

where ${\sigma }_n$ is the standard deviation of the additive noise n in (2), ${\mathbf{I}}_{3N}$ is a $3N\times 3N$ identity matrix, and $diag\{.\}$ is a diagonal operator. Moreover, ${}^{†}$ and E are Hermitian and expectation operators, respectively. Multi-look SAR tomography is generally concerned with reconstructing the backscattering distribution ($\mathbf{g}= |{\mathbf{s}}_M|$) for each resolution cell using the covariance matrix estimated by a set of L independent and homogeneous looks. Here, the multi-look operation is used, since estimating the covariance matrix of the data with a single look gives a poor estimate.

The estimation of $\mathbf{g}$ can be approached in the framework of spectral estimation theory [24], and beamforming is the simple method that recovers the backscattering distribution from a given scatterer while limiting the other interfering spectral components. Therefore, the recovered power of the filtered signal from the mth scatterer is given by [25]:

$$\mathbf{g}\left({\mathbf{\theta }}_m\right)=\underset{{\mathbf{k}}_m}{max}\mathbf{a}{({\mathbf{\theta }}_m,{\mathit{k}}_m)}^{†}\mathbf{R}\mathbf{a}({\mathbf{\theta }}_m,{\mathit{k}}_m).$$

(5)

Since different realizations of an image are usually not available, the actual covariance matrix R is unknown. It is straightforward, however, to approximate R by the sample covariance matrix $\widehat{\mathbf{R}}$ by spatially averaging L similar pixels in the images that share statistical similarity.

$$\widehat{\mathbf{R}}=\frac{{\sum }_{l=1}^{L}w\left(l\right)\mathbf{x}\left(l\right){\mathbf{x}}^{†}\left(l\right)}{{\sum }_{l=1}^{L}w\left(l\right)},$$

(6)

where $w\left(l\right)$ indicates the relative importance of the lth pixel during averaging. Currently, the developed tomographic reconstruction methods are affected by decorrelating sources and the temporal decoherence of the multi-temporal data. This indeed leads to defocusing and blurring effects in the reconstruction. In such a situation, it is important to distinguish the reliable scatterers and false alarms, e.g., using detection methods such as GLRT. The next section gives an overview of the scattering-based support GLRT.

3. Overview of Scattering-Based Support GLRT

To solve the problem of detecting reliable and genuine scatterers, scattering support GLRT is among the most competitive techniques that allow for the identification of scatterers by controlling the false alarm rate [13,19]. According to the signal model in (2), the multiple hypothesis tests are given as [19]:

$$\begin{array}{c}{\mathcal{H}}_0:\mathbf{x}\left(l\right)=\mathbf{n}\left(l\right)\hfill \\ {\mathcal{H}}_1:\mathbf{x}\left(l\right)=\mathbf{A}({\Theta }_1,{\mathbf{K}}_1){\mathbf{s}}_1\left(l\right)+\mathbf{n}\left(l\right)\hfill \\ {\mathcal{H}}_2:\mathbf{x}\left(l\right)=\mathbf{A}({\Theta }_2,{\mathbf{K}}_2){\mathbf{s}}_2\left(l\right)+\mathbf{n}\left(l\right)\hfill \\ ⋮\hfill \\ {\mathcal{H}}_m:\mathbf{x}\left(l\right)=\mathbf{A}({\Theta }_m,{\mathbf{K}}_m){\mathbf{s}}_m\left(l\right)+\mathbf{n}\left(l\right)\hfill \end{array}$$

(7)

where the hypothesis ${\mathcal{H}}_0$ assumes the absence of a scatterer, i.e., the observed data $\mathbf{x}\left(l\right)$ are a circularly symmetric Gaussian random vector with zero mean and covariance matrix ${\sigma }_n^2\mathbf{I}$, while ${\mathcal{H}}_m$ assumes the presence of m scatterers whose parameters are in the unknown matrix ${\Theta }_m=[{\mathbf{\theta }}_1\dots {\mathbf{\theta }}_m]$, while ${\mathbf{s}}_m$ contains the reflectivity of the m scatterers and ${\mathbf{K}}_m=[{\mathbf{k}}_1\dots {\mathbf{k}}_m]$. Moreover, $\mathbf{A}\left({\Theta }_m\right)$ is formed by their m corresponding steering vectors.

If ${\mathcal{H}}_0$ is rejected in favor of ${\mathcal{H}}_m$, then the observed data can be assumed to follow a circularly symmetric Gaussian random vector with zero mean and covariance matrix R given in (4). Here, we assume the presence of up to two scatterers per resolution cell, i.e., $m=2$. Although the identification of m scatterers is generally and theoretically possible, detection of up to two scatterers is challenging in practice due to the limited number and distribution of baselines and may not be feasible, especially when dealing with areas with short buildings and no skyscrapers, such as our study area. With this assumption, the likelihood ratio for the binary test (${\mathcal{H}}_0$, $\overline{{\mathcal{H}}_0}$) equals:

$$\frac{\underset{{\sigma }_n^2}{max}f(\mathbf{x}\left(1\right),\dots ,\mathbf{x}\left(l\right),{\sigma }_n^2|{\mathcal{H}}_0)}{\underset{{\Theta }_2{\mathbf{K}}_2,{\mathbf{s}}_2,{\sigma }_n^2}{max}f(\mathbf{x}\left(1\right),\dots ,\mathbf{x}\left(l\right),{\Theta }_2,{\mathbf{K}}_2,{\mathbf{s}}_2,{\sigma }_n^2|\overline{{\mathcal{H}}_0})}\stackrel{{\mathcal{H}}_0}{\underset{\overline{{\mathcal{H}}_0}}{\gtrless }}th,$$

(8)

where f is the joint probability density function of the looks $\mathbf{x}\left(l\right)$, e.g., a Gaussian zero-mean model. In (8), $\overline{{\mathcal{H}}_0}$ means that the hypothesis of the absence of the scatterer (${\mathcal{H}}_0$) is rejected and the data stack consists of either one (${\mathcal{H}}_1$) or two (${\mathcal{H}}_2$) scatterers. By substituting the maximum likelihood estimation (MLE) of the unknowns (${\mathbf{s}}_2,{\sigma }_n^2$) into (8), the first stage of the scattering property-based support GLRT [19], which serves as a decision rule for the presence or absence of scatterers, can be written as follows:

$$\frac{\underset{{\Theta }_2,{\mathbf{K}}_2}{max}trace\{{\mathbf{D}}^{†}({\Theta }_{\mathbf{2}},{\mathbf{K}}_2)\widehat{\mathbf{R}}\mathbf{D}({\Theta }_{\mathbf{2}},{\mathbf{K}}_2\}}{trace\{\widehat{\mathbf{R}}\}}\stackrel{{\overline{\mathcal{H}}}_0}{\underset{{\mathcal{H}}_0}{\gtrless }}th,$$

(9)

where $\mathbf{D}({\Theta }_2,{\mathbf{K}}_2)=\mathbf{A}({\Theta }_2,{\mathbf{K}}_2){\mathbf{V}}_d$, in which the matrix ${\mathbf{V}}_d$ is satisfying the equality ${\mathbf{V}}_d{\mathbf{V}}_d^{†}={\left({\mathbf{A}}^{†}({\Theta }_2,{\mathbf{K}}_2)\mathbf{A}({\Theta }_2,{\mathbf{K}}_2)\right)}^{-1}$, and obtained by Cholesky decomposition. The threshold $th$ in the above equations can be defined using the FAR approach [13]. It can be shown that the solution of the maximization in (9) with respect to ${\mathbf{K}}_2$ is equivalent to the maximum eigenvalue of the following matrix:

$$\frac{{\lambda }_{max}({\mathbf{B}}^{†}\left({\Theta }_{\mathbf{2}}\right)\widehat{\mathbf{R}}\mathbf{B}\left({\Theta }_{\mathbf{2}}\right))}{trace\{\widehat{\mathbf{R}}\}}\stackrel{{\overline{\mathcal{H}}}_0}{\underset{{\mathcal{H}}_0}{\gtrless }}th,$$

(10)

where ${\lambda }_{max}\left(\mathbf{X}\right)$ and $trace\left(\mathbf{X}\right)$ are, respectively, the maximum eigenvalue and trace of matrix X, $\mathbf{B}\left({\Theta }_{\mathbf{2}}\right)=[{\mathbf{I}}_3\otimes \mathbf{a}\left({\mathbf{\theta }}_{\mathbf{1}}\right){\mathbf{I}}_3\otimes \mathbf{a}\left({\mathbf{\theta }}_{\mathbf{2}}\right)]{\mathbf{V}}_b={\mathbf{U}}_b{\mathbf{V}}_b$, and ${\mathbf{V}}_b{\mathbf{V}}_b^{†}={\left({\mathbf{U}}_b^{†}{\mathbf{U}}_b\right)}^{-1}$, and ${\mathbf{I}}_3$ is a $3\times 3$ identity matrix.

The second stage of the scattering support GLRT, instead, operates the decision according to the rule of the presence of single or double scatterers $({\mathcal{H}}_1,{\mathcal{H}}_2)$ [19]:

$$\frac{1-\frac{{\lambda }_{max}({\mathbf{B}}^{†}\left({\Theta }_{\mathbf{2}}\right)\widehat{\mathbf{R}}\mathbf{B}\left({\Theta }_{\mathbf{2}}\right))}{trace\{\widehat{\mathbf{R}}\}}}{1-\frac{{\lambda }_{max}({\mathbf{B}}^{†}\left({\Theta }_{\mathbf{1}}\right)\widehat{\mathbf{R}}\mathbf{B}\left({\Theta }_{\mathbf{1}}\right))}{trace\{\widehat{\mathbf{R}}\}}}\stackrel{{\mathcal{H}}_1}{\underset{{\mathcal{H}}_2}{\gtrless }}th.$$

(11)

In analogy to the first stage of scattering support GLRT, Equation (11) is derived by substituting the unknown parameters into the likelihood ratio of the binary test $({\mathcal{H}}_1,{\mathcal{H}}_2)$. The parameters ${\Theta }_1=\left[{\mathbf{\theta }}_{\mathbf{1}}\right]$ and ${\Theta }_2=\left[{\mathbf{\theta }}_{\mathbf{1}}{\mathbf{\theta }}_2\right]$ can be obtained by the MLE approach, and the fast implementation allows decoupled MLEs as:

$$\begin{array}{c}{{\displaystyle \hat{\mathbf{\theta }}}}_1=\underset{{\mathbf{\theta }}_1}{max}{\lambda }_{max}\left({\mathbf{B}}^{†}\left({\mathbf{\theta }}_{\mathbf{1}}\right){\displaystyle \hat{\mathbf{R}}}\mathbf{B}\left({\mathbf{\theta }}_{\mathbf{1}}\right)\right)\hfill \\ {{\displaystyle \hat{\mathbf{\theta }}}}_2=\underset{{\mathbf{\theta }}_2}{max}{\lambda }_{max}\left({\mathbf{B}}^{†}({{\displaystyle \hat{\mathbf{\theta }}}}_1,{\mathbf{\theta }}_2)\widehat{\mathbf{R}}\mathbf{B}({{\displaystyle \hat{\mathbf{\theta }}}}_1,{\mathbf{\theta }}_2)\right).\hfill \end{array}$$

(12)

In [19], it is shown that the scattering support GLRT (10) and (11) outperforms the conventional support GLRT [18] in the detection of both single and double scatterers. The attempt of the current work is to further improve the detection performance by introducing an extension to the scattering support GLRT method.

4. The Proposed Method

Polarimetric radar systems typically measure complex backscattering of any target on a linear horizontal and vertical orthogonal polarization basis $\{h,v\}$. This is commonly equivalent to the acquisition of images in the polarization channels hh, hv, vh, and vv. Therefore, in the case of multi-baseline polarimetric data, the response of a reciprocal medium can be represented by the vector in Equation (2) $\mathbf{x}={\left[{\mathbf{x}}_{hh}{\mathbf{x}}_{hv}{\mathbf{x}}_{vv}\right]}^T\in {\mathbb{C}}^{3N\times 1}$ that has $3N$ unique elements, where ${\mathbf{x}}_{hv}\in {\mathbb{C}}^{N\times 1}$ is the complex backscattering vector from the vertical polarized transmitted signal and the horizontal polarized return in all baselines.

Since any orthogonal set of elliptically polarized states can form a polarization basis, the polarimetric response vector can be represented in any arbitrary orthogonal elliptic basis $\{p,q\}$, where q is the orthogonal complement of p. In this way, a complete characterization of the backscattering of targets can be obtained using the polarization synthesis theorem [21]. Thus, the multi-baseline polarimetric response vector in Equation (2) can be transformed from the linear basis $\{h,v\}$ (or polarization channels {hh, hv, vv}) to a desired $\{p,q\}$ basis (or polarization channels {pp, pq, qq}) as [26]:

$$\mathbf{x}{\left(l\right)}^{\{p,q\}}={\left[{\mathbf{x}}_{pp}{\mathbf{x}}_{pq}{\mathbf{x}}_{qq}\right]}^T=\mathbf{W}\mathbf{A}({\Theta }_m,{\mathbf{K}}_m){\mathbf{s}}_m\left(l\right)+\mathbf{n}\left(l\right),$$

(13)

where

$$\mathbf{W}=\frac{1}{1+\rho {\rho }^{†}}\left[\begin{array}{ccc}1& \sqrt{2}\rho & {\rho }^2\\ -\sqrt{2}{\rho }^{†}& 1-\rho {\rho }^{†}& \sqrt{2}\rho \\ {\rho }^{†}{}^2& -\sqrt{2}{\rho }^{†}& 1\end{array}\right]\otimes {\mathbf{I}}_N.$$

(14)

In (14), $\rho$ is a complex polarization ratio that represents the polarization state of an electromagnetic wave with a specific orientation ($\chi$) and ellipticity ($\tau$) as:

$$\rho =\frac{cos\left(2\chi \right)sin\left(2\tau \right)+isin\left(2\chi \right)}{1+cos\left(2\chi \right)cos\left(2\tau \right)}.$$

(15)

Accordingly, the first stage of the scattering support GLRT in (10), for any specific polarization basis, can be revised as follows:

$$\frac{{\lambda }_{max}({\mathbf{B}}^{†}\left({\Theta }_{\mathbf{2}}\right)\mathbf{W}\widehat{\mathbf{R}}{\mathbf{W}}^{†}\mathbf{B}\left({\Theta }_{\mathbf{2}}\right))}{trace\{\mathbf{W}\widehat{\mathbf{R}}{\mathbf{W}}^{†}\}}\stackrel{{\overline{\mathcal{H}}}_0}{\underset{{\mathcal{H}}_0}{\gtrless }}th,$$

(16)

while the second stage of the test is re-written as:

$$\frac{1-\frac{{\lambda }_{max}({\mathbf{B}}^{†}\left({\Theta }_{\mathbf{2}}\right)\mathbf{W}\widehat{\mathbf{R}}{\mathbf{W}}^{†}\mathbf{B}\left({\Theta }_{\mathbf{2}}\right))}{trace\{\mathbf{W}\widehat{\mathbf{R}}{\mathbf{W}}^{†}\}}}{1-\frac{{\lambda }_{max}({\mathbf{B}}^{†}\left({\Theta }_{\mathbf{1}}\right)\mathbf{W}\widehat{\mathbf{R}}{\mathbf{W}}^{†}\mathbf{B}\left({\Theta }_{\mathbf{1}}\right))}{trace\{\mathbf{W}\widehat{\mathbf{R}}{\mathbf{W}}^{†}\}}}\stackrel{{\mathcal{H}}_1}{\underset{{\mathcal{H}}_2}{\gtrless }}th.$$

(17)

In analogy to the scattering-based support GLRT, the threshold $th$ in the proposed method can be calculated using the FAR approach (see [13]). Note that the sample covariance matrix $\widehat{\mathbf{R}}$ is still computed on the basis ($\{h,v\}$) of the original data set, while the matrix W is in charge of changing the polarization basis of the original data to any arbitrary transmitted/received wave by radar antenna. The detection tests in (16) and (17) are a function of $\Theta$ and the geometrical parameters $\chi$ and $\tau$. Therefore, ${\mathbf{\theta }}_{\mathbf{1}}$ and ${\mathbf{\theta }}_2$ are estimated with (12), while the MLEs of $\chi$ and $\tau$ are represented as:

$$\hat{\chi },\widehat{\tau }=\underset{\chi ,\tau }{max}{\lambda }_{max}({\mathbf{B}}^{†}(\hat{\Theta })\mathbf{W}\widehat{\mathbf{R}}{\mathbf{W}}^{†}\mathbf{B}({\displaystyle \hat{\Theta }})).$$

(18)

5. Results

5.1. Numerical Examples

This subsection is devoted to evaluating the detection performance of the proposed method using numerical experiments. Thus, a tomographic data stack is simulated by using the system parameters of a dual-baseline $N=3$ polarimetric ESAR sensor with baseline distributions and system parameters that correspond to our real data sets, which are explained in the next subsection. The data are simulated in the presence of a single scatterer with parameter vector ${\mathbf{\theta }}_{\mathbf{1}}={[z_1,v_1,{\gamma }_1]}^T={[0,0,0]}^T$, while the additional scatterer is given by ${\mathbf{\theta }}_{\mathbf{2}}={[z_2,v_2,{\gamma }_2]}^T={[1.5ru,0,0]}^T$, where $ru$ is the vertical Rayleigh resolution. The total signal power is defined by ${\sigma }_s={\sigma }_{s_1}+{\sigma }_{s_2}$, where ${\sigma }_{s_i}$ is the polarimetric signal power of the ith scatterer and ${\sigma }_{s_2}/{\sigma }_{s_1}=0.8$. The probabilities of detection (PD) and false alarm (PFA) are evaluated using Monte Carlo techniques over ${10}^5$ realizations of simulated data. Furthermore, the parameter spaces $\mathbf{\theta }$ are scanned in a uniformly discretized interval of $[-4ru,4ru]$ in all domains, and the geometrical parameters ($\chi$ and $\tau$) are respectively scanned in uniformly discretized intervals of $\left[0180\right]$ and $[-4545]$ with increments of 1 degree. Note that the detectors are applied with $L=16$, where the number of looks forms a $4\times 4$ average window for sample covariance matrix estimation. In this way, the relative importance of all samples is set to be comparable, i.e., $w\left(l\right)=1,l=1,2,\dots ,L$. In the following, the detection performance of the proposed approach (Equations (16) and (17)) is compared with scattering support GLRT [19] (Equations (10) and (11)). The comparison is made using the evaluation of the probability of scatterer detection when the data stack is simulated with different signal-to-noise ratios.

In the first experiment shown in Figure 1, the detection performance of the proposed method is evaluated for both single and double scatterers with simulated data. The performance of the proposed tests in (16) and (17) is compared with the existing scattering support GLRT given by (10) and (11). Two action lines are followed below to evaluate the results in Figure 1.

5.1.1. Performance of the First Stage of the Tests

The probabilities of detection and false alarm, as the measures of the performance assessment, indicate the probability of deciding for the specific hypothesis when it is true or it is false, respectively [5]. In order to evaluate the capability of the first step of the proposed test in the hypothesis $({\mathcal{H}}_0,{\mathcal{H}}_2)$, the Receiver Operating Characteristic (ROC) curve of PD and PFA, as well as the plot of PD versus SNR, are computed. In Figure 1a, two ROC curves, both for the same data set with a fixed SNR = 1 dB, are given. As can be seen, the red curves correspond to the proposed approach and outperform the blue curves corresponding to existing scattering support GLRT. For more experiments, the simulations are extended by varying the SNR up to 16 dB, and the probability of scatterer detection using (16) and (10) are computed and plotted in Figure 1b. Again, the superiority of the proposed detection approach over the scattering support GLRT is evident for different signal to noise powers.

5.1.2. Performance of the Second Stage of the Tests

The performance of the second stage of the tests in (17) and (11) is evaluated in the same way as with PD and PFA. The ROC curves in Figure 1c (when data SNR = 1 dB) show the superior capability of the proposed detector compared to scattering support GLRT in discriminating the two superimposed scatterers. Moreover, the extensive experiments are reported by giving the PD of the detectors with respect to different SNR values in Figure 1d. As expected, by increasing the SNR, the probability of detection is improved for both detectors. Thus, in this case as well, the proposed detector (17) significantly outperforms the scattering support GLRT (11).

5.2. Real Polarimetric Data

5.2.1. Experiment Using the ESAR Data Set

Next, the performance of our proposed method is assessed using dual-baseline full-polarimetric L-band real data acquired by the DLR ESAR sensor over the city of Dresden in Germany, on 1 August 2000. The image has a spatial resolution equal to $3\times 2.2$ m in azimuth and range directions, respectively (see

Table 1. Characteristics of the data set used.

	ESAR	UAVSAR
Number of acquisitions	3	9
Acquisition dates	1 August 2000	6 March 2016
Study area	Dresden, Germany	L’Amitie stadium in Mondah, Gabon
Coordinates of the center of SAR images	(51°02′25.21″ N, 13°46′52.53″ E)	(0°31′19.72″ N, 9°23′35.15″ E)
Mid-range incidence angle	${40}^{\circ }$	${45}^{\circ }$
Range resolution	2.2 m	1.67 m
Azimuth resolution	3 m	0.60 m
Wavelength	23 cm	23 cm
Polarization channels	HH, HV, VV	HH, HV, VV

for further information about the polarimetric data set). A subset of data with a size of $1600\times 455$ pixels was selected, with the corresponding Pauli image of the polarimetric master data shown in Figure 2a. To give the reader more information about the study area, the location and the optical image of the test area are also shown in Figure 2b. Note that the optical image is only presented as additional information and is not used at all in the experiments. The overall temporal baseline is less than one hour, and the spatial baseline configuration is non-uniform, spanning from 10 m to about 40 m with a mid-range incidence angle and a range distance of ${40}^{\circ }$ and 4486 m, respectively, resulting in a Fourier height resolution of 15 m in the mid-range. The dual-baseline polarimetric sample covariance matrix $\widehat{\mathbf{R}}$ of the data set was estimated by $L=25$ looks, obtained by a $5\times 5$ moving average window, and the reconstruction processes were performed in 3D space.

Figure 3a,b show the maps of the detected single scatterers using the existing scattering support GLRT (10) and the proposed method (16), respectively. Moreover, the proposed detector (17) and scattering support GLRT (11) tailored to detect the double scatterers were also applied to the dual-baseline ESAR images. The elevation maps of the detected double scatterers are shown in Figure 3c,d.

From the detected single scatterers shown in the first row of Figure 3, and as expected from the simulated experiments, the proposed method is superior to the scattering support GLRT in terms of detection performance. As can be noticed in Figure 3a,b, the number of detected single scatterers is increased from 265,231 with the scattering support GLRT to 339,840 with the proposed method. This corresponds to an increase in the average density of detected single scatterers from 2379 to 3056 points per square kilometer with the existing and proposed methods, respectively. From the maps, most of the scatterers detected by both approaches belong to the non-vegetation areas. This is due to the fact that a hard threshold was set for both approaches, which was achieved by fixing a very low probability of false alarms.

The second row of the figure shows the double scatterers detected by both approaches. As can be seen from the results of both methods, particularly in the upper right part of the images, mainly the ground scatterers are detected as double scatterers (Figure 3c,d), while the roofs of the residential buildings are detected as single scatterers (see Figure 3a,b). In agreement with the results of the numerical experiments, the comparison of the detected double scatterers given in Figure 3c,d confirms the superior performance of the proposed approach over the existing scattering support GLRT. In particular, the proposed method detects a larger number of double scatterers, which are mainly found on vertically extending structures, while the existing scattering support GLRT failed to detect them. The number of detected double scatterers is 41,556 with the scattering support GLRT, and it increases to 45,294 with the proposed method.

In addition, two representative pixels are selected for detailed analysis with locations shown as green circles in Figure 2. These cells were chosen to observe how Equation (18) finds the parameters of optimal polarization for detecting PS. Figure 4 shows a three-dimensional plot of ${\mathrm{\Lambda }}_{\rho }={\lambda }_{max}({\mathbf{B}}^{†}(\widehat{\Theta })\mathbf{W}\widehat{\mathbf{R}}{\mathbf{W}}^{†}\mathbf{B}(\widehat{\Theta }))$ over the geometrical parameters $\chi$ and $\tau$. The optimal polarization basis shows up with highest power in the plots. The plots show the estimated values of ${\mathrm{\Lambda }}_{\rho }$ on all polarization bases derived by changing the orientation and the ellipticity angles in their interval ranges. From the 3D power plot of pixel A in Figure 4a, it can be seen that the most optimal polarization basis is given when $\chi ={24}^{\circ }$ and $\tau =-4^{\circ }$. For pixel B, on the other hand, and according to Figure 4b, the best polarization basis is represented by $\chi ={179}^{\circ }$ and $\tau =-{21}^{\circ }$. In general, the shape of the 3D plots and the amount of the power ${\mathrm{\Lambda }}_{\rho }$ on different bases depend on the characteristics of the sensed target. It is important, however, to show that optimization of the polarization basis with the proposed method can identify the presence of additional weaker scatterers that cannot be detected with the conventional polarization basis ($\chi =0^{\circ }$ and $\tau =0^{\circ }$). This strategy obviously increases the number of detectable scatterers, which leads to an increase in the density of PSs, as confirmed by the results in Figure 3.

To evaluate the results represented in Figure 3 and to identify the optimal polarization bases for detecting permanent scatterers in the tomographic frame, the histograms of the polarizations characterized in (18) are plotted with respect to the geometrical parameters, orientation, and ellipticity angles, as seen in Figure 5. The histograms provide detailed insight into the discovered polarization bases, which increase the probability of PS detection. Considering the mean value of the geometric parameters in the histograms, the polarization bases with parameters of ($\chi ={20}^{\circ }$ and $\tau =0^{\circ }$) and ($\chi ={110}^{\circ }$ and $\tau =0^{\circ }$) are recognized as the optimal polarization bases for detecting the permanent scatterers in the study area.

5.2.2. Experiment Using the UAVSAR Data Set

Additionally, the experiments are extended to a UAVSAR data set over l’Amitie stadium in Mondah, Gabon. The data set was acquired by the National Aeronautics and Space Administration (NASA), during the mission of the AfriSAR, on 6 March 2016. The data include a stack of 9 polarimetric L-band SAR images with a spatial resolution equal to $0.6\times 1.67$ m in azimuth and slant range directions, respectively (See Table 1 for further information). In contrast to the ESAR data set, the spatial baseline distribution is uniform, and the flight lines were in a vertical plane, with fixed offsets of 20 m with respect to the reference flight line. The incidence angle changes from 25 degrees at near range to 65 degrees at far range. The subset of the data with a size of $970\times 620$ pixels covering l’Amitie stadium was selected, and the corresponding Pauli image of the master data is shown in Figure 6a. For better visual inspection of the study area, the location of the test site and the corresponding optical image are also shown in Figure 6b.

The elevation maps of the detected single and double scatterers using the existing scattering-based support GLRT and the proposed method on the UAVSAR data set are shown in Figure 7. The comparison of the detected single scatterers shows that the proposed method provides denser and more valid PS points in analogy with the simulated and real ESAR data experiments. In particular, 115,850 points are detected as single scatterers with the scattering-based support GLRT in Figure 7a, while the proposed method was able to detect 41,340 additional scatterers and increase the number of single scatterers to 157,190 points in Figure 7b. The second line of Figure 7 belongs to the double scatterers detected by both methods. As can be seen from the results, the detected double scatterers mainly belong to buildings and the stadium. The results show that the double scatterers are due to either the ground–wall interaction of the SAR signal or the direct backscattering from the upper layers of the buildings and roofs. From the comparison of the results, it is also seen that the proposed method is able to increase the number of double scatterers from 5180 points with the scattering-based support GLRT (Figure 7c) to 6987 PS points (Figure 7d). The comparison using the UAVSAR data set, thus, reconfirms the performance of the proposed method over the existing scattering-based support GLRT. In general, it can be observed that the proposed method detects a larger number of additional scatterers that the existing method could not identify.

Further evaluation and assessment of the results yielded by the proposed method are reported in Figure 8. The figure shows the histograms of the identified optimal polarization bases using Equation (18) in the detection process. In analogy with the experiments from the ESAR data set (Figure 5), the discovered optimal polarization bases in Figure 8, when plotted as a function of the geometrical parameters (orientation and ellipticity angles), show that the polarization bases with parameters of ($\chi ={40}^{\circ }$ and $\tau =0^{\circ }$) and ($\chi ={140}^{\circ }$ and $\tau =0^{\circ }$) are mainly identified as the optimal polarization bases for the detection of the permanent scatterers over the selected test site in Mondah, Gabon.

6. Discussion

The results and analyses presented in the previous section are discussed here. Based upon the numerical experiments, the proposed method increases the probability of detection of both single and double scatterers compared to the existing scattering support GLRT under the same conditions. The improvement in detection at both stages of the test is mainly related to the optimization of the polarization bases by the proposed method. The detection tests, in general, depend indirectly on the polarization basis of the data stack via the covariance matrix. Our method uses the best combination of transmit and receive polarizations to achieve the best detection performance at a constant false alarm rate.

In analogy with the simulated experiments, the results with real data presented in Figure 3 confirm the performance of the proposed method over the existing scattering support GLRT. The validity and reliability of the GLRT-based techniques, especially the support GLRT and the scattering support GLRT, have been extensively evaluated and confirmed in the literature [7,13,14,15,16,17,18,19]. According to the results, and in the context of the multi-look detection process with polarimetric data, it could be more efficient to design the detection test based upon polarization basis optimization rather than relying on the use of the conventional horizontal and vertical transmit and receive polarizations in the radar data. The efficiency of such an optimization is that it provides a higher density of detected PSs. The proposed method provides the possibility of revealing the presence of weaker scatterers by using the optimal polarization basis; these scatterers could not be detected using the detection tests relying only on the use of the conventional horizontal and vertical polarization basis. This is in line with the point that detection performance depends upon the polarization basis used, as shown in Figure 4 for two sample resolution cells. Typically, the best polarization basis for the detection process will vary from pixel to pixel, and the optimal basis depends mainly upon the sensed target. The statistics of the optimal polarization over different detected PSs, however, can provide important information about the effects of the polarization design of a radar antenna on the probability of scatterer identification. The analysis of the reported histograms using the ESAR data set in Figure 5 provides an important indication: the linear polarization basis ($\tau \approx 0$) with orientation angles of $\chi \approx {20}^{\circ }$ and $\chi \approx {110}^{\circ }$ is the optimal polarization for detecting permanent scatterers in the Dresden test area. Similarly, the histograms obtained from the UAVSAR data set in Figure 8 confirm the capability of the linear polarization basis ($\tau \approx 0$) in the PS detection process. At the Mondah test site, however, linear polarization with orientation angles of $\chi \approx {40}^{\circ }$ and $\chi \approx {140}^{\circ }$ is identified as the most optimal polarization basis. Recall that the change in optimal polarization from pixel to pixel and from study area to study area is due to the fact that the optimal polarization basis is a function of the backscattering mechanism of the sensed targets. Therefore, each study area may have a different optimal polarization, depending on the targets presented.

Generally, the use of such polarization bases in both test sites made it possible to detect the presence of additional weaker scatterers that were not detectable with the conventional horizontal-vertical polarization basis.

The efficiency of the implemented methods for processing our real data stack in Figure 3 is given in terms of processing time. The experiments were performed in the MATLAB environment on a desktop system with a configuration of Intel(R) Core(TM) i7-4510U CPU @ 2.00 GHz and 12 GB RAM. The total runtime for the existing scattering support GLRT is about $7.6$ h, while the runtime for the proposed method increases to $16.4$ h. Although this increase in runtime may limit the application of the proposed framework in the tomographic processing of large data sets, some tips can be considered to reduce the processing time. For example, the geometrical parameters ($\chi$ and $\tau$) can be scanned at widely discretized intervals. In addition, the use of Matlab Mex files or parallel computing and Graphics Progressing Units (GPUs) can also significantly reduce running times. Although online processing is an interesting strategy to extend the application of the proposed method, the requirement of the method for large amounts of data to be uploaded as input may limit online processing tools.

Although not examined in this study, the proposed method, like other existing tomographic GLRT-based methods, is expected to be applicable to detecting permanent scatterers in built-up regions with high-rise structures and also for mapping elevations and displacements of the earth’s surface in non-urban regions by detecting PSs from coherent points such as rocks, etc. The performance of the method in such scenarios has yet to be evaluated. In addition, the detection of PSs from single and double scatterers can be extended to additional scatterers (e.g., triple and quadruple scatterers). This calls for the extension of testing to multiple stages. Testing and evaluating the effectiveness of the proposed method with multiple stages in a study area with high-rise skyscrapers will be a future line of research after this study.

7. Conclusions

This paper presented a new method for detecting permanent scatterers using tomographic polarimetric data sets. The proposed method extends the existing scattering support GLRT method and takes advantage of multi-look and polarimetric data to enhance the performance of the detection scheme. The main insight of the paper lies in a polarization synthesis theorem that allows the full characterization of scatterers in any desired polarization basis. The proposed method optimizes the detection performance with respect to the polarization basis. Experiments on one simulated and two different real data sets showed that the proposed method outperforms the existing scattering support GLRT method for single and double scatterer detection. The results showed that linear polarizations with non-zero orientation angles commonly provide the possibility of detecting additional weaker scatterers compared to the use of the conventional horizontal-vertical polarization basis. Improvement in detection performance is achieved at the cost of higher computational processing requirements. The existing GLRT-based tests are based upon optimization with respect to the target parameters ($\Theta$), while the proposed method performs optimization with respect to both the target parameters and the polarization basis, resulting in a higher computational cost. However, strategies such as those listed in the previous section can improve the proposed method in terms of computational cost.

In conclusion, the outcomes of the proposed PS detection method are particularly useful for optimal data selection for PS detection-based applications and for the design of new spaceborne/airborne missions based on repeat-pass processing. Municipalities and city managers, in particular, will benefit from the high density PS extracted by this study, either to update their digital surface models or to accurately model and inspect the urban environment at a fine level of detail so that they can better predict risk situations when the detection process is extended to 4D and 5D cases.