A MIMO-SAR Tomography Algorithm Based on Fully-Polarimetric Data

Abstract

A fully-polarimetric unitary multiple signal classification (UMUSIC) tomography algorithm is proposed, which can be used for acquiring high-resolution three-dimensional (3D) imagery, in a polarimetric multiple-input multiple-output synthetic aperture radar (MIMO-SAR) with a small number of baselines. In terms of the elevation resolution, UMUSIC provides an improvement over standard MUSIC by utilizing the conjugate of the complex sample data and converting the complex covariance matrix into a real matrix. The combination of UMUSIC and fully-polarimetric data permits a further reduction of the noise of the sample covariance matrix, which is obtained through pixel averaging of multiple two-dimensional (2D) images. Considering the consistency of four polarizations, this algorithm not only makes scattering centers have the same estimated height in four polarizations, but it also improves the estimation accuracy. Simulation results show that this algorithm outperforms the popular distributed compressed sensing (DCS). Image processing of measured data of an aircraft model using a multiple-input multiple-output synthetic aperture radar (MIMO-SAR) with six baselines is presented to validate the proposed algorithm.

1. Introduction

Multiple-input multiple-output synthetic aperture radar (MIMO-SAR) is an enabling technique capable of imaging a target [1,2,3,4,5,6,7], which is different from a rail synthetic aperture radar (SAR) and a turntable inverse synthetic aperture radar (ISAR). Two-dimensional (2D) virtual apertures can be synthesized through different combinations of transceiver antenna elements. A large number of virtual apertures in the cross-range and elevation directions are beneficial to obtain high-resolution three-dimensional (3D) radar images [4,5,6], but they also result in a high cost and large size of radar systems due to the increase of the number of antenna elements. When the measured target has several scattering centers in an elevation direction, such as airplanes, an affordable array strategy can be adopted: the priority is given to ensuring adequate cross-range virtual apertures for high-resolution two-dimensional (2D) radar images [7], and then a small number of elevation virtual apertures ensure high-resolution 3D radar images by SAR tomography [8,9,10,11,12,13,14,15,16,17,18].

The reconstruction quality of SAR tomography depends on the product of the number of baselines and the signal-to-noise ratio (SNR) [8]. When the number of baselines is limited, the SNR of radar images can be equivalently improved by filtering [8,9], auxiliary information [10,11,12], and polarization [13,14,15,16,17]. The SNR can be improved by integrating nonlocal filtering into the compressed sensing (CS) algorithm and a reasonable reconstruction of buildings from only seven baselines is feasible [8]. In addition, [9] investigates the possibility related to the use of a multi-looking approach for fine resolution analysis of ground structures that combines SAR tomography. Filtering consisting of averaging pixels is bound to reduce the range-azimuth resolution and therefore is not suitable for high-resolution radar images of artificial targets. The auxiliary information added to the standard Capon and multiple signal classification (MUSIC) algorithms can be exploited to reduce the ambiguity and resolve the superimposition of the scatterers in the case of a limited number of radar images [10]. Atmospheric phases for SAR tomography in mountainous regions are regressed against the spatial coordinates in map geometry at persistent scatterers locations [11]. The high-resolution 3D positions of a large amount of natural scatterers are obtained by a geodetic SAR tomography framework that fuses SAR tomography and SAR image geodesy compensating SAR measurement error [12]. The auxiliary information can effectively improve the image quality, but it needs to be obtained by other technologies which increases the complexity of the algorithm and the cost of the imaging. A distributed compressed sensing (DCS) algorithm based on fully-polarimetric data is proposed in [13,14,15,16] to improve the accuracy of the estimation. However, the CS algorithm suffers from a high computational expense and is hard to extend to fast practice [18]. In [17], a comparison among tomograms obtained in different polarizations is made to analyze how polarimetry can enhance target signatures.

To address these problems, the combination of spectral analysis and full polarization is an attractive way to improve resolution and processing speed for a small number of baselines. This paper explores a fully-polarimetric unitary multiple signal classification (UMUSIC) technique for polarimetric MIMO-SAR tomography [7]. The remainder of the paper is organized as follows. Section 2 describes a signal model based on fully-polarimetric data. A fast and high-resolution UMUSIC algorithm is developed in Section 3. In Section 4, two algorithms are compared through simulation of different point scatterers. Finally, Section 5 contains measured tomography results of an aircraft model to validate the proposed algorithm.

2. Polarization Signal Model

SAR tomography allows us to obtain 3D imagery to describe the electromagnetic property of illuminated objects. The geometry of MIMO-SAR tomography is shown in Figure 1, where x,y,z denote coordinates originating from the center of the imagery scene, R_n represents the distance from the target to the nth baseline, and R₀ is the projection distance from the radar to the center of the imagery scene on the y-axis. The orange triangles and blue circles denote receivers and transmitters, respectively. Each baseline represents a linear array, where transmitters are at both ends of the array and receivers are in the middle of the array. The 2D image for the nth baseline is represented by the following form [19].

$$g_n(x,y,z_n)={\displaystyle {\int }_{-h/2}^{h/2}\sqrt{\stackrel{⌢}{\sigma }(x,y,z)}{e}^{-j\frac{4\pi }{\lambda }{R}_n}}dz$$

(1)

where $\sqrt{\stackrel{⌢}{\sigma }(x,y,z)}$ denotes the target scattering function that needs to be solved, h is the height of the imagery scene, λ represents the wavelength and z_n refers to the z-coordinates of the nth baseline. Under the Born weak-scattering approximation, $R_n$ representing the distance from the point target at (x, y, z) to the nth baseline, is approximated as:

$$R_n=\sqrt{{(y+R_0)}^2+{(z-z_n)}^2}\approx y+R_0+\frac{z_n^2+z^2-2z_{n}z}{2R_0}$$

(2)

The first three terms in (2) are irrelevant to z, the fourth term is the residual phase term that can be merged into $\sqrt{\stackrel{⌢}{\sigma }(x,y,z)}$, and the fifth term is the phase term for imaging in the elevation direction. We can choose one of N baselines as a reference baseline as:

$$g_{ref}(x,y,z_{ref}){=\mathrm{e}}^{-j\frac{4\pi }{\lambda }(y+R_0+\frac{z_n^2}{2R_0})}$$

(3)

After phase compensation, the N 2D images $g_n^{\prime }(x,y,z_n)$ and ${\stackrel{⌢}{\sigma }}^{\prime }(x,y,z)$ are in a Fourier transform relation.

$$g_n^{\prime }(x,y,z_n)=g_n(x,y,z_n)g_{ref}^{*}(x,y,z_{ref})={\displaystyle {\int }_{-h/2}^{h/2}{\stackrel{⌢}{\sigma }}^{\prime }(x,y,z)e^{-j2\pi w_{n}z}dz}$$

(4)

with

$$w_n=\frac{2z_n}{\lambda R_0}$$

(5)

$${\stackrel{⌢}{\sigma }}^{\prime }(x,y,z)=\sqrt{\stackrel{⌢}{\sigma }(x,y,z)}{e}^{-j\frac{2\pi z^2}{\lambda R_0}}$$

(6)

where ${\stackrel{⌢}{\sigma }}^{\prime }(x,y,z)$ integrated with the residual phase term still has the same amplitude as $\sqrt{\stackrel{⌢}{\sigma }(x,y,z)}$, which has no effect on the 3D imagery. $g_n^{\prime }(x,y,z_n)$ can be discretized as the multiplication of two matrices.

$$g_n^{\prime }(x,y,z_n)={\displaystyle \sum _{l=1}^L{\stackrel{⌢}{\sigma }}^{\prime }(x,y,z_l){\mathrm{e}}^{-j2\pi w_n{z}_l}}=a_{n}s$$

(7)

with

$$a_n=[{\mathrm{e}}^{-j2\pi w_n{z}_1},{\mathrm{e}}^{-j2\pi w_n{z}_2},\dots ,{\mathrm{e}}^{-j2\pi w_n{z}_L}]$$

(8)

$$s={[{\stackrel{⌢}{\sigma }}^{\prime }(x,y,z_1),{\stackrel{⌢}{\sigma }}^{\prime }(x,y,z_2),\dots ,{\stackrel{⌢}{\sigma }}^{\prime }(x,y,z_L)]}^T$$

(9)

where L represents the number of scattering centers in the elevation direction. Before imaging, the imagery scene in the elevation direction needs to be divided into many discrete points to represent the range of L. As a consequence, the matrix $a_n$ is very large, which is the reason why the imaging algorithm takes a long time to locate scattering centers and determine L in the simulation and measurement. In combination with N 2D images, the polarimetric tomography model is given by:

$$g=\left[\begin{array}{c}{g}_1^{\prime }(x,y,z_1)\\ g_2^{\prime }(x,y,z_2)\\ ⋮\\ g_N^{\prime }(x,y,z_N)\end{array}\right]=\left[\begin{array}{cccc}{\mathrm{e}}^{-j2\pi w_1{z}_1}& {\mathrm{e}}^{-j2\pi w_1{z}_2}& \cdots & {\mathrm{e}}^{-j2\pi w_1{z}_L}\\ {\mathrm{e}}^{-j2\pi w_2{z}_1}& {\mathrm{e}}^{-j2\pi w_2{z}_2}& \cdots & {\mathrm{e}}^{-j2\pi w_2{z}_L}\\ ⋮& ⋮& \ddots & ⋮\\ {\mathrm{e}}^{-j2\pi w_N{z}_1}& {\mathrm{e}}^{-j2\pi w_N{z}_2}& \cdots & {\mathrm{e}}^{-j2\pi w_N{z}_L}\end{array}\right]s=As$$

(10)

(10) can be further developed by merging with the fully-polarimetric data.

$$G=[\begin{array}{cccc}{g}_{HH}& g_{HV}& g_{VH}& g_{VV}\end{array}]={[a_1,a_2,\dots ,a_N]}^T[\begin{array}{cccc}{s}_{HH}& s_{HV}& s_{VH}& s_{VV}\end{array}]=AS$$

(11)

where G denotes the N×4 2D imagery matrix, S refers to the L×4 3D imagery matrix, and A is the N × L transformation matrix.

3. Tomography Algorithm

Multiple signal classification (MUSIC) is a spectral analysis algorithm based on the eigen decomposition of the sample covariance matrix. It is necessary to reduce the matrix noise by some techniques, including snapshot in direction-of-arrival and multi-looking in SAR tomography that also leads to a decrease in range-azimuth resolution [13]. In order to dispel the influence, we employ fully-polarimetric data and their conjugation to obtain the matrix.

The polarization tomography model (11) including Gaussian white noise matrix $W$ is rewritten as:

$$G=AS+W$$

(12)

The covariance matrix $R$ is given by:

$$\begin{array}{l}R=E\left\{G{G}^H\right\}=E\left\{\left[AS+W\right]\left[S^H{A}^H+W^H\right]\right\}\\ =AE\left\{S{S}^H\right\}{A}^H+E\left\{W{W}^H\right\}=AP{A}^H+{\sigma }^{2}I\end{array}$$

(13)

where ${\sigma }^{2}I=diag\left\{{\left|{\sigma }_1\right|}^2\cdots {\left|{\sigma }_N\right|}^2\right\}$. The Hermite matrix $AP{A}^H$ composed of the positive definite diagonal matrix $P$ and the column full rank matrix $A$ can be eigen decomposed to obtain the noise subspace.

Multi-sample data is critical for the MUSIC algorithm to ensure a high-quality covariance matrix. In SAR tomography, the average of pixels with the same scattering characteristics is called multi-looking. However, the number of these pixels is limited, and the average processing also reduces the range-azimuth resolution. The combination of observed data and their conjugation can be equivalent to doubling the number of these pixels, which not only improves the estimation accuracy but also solves the problem of coherent signal estimation. To force the Hermite property of $AP{A}^H$, the average of the covariance matrices is computed from forward and backward data samples [20].

$$R_M=\frac{1}{2}(R+J_N{R}^{*}{J}_N)=A\tilde{P}{A}^H+{\sigma }^{2}I$$

(14)

with

$$\tilde{P}=\frac{1}{2}(P+D{P}^{*}{D}^H)$$

(15)

$$D=diag\left\{{\mathrm{e}}^{j2\pi (w_N+w_1)z_1}\dots {\mathrm{e}}^{j2\pi (w_N+w_1)z_L}\right\}$$

(16)

where superscript ^* represents the conjugation, $J_N$ denotes the N×N exchange matrix with ones on its antidiagonal and zeros elsewhere. In order to reduce the computational complexity, $R_M$ can be transformed into a real covariance matrix by unitary transformation.

$$R_U=Q_N^H{R}_M{Q}_N^{}=Q_N^{H}A\tilde{P}{A}^H{Q}_N^{}+{\sigma }^{2}I$$

(17)

where $Q_N^{}$ is any N×N unitary matrix to satisfy column conjugate symmetry. A simple form can be chosen as [21]

$$Q_{N:\mathrm{even}}=\frac{1}{\sqrt{2}}\left[\begin{array}{cc}{I}_{N/2}& \mathrm{j}{I}_{N/2}\\ J_{N/2}& -\mathrm{j}{J}_{N/2}\end{array}\right]$$

(18)

$$Q_{N:odd}=\frac{1}{\sqrt{2}}\left[\begin{array}{ccc}{I}_{(N-1)/2}& 0& \mathrm{j}{I}_{(N-1)/2}\\ 0& \sqrt{2}& 0\\ J_{(N-1)/2}& 0& -\mathrm{j}{J}_{(N-1)/2}\end{array}\right]$$

(19)

As mentioned in (14) and (17), the real covariance matrix $R_U$ can be rewritten as:

$$\begin{array}{l}{R}_U=\frac{1}{2}(Q_N^{H}R{Q}_N^{}+Q_N^H{J}_N{R}^{*}{J}_N{Q}_N^{})\\ =\frac{1}{2}(Q_N^{H}R{Q}_N^{}+(Q_N^H{)}^{*}{R}^{*}{J}_N{Q}_N^{*})=\mathrm{Re}\{Q_N^{H}R{Q}_N^{}\}\end{array}$$

(20)

The real matrix is eigen decomposed as:

$$\mathrm{Re}\{Q_N^{H}R{Q}_N^{}\}={\displaystyle \sum _{i=1}^N{\lambda }_i{u}_i{u}_i^H}+{\sigma }^2{\displaystyle \sum _{i=1}^N{u}_i{u}_i^H}={\displaystyle \sum _{i=1}^L{\lambda }_i{u}_i{u}_i^H}+{\sigma }^2{\displaystyle \sum _{i=1}^N{u}_i{u}_i^H}$$

(21)

where ${\lambda }_1$,…,${\lambda }_N$ are eigenvalues and $u_1$,…,$u_N$ represent corresponding orthogonal normalized eigenvectors. Among N eigenvectors of the $R_U$, L eigenvalues are related to the signal, and N-L eigenvalues are related to the noise. By using the noise subspace $E_N=span\{u_{N-L},\dots ,u_N\}$, the fully-polarimetric pseudo-spectrum is expressed as:

$$P_{MUSIC}^{FP}(w)=\frac{1}{A^H(w)E_N{E}_N^{H}A(w)}$$

(22)

We can find out the peaks of the spectrums to locate different scattering centers in the elevation dimension and estimate the scattering intensity by the least square method (LSM) in four polarizations.

$$S={(A^{H}A)}^{-1}{A}^{H}G$$

(23)

where $A$ is updated according to positions of the scattering centers.

4. Simulation

We adopt the fully-polarimetric DCS as the comparison item. Considering crosstalk, noise, and dispersion, point scatterers with typical polarimetric scattering matrices (PSMs) are simulated to generate return signals. The simulation scene is shown in Figure 2, where scatterers with different heights are located in the coordinate origin. The working frequency is 8 GHz–12 GHz. The down-range, cross-range and elevation Rayleigh limits of the MIMO radar are 0.037 m, 0.047 m, and 0.188 m, respectively. As shown in Figure 3, we simulate three cases to compare the two algorithms.

4.1. Case 1: Two Point Scatterers with a Spacing of 0.18 m

A cylinder and a 90° rotated dihedral reflector are located at −0.09 m and 0.09 m in the elevation dimension, respectively. When the two scatterers spacing is 0.18 m (close to the Rayleigh limit), the simulation results of the two algorithms are seen in Figure 4, where lines represent pseudo-spectrums and points denote estimated results including height and scattering intensity of scatterers.

Two fully-polarimetric algorithms make scatterers have the same estimated height in four polarizations. The specific estimated results are listed in

Table 1. Estimated results of two point scatterers with a spacing of 0.18 m.

			Cylinder (–0.09 m)	90° Rotated Dihedral Reflector (0.09 m)
Fully-polarimetric DCS	Height		−0.091 m	0.091 m
	Scattering intensity	HH	−1.8 dB	0 dB
		HV	−18.5 dB	−39.6 dB
		VH	−18.5 dB	−36.5 dB
		VV	−1.8 dB	0 dB
Fully-polarimetric UMUSIC	Height		−0.090 m	0.090 m
	Scattering intensity	HH	−1.8 dB	0 dB
		HV	−18.5 dB	−39.1 dB
		VH	−18.5 dB	−36.2 dB
		VV	−1.8 dB	0 dB

. The PSMs of the two scatterers are estimated, where the scattering intensity of the cylinder return signal at −0.09 m is inconsistent with the truth because of polarimetric distortion. In Figure 4a, it is noteworthy that the pseudo-spectrum of CS is leaked to form false scattering points. There are two main reasons for signal leakage [13]: on the one hand, if the regularization parameter is too small in the optimization model, it can lead to over-fitting of data; on the other hand, the observed data does not satisfy the sparsity in the unit orthogonal basis. Therefore, it is necessary to use a sliding window to suppress signal leakage.

4.2. Case 2: Two Point Scatterers with a Spacing of 0.06 m

When the spacing is reduced to 0.06 m (one-third of elevation Rayleigh limit), the estimated results of the two algorithms are displayed in Figure 5 and

Table 2. Estimated results of two point scatterers with a spacing of 0.06 m.

			Cylinder (−0.06 m)	90° Rotated Dihedral Reflector (0 m)
Fully-polarimetric DCS	Height		−0.062 m	0.004 m
	Scattering intensity	HH	−3.0 dB	−1.6 dB
		HV	−17.6 dB	−29.5 dB
		VH	−17.6 dB	−37.5 dB
		VV	−3.2 dB	0 dB
Fully-polarimetric UMUSIC	Height		−0.060 m	0.004 m
	Scattering intensity	HH	−2.9 dB	−0.7 dB
		HV	−17.6 dB	−28.7 dB
		VH	−17.6 dB	−36.4 dB
		VV	−3.1 dB	0 dB

. Two fully-polarimetric algorithms still have high-resolution. Furthermore, polarimetric distortion of the cylinder return signals becomes more severe as the spacing decreases. Consequently, polarimetric calibration is necessary for a fully-polarimetric radar system.

4.3. Case 3: Four Point Scatterers with a Spacing of 0.09 m

The four point scatterers are a cylinder, a 67.5° rotated dihedral reflector, a 90° rotated dihedral reflector, and a plate and their PSMs are listed in

Table 3. Polarimetric scattering matrix (PSM) of four point scatterers

Parameters	Cylinder	67.5° Rotated Dihedral Reflector	90° Rotated Dihedral Reflector	Plate
HH	−1	$\sqrt{2}/2$	1	−1
HV	0	$\sqrt{2}/2$	0	0
VH	0	$\sqrt{2}/2$	0	0
VV	−1	−$\sqrt{2}/2$	−1	−1

. According to the MIMO configuration, the bistatic angles of all transceiver channels are less than 10°. To simplify the simulation, we assume that the PSMs listed in Table 3 are applicable to all transceiver channels. It can be seen from Figure 6 that the pseudo-spectrums of two fully-polarimetric algorithms are not affected when the number of scatterers increase. We summarize the estimation results in

Table 4. Estimated results of two point scatterers with a spacing of 0.09 m.

			Cylinder (−0.13 m)	67.5° Rotated Dihedral Reflector (−0.04 m)	90° Rotated Dihedral Reflector (0.05 m)	Plate (0.14 m)
Fully-polarimetric DCS	Height		−0.141 m	−0.030 m	0.041 m	0.141 m
	Scattering intensity	HH	−2.4 dB	−2.5 dB	−2.4 dB	−2.2 dB
		HV	−16.2 dB	−1.6 dB	−24.5 dB	−15.1 dB
		VH	−16.2 dB	−1.6 dB	−24.1 dB	−15.1 dB
		VV	−1.2 dB	−2.4 dB	−2.5 dB	0 dB
Fully-polarimetric UMUSIC	Height		−0.131 m	−0.039 m	0.057 m	0.140 m
	Scattering intensity	HH	−0.3 dB	0 dB	−1.7 dB	−0.5 dB
		HV	−14.6 dB	−2.4 dB	−36.8 dB	−15.1 dB
		VH	−14.6 dB	−2.4 dB	−37.2 dB	−15.1 dB
		VV	−0.3 dB	−4.0 dB	−3.0 dB	0 dB

, which demonstrates the estimation accuracy of fully-polarimetric UMUSIC is higher than that of the fully-polarimetric DCS. The CS, which is essentially an optimization problem, needs to be solved iteratively, therefore, its processing speed is bound to be limited by the number of iterations. The simulation results show that for a pixel, the processing speed of the fully-polarimetric UMUSIC is more than five times faster than that of the fully-polarimetric DCS in the same computing condition.

5. Experiment

An experimental polarimetric MIMO array has been upgraded based on the radar system in [7], and baselines with different heights are controlled by an elevator. It can be seen from Figure 7 that the polarimetric MIMO array consists of 20 receive elements and 6 transmit elements, where the combinations among them synthesize 80 transceiver channels. The measured target is an aircraft model with an elevation angle of 16 degrees on a foam support, as shown in Figure 8. M1, M2, and M3 represent three missile models mounted on the wing, respectively. To avoid complex scattering properties of cavity structures, the inlet of the aircraft model is sealed with copper foils. The measurement parameters are the same as the simulation parameters in Section 4.

It can be seen from Figure 9 that the scattering mechanisms of the aircraft model are different in four polarizations. In the HH image, there are three strong scattering centers including two parts that are not distinguished (see Figure 9a). Compared with the HH image, more components can be distinguished from the VV image. The scattering intensity of the two cross-polarization images is low. Figure 10, Figure 11, Figure 12 and Figure 13 illustrate the 3D point cloud maps obtained from 24 2D images. The top views are similar to the 2D image, which proves that the 3D scattering intensity can be estimated by LSM. It can be seen from the bottom and side views that scattering centers with different heights are basically consistent with the aircraft model. In addition, the scattering intensity in front of the fuselage is higher than that of the fuselage tail due to the shielding of the supporting foam. By comparing the 3D point cloud maps in different polarizations, we can analyze its scattering mechanism.

Figure 14 illustrates tomographic image slices along the down range for HH (Figure 10), HV (Figure 11), VH (Figure 12), and VV (Figure 13). It can be seen from the figures that scattering of the model shows obvious variety with heights. We summarize components of the model in

Table 5. Components of the aircraft in tomographic image slices.

Down Range	−0.4 m	−0.3 m	−0.2 m	−0.1 m	0 m	0.1 m	0.2 m
Components	Aircraft head	Front wheel	Inlet	M1 head	M1 tail+ M2 head	M2 tail+ Rear wheel	M3 tail
HH	x	x	x	x	x	x	x
HV				x	x	x
VH			x	x	x	x
VV	x	x	x	x	x	x	x

, where M1 tail and M2 head cannot be distinguished because they have the same height, so do M2 tail and rear wheel.

6. Conclusions

This paper proposes a fully-polarimetric UMUSIC tomography algorithm to acquire high-resolution 3D radar imagery for a MIMO-SAR with a small number of baselines. In order to mitigate the effect of multi-looking on the range-azimuth resolution, we employ fully-polarimetric data and their conjugation to obtain the sample covariance matrix. Two algorithms including the fully-polarimetric DCS and the fully-polarimetric UMUSIC, are compared through numeric simulation of different point scatterers. Simulation results demonstrate that the fully-polarimetric UMUSIC outperforms the popular fully-polarimetric DCS in processing speed and estimation accuracy. Measurements for an aircraft model are conducted using an X-band experimental polarimetric MIMO-SAR which was upgraded from a previous system [7]. The resulting 3D images using six baselines demonstrate the usefulness of the algorithm for 3D imagery of complex radar targets.