A Fast IAA–Based SR–STAP Method for Airborne Radar

Abstract

Space–time adaptive processing (STAP) is an effective technology in clutter suppression and moving target detection for airborne radar. When working in the heterogeneous environment, the number of training samples that satisfy independent and identically distributed (IID) conditions is insufficient, making it difficult to ensure the estimation accuracy of the clutter plus noise covariance matrix for traditional STAP methods. Sparse recovery–based STAP (SR–STAP) methods have received widespread attention in the past few years. The accurate estimation of the clutter plus noise covariance matrix can be achieved using only a few training samples. The iterative adaptive approach (IAA) can quickly and accurately estimate the power spectrum, but applying this method directly to the STAP method cannot produce good performance. In this paper, a fast IAA–based SR–STAP method is proposed. Based on the weighted $l_1$ problem, the IAA spectrum is used as a weighted term to obtain a good approximation. In order to obtain an analytical solution, we use the weighted $l_2$ norm to approximate the weighted $l_1$ norm without loss of performance. Compared with the IAA–STAP method, the proposed method is more robust to errors. Moreover, the proposed method has a fast computational speed. The effectiveness of the proposed method is demonstrated by simulations.

1. Introduction

Space–time adaptive processing (STAP) [1,2] has been widely used in airborne early warning (AEW) radar since the 1970s and it is a key technology for clutter suppression. After nearly 40 years of development, the application of the STAP method has also expanded to many aeras, such as spaceborne battlefield surveillance radar [3,4] and radar imaging [5,6]. Detecting a weak and slow–moving target that is submerged by ground clutter is a challenging problem, and the STAP method plays an important role in suppressing strong clutter. The core of the STAP method lies in the accurate estimation of the clutter plus noise covariance matrix (CNCM), but accurately estimating the CNCM requires a large number of training samples that satisfy independent and identically distributed (IID) conditions. Unfortunately, it is very difficult to obtain sufficient training samples in a nonhomogeneous environment. Therefore, the estimation accuracy of the CNCM cannot be guaranteed and the performance of the traditional STAP method is severely degraded.

To solve this problem, researchers have made improvements in many ways, and several methods have been proposed. Airborne radar systems typically have a high degree of freedom, and reducing the dimensionality can effectively reduce the demand for training samples. Many reduced–dimension and reduced–rank methods [7,8,9,10] have been proposed. The reduced–dimension methods greatly reduce the requirements for IID training samples while ensuring the clutter suppression performance. However, these methods still need many training samples, especially for systems with a large DOF. The reduced–rank methods have the advantage of small calculation but require the correct selection of the clutter rank. If its value is too large or too small, it seriously affects the performance. Both reduced–dimension methods and reduced–rank methods cannot reduce the dimensionality without limitations and still need to ensure the clutter suppression performance. Thus, the problem of insufficient training samples remains serious. The direct data domain (DDD) method [11,12] further reduces the requirements for training samples. This method only utilizes the cell under test (CUT) itself to achieve good performance in clutter suppression, but aperture loss is also inevitable. The knowledge–aided (KA) method [13,14] has also been proposed. The KA method fully utilizes the available prior information to achieve good performance in clutter suppression, but the performance of these methods is affected by the accuracy of the prior information. When there is a deviation between the prior information and true information, the performance of these methods rapidly decreases. Furthermore, accurate information is difficult to obtain in practice.

Over the past twenty years, researchers have combined the sparse recovery framework with the STAP framework, and many SR–STAP methods have been proposed. Through utilizing the sparsity of clutter in the angle–Doppler domain, these SR–STAP methods can achieve the accurate estimation of the CNCM with only a few training samples. However, the original $l_0$ norm problem is NP–hard; it cannot be solved directly. Such problems are solved by seeking approximations. Researchers have proposed many greedy methods [15,16] to solve such problems, and the orthogonal matching pursuit method (OMP) [16] is one example. This method is widely used because of its easy implementation and high solving efficiency. However, it cannot guarantee that the solution is the global optimal solution, which affects the sparse recovery performance. Another way to solve the $l_0$ norm problem is to replace the $l_0$ norm penalty term with the $l_1$ norm penalty term [17,18]. Although the $l_1$ norm is the tightest convex hull of the $l_0$ norm, the $l_1$ norm penalty term takes into account the magnitude of the non–zero entries of the sparse coefficient, which may result in an “over–penalty” and exaggerate the role of large non–zero elements of the sparse coefficient. Thus, it will reduce the estimation accuracy of the CNCM. The weighted $l_1$ norm method [19] can overcome this problem. There are also some SR–STAP methods based on the sparse Bayesian learning (SBL) framework [20,21,22], which can also achieve good clutter suppression performance, but these methods require a huge amount of computation, which limits the application of these methods. Assigning other sparse priors such as the generalized double Pareto (GDP) prior to the SBL method can reduce the computational time compared with the SBL method [23,24].

In this paper, a fast IAA–based SR–STAP method for airborne radar is proposed, which can achieve a balance between the clutter suppression performance and computing speed. We design a weighted $l_1$ norm penalty term instead of using the $l_1$ norm penalty term directly, making full use of the IAA spectrum and using it as the weighting term. Since there is no analytical solution for the weighted $l_1$ norm, we replace the weighted $l_1$ norm with the weighted $l_2$ norm, without affecting the performance. The main contributions are as follows.

• (1). We combine the IAA spectrum method with the weighted $l_1$ norm method, and a fast IAA–based SR–STAP method is proposed. Compared with the STAP method, which directly uses the IAA method, the proposed method can estimate the CNCM more accurately. Compared with the weighted $l_1$ norm method, the proposed method has an analytical solution.

• (2). The proposed method has fast convergence performance, a shorter running time, and good clutter suppression performance.

• (3). Through a comparison with other STAP methods, simulation results and a performance analysis are given to demonstrate the effectiveness of the proposed method.

The work is organized as follows. Section 2 introduces the data model and gives the details of the proposed method. Section 3 provides the simulation results and performance analyses of the proposed method. Conclusions are drawn in Section 4.

Notations: ${[\cdot ]}^{-1}$, ${[\cdot ]}^T$, and ${[\cdot ]}^H$ represent the matrix inverse, transpose, and conjugate transpose, respectively. ${‖\cdot ‖}_0$, ${‖\cdot ‖}_1$, and ${‖\cdot ‖}_2$ represent the $l_0$ norm, $l_1$ norm, and $l_2$ norm, respectively. The vector is represented in boldface lowercase and the matrix is represented in boldface uppercase. $\odot$ and $\otimes$ denote the Hadamard product and Kronecker product, respectively. $I_N$ is an $N\times N$ identity matrix.

The acronyms employed throughout the main text are shown in Table 1.

2. Signal Model and Problem Formulation

2.1. Signal Model

This paper considers an airborne radar system equipped with a uniform linear array (ULA) of $N$ elements and the element spacing between two elements is set to a half–wavelength. The platform geometry of the airborne radar that is used in this paper is shown in Figure 1.

This radar system transmits $M$ pulses at a constant pulse repetition frequency (PRF) $f_r$. The height of the radar platform is $H$. The azimuth and elevation angles are denoted as $\phi$ and $\theta$, respectively. The velocity of the airborne radar is $v$. By uniformly dividing the distance ring into $Nc$ patches and neglecting the effect of range ambiguity, we obtain the following signal model:

(1)$\begin{array}{ll}\hfill x& =c+n\hfill \\ & ={\displaystyle \sum _{n=1}^{N\mathrm{c}}{a}_{n}s(f_{d,n},f_{s,n})+n}\hfill \\ & ={\displaystyle \sum _{n=1}^{N\mathrm{c}}{a}_n}(s(f_{d,n})\otimes s(f_{s,n}))+n\hfill \end{array}$

where $c$ and $n$ denote the clutter component and white Gaussian noise component, respectively. $a_n$ and $s_n$ are the complex amplitude and the space–time steering vector of the $n-th$ clutter patch, respectively. $s(f_{d,n})$ and $s(f_{s,n})$ are the temporal steering vector and spatial steering vector corresponding to the $n-th$ clutter patch, respectively, which can be expressed as follows:

(2)$\begin{array}{l}s(f_{d,n})={\left[1, \exp (j2\pi f_{d,n}), \cdots , \exp (j2\pi (M-1)f_{d,n})\right]}^T\\ s(f_{s,n})={\left[1, \exp (j2\pi f_{s,n}), \cdots , \exp (j2\pi (N-1)f_{s,n})\right]}^T\end{array}$

where $f_{d,n}$ and $f_{s,n}$ denote the normalized Doppler frequency and normalized spatial frequency of the $n-th$ clutter patch, respectively, which can be expressed as follows:

(3)$\begin{array}{l}{f}_{d,n}=2v\cos {\theta }_n\cos {\phi }_n/(\lambda f_r)\\ f_{s,n}=d\cos {\theta }_n\cos {\phi }_n/\lambda \end{array}$

Assuming that the clutter component and noise component are mutually uncorrelated, the clutter plus noise covariance matrix (CNCM) can be written as

(4)$R=\mathrm{E}[x{x}^{\mathrm{H}}]=R_c+R_n$

where $\mathrm{E}[\cdot ]$ stands for the expectation operation. $R_c$ and $R_n$ denote the clutter covariance matrix and noise covariance matrix, respectively, which have the following form:

(5)$\begin{array}{l}{R}_c=\mathrm{E}[x_c{x}_c^{\mathrm{H}}]={\displaystyle \sum _{n=1}^{N\mathrm{c}}{\left|a_n\right|}^{2}s(f_{d,n},f_{s,n})s^{\mathrm{H}}(f_{d,n},f_{s,n})}\\ R_n=\mathrm{E}[n{n}^{\mathrm{H}}]={\sigma }^2{I}_{MN}\end{array}$

where ${\sigma }^2$ is the noise power. $I_{NM}$ is the $MN\times MN$ identity matrix.

In practice, the ideal CNCM is not easy to obtain directly. Taking the adjacent samples of the cell under test as the training samples and using the maximum likelihood method, we can estimate the correct CNCM using Formula (6), and the accuracy of CNCM estimation is related to the number of training samples. The number of training samples should be greater than twice the system’s degrees of freedom according to the RMB criterion.

(6)$\widehat{R}=\frac{1}{L}{\displaystyle \sum _{l=1}^L{x}_l{x}_l^{\mathrm{H}}}$

where $L$ denotes number of training samples.

Under the criterion of the maximum SINR, the optimal STAP weight vector can be obtained by solving the following problem:

(7)$\left\{\begin{array}{l}\underset{w}{\min }\hspace{1em}{w}^{\mathrm{H}}\widehat{R}w\\ \mathrm{s}.\mathrm{t}.w^{\mathrm{H}}s(f_{d,t},f_{s,t})=1\end{array}\right.$

Solving the problem in Formula (7), the optimal STAP weight vector of the filter is given by

(8)$w=\frac{{\widehat{R}}^{-1}s(f_{d,t},f_{s,t})}{s^{\mathrm{H}}(f_{d,t},f_{s,t}){\widehat{R}}^{-1}s(f_{d,t},f_{s,t})}$

where $s(f_{d,t},f_{s,t})$ is the space–time steering vector of the target.

Combining the sparse recovery techniques with STAP to make full use of the sparse characteristics of clutter data, the signal model can be represented as

(9)$\mathbf{x}=\mathbf{B}\mathbf{a}+\mathbf{n}$

where ${\rm B}=[b_1,\dots ,b_K]\in {ℂ}^{MN\times K}$ is a dictionary matrix, which can be obtained by uniformly dividing the angle–angular two–dimensional plane into $K=N_s{M}_d$ grid points. $N_s$ denotes the number of points in the spatial domain and $M_d$ denotes the number of points in the temporal domain. Usually, $N_s$ is set to five times the number of array elements and $M_d$ is set to five times the number of pulses. Each point corresponds to a space–time steering vector $b_k$ and we call it an atom. $a\in {ℂ}^{K\times 1}$ is a sparse coefficient and most elements in $a$ are zero, which correspond to signal components. $n\in {ℂ}^{MN\times 1}$ is the zero–mean Gaussian noise vector.

From the theory of sparse recovery, if we wish to solve the problem of representing $x$, we need to estimate $a$ by solving the following objective function, which can be expressed as

(10)$\begin{array}{l}\underset{a}{\min }\hspace{1em}{‖a‖}_0\\ \mathrm{s}.\mathrm{t}.\hspace{1em}{‖x-Ba‖}_2^2\le \epsilon \end{array}$

where $\epsilon$ denotes the fitting error tolerance.

However, solving this problem is difficult because the original $l_0$ norm problem is NP–hard. We need to replace the $l_0$ norm penalty term with the $l_1$ norm penalty term to obtain the solution:

(11)$\begin{array}{l}\underset{a}{\min }\hspace{1em}{‖a‖}_1\\ \mathrm{s}.\mathrm{t}.\hspace{1em}{‖x-Ba‖}_2^2\le \epsilon \end{array}$

In order to obtain the solution of sparse coefficient $a$, we must solve the objective function Formula (11) is equivalent to solving an underdetermined linear equation, which has the following form

(12)$a=\arg \underset{a}{\min }{‖x-Ba‖}_2^2+\lambda {‖a‖}_1$

where $\lambda$ denotes the regularization parameter.

2.2. Review of IAA

The classical IAA method aims to solve the weighted least squares problem, which is given by

(13)$\underset{v_k}{\min }{[x_l-v_k{b}_k]}^{\mathrm{H}}{R}_k^{-1}[x_l-v_k{b}_k]$

where $R_k=R-v_k^2{b}_k{b}_k^H$ represents the CNCM that does not contain the $b_k$ component. Paper [25] extends the conventional IAA algorithm to the MMV case, and the MMV case of the signal can be expressed as

(14)$\mathbf{X}=\mathbf{B}\mathbf{A}+\mathbf{N}$

where $X=[x_1,\dots ,x_L]\in {ℂ}^{MN\times L}$ is the clutter data containing $L$ training samples. The sparse coefficient matrix is denoted by $A=[a_1,\dots ,a_L]\in {ℂ}^{K\times L}$. $N=[n_1,\dots ,n_L]\in {ℂ}^{MN\times L}$ is the zero–mean Gaussian noise matrix.

The joint likelihood function with respect to the clutter data matrix $X$ can be expressed as

(15)$p(X|R)=\frac{1}{{\pi }^{NML}|R{|}^L}{e}^{-\mathrm{Tr}(X^{\mathrm{H}}{R}^{-1}X)}$

The maximum likelihood estimation (MLE) of the $v_k^2$ is equivalent to maximizing the likelihood function with respect to the clutter data matrix $X$. In order to simplify the calculations, the joint likelihood function with respect to the clutter data matrix $X$ will be written in the corresponding logarithmic form; then, the original problem will be transformed into the following form:

(16)$\begin{array}{cc}\hfill {\widehat{v}}_{m,n}^{MLE}& =\arg \max (p(X|R))\hfill \\ & =\arg \min (-\ln p(X|R))p(X|R)\hfill \\ & =\arg \min (L\ln |R|+\mathrm{Tr}(X^{\mathrm{H}}{R}^{-1}X))\hfill \end{array}$

where $\left|R\right|$ and $R^{-1}$ have the following form:

(17)$\begin{array}{cc}\hfill \left|R\right|& =|R_k+v_k^2{b}_k{b}_k^H|\hfill \\ & =|R_k\left|\right(1+v_k^2{b}_k^H{R}_k^{-1}{b}_k)\hfill \end{array}$

(18)$R^{-1}=R_k^{-1}-\frac{v_k^2{R}_k^{-1}{b}_k{b}_k^H{R}_k^{-1}}{1+v_k^2{b}_k^H{R}_k^{-1}{b}_k}$

Substituting (17) and (18) into (16), we have

(19)$\begin{array}{l}{\widehat{v}}_k^{MLE}=\arg \min (L\ln |R_k\left|\right(1+v_k^2{b}_k^H{R}_k^{-1}{b}_k)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} +\mathrm{Tr}(X^{\mathrm{H}}{R}_k^{-1}X-X^{\mathrm{H}}\frac{v_k^2{R}_k^{-1}{b}_k{b}_k^H{R}_k^{-1}}{1+v_k^2{b}_k^H{R}_k^{-1}{b}_k}X))\\ \hspace{1em}\hspace{1em} =\arg \min (L\ln (1+v_k^2{b}_k^H{R}_k^{-1}{b}_k)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} +L\ln \left|R_k\right|+\mathrm{Tr}(X^{\mathrm{H}}{R}_k^{-1}X)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} -\mathrm{Tr}(X^{\mathrm{H}}{v}_k^2{b}_k^H{R}_k^{-1}{b}_{k}X))\end{array}$

The derivative of Equation (19) with respect to $v_k^2$ is as follows:

(20)${\widehat{v}}_k^{MLE}=\frac{b_k^H{R}_k^{-1}({\widehat{R}}^{MLE}-R_k)R_k^{-1}{b}_k}{{(b_k^H{R}_k^{-1}{b}_k)}^2}$

where ${\widehat{R}}^{MLE}=\frac{1}{L}X{X}^{\mathrm{H}}$ represents the CNCM estimated from the IID training samples. Applying the matrix inversion lemma to obtain $R_k^{-1}$ and substituting it into Equation (20), we have

(21)$\begin{array}{cc}\hfill {\widehat{v}}_k^{MLE}& ={\widehat{v}}_k^{IAA}+v_k-{\widehat{v}}_k^{CAPON}\hfill \\ & =\frac{b_k^H{R}^{-1}{\widehat{R}}^{MLE}{R}^{-1}{b}_k}{{(b_k^H{R}^{-1}{b}_k)}^2}+v_k-\frac{1}{b_k^H{R}^{-1}{b}_k}\hfill \end{array}$

From Equation (21), we know that if $R$ is known, then $v_k\approx {\widehat{v}}_k^{CAPON}$. Thus, we can obtain ${\widehat{v}}_k^{MLE}\approx {\widehat{v}}_k^{IAA}$. It can be seen from Equation (21) that $R$ is required in the calculation of the ${\widehat{v}}_k^{IAA}$, so the solving process needs to be carried out in an iterative way.

2.3. Proposed Method

In this subsection, a fast IAA–based SR–STAP method for airborne radar is presented.

To begin, we extend the objective function in Equation (11) to the MMV case, and it is equivalent to the following objective function:

(22)$A^{*}=\underset{A}{\arg \min }{‖X-BA‖}_F^2+\lambda {‖A‖}_{2,0}$

where ${‖\cdot ‖}_{2,0}$ denotes the $l_{2,0}$ norm and it is calculated in two steps: first, we calculate the $l_2$ norm of each row of matrix A to form a new vector; second, we calculate the $l_0$ norm of each row of the new vector.

In order to solve the objective function (22), we define a new variable as follows:

(23)$\begin{array}{l}d={[d_1,\dots ,d_K]}^{\mathrm{T}}\in {ℂ}^{K\times 1}\\ d_k=\sqrt{\frac{1}{L}{\displaystyle \sum _{l=1}^L|a_{k,l}{|}^2}}\end{array}$

Then, Equation (22) can be written in the following form:

(24)$A^{*}=\underset{A}{\arg \min }{‖X-BA‖}_F^2+\lambda {‖d‖}_0$

As mentioned above, the $l_0$ norm problem is an NP–hard problem and the existing methods cannot directly solve this problem. Thus, the $l_0$ norm problem is relaxed to the $l_1$ norm problem, which can be written as

(25)$A^{*}=\underset{A}{\arg \min }{‖X-BA‖}_F^2+\lambda {‖d‖}_1$

However, in clutter signal recovery, the $l_1$ norm problem may not be a good approximation of the $l_0$ norm problem, because it also takes into account the amplitude of the sparse coefficient, which may ignore the role of smaller non–zero values of the sparse coefficient. This may lead to inaccurate results. According to the paper [19], we know that if we can find the upper bound of $d={[d_1,\dots ,d_K]}^{\mathrm{T}}\in {ℂ}^{K\times 1}$, then Equation (24) can be written as

(26)$\begin{array}{l}{A}^{*}=\underset{A}{\arg \min }{‖X-BA‖}_F^2+\lambda 1^{\mathrm{T}}z\\ s.t. d_k\le u_k{z}_k,z_k\in \left\{0,1\right\},k=1,\cdots ,K\end{array}$

where $u_k$ is the upper bound of $d_k$.

Further, we relax the constraint of Equation (26), and we can obtain

(27)$\begin{array}{l}{A}^{*}=\underset{A}{\arg \min }{‖X-BA‖}_F^2+\lambda 1^{\mathrm{T}}z\\ s.t. d_k\le u_k{z}_k,0\le z_k\le 1,k=1,\cdots ,K\end{array}$

Setting $h_k=1/u_k$, Equation (27) can be written as

(28)$A^{*}=\underset{A}{\arg \min }{‖X-BA‖}_F^2+\lambda {\displaystyle \sum _{k=1}^K{h}_k{d}_k}$

The objective function in Equation (28) is actually the standard form of the weighted $l_1$ norm directly derived from the canonical $l_0$ norm, which is better than the traditional $l_1$ norm [15].

From Equation (28), we know that $h_k$ is in fact the weight of the sparse coefficients. We need to choose the appropriate weights to meet not only the performance requirements but also the constraints mentioned above. Moreover, the setting of the weights is preferably related to the clutter environment. Based on this, we set the value of the weights to the IAA spectrum described above, which is

(29)$h_k=\sqrt{\frac{{(b_k^H{R}^{-1}{b}_k)}^2}{b_k^H{R}^{-1}{\widehat{R}}^{MLE}{R}^{-1}{b}_k}}=\sqrt{\frac{1}{{\widehat{v}}_k^{IAA}}}$

where ${\widehat{R}}^{MLE}$ is the CNCM estimated from the training samples, which is shown in (30). $R$ is the CNCM that is estimated using dictionary matrix $d_k^2$ and it is shown in (31).

(30)${\widehat{R}}^{MLE}=\frac{1}{L}X{X}^{\mathrm{H}}$

(31)$\mathbf{R}=\frac{1}{L}\sum _{l=1}^L\sum _{k=1}^K\left({\left|a_{k,l}\right|}^2\right){\mathbf{b}}_k{\mathbf{b}}_k^{\mathrm{H}}+{\sigma }^2\mathbf{I}$

According to the previous definition of $d_k$, is the average clutter power of the $k-th$ atom and ${\widehat{v}}_k^{IAA}$ represents the clutter plus noise power of the $k-th$ atom. Such a setting satisfies the condition $d_k\le u_k{z}_k$.

After choosing the appropriate weights, Equation (28) still presents a problem: solving the problem of the weighted $l_1$ norm requires expensive computation and it does not have an analytic solution. In practice, the clutter power is much greater than the noise power. The value of $d_k{h}_k$ is approximately equal to ${(d_k{h}_k)}^2$, which means that we use the weighted $l_2$ norm instead of the weighted $l_1$ norm. The advantage of this transformation is obvious: the following objective function has an analytic solution.

(32)$A^{*}=\underset{A}{\arg \min }{‖X-BA‖}_F^2+\lambda {\displaystyle \sum _{k=1}^K{h}_k^2{b}_k^2}$

The analytic solution of the objective function in Equation (32) is

(33)$A=W{B}^{\mathrm{H}}{\left(\lambda I+BW{B}^{\mathrm{H}}\right)}^{-1}X$

(34)$W=\left[\begin{array}{ccc}{\left(1/h_1\right)}^2& & \\ & \ddots & \\ & & {\left(1/h_K\right)}^2\end{array}\right]$

As we can see from Equations (29), (33), and (34), calculating $A$ requires the information of $h_k=1, k=1,\dots ,K$ and calculating $h_k=1, k=1,\dots ,K$ requires the information of $A$. Therefore, the proposed method must work in an iterative way.

Moreover, we need to estimate noise power ${\sigma }^2$ and choose the proper value of regularization parameter $\lambda$. According to paper [26], the value of $\lambda$ can be set equal to ${\sigma }^2$. The noise is updated by

(35)${\sigma }^2=\frac{{‖X-BA‖}_F^2}{MNL}$

The pseudo–code for the proposed method is given in Table 2.

Next, the computational complexity of the method proposed in this section is explained. Computational complexity is measured by the number of complex multiplications in a single iteration. The calculation amount of the method proposed in this section mainly comes from the calculation of each item in the sparse coefficient matrix and the weighted coefficient, which correspond to Formulas (33) and (34), respectively. In an iterative process, the computational complexity of the sparse coefficient matrix is $o\left({(MN)}^3+2MN{K}^2+2K{(MN)}^2+MNKL\right)$. The computational complexity of the weighted coefficient is $o\left(4{(MN)}^2+2MN\right)$.

3. Performance Assessment

In this section, simulation results will be given to show the effectiveness of the proposed method. This paper considers the uniform linear array and the clutter simulation is based on Ward’s clutter model. We choose the 600th range gate as CUT and six IID training samples are chosen. The dictionary matrix $B$ is obtained by dividing the spatial domain and Doppler domain evenly into $4N$ and $4M$ parts, respectively. The main parameters of a side–looking airborne radar system are listed in Table 3.

In this section, the performance of the proposed method is measured with the metrics of the improvement factor (IF), the CAPON spectrum, the average SINR loss versus the number of training samples and the number of iterations, the target detection performance, and the STAP output power against the range cell. At the end of this section, we will also give the average running time of the proposed method and the method used for comparison.

First, we briefly introduce the concepts of the improvement factor and SINR loss. The improvement factor is defined as the ratio of the output SINR to the input SINR and its value reflects the improvement in the detection performance of the airborne radar. The IF is defined as follows:

(36)$\mathrm{IF}=(CNR+1){\sigma }^2\frac{\left|w^{H}s\right|}{w^{H}Rw}$

The SINR loss has a linear relationship with the improvement factor and both of them can reflect the clutter suppression performance of the STAP method. The SINR loss is defined as follows:

(37)${\mathrm{SIN}\mathrm{R}}_{loss}=\frac{{\sigma }^2}{MN}\frac{\left|w^{H}s\right|}{w^{H}Rw}$

The proposed method is compared with the loaded sample matrix inversion (LSMI) method, the multiple orthogonal matching pursuit (M–OMP) method, the multiple sparse Bayesian learning with generalized double Pareto (GDP) prior (M–GDP) method, and the SR–STAP algorithm with the log summation penalty function (M–LOG). It should be noted that the difference between the M–GDP method and the M–SBL method lies in the difference in the sparse priors. The whole experiment is divided into two subsections: one subsection corresponds to the ideal case and the other subsection corresponds to the non–ideal case. For the non–ideal case, the amplitude error (standard deviation 0.03) and phase random error (standard deviation 2°) are added.

3.1. Ideal Case

In the first experiment, we plot the CAPON spectra of different STAP methods, as in Figure 2a–f, to illustrate the performance. From Figure 2a–f, we can see that the CAPON spectra plotted by M–LOG, M–GDP, and the proposed method are closest to the ideal CAPON spectrum, either in power or position. M–OMP can barely recover the correct CAPON spectrum and the performance of M–OMP is not acceptable. The M–OMP method is a greedy algorithm. Although the operation speed is fast, it is easy to regard the local optimal solution as the global optimal solution, so the estimated CAPON spectrum is not accurate. Due to the insufficient number of training samples, the LSMI method cannot estimate the correct CAPON spectrum.

In the second experiment, we plot the IF versus normalized Doppler frequency curves, as in Figure 3. Similar to the previous experimental conclusions, the M–LOG method, the M–GDP method, and the proposed method all can achieve near–optimal performance, while the M–OMP method and LSMI method achieve poor performance in both the main lobe and sidelobe regions.

In the third experiment, we plot the average SINR loss against the number of training samples. From the experimental results shown in Figure 4, it can be seen that under the ideal conditions, the proposed method in this paper has the best performance with only three training samples. When the number of training samples is less than three, the proposed method also has the minimum average SINR loss. The M–GDP method and the M–LOG method also can achieve good performance with three training samples, but the performance of these methods is not as high as that of the proposed method when the number of training samples is less than three. The average SINR loss of the M–OMP method and the LSMI method is poorer than that of the proposed method, the convergence speed is slow, and the average SINR loss is quite large.

In the fourth experiment, we evaluate the target detection performance of the different methods. We plot the probability of detection (PD) versus signal–to–noise (SNR) curves to show the target detection performance. Although there are many excellent detectors [27,28], the focus of our simulation is to observe the relative relationships between different methods, and there are no strict restrictions on the choice of detectors. In this section of the simulation experiment, the cell average constant false alarm rate (CA–CFAR) detector is adopted, and the probability false alarm rate is set as ${10}^{-4}$. We add targets from the 550th range gate to the 650th range gate, and the targets are located in the 77th–89th Doppler bins. As shown in Figure 5, except for the M–OMP method, all methods achieve better target detection performance.

In the last experiment in this subsection, to further illustrate the target detection performance, we compare the output power against the range gates for different methods. We set the SNR as 10 dB to demonstrate the performance. Targets are added in range gate 574 and two cases with different velocities are considered, where the normalized Doppler frequency is set to 0.1 for slow–moving targets and 0.3 for fast–moving targets. The experimental results are normalized for ease of presentation. It can be seen from Figure 6 that under the ideal conditions, the proposed method, the M–GDP method, and the M–LOG method can obtain better output power at the position of the target range gate. The M–OMP method has poor performance and it cannot form the output power at the position of the target range gate, which means that it is difficult to detect the target.

3.2. Non–Ideal Case

In the first experiment, we plot the CAPON spectra of the different STAP methods, as in Figure 7a–f, to illustrate the performance. Due to the presence of the amplitude error and phase error, the performance of all methods degrades to some extent. The LSMI method still faces the problem of insufficient training samples and it cannot recover the correct CAPON spectrum. The CAPON spectra of the M–OMP method and M–LOG method are spread out; the situation for the M–LOG method is more serious because this method is more sensitive to the error. M–GDP and the proposed method are closest to the ideal CAPON spectrum.

In the second experiment, we plot the IF versus normalized Doppler frequency curves, as in Figure 8. It can be clearly seen that the M–GDP method and the proposed method can achieve good performance compared with the other methods. The M–LOG method is slightly sensitive to the error and its performance is degraded severely in the main lobe region. The performance of the M–OMP method and LSMI method is not acceptable.

In the third experiment, we plot the average SINR loss against the number of training samples. From the experimental results shown in Figure 9, it can be seen that due to the presence of errors, the average SINR loss is increased for all methods, and the proposed method in this paper still has the best performance. When the number of training samples is more than three, the M–GDP method and M–LOG method can achieve the same performance as the proposed method. Otherwise, the performance of the proposed method is better. The average SINR loss of the M–OMP method and the LSMI method increases.

In the fourth experiment, we evaluate the target detection performance of the different methods. We plot the probability of detection (PD) versus signal–to–noise (SNR) curves to show the target detection performance. The other settings are consistent with the corresponding experiments in the previous subsection. As shown in Figure 10, we can see that our proposed method is slightly better than the M–GDP method. The target detection performance of M–LOG is affected by the error, which is mainly reflected in the poor detection ability for slow–moving targets.

In the last experiment in this subsection, to further illustrate the target detection performance, we compare the output power against the range gate for the different methods. Targets are added in range gate 574 and two cases with different velocities are considered, where the normalized Doppler frequency is set to 0.1 for slow–moving targets and 0.3 for fast–moving targets. The SNR of the target is set to 10 dB. It can be seen from Figure 11 that the proposed method can obtain better output power at the position of the target range gate. When dealing with slow–moving targets, the M–LOG method cannot form the highest output power at the position of the target range gate, which means that this method is not successful in detecting the slow–moving target. The M–OMP method has poor performance too. When dealing with fast–moving targets, all methods except the M–OMP method can form the output power at the position of the target range gate.

We illustrate the convergence performance of the different methods by plotting the average SINR loss against the number of iterations. Since only the M–GDP method and the M–LOG method are related to the number of iterations among the comparison methods used in this paper, the experimental results only refer to these methods. As can be seen from Figure 12a, the method proposed in this paper has an extremely high convergence speed, which is better than that of the other two methods in ideal conditions. In non–ideal conditions, the proposed method still has a fast convergence speed and the minimum average SINR loss.

At the end of this section, we compare the average running times of the different methods and the results are as follows. One hundred independent Monte Carlo trials are performed to obtain the average running time. The simulation platform is a notebook computer and the processor is an AMD Ryzen7 5800H with 8 cores and 16 threads. The MATLAB version is 2021B. It is clear from Table 4 that the proposed method in this paper runs significantly faster than the M–GDP method. The M–OMP method and M–LOG method run slightly faster than our proposed method, but the performance of these two methods is far weaker than that of the method proposed in this paper.

At the end of the simulation experiment, we also give the CAPON spectrum estimated by directly applying the IAA method to the STAP method under ideal conditions and error conditions. From Figure 13, we can see that the CAPON spectrum of this method is not acceptable under these conditions. The results are all far from the ideal CAPON spectrum, and it can be easily seen that the existence of errors causes the large spread of the spectrum.

3.3. Measured Data

In this section, we evaluate the performance of the proposed method via a public data set, the Mountain–Top data set (t38pre01v1CPI6). These data were measured by the Lincoln Laboratory, and the numbers of array elements and pulses of the radar system are 14 and 16, respectively. This data set contains 403 range gates and the target is at the 147th range gate. Using all 403 range gates, the clutter CAPON spectrum can be estimated, which is shown in Figure 14.

The STAP output power against the 130th to 165th range gates is depicted in Figure 15. In this experiment, the results are normalized. Due to the poor performance of the M–OMP method, which has been verified in the previous simulation experiments, the experimental results of this method are not shown. It can be seen from the figure that the proposed method, the M–GDP method, and the M–LOG method can all obtain better output power at the range gate where the target is located.

To further illustrate the performance of the proposed method, we present the average background noise power of the different methods in Table 5. The average background noise power is calculated from the output power of the remaining range gates after excluding the data of the range gate where the target is located and the data of the left and right protection range gates. It can be seen from the table that the method proposed in this paper is slightly better than the comparison methods.

4. Discussion

We present the performance of the proposed method in various experiments and thus we have a preliminary understanding of this method. In the ideal case, the performance of the proposed method is close to that of the M–GDP method and M–LOG method, and both of them can achieve excellent clutter suppression performance. In the non–ideal case, the performance of the M–LOG method is degraded. Only our proposed method and the M–GDP method can achieve better clutter suppression in this case. It is worth noting that the proposed method can achieve the same performance as the M–GDP method in a shorter running time, which is also the greatest advantage of the proposed method over the M–GDP method. According to previous research [24], the running time of the M–GDP method is much shorter than that of the M–SBL method. Meanwhile, the proposed method has a high convergence rate compared with the M–GDP method. From the results of the measured data, it can be seen that the proposed method also shows good clutter suppression performance.

5. Conclusions

In this paper, a fast IAA–based SR–STAP method for airborne radar was proposed. The proposed method combined the IAA method and weighted $l_1$ norm method. Compared with the STAP method, which directly uses the IAA method, the proposed method could estimate the CNCM more accurately. Compared with the weighted $l_1$ norm method, the proposed method had an analytical solution and could estimate the CNCM quickly and accurately. At the end of this paper, simulation results and a performance analysis were given to show the effectiveness of the proposed method.

Although the proposed method can estimate the clutter plus noise covariance matrix quickly and accurately, there are still many challenges to be addressed. In future research, it is necessary to propose a method to alleviate the influence of the presence of amplitude and phase errors. Furthermore, grid mismatch will affect the performance of the STAP method, and we also need to improve the performance of the method in such cases.

Author Contributions

Conceptualization, S.Z.; methodology, S.Z.; software, S.Z.; validation, S.Z., C.L. and B.R.; formal analysis, S.Z.; investigation, C.L.; resources, B.R.; writing—original draft preparation, S.Z.; writing—review and editing, S.Z.; visualization, S.Z.; supervision, T.W.; project administration, T.W.; funding acquisition, T.W. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors on request.

Conflicts of Interest

There are no conflicts of interest to be declared.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Platform geometry of airborne radar.

Figure 2. CAPON spectra of different methods. (a) The ideal CAPON spectrum; (b) the CAPON spectrum of LSMI; (c) the CAPON spectrum of M–OMP; (d) the CAPON spectrum of M–GDP; (e) the CAPON spectrum of M–LOG; (f) the CAPON spectrum of the proposed method.

Figure 2. CAPON spectra of different methods. (a) The ideal CAPON spectrum; (b) the CAPON spectrum of LSMI; (c) the CAPON spectrum of M–OMP; (d) the CAPON spectrum of M–GDP; (e) the CAPON spectrum of M–LOG; (f) the CAPON spectrum of the proposed method.

Figure 3. IF curves.

Figure 4. The average SINR loss against the number of training samples.

Figure 5. IF curves of probability of detection (PD) versus signal–to–noise (SNR).

Figure 6. STAP output power against the range gate for different methods. (a) Slow–moving target with SNR of 10 dB; (b) fast–moving target with SNR of 10 dB.

Figure 7. CAPON spectra of different methods. (a) The ideal CAPON spectrum; (b) the CAPON spectrum of LSMI; (c) the CAPON spectrum of M–OMP; (d) the CAPON spectrum of M–GDP; (e) the CAPON spectrum of M–LOG; (f) the CAPON spectrum of the proposed method.

Figure 7. CAPON spectra of different methods. (a) The ideal CAPON spectrum; (b) the CAPON spectrum of LSMI; (c) the CAPON spectrum of M–OMP; (d) the CAPON spectrum of M–GDP; (e) the CAPON spectrum of M–LOG; (f) the CAPON spectrum of the proposed method.

Figure 8. IF curves.

Figure 9. The average SINR loss against the number of training samples.

Figure 10. Probability of detection (PD) versus signal–to–noise (SNR) curves.

Figure 11. STAP output power against the range gate for different methods. (a) Slow–moving target with SNR of 10 dB; (b) fast–moving target with SNR of 10 dB.

Figure 12. The average SINR loss against the number of iterations: (a) ideal case; (b) non–ideal case.

Figure 13. CAPON spectrum: (a) ideal case; (b) non–deal case.

Figure 14. Estimated clutter CAPON spectrum with all 403 samples.

Figure 15. STAP output power against the range gate for different methods.

Table 1

Acronyms and definitions.

Acronyms	Full Name
AEW	Airborne early warning
CNCM	Clutter plus noise covariance matrix
CUT	Cell under test
DDD	Direct data domain
IAA	Iterative adaptive approach
IID	Independent and identically distributed
KA	Knowledge–aided
SBL	Sparse Bayesian learning
STAP	Space–time adaptive processing

Table 2

Pseudocode for Proposed method.

Step 1	Input data $X$ and dictionary matrix $B$
Step 2	Initialize ${\sigma }^2=1$ and $h_k=1, k=1,\dots ,K$
Step 3	Calculate $h_k=1, k=1,\dots ,K$ using (29)
Step 4	Calculate $W$ using (34)
Step 5	Calculate $A$ using (33)
Step 6	Calculate ${\sigma }^2$ using (35)
Step7	Repeat step 3, step 4, step 5, and step 6 until a stopping criterion is satisfied
Step 8	Calculate CNCM $R=\frac{1}{L}{\displaystyle \sum _{l=1}^L{\displaystyle \sum _{k=1}^K\left({\left|a_{k,l}\right|}^2\right)b_k{b}_k^{\mathrm{H}}}}+{\sigma }^2{\rm I}$
Step 9	Compute STAP weight $w$

Table 3

Simulation parameters of the radar system.

Parameter	Value
Number of elements in array	8
Number of pulses per CPI	8
Pulse repetition frequency	2000 Hz
Receiver bandwidth	2.5 MHz
Platform height	9000 m
Wavelength	0.3 m
Platform velocity	150 m/s
Clutter–to noise ratio	40 dB
Operation frequency	1 Ghz

Table 4

The average running time.

Algorithm	Average Running Time (s)
M–OMP	0.03
M–GDP	7.26
M–LOG	0.45
Proposed Method	0.22

Table 5

Average background noise power.

Algorithm	Power (dB)
Proposed Method	−12.0662
M–GDP	−10.5169
M–LOG	−10.1075

References

1. Ward, J. Space–Time Adaptive Processing for Airborne Radar; Lincoln Laboratory: Lexington, MA, USA, 1998.

2. Brennan, L.E.; Reed, I.S. Theory of adaptive radar. IEEE Trans. Aerosp. Electron. Syst.; 1973; 9, pp. 237-252. [DOI: https://dx.doi.org/10.1109/TAES.1973.309792]

3. Sedwick, R.J.; Hacker, T.L.; Marais, K. Performance analysis for an interferometric space–based GMTI system. Proceedings of the Record of the IEEE 2000 International Radar Conference [Cat. No. 00CH37037]; Alexandria, VA, USA, 12 May 2000; pp. 689-694.

4. Maher, J.; Callahan, M.; Lynch, D. Effects of clutter modeling in evaluating STAP processing for space–based radars. Proceedings of the Record of the IEEE 2000 International Radar Conference [Cat. No. 00CH37037]; Alexandria, VA, USA, 12 May 2000; pp. 565-570.

5. Devaney, A.J. Time reversal imaging of obscured targets from multistatic data. IEEE Trans. Antennas Propag.; 2005; 53, pp. 1600-1610. [DOI: https://dx.doi.org/10.1109/TAP.2005.846723]

6. Ciuonzo, D. On time–reversal imaging by statistical testing. IEEE Signal Process. Lett.; 2017; 24, pp. 1024-1028. [DOI: https://dx.doi.org/10.1109/LSP.2017.2704612]

7. Goldstern, J.S.; Reed, I.S.; Zulch, P.A. Multistage partially adaptive STAP CFAR detection algorithm. IEEE Trans. Aerosp. Electron. Syst.; 1999; 35, pp. 645-661. [DOI: https://dx.doi.org/10.1109/7.766945]

8. Tong, Y.L.; Wang, T.; Wu, J.X. Improving EFA–STAP performance using persymmetric covariance matrix estimation. IEEE Trans. Aerosp. Electron. Syst.; 2015; 51, pp. 924-936. [DOI: https://dx.doi.org/10.1109/TAES.2015.130264]

9. Goldstein, J.S.; Reed, I.S. Reduced rank adaptive filtering. IEEE Trans. Signal Process.; 1997; 45, pp. 493-496. [DOI: https://dx.doi.org/10.1109/78.554317]

10. Zhang, W.; He, Z.S.; Li, H.Y.; Li, J.; Duan, X. Beam–space reduced–dimension space–time adaptive processing for airborne radar in sample starved heterogeneous environments. IET Radar Sonar Navig.; 2016; 10, pp. 1627-1634.

11. Sarkar, T.K.; Wang, H.; Park, S.; Adve, R.; Koh, J.; Kim, K.; Zhang, Y.; Wicks, M.C.; Brown, R.D. A deterministic least squares approach to space time adaptive processing (STAP). IEEE Trans. Aerosp. Electron. Syst.; 2001; 49, pp. 91-103.

12. Cristallini, D.; Burger, W. A Robust Direct Data Domain Approach for STAP. IEEE Trans. Signal Process.; 2012; 60, pp. 1283-1294. [DOI: https://dx.doi.org/10.1109/TSP.2011.2176335]

13. Melvin, W.L.; Showman, G.A. An approach to knowledge–aided covariance estimation. IEEE Trans. Aerosp. Electron. Syst.; 2006; 42, pp. 1021-1042. [DOI: https://dx.doi.org/10.1109/TAES.2006.248216]

14. Stoica, P.; Li, J.; Zhu, X. On Using a priori Knowledge in Space–Time Adaptive Processing. IEEE Trans. Signal. Process.; 2008; 56, pp. 2598-2602. [DOI: https://dx.doi.org/10.1109/TSP.2007.914347]

15. Mallet, S.; Zhang, Z. Matching pursuits with time frequency dictionaries. IEEE Trans. Signal Process.; 1993; 41, pp. 3397-3415. [DOI: https://dx.doi.org/10.1109/78.258082]

16. Pati, Y.C.; Rezaiifar, R.; Krishnaprasad, P.S. Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition. Proceedings of the 27th Asilomar Conference on Signal, Systems and Computers; Pacific Grove, CA, USA, 1–3 November 1993; pp. 40-44.

17. Candes, E.; Wakin, M.; Boyd, S. Enhancing sparsity by reweighted minimization. J. Fourier Anal. Appl.; 2008; 14, pp. 877-905. [DOI: https://dx.doi.org/10.1007/s00041-008-9045-x]

18. Chartrand, R.; Yin, W. Iteratively reweighted algorithms for compressive sensing. Proceedings of the 2008 IEEE International Conference on Acoustics, Speech and Signal Processing; Las Vegas, NV, USA, 31 March–4 April 2008; pp. 3869-3872.

19. Xu, X.; Wei, X.H.; Ye, Z.F. DOA estimation based on sparse signal recovery utilizing weighted penalty. IEEE Signal Process. Lett.; 2012; 19, pp. 155-158. [DOI: https://dx.doi.org/10.1109/LSP.2012.2183592]

20. Tipping, M.E. Spars Bayesian learning and the relevance vector machine. J. Math. Learn.; 2001; 1, pp. 211-244.

21. Wipf, D.P.; Rao, B.D. An empirical Bayesian strategy for solving the simultaneous sparse approximation problem. IEEE Trans. Signal Process.; 2007; 55, pp. 3704-3716. [DOI: https://dx.doi.org/10.1109/TSP.2007.894265]

22. Duan, K.Q.; Wang, Z.T.; Xie, W.C.; Chen, H.; Wang, Y.L. Sparsity–based STAP algorithm with multiple measurement vectors via sparse Bayesian learning strategy for airborne radar. IET Signal Process.; 2017; 11, pp. 544-553. [DOI: https://dx.doi.org/10.1049/iet-spr.2016.0183]

23. Wang, Q.; Yu, H.; Li, J. Sparse Bayesian Learning Using Generalized Double Pareto Prior for DOA Estimation. IEEE Signal Process. Lett.; 2021; 28, pp. 1744-1748. [DOI: https://dx.doi.org/10.1109/LSP.2021.3104503]

24. Zhang, S.; Wang, T.; Liu, C.; Wang, D. A Space–Time Adaptive Processing Method Based on Sparse Bayesian Learning for Maneuvering Airborne Radar. Sensors; 2022; 22, 5479. [DOI: https://dx.doi.org/10.3390/s22155479] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35897983]

25. Yang, Z.C.; Li, X.; Wang, H.; Jiang, W. Adaptive clutter suppression based on iterative adaptive approach for airborne radar. Signal Process.; 2013; 93, pp. 3567-3577. [DOI: https://dx.doi.org/10.1016/j.sigpro.2013.03.033]

26. Fang, J.; Wang, F.Y.; Shen, Y.N.; Li, H.B.; Rick, S.B. Super–resolution compressed sensing for line spectral estimation: An iterative reweighted approach. IEEE Trans. Signal Process.; 2016; 64, pp. 4649-4662. [DOI: https://dx.doi.org/10.1109/TSP.2016.2572041]

27. Bandiera, F.; Miao, A.D.; Riccie, G. Adaptive CFAR radar detection with conic rejection. IEEE Trans. Signal Process.; 2007; 55, pp. 2533-2541. [DOI: https://dx.doi.org/10.1109/TSP.2007.893763]

28. Liu, W.J.; Liu, J.; Hao, C.P.; Gao, Y.C.; Wang, Y.L. Multichannel adaptive signal detection: Basic theory and literature review. Sci. China Inf. Sci.; 2022; 65, 121301. [DOI: https://dx.doi.org/10.1007/s11432-020-3211-8]