Clutter-Sensing-Driven Space-Time Adaptive Processing Approach for Airborne Sub-Array-Level Digital Array

Abstract

Sub-array-level digital arrays effectively diminish the computational complexity and sample demand of space-time adaptive processing (STAP), thus finding extensive applications in many airborne platforms. Nonetheless, airborne sub-array-level digital array radar still encounters pronounced performance deterioration in highly heterogeneous clutter environments due to inadequate training samples. To address this issue, a clutter-sensing-driven STAP approach for airborne sub-array-level digital arrays is proposed in this paper. Firstly, we derive a signal model of sub-array-level clutter sensing in detail and then further analyze the influence of the sidelobe characteristics of the conventional sub-array joint beam on clutter sensing. Secondly, a sub-array joint beam optimization model is proposed, which optimizes the sub-array joint beam into a wide beam with flat-top characteristics to improve the clutter-sensing performance in the beam sidelobe region. Finally, we decompose the complex optimization problem into two subproblems and then relax them into the low sidelobe-shaped beam pattern synthesisproblem and second-order cone programming problem, which can be effectively solved. The effectiveness of the proposed approach is validated in a real clutter environment through numerical experiments.

1. Introduction

As an airborne surveillance sensor, airborne radar has the characteristics of long sight distance, a wide observation field, and flexible deployment. However, airborne radar usually works in down-looking mode, and its echo contains high-intensity scattered echoes of ground objects, which leads to the deterioration of its detection performance of slow small targets. Space-time adaptive processing (STAP) is an effective clutter suppression method, which can improve the moving target detection performance of airborne radar [1]. By combining spatial and temporal degrees of freedom, STAP constructs an optimal adaptive filter in the space-time two-dimensional plane, so it can effectively suppress strong clutter coupled with angle and Doppler [2,3]. The performance of a space-time adaptive filter depends on the estimation accuracy of the clutter plus the noise covariance matrix (CNCM) of the cell under test (CUT). According to the well-known Reed–Mallett–Brennan (RMB) rule [4], STAP requires independent and identically distributed (IID) training samples with at least twice the system degrees of freedom (DoFs) to ensure that the signal-to-clutter-noise ratio (SCNR) loss is less than 3 dB. However, in most practical applications, due to the influence of heterogeneous clutter, it is difficult for space-time adaptive filters to obtain sufficient training samples. Moreover, in some extreme cases, each adjacent range ring has a completely different clutter distribution. Therefore, determining how to improve the performance of space-time adaptive filters under limited-sample conditions is essential.

To tackle the aforementioned issues, several methods have been proposed. The existing heterogeneous clutter suppression methods are mainly divided into five classifications: the power heterogeneous suppression, non-homogeneity detector (NHD), direct data domain STAP (DDD-STAP), sparse recovery STAP (SR-STAP), and knowledge-aided STAP (KA-STAP) methods.

The power heterogeneous suppression methods assume that the non-homogeneity of clutter primarily originates from significant power disparities among training samples from different range cells [5,6]. The primary idea of these methods is to select samples that are strong enough compared to the power level of the sample from the CUT. Rabideau et al. [5] tackled this issue by measuring the sample power of the CUT and then selecting samples with sufficiently high power as training data. Kogon et al. [6] select samples with both sufficiently high power and a phase distribution similar to that of clutter as training samples. These two methods effectively solve the issue of insufficient nulls in space-time filters due to the difference in sample power.

The NHD methods primarily address the target cancelation issue resulting from target-like interference within the training samples. These methods aim to detect and remove samples containing target-like interference by constructing appropriate statistical measures [7,8]. Melvin et al. [7] initially proposed using the generalized inner product (GIP) as a statistical measure to detect and eliminate training samples containing target-like interference, validating this method’s performance using real MCARM data. Guo et al. [8] employed the GIP method to filter samples, followed by utilizing the volume cross correlation function as a similarity measure between the training sample and the CUT sample, to select appropriate samples. This approach enabled more comprehensive extraction of reliable samples from the sample set.

The DDD-STAP methods solely utilize samples from the target cell to mitigate the impact of heterogeneous clutter [9]. Sarkar et al. [9] derived training samples by applying spatial or temporal smoothing to the data from the CUT. These methods solely rely on samples from the CUT; theoretically, they can completely address the issue of heterogeneous clutter. However, the smoothing operations in this approach reduce the system’s DoFs and limit its only applicability to uniformly spaced arrays.

The SR-STAP methods have been widely researched in recent years. The existing research results show that SR-STAP still has excellent clutter suppression performance under small-sample conditions [10,11,12,13,14,15]. The existing SR-STAP methods are based on the sparse characteristics of the clutter space-time spectrum [16]. In such methods, the space-time plane is divided into several grids. For the CUT echo data, the clutter ridge only occupies a small part of the grid and has a large space-time spectrum coefficient. The noise signal occupies the remaining grids, and the space-time spectral coefficient is relatively small. The SR-STAP method first uses the sparse recovery (SR) technique to estimate the space-time spectrum of the CUT, and then uses the space-time spectrum to estimate the CNCM of the CUT. However, SR-STAP methods often require high computational resources and are difficult to directly apply to actual airborne radar systems [17].

The KA-STAP methods enhance STAP’s performance to suppress heterogeneous clutter by introducing prior knowledge of clutter [18,19,20,21,22,23,24,25]. Generally speaking, KA-STAP methods can be divided into two categories [26]: (i) The first category involves the KA-STAP method indirectly using knowledge, using prior knowledge to guide the selection of training strategies and training samples, and then estimate the CNCM for adaptive filtering. (ii) The second category involves the KA-STAP method using prior knowledge to construct the prior CNCM directly and fusing it with the sampling covariance matrix (SCM) to generate the CNCM for the adaptive filter. At present, the knowledge used by KA-STAP mainly includes terrain data (digital terrain elevation map, land cover, land use, etc.) or other sensor data (synthetic aperture radar (SAR) image, satellite optical image, etc.). The matching degree between prior knowledge and the actual clutter environment directly determines the performance of KA-STAP [26]. However, when the external clutter environment changes, the poor timeliness of terrain data and sensor data would degrade clutter suppression performance.

Following the concept of cognitive radars, a STAP method based on dynamic clutter sensing is proposed in [27,28], in which, before target detection, the radar transmits a set of orthogonal signals to sense the clutter information and then reconstruct the RCS of clutter patches from the clutter echo as the prior clutter information. In the target detection stage, prior clutter information is used to assist STAP in heterogeneous clutter suppression. The experimental results show that this method can effectively improve the heterogeneous clutter suppression performance of STAP. However, this method will face the following problems in practical applications: (i) only some of the array elements are used as transmitting array elements, and the transmitting power of the array is not fully utilized. (ii) implementing processing at the element level for airborne early warning radar systems is almost impossible due to the difficulty of engineering implementation, the dramatic increase in algorithmic complexity due to too many DoFs, and the expensive cost. So far, almost all existing large airborne early warning (AEW) radars use sub-array structures [29]. However, the signal model in [27,28] adopts an element-level array and it is not applicable for actual airborne radar systems. Therefore, it is necessary to study a clutter-sensing approach suitable for sub-array-level digital array radar in practice.

To address the issues mentioned above, a clutter-sensing-driven STAP approach for an airborne sub-array-level digital array is proposed in this paper. The main contributions of this paper are summarized as follows:

• This work focuses on the sub-array architecture adopted in many airborne early warning radar systems. We aim to effectively resolve the issue of the existing clutter-sensing STAP method’s inability to be directly applied to sub-array-level digital arrays. To address the degradation of clutter-sensing performance in the sidelobe region of sub-array-level digital arrays caused by the low sidelobe gain of the static sub-array joint beam, we propose a sub-array joint beam optimization method. This optimization method enhances clutter-sensing performance in the sidelobe region by configuring the sub-array joint beam into a wide beam with specific flat-top characteristics.

• In order to solve the complex optimization problem effectively, we decompose the original optimization problem into two subproblems. After decomposition, the transmitting sub-array weight design and receiving sub-array weight design are decoupled, which greatly reduces the complexity of the optimization problem.

• To address the beam optimization problem of transmitting sub-array, we convert the optimization problem into a low sidelobe-shaped beam pattern synthesis (SBPS) problem, and use the existing penalty dual decomposition (PDD) method [30] to solve the problem effectively. To resolve the receiving sub-array beam optimization problem, we relax the non-convex optimization problem into a convex second-order cone programming problem, which can be solved effectively by using the existing interior-point method [31].

The rest of this article is organized as follows. In Section 2, the echo signal model of an airborne sub-array-level digital array in the sensing stage is given, and the corresponding SR model is constructed. In addition, we further analyze the adverse effect of the static sub-array joint beam on clutter sensing. In Section 3, we give the corresponding optimization method of the sub-array joint beam. Section 4 gives the method of prior CNCM estimation and the design method of the space-time adaptive filter. In Section 5, numerical simulations are given. Finally, the conclusions are drawn in Section 6.

Notation: Vectors are denoted by lowercase bold letters, and matrixes are denoted by uppercase bold letters. ⊗ and ⊙, respectively, represent the Kronecker and Hadamard products. * denotes the convolution operation. ${\mathbb{C}}^{N\times 1}$ and ${\mathbb{C}}^{N\times N}$ represent an N-dimensional complex vector and an $N\times N$-dimensional complex matrix, respectively. ${\mathbb{R}}^{N\times 1}$ represents an N-dimensional real vector. The matrix transpose, complex conjugate, and conjugate transpose are denoted by ${\left(·\right)}^{\mathrm{T}}$, ${\left(·\right)}^{\ast }$, and ${\left(·\right)}^{\mathrm{H}}$, respectively. $\left|\mathbf{x}\right|$ is the modulus of $\mathbf{x}$. ${∥\mathbf{x}∥}_0$ and ${∥\mathbf{x}∥}_2$ represent the ${\ell }_0$-norm and the ${\ell }_2$-norm of $\mathbf{x}$, respectively. ${\left(·\right)}^2$ is defined as a square operation. $\mathrm{tr}\left(·\right)$ denotes the trace of the matrix.

2. Signal Model in Sensing Stage

In this paper, we consider a side-looking airborne sub-array-level digital array radar system with a uniform linear array (ULA) composed of N elements, and the inter-element spacing is d, as shown in Figure 1. The receiving antenna and transmitting antenna of the radar are co-located, and the sub-array division diagram is shown in Figure 2. The transmitting array is divided into $N_1$ uniform non-overlapping sub-arrays, and the number of antennas in each sub-array is equal to $N_{sub1}$, i.e., $N_1{N}_{sub1}=N$. Therefore, the spacing between the first antennas of the adjacent transmitting sub-array is $d_1=N_{sub1}d$. The receiving array is divided into $N_2$ uniform overlapping sub-arrays; the number of antennas in each sub-array is $N_{sub2}$, and the spacing between the first antennas of the adjacent receiving sub-arrays is $d_2=d$, i.e., $N_2+N_{sub2}-1=N$. Within a coherent processing interval (CPI), the radar transmits $K_s$ pulses with a fixed pulse repetition frequency (PRF) $f_r=1/T_r$, where $T_r$ is the pulse repetition interval (PRI). The altitude and velocity of the platform are h and $v_0$, respectively. Figure 3 depicts the signal processing flow of the radar. Before the target detection, the radar transmits signals to actively sense the external clutter environment in what is referred to as the sensing stage.

As shown in Figure 4, the sensing stage and the detection stage are arranged at different CPIs, that is, they are serial in time. Specifically, in the sensing stage (sensing CPI), the radar acquires and stores prior knowledge; in the detection stage, the radar uses prior knowledge to assist in target detection. The prior knowledge acquired by a sensing CPI can be used in several subsequent detection CPIs.

2.1. Echo Signal Model in Sensing Stage

During the sensing stage, a set of orthogonal waveforms are generated as the transmitting signals of the different sub-arrays, which can be expressed as $\mathbf{p}\left(t\right)={\left[p_1\left(t\right),\dots ,p_{N_1}\left(t\right)\right]}^{\mathrm{T}}$. All antennas within the $n_1$-th sub-array coherently emit the signal $p_{n_1}\left(t\right)$, $n_1=1,\dots ,N_1$. Let us introduce the $N\times 1$ vector ${\mathbf{v}}_{n_1}\in {\mathbb{R}}^{N\times 1}$, whose elements contain only 0 and 1. Here, 1 indicates that the array antenna corresponding to its index is located in the $n_1$-th transmitting sub-array, while 0 indicates that the array antenna corresponding to its index does not belong to the $n_1$-th transmitting sub-array. Then, the spatial steering vector associated with the $n_1$-th transmitting sub-array ${\tilde{\mathbf{a}}}_{t,n_1}\left(\phi \right)\in {\mathbb{C}}^{N\times 1}$ can be expressed as

(1)${\tilde{\mathbf{a}}}_{t,n_1}\left(\phi \right)\triangleq {\mathbf{v}}_{n_1}\odot {\widehat{\mathbf{a}}}_t\left(\phi \right),n_1=1,\cdots ,N_1,$

where ${\widehat{\mathbf{a}}}_t\left(\phi \right)\in {\mathbb{C}}^{N\times 1}$ is the spatial steering vector for the entire transmitting array, which can be expressed as

(2)${\widehat{\mathbf{a}}}_t\left(\phi \right)={\left[1,e^{j{\displaystyle \frac{2\pi dsin\phi }{\lambda }}},\cdots ,e^{j{\displaystyle \frac{2\pi \left(N-1\right)dsin\phi }{\lambda }}}\right]}^{\mathrm{T}}.$

The signal emitted by the $n_1$-th transmitting sub-array can be expressed as

(3)${\mathbf{S}}_{n_1}\left(t\right)={\tilde{\mathbf{w}}}_{t,n_1}^{\ast }{p}_{n_1}\left(t\right),n_1=1,2,\dots ,N_1$

where ${\tilde{\mathbf{w}}}_{t,n_1}\in {\mathbb{C}}^{N\times 1}$ denotes the transmitting beam-forming weight vector which contains $N_{sub1}$ non-zero elements corresponding to the contained antenna of the $n_t$-th transmitting sub-array.

Suppose there is a clutter patch, clutter patch C, whose azimuth angle, pitch angle, and slant range are denoted as $\theta$, $\varphi$, and r, respectively. The baseband component of the reflected signal of clutter patch C can be expressed as

(4)$\begin{array}{cc}\hfill & r\left(t\right)=\gamma ·\overline{\alpha }\sum _{n_1=1}^{N_1}{\tilde{\mathbf{w}}}_{t,n_1}^{\mathrm{H}}{\tilde{\mathbf{a}}}_{t,n_1}\left(\phi \right)p_{n_1}\left(t-{\displaystyle \frac{\tau }{2}}\right)\hfill \\ \hfill & =\gamma ·\overline{\alpha }\sum _{n_1=1}^{N_1}\left[{\mathbf{w}}_{t,n_1}^{\mathrm{H}}{\mathbf{a}}_{t,n_1}\left(\phi \right)e^{-j{\zeta }_{t,n_1}\left(\phi \right)}{p}_{n_1}\left(t-{\displaystyle \frac{\tau }{2}}\right)\right]\hfill \end{array},$

where $\gamma$ is the RCS of clutter patch C, and $\overline{\alpha }=\sqrt{{\displaystyle \frac{P_t}{4\pi r^2}}}{e}^{-j{\displaystyle \frac{2\pi f_{0}r}{c}}}$. $P_t$ is the transmitting power of a single transmitting antenna, $f_0$ is the carrier frequency, c represents the speed of light, $\phi$ is the cone angle of clutter patch C, with $sin\phi =sin\theta cos\varphi$, and $\tau =\raisebox{1ex}{$2r$}\!\left/ \!\raisebox{-1ex}{$c$}\right.$ is round-trip propagation delay of the signal. The phase delay ${\zeta }_{t,n_1}\left(\phi \right)={\displaystyle \frac{2\pi \left(n_1-1\right)d_{1}sin\phi }{\lambda }}$, $n_1=1,2,\dots ,N_1$ is caused by the spacing between the first antenna of the entire transmitting array and the first antenna of the $n_1$-th transmitting sub-array. ${\mathbf{w}}_{t,n_1}\in {\mathbb{C}}^{N_{sub1}\times 1}$ and ${\mathbf{a}}_{t,n_1}\left(\phi \right)\in {\mathbb{C}}^{N_{sub1}\times 1}$ denote the beam-forming weight vector and the spatial steering vector corresponding to the antennas contained in the $n_1$-th transmitting sub-array, respectively. Suppose the reference antenna of ${\mathbf{a}}_{t,n_1}\left(\phi \right)$ is the first antenna of the $n_1$-th transmitting sub-array; then, ${\mathbf{a}}_{t,n_1}\left(\phi \right)$ can be expressed as

(5)${\mathbf{a}}_{t,n_1}\left(\phi \right)={\left[1,e^{j{\displaystyle \frac{2\pi dsin\phi }{\lambda }}},\cdots ,e^{j{\displaystyle \frac{2\pi \left(N_{sub1}-1\right)dsin\phi }{\lambda }}}\right]}^{\mathrm{T}},n_1=1,2,\dots ,N_1$

Let us introduce the transmitting sub-array beam gain vector ${\mathbf{g}}_t\left(\phi \right)\in {\mathbb{C}}^{N_1\times 1}$ and waveform diversity vector ${\mathbf{a}}_{\mathrm{T}}\left(\phi \right)\in {\mathbb{C}}^{N_1\times 1}$, i.e.,

(6)${\mathbf{g}}_t\left(\phi \right)={\left[{\mathbf{w}}_{t,1}^{\mathrm{H}}{\mathbf{a}}_{t,1}\left(\phi \right),\cdots ,{\mathbf{w}}_{t,N_1}^{\mathrm{H}}{\mathbf{a}}_{t,N_1}\left(\phi \right)\right]}^{\mathrm{T}}$

and

(7)${\mathbf{a}}_{\mathrm{T}}\left(\phi \right)={\left[e^{-{\zeta }_{t,1}\left(\phi \right)},\cdots ,e^{-{\zeta }_{t,N_1}\left(\phi \right)}\right]}^{\mathrm{T}},$

Then, the space-time echo signal from clutter patch C can be expressed as

(8)$\begin{array}{cc}\hfill \mathbf{x}& \left(t\right)=\alpha ·\gamma ·\left[{\mathbf{a}}_d\left(\phi \right)\otimes {\widehat{\mathbf{a}}}_r\left(\phi \right)\right]·\left[{\left({\mathbf{a}}_{\mathrm{T}}\left(\phi \right)\odot {\mathbf{g}}_t\left(\phi \right)\right)}^{\mathrm{T}}\mathbf{p}\left(t-\tau \right)\right]\hfill \end{array}$

where $\alpha ={\displaystyle \frac{\sqrt{P_t}}{4\pi r^2}}{e}^{-j{\displaystyle \frac{4\pi f_{0}r}{c}}}$. ${\mathbf{a}}_d\left(\phi \right)\in {\mathbb{C}}^{K_s\times 1}$ is the temporal steering vector, which can be expressed as

(9)${\mathbf{a}}_d\left(\phi \right)={\left[1,e^{j{\displaystyle \frac{4\pi v_{0}sin\phi }{\lambda f_r}}},\cdots ,e^{j{\displaystyle \frac{4\pi \left(K_s-1\right)v_{0}sin\phi }{\lambda f_r}}}\right]}^{\mathrm{T}}.$

${\widehat{\mathbf{a}}}_r\left(\phi \right)\in {\mathbb{C}}^{N\times 1}$ is the spatial steering vector for the entire receiving array.

Similar to the definitions of ${\mathbf{w}}_{t,n_1}$ and ${\mathbf{a}}_{t,n_1}\left(\phi \right)$, we define ${\mathbf{w}}_{r,n_2}\in {\mathbb{C}}^{N_{sub2}\times 1}$ and ${\mathbf{a}}_{r,n_2}\left(\phi \right)\in {\mathbb{C}}^{N_{sub2}\times 1}$ as the beam-forming weight vector and the spatial steering vector of the $n_2$-th receiving sub-array, respectively, where $n_2=1,\dots ,N_2$. The specific form of ${\mathbf{a}}_{r,n_2}\left(\phi \right)\in {\mathbb{C}}^{N_{sub2}\times 1}$ can be expressed as

(10)${\mathbf{a}}_{r,n_2}\left(\phi \right)={\left[1,e^{j{\displaystyle \frac{2\pi dsin\phi }{\lambda }}},\cdots ,e^{j{\displaystyle \frac{2\pi \left(N_{sub2}-1\right)dsin\phi }{\lambda }}}\right]}^{\mathrm{T}},n_2=1,2,\dots ,N_2$

Then, after receiving sub-array beam-forming, the space-time echo signal of clutter patch C can be expressed as

(11)$\begin{array}{cc}\hfill {\mathbf{x}}^{\mathrm{dbf}}\left(t\right)=& \alpha ·\gamma ·\left\{{\mathbf{a}}_d\left(\phi \right)\otimes {\left[{\mathbf{a}}_R\left(\phi \right)\odot {\mathbf{g}}_r\left(\phi \right)\right]}^{\mathrm{T}}\right\}·\left\{{\left[{\mathbf{a}}_{\mathrm{T}}\left(\phi \right)\odot {\mathbf{g}}_t\left(\phi \right)\right]}^{\mathrm{T}}\mathbf{p}\left(t-\tau \right)\right\}\hfill \end{array}.$

The specific forms of ${\mathbf{g}}_r\left(\phi \right)\in {\mathbb{C}}^{N_2\times 1}$ and ${\mathbf{a}}_R\left(\phi \right)\in {\mathbb{C}}^{N_2\times 1}$ are as follows

(12)${\mathbf{g}}_r\left(\phi \right)={\left[{\mathbf{w}}_{r,1}^{\mathrm{H}}{\mathbf{a}}_{r,1}\left(\phi \right),\cdots ,{\mathbf{w}}_{r,N_2}^{\mathrm{H}}{\mathbf{a}}_{r,N_2}\left(\phi \right)\right]}^{\mathrm{T}},$

(13)${\mathbf{a}}_{\mathrm{R}}\left(\phi \right)={\left[e^{-{\zeta }_{r,1}\left(\phi \right)},\cdots ,e^{-{\zeta }_{r,N_2}\left(\phi \right)}\right]}^{\mathrm{T}},$

where ${\zeta }_{r,n_2}\left(\phi \right)={\displaystyle \frac{2\pi \left(n_2-1\right)d_{2}sin\phi }{\lambda }},n_2=1,2,\dots ,N_2$ is the phase delay caused by the spacing between the first antenna of the entire receiving array and the first antenna of the $n_2$-th receiving sub-array.

The transmitting and receiving arrays are uniformly divided into different sub-arrays as described before. Suppose that the same beam-forming weight vectors are used for each sub-array, i.e., ${\mathbf{w}}_t={\mathbf{w}}_{t,1}=\cdots ={\mathbf{w}}_{t,N_1}$ and ${\mathbf{w}}_r={\mathbf{w}}_{r,1}=\cdots ={\mathbf{w}}_{r,N_2}$. Thus, ${\mathbf{g}}_t\left(\phi \right)$ and ${\mathbf{g}}_r\left(\phi \right)$ can be expressed as

(14)$\begin{array}{cc}\hfill & {\mathbf{g}}_t\left(\phi \right)={\left[{\mathbf{w}}_{t,1}^{\mathrm{H}}{\mathbf{a}}_{t,1}\left(\phi \right),\cdots ,{\mathbf{w}}_{t,N_1}^{\mathrm{H}}{\mathbf{a}}_{t,N_1}\left(\phi \right)\right]}^{\mathrm{T}}\hfill \\ \hfill & ={\mathbf{w}}_t^{\mathrm{H}}{\mathbf{a}}_t\left(\phi \right)·{\mathbf{1}}_{N_1\times 1}\hfill \\ \hfill & ={\mathrm{g}}_t\left(\phi \right)·{\mathbf{1}}_{N_1\times 1}\hfill \end{array},$

and

(15)$\begin{array}{cc}\hfill & {\mathbf{g}}_r\left(\phi \right)={\left[{\mathbf{w}}_{r,1}^{\mathrm{H}}{\mathbf{a}}_{r,1}\left(\phi \right),\cdots ,{\mathbf{w}}_{r,N_2}^{\mathrm{H}}{\mathbf{a}}_{r,N_2}\left(\phi \right)\right]}^{\mathrm{T}}\hfill \\ \hfill & ={\mathbf{w}}_r^{\mathrm{H}}{\mathbf{a}}_r\left(\phi \right)·{\mathbf{1}}_{N_2\times 1}\hfill \\ \hfill & ={\mathrm{g}}_r\left(\phi \right)·{\mathbf{1}}_{N_2\times 1}\hfill \end{array},$

where ${\mathbf{1}}_{N_1\times 1}$ and ${\mathbf{1}}_{N_2\times 1}$ denote the column vector with all elements being 1.

After matching filters, the echo signal of clutter patch C can be expressed as

(16)$\begin{array}{c}\hfill \mathbf{y}\left(t\right)=\alpha ·\gamma ·g_t\left(\phi \right)·g_r\left(\phi \right)·{\mathbf{a}}_d\left(\phi \right)\otimes {\mathbf{a}}_R\left(\phi \right)\otimes \mathbf{R}\left(t\right){\mathbf{a}}_{\mathrm{T}}\left(\phi \right)\end{array}.$

According to the Ward clutter model [1], the echo signal of a single range ring can be expressed as the sum of the echo signals of all clutter patches on the corresponding range ring. Therefore, the space-time snapshot vector of the i-th range cell can be expressed as

(17)$\begin{array}{cc}\hfill {\mathbf{y}}_i=& \sum _{j=1}^{N_c}{\mathbf{y}}_{ij}\left(t_i\right)+\mathbf{n}\left(t_i\right)\hfill \\ \hfill =& \sum _{j=1}^{N_c}\{{\alpha }_i·{\gamma }_{ij}·g_t\left({\phi }_{ij}\right)·g_r\left({\phi }_{ij}\right)·{\mathbf{a}}_d\left({\phi }_{ij}\right)\otimes {\mathbf{a}}_R\left({\phi }_{ij}\right)\otimes {\mathbf{R}}_i\left(t_i\right){\mathbf{a}}_{\mathrm{T}}\left({\phi }_{ij}\right)\}+\mathbf{n}\left(t_i\right)\hfill \end{array},$

where $t_i$ represents the sampling time corresponding to the i-th range cell, $N_c$ is the number of clutter patches in a single range ring, $\mathbf{n}\left(t\right)$ stands for zero-mean Gaussian white noise, and ${\mathbf{R}}_i\left(t\right)$ is the matched filter coefficient matrix which is given in (18).

(18)${\mathbf{R}}_i\left(t\right)=\left[\begin{array}{cccc}{p}_1\left(t-{\tau }_i\right)\ast p_1^{\ast }\left(-t\right)& p_2\left(t-{\tau }_i\right)\ast p_1^{\ast }\left(-t\right)& \cdots & p_{N_1}\left(t-{\tau }_i\right)\ast p_1^{\ast }\left(-t\right)\\ p_1\left(t-{\tau }_i\right)\ast p_2^{\ast }\left(-t\right)& p_2\left(t-{\tau }_i\right)\ast p_2^{\ast }\left(-t\right)& \cdots & p_{N_1}\left(t-{\tau }_i\right)\ast p_2^{\ast }\left(-t\right)\\ ⋮& ⋮& \ddots & ⋮\\ p_1\left(t-{\tau }_i\right)\ast p_{N_1}^{\ast }\left(-t\right)& p_2\left(t-{\tau }_i\right)\ast p_{N_1}^{\ast }\left(-t\right)& \cdots & p_{N_1}\left(t-{\tau }_i\right)\ast p_{N_1}^{\ast }\left(-t\right)\end{array}\right]$

2.2. Sparse Recovery Model

Compactly, Equation (17) can be expressed as the following linear equation:

(19)${\mathbf{y}}_i={\mathbf{H}}_i{\gamma }_i+\mathbf{n}\left(t_i\right)$

where ${\gamma }_i={\left[{\gamma }_{i1},{\gamma }_{i2},\cdots ,{\gamma }_{i{N}_c}\right]}^{\mathrm{T}}$ is the RCS vector corresponding to all clutter patches in i-th range ring. ${\mathbf{H}}_i\in {\mathbb{C}}^{N_1{N}_2{M}_s\times N_c}$ is an over-complete dictionary, which can be express as

(20)$\begin{array}{cc}\hfill & {\mathbf{H}}_i={\alpha }_i\left[\begin{array}{c}g\left({\phi }_{i1}\right)·{\mathbf{a}}_d\left({\phi }_{i1}\right)\otimes {\mathbf{a}}_R\left({\phi }_{i1}\right)\otimes {\mathbf{R}}_i\left(t_i\right){\mathbf{a}}_T\left({\phi }_{i1}\right),\\ g\left({\phi }_{i2}\right)·{\mathbf{a}}_d\left({\phi }_{i2}\right)\otimes {\mathbf{a}}_R\left({\phi }_{i2}\right)\otimes {\mathbf{R}}_i\left(t_i\right){\mathbf{a}}_T\left({\phi }_{i2}\right),\\ ⋮\\ g\left({\phi }_{i{N}_c}\right)·{\mathbf{a}}_d\left({\phi }_{i{N}_c}\right)\otimes {\mathbf{a}}_R\left({\phi }_{i{N}_c}\right)\otimes {\mathbf{R}}_i\left(t_i\right){\mathbf{a}}_T\left({\phi }_{i{N}_c}\right)\end{array}\right],\hfill \end{array}$

where $g\left(\phi \right)$ in (20) can be expressed as $g\left(\phi \right)=g_t\left(\phi \right)g_r\left(\phi \right)$.

Since the heterogeneity of the clutter scene comes from discrete and strong patches, generally speaking, the RCS of this kind of clutter patch is much larger than that of other clutter patches in the clutter scene, so only some of the elements in the RCS vector have a relatively large modal value, and the modal values of the rest of the elements are relatively small or close to zero [32]. It should be noted that the RCS vector ${\gamma }_i$ of the i-th range loop can be obtained by solving linear Equation (19), and the row dimension of matrix ${\mathbf{H}}_i$ is much less than the column dimension, i.e., $N_1{N}_2{K}_s<N_c$; hence, Equation (19) is heavily ill posed. According to [33], under such conditions on matrix ${\mathbf{H}}_i$ and vector ${\gamma }_i$, Equation (19) for its solution ${\gamma }_i$ can be formulated as

(21)${\widehat{\gamma }}_i=\underset{{\gamma }_i}{argmin}{∥{\gamma }_i∥}_0\mathrm{s}.\mathrm{t}.{∥{\mathbf{y}}_i-{\mathbf{H}}_i{\gamma }_i∥}_2\le \epsilon ,$

where $\epsilon$ represents the error tolerance, which is generally set as the power value of radar system noise. Several algorithms have been proposed to solve (21), such as orthogonal matching pursuit (OMP) [34], basis pursuit denoising (BPDN) [35], and so on. Among them, the BPDN algorithm is widely used to solve sparse recovery problems because of its computational efficiency. This algorithm uses $l_1$-norm instead of $l_0$-norm to solve the NP-hard optimization problem in (21), which can be expressed as

(22)${\widehat{\gamma }}_i=\underset{{\gamma }_i}{argmin}{∥{\gamma }_i∥}_1\mathrm{s}.\mathrm{t}.{∥{\mathbf{y}}_i-{\mathbf{H}}_i{\gamma }_i∥}_2\le \epsilon ,$

Solving the constrained optimization problem in the form of (22) can transform the constraint conditions into penalty terms, thus obtaining the following unconstrained optimization problem:

(23)${\widehat{\gamma }}_i=\underset{{\gamma }_i}{argmin}{\displaystyle \frac{1}{2}}{∥{\mathbf{y}}_i-{\mathbf{H}}_i{\gamma }_i∥}_2^2+\kappa {∥{\gamma }_i∥}_1,$

where $\kappa$ is a user-defined parameter to balance the precision and sparsity of the final solution. Problem (23) is a second-order cone programming problem, which can be solved by the interior-point method [36]. It is worth noting that the observation signal ${\mathbf{y}}_i$ contains both the expected signal ${\mathbf{H}}_i{\gamma }_i$ and useless noise signal $\mathbf{n}\left(t_i\right)$. Since the noise level is fixed, the value of each element in ${\mathbf{H}}_i$ heavily effects the estimation accuracy of ${\gamma }_i$, which in fact depends on beam patterns $g_t\left(\phi \right)$ and $g_r\left(\phi \right)$ (c.f. Equation (20)). Intuitively, this observation is straightforward, as the shape of the radar beam determines the intensity of the echo in different directions, which, in turn, affects the sensing accuracy of clutter patches. Later, we will employ a simple calculation to demonstrate that if the conventional beam-forming weight vectors are adopted, specifically ${\mathbf{w}}_t={\mathbf{a}}_t\left({\phi }_0\right)$ and ${\mathbf{w}}_r={\mathbf{a}}_r\left({\phi }_0\right)$, where ${\phi }_0$ represents the beam direction of the sub-array, then the variation in $g\left(\phi \right)=g_t\left(\phi \right)g_r\left(\phi \right)$ with respect to different angles $\phi$ would exceed 40 dB. Such a significant difference profoundly affects the estimation accuracy of ${\gamma }_i$.

Recall $g_t\left(\phi \right)={\mathbf{w}}_t^{\mathrm{H}}{\mathbf{a}}_t\left(\phi \right)$ and $g_r\left(\phi \right)={\mathbf{w}}_r^{\mathrm{H}}{\mathbf{a}}_r\left(\phi \right)$ in (14) and (15), which represent the beam gain of the transmitting sub-array and receiving sub-array, respectively. Then, the joint beam gain $g\left(\phi \right)$ of the sub-array for the clutter patches at different angles can be expressed as

(24)$\begin{array}{cc}\hfill & \left|g\left(\phi \right)\right|=\left|g_t\left(\phi \right)·g_r\left(\phi \right)\right|\hfill \\ \hfill =& \left|\sum _{n_{sub1}=1}^{N_{sub1}}{e}^{{\displaystyle \frac{-j2\pi n_{sub1}d\left(sin\phi -sin{\phi }_0\right)}{\lambda }}}·\sum _{n_{sub2}=1}^{N_{sub2}}{e}^{{\displaystyle \frac{-j2\pi n_{sub2}d\left(sin\phi -sin{\phi }_0\right)}{\lambda }}}\right|\hfill \\ \hfill =& \left|{\displaystyle \frac{sin\left[\left(N_{sub1}d\pi /\lambda \right)\left(sin\phi -sin{\phi }_0\right)\right]}{sin\left[\left(d\pi /\lambda \right)\left(sin\phi -sin{\phi }_0\right)\right]}}\right|·\left|{\displaystyle \frac{sin\left[\left(N_{sub2}d\pi /\lambda \right)\left(sin\phi -sin{\phi }_0\right)\right]}{sin\left[\left(d\pi /\lambda \right)\left(sin\phi -sin{\phi }_0\right)\right]}}\right|\hfill \end{array}.$

Equation (24) shows that the sub-array joint beam gain can be regarded as the product of the beam gain of a single transmitting sub-array and the beam gain of a single receiving sub-array. Figure 5 shows the sub-array joint beam gain at different angles. Obviously, the static joint beam gain has the highest gain near the beam direction ${\phi }_0$; however, the sidelobe gain is much lower than the mainlobe gain, i.e.,

(25)$20·{log}_{10}{\displaystyle \frac{g\left(\phi \right)}{g\left({\phi }_0\right)}}<-40\mathrm{dB},\mathrm{if}\left|\phi -{\phi }_0\right|>{11.6}^{\circ }.$

Lower gain in the sidelobes suggests that clutter patches falling within the sidelobe region will exhibit reduced echo, which could result in greater estimation errors for the vector ${\gamma }_i$.

3. Joint Optimization of Sub-Array Beams

To address the problem of low sidelobe gain associated with a static sub-array joint beam, one straightforward approach is to boost the average transmission power of the radar. However, this is constrained in practice by the radar transmitter’s peak power capacity and duty cycle, preventing arbitrary increases in power. This paper proposes a solution that focuses on optimizing the beam-forming vectors for both the transmitting and receiving sub-arrays. The objective is to create a joint beam with an expanded mainlobe, which helps to minimize the disparity in clutter patch reconstruction across various spatial domains.

Let the desired range of the mainlobe for the sub-array joint beam be defined as ${\Theta }_{\mathrm{m}}\triangleq \left[{\phi }_l,{\phi }_r\right]$, and the sidelobe range be defined as ${\Theta }_{\mathrm{s}}\triangleq \left[-{90}^{\circ },{\phi }_l\right]\cup \left[{\phi }_r,{90}^{\circ }\right]$. To achieve the characteristics of a flat-top beam within the sub-array joint beam, the optimized joint beam gain should be confined to a certain range within the mainlobe region, while simultaneously minimizing sidelobe gain. Furthermore, to ensure that each transmitting element transmits at full power, a constant modulus constraint is applied to the weight vector of the transmitting sub-array. For simplicity, we can normalize the transmission power of each element without any loss of generality, and the proposed optimization problem is

$\begin{array}{cc}\hfill \underset{\eta ,{\mathbf{w}}_t,{\mathbf{w}}_r}{min}& \eta \hfill \end{array}$

(26a)$\begin{array}{cc}\hfill \mathrm{s}.\mathrm{t}.& L_{\mathrm{m}}\le {\left|g_t\left({\phi }_{\mathrm{m}}\right)·g_r\left({\phi }_{\mathrm{m}}\right)\right|}^2\le U_{\mathrm{m}},{\phi }_{\mathrm{m}}\in {\Theta }_{\mathrm{m}}\hfill \end{array}$

(26b)$\begin{array}{c}{\left|g_t\left({\phi }_{\mathrm{s}}\right)·g_r\left({\phi }_{\mathrm{s}}\right)\right|}^2\le \eta ,{\phi }_{\mathrm{s}}\in {\Theta }_{\mathrm{s}}\hfill \end{array}$

(26c)$\begin{array}{c}{g}_t\left({\phi }_{\mathrm{m}}\right)={\mathbf{w}}_t^{\mathrm{H}}{\mathbf{a}}_t\left({\phi }_{\mathrm{m}}\right),{\phi }_{\mathrm{m}}\in {\Theta }_{\mathrm{m}}\hfill \end{array}$

(26d)$\begin{array}{c}{g}_r\left({\phi }_{\mathrm{m}}\right)={\mathbf{w}}_r^{\mathrm{H}}{\mathbf{a}}_r\left({\phi }_{\mathrm{m}}\right),{\phi }_{\mathrm{m}}\in {\Theta }_{\mathrm{m}}\hfill \end{array}$

(26e)$\begin{array}{c}{g}_t\left({\phi }_{\mathrm{s}}\right)={\mathbf{w}}_t^{\mathrm{H}}{\mathbf{a}}_t\left({\phi }_{\mathrm{s}}\right),{\phi }_{\mathrm{s}}\in {\Theta }_{\mathrm{s}}\hfill \end{array}$

(26f)$\begin{array}{c}{g}_r\left({\phi }_{\mathrm{s}}\right)={\mathbf{w}}_r^{\mathrm{H}}{\mathbf{a}}_r\left({\phi }_{\mathrm{s}}\right),{\phi }_{\mathrm{s}}\in {\Theta }_{\mathrm{s}}\hfill \end{array}$

(26g)$\begin{array}{c}\left|w_{t,n_{sub1}}\right|=\sqrt{{\displaystyle \frac{1}{N_{sub1}}}},n_{sub1}=1,\dots ,N_{sub1}\hfill \end{array}$

In (26), $L_{\mathrm{m}}$ and $U_{\mathrm{m}}$ are the lower and upper bounds in the desired broad mainlobe, $\eta$ is the sidelobe level to be minimized, ${\phi }_{\mathrm{m}}$ denotes any sampling angle within the mainlobe region, and ${\phi }_{\mathrm{s}}$ denotes any sampling angle within the sidelobe region. It is noteworthy that we did not restrict the absolute value of the joint beam gain in (26a). This is because the amplitudes of elements in the receiving sub-array beam-forming vector were not constrained. Therefore, the optimized receiving sub-array weight vector naturally converges to a solution that satisfies (26a).

Constraints (26a), (26b), and (26g) make problem (26) non-convex. The optimization of ${\mathbf{w}}_t$ must be phase-only, while the optimization of ${\mathbf{w}}_r$ does not need to restrict the magnitude of the elements to be equal. To tackle problem (26), we decompose it into the following two subproblems, and a solution can be attained by solving these two subproblems.

Subproblem for transmitting beam optimization:

$\begin{array}{cc}\hfill \underset{\eta ,{\mathbf{w}}_t}{min}& \eta \hfill \end{array}$

(27a) $\begin{array}{cc}\hfill \mathrm{s}.\mathrm{t}.& L_{\mathrm{m}}\le {\left|{\mathbf{w}}_t^{\mathrm{H}}{\mathbf{a}}_t\left({\phi }_{\mathrm{m}}\right)\right|}^2\le U_{\mathrm{m}},{\phi }_{\mathrm{m}}\in {\Theta }_{\mathrm{m}}\hfill \end{array}$

(27b) $\begin{array}{c}{\left|{\mathbf{w}}_t^{\mathrm{H}}{\mathbf{a}}_t\left({\phi }_{\mathrm{s}}\right)\right|}^2\le \eta ,{\phi }_{\mathrm{s}}\in {\Theta }_{\mathrm{s}}\hfill \end{array}$

(27c) $\begin{array}{c}\left|w_{t,n_{sub1}}\right|=\sqrt{{\displaystyle \frac{1}{N_{sub1}}}},n_{sub1}=1,\dots ,N_{sub1}\hfill \end{array}$

Subproblem for receiving beam optimization:

$\begin{array}{cc}\hfill \underset{\eta ,{\mathbf{w}}_r}{min}& \eta \hfill \end{array}$

(28a) $\begin{array}{cc}\hfill \mathrm{s}.\mathrm{t}.& L_{\mathrm{m}}\le {\left|{\mathbf{w}}_r^{\mathrm{H}}{\mathbf{a}}_r\left({\phi }_{\mathrm{m}}\right)\right|}^2\le U_{\mathrm{m}},{\phi }_{\mathrm{m}}\in {\Theta }_{\mathrm{m}}\hfill \end{array}$

(28b) $\begin{array}{c}{\left|{\mathbf{w}}_r^{\mathrm{H}}{\mathbf{a}}_r\left({\phi }_{\mathrm{s}}\right)\right|}^2\le \eta ,{\phi }_{\mathrm{s}}\in {\Theta }_{\mathrm{s}}\hfill \end{array}$

Problem (26) is divided into two separate optimization problems: the first pertains to the transmitting sub-array weight vector, and the second concerns the receiving sub-array weight vector. This bifurcation successfully reduces the issue of mutual coupling between the weight vectors ${\mathbf{w}}_t$ and ${\mathbf{w}}_r$.

Next, we develop algorithms to solve the two subproblems.

• (a). Solving subproblem (27)

For the sake of discussion, we rewrite subproblem (27) as follows:

$\begin{array}{cc}\hfill \underset{\eta ,{\mathbf{w}}_t}{min}& \eta \hfill \end{array}$

(29a)$\begin{array}{cc}\hfill \mathrm{s}.\mathrm{t}.& L_{\mathrm{m}}\le {\left|{\mathbf{w}}_t^{\mathrm{H}}{\mathbf{a}}_t\left({\phi }_{\mathrm{m}}\right)\right|}^2\le U_{\mathrm{m}},{\phi }_{\mathrm{m}}\in {\Theta }_{\mathrm{m}}\hfill \end{array}$

(29b)$\begin{array}{c}{\left|{\mathbf{w}}_t^{\mathrm{H}}{\mathbf{a}}_t\left({\phi }_{\mathrm{s}}\right)\right|}^2\le \eta ,{\phi }_{\mathrm{s}}\in {\Theta }_{\mathrm{s}}\hfill \end{array}$

(29c)$\begin{array}{c}\left|w_{t,n_{sub1}}\right|=\sqrt{{\displaystyle \frac{1}{N_{sub1}}}},n_{sub1}=1,\dots ,N_{sub1}\hfill \end{array}$

For given amplitude constraints (29a) and (29c), according to the Cauchy–Schwarz inequality, there must be such an upper bound $U_b$ that ${\left|{\mathbf{w}}_t^{\mathrm{H}}{\mathbf{a}}_t\left({\phi }_p\right)\right|}^2\le U_b$. If the lower bound of the given mainlobe constraint satisfies $L_{\mathrm{m}}>U_b$, then constraints (29a) and (29c) cannot be satisfied at the same time, which leads to no feasible solution for problem (29).

In order to avoid the existence of the above scale uncertainty problem, we need to introduce the scale scaling factor s, and problem (29) can be relaxed as follows:

$\begin{array}{cc}\hfill \underset{s,\eta ,{\mathbf{w}}_t}{min}& \eta \hfill \end{array}$

(30a)$\begin{array}{cc}\hfill \mathrm{s}.\mathrm{t}.& L_{\mathrm{m}}\le {\left|s·{\mathbf{w}}_t^{\mathrm{H}}{\mathbf{a}}_t\left({\phi }_{\mathrm{m}}\right)\right|}^2\le U_{\mathrm{m}},{\phi }_{\mathrm{m}}\in {\Theta }_{\mathrm{m}}\hfill \end{array}$

(30b)$\begin{array}{c}{\left|s·{\mathbf{w}}_t^{\mathrm{H}}{\mathbf{a}}_t\left({\phi }_{\mathrm{s}}\right)\right|}^2\le \eta ,{\phi }_{\mathrm{s}}\in {\Theta }_{\mathrm{s}}\hfill \end{array}$

(30c)$\begin{array}{c}\left|w_{t,n_{sub1}}\right|=\sqrt{{\displaystyle \frac{1}{N_{sub1}}}},n_{sub1}=1,\dots ,N_{sub1}\hfill \end{array}$

Problem (30) can be regarded as an SBPS problem with amplitude constraints. Solving the aforementioned optimization problem is challenging due to two main factors: (1) constraint (30a) is non-convex, and (2) the variables s and ${\mathbf{w}}_t$ are interdependent, impacting each other throughout the solution process. To address problem (30), the authors in [30] present an efficient algorithm based on PDD, the details of which are not included here.

• (b). Solving subproblem (28)

Similarly, for the sake of discussion, we rewrite problem (28) as follows:

$\begin{array}{cc}\hfill \underset{\eta ,{\mathbf{w}}_r}{min}& \eta \hfill \end{array}$

(31a)$\begin{array}{cc}\hfill \mathrm{s}.\mathrm{t}.& L_{\mathrm{m}}\le {\left|{\mathbf{w}}_r^{\mathrm{H}}{\mathbf{a}}_r\left({\phi }_{\mathrm{m}}\right)\right|}^2\le U_{\mathrm{m}},{\phi }_{\mathrm{m}}\in {\Theta }_{\mathrm{m}}\hfill \end{array}$

(31b)$\begin{array}{c}{\left|{\mathbf{w}}_r^{\mathrm{H}}{\mathbf{a}}_r\left({\phi }_{\mathrm{s}}\right)\right|}^2\le \eta ,{\phi }_{\mathrm{s}}\in {\Theta }_{\mathrm{s}}\hfill \end{array}$

Because of the absolute operation in the constraints (31a) and (31b), problem (31) remains non-convex. Therefore, we relax the constraints with absolute operations into norm constraints.

First, we relax the mainlobe gain constraint to

(32)${∥{\mathbf{w}}_r^{\mathrm{H}}{\mathbf{A}}_m-\mathbf{1}∥}_2\le {\epsilon }_m.$

In (32), ${\mathbf{A}}_m$ is the steering vector matrix corresponding to the mainlobe of the flat-top broad beam, which can be expressed as

(33)${\mathbf{A}}_m=\left[{\mathbf{a}}_r\left({\phi }_{m,1}\right),{\mathbf{a}}_r\left({\phi }_{m,2}\right),\cdots ,{\mathbf{a}}_r\left({\phi }_{m,P_m}\right)\right],$

where ${\phi }_{m,1}<{\phi }_{m,2}<\cdots <{\phi }_{m,P_m}\in {\Theta }_m$, $P_m$ represents the sampling number of angles in the mainlobe region. ${\epsilon }_m$ is a pre-defined toleance threshold used to control the consistency of beam gain at various angles within the mainlobe. Inequality (32) can be understood to mean that the norm value of the error between the beam gain at different angles in the constrained mainlobe and the given value 1 is less than ${\epsilon }_m$. Then, we similarly give the sidelobe constraints as follows:

(34)${∥{\mathbf{w}}_r^{\mathrm{H}}{\mathbf{A}}_s∥}_2\le \eta .$

In (34), ${\mathbf{A}}_s$ is the steering vector matrix corresponding to the sidelobe of the flat-top beam, which can be expressed as

(35)${\mathbf{A}}_s=\left[{\mathbf{a}}_r\left({\phi }_{s,1}\right),{\mathbf{a}}_r\left({\phi }_{s,2}\right),\cdots ,{\mathbf{a}}_r\left({\phi }_{s,P_s}\right)\right],$

where ${\phi }_{s,1}<{\phi }_{s,2}<\cdots <{\phi }_{s,P_s}\in {\Theta }_s$. $P_s$ represents the sampling number of angles in the sidelobe region, and $\eta$ is a tolerance to be optimized to control the sidelobe gain of the receiving sub-array.

By combining (32) and (34), we have the following receiving sub-array optimization model:

$\begin{array}{cc}\hfill \underset{\eta ,{\mathbf{w}}_r}{min}& \eta \hfill \end{array}$

(36a)$\begin{array}{cc}\hfill \mathrm{s}.\mathrm{t}.& {∥{\mathbf{w}}_r^{\mathrm{H}}{\mathbf{A}}_m-\mathbf{1}∥}_2\le {\epsilon }_m\hfill \end{array}$

(36b)$\begin{array}{c}{∥{\mathbf{w}}_r^{\mathrm{H}}{\mathbf{A}}_s∥}_2\le \eta \hfill \end{array}$

The constraints (36a) and (36b) can be considered a special form of second-order cone constraint, so problem (36) is a second-order cone program (SOCP), which can be easily solved using the interior-point method [31].

According to the literature [30], the PDD algorithm first needs to calculate the eigenvalue of the matrix once outside the algorithm loop, with a complexity of $O\left(N^3\right)$, and the single iteration complexity of the weight vector is $O\left(max\left(NL,N^2\right)\right)$. Therefore, the overall complexity of the algorithm is $O\left(N^3+I_{in}·I_{out}max\left(NL,N^2\right)\right)$, where $I_{in}$ and $I_{out}$ represent the maximum number of inner and outer loops of the algorithm, respectively, N represents the dimension of the optimization weight vector, and L represents the number of inequality constraints. Assume that when solving problem (30), the number of angular sampling points in the main lobe region and sidelobe region are $L_m$ and $L_s$, respectively. For the optimization problem (30), the dimension of the optimization weight vector ${\mathbf{w}}_t$ is $N_{sub1}\times 1$, and the number of inequality constraints is $L_m+L_s$. Therefore, the computational complexity of using the PDD algorithm to solve problem (30) is $O\left(N_{sub1}^3+I_{in}·I_{out}max\left(N_{sub1}\left(L_m+L_s\right),N_{sub1}^2\right)\right)$.

According to the literature [37], the number of iterations of the interior-point method in the worst case is $O\left(max{\left\{L,N\right\}}^4{N}^{0.5}log\left({\epsilon }^{-1}\right)\right)$, where L represents the number of constraints, N represents the number of variables, and $\epsilon >0$ represents the given solution accuracy. Therefore, the number of iterations in the worst case when the interior-point method is used to solve problem (36) is $O\left(N_{sub2}^{4.5}log\left({\epsilon }^{-1}\right)\right)$. For a single iteration, the computational complexity is $O\left(\left(N_{sub2}+1\right)\left(L_m+L_s\right)\right)$, so the total algorithm complexity of the interior-point method for problem (36) is $O\left(\left(N_{sub2}+1\right)\left(L_m+L_s\right)N_{sub2}^{4.5}log\left({\epsilon }^{-1}\right)\right)$.

4. Prior CNCM Estimation and Space-Time Filter Design

Upon finishing the sensing stage, the radar transitions to the detection stage. It is a recognized fact in an AEW radar scenario that a specific ground area will contribute to clutter echoes across multiple coherent processing intervals (CPIs) [38]. This indicates that the scattering properties of the clutter in consecutive CPIs evolve gradually. Consequently, we can leverage the clutter’s prior information, acquired during the sensing phase, to estimate the prior Clutter-to-Noise Covariance Matrix (CNCM), which aids in clutter mitigation during the detection phase.

4.1. Estimation of Prior CNCM

Due to the motion of the platform, the observation view of the same clutter patch differs between the sensing stage and the detection stage. Therefore, it is necessary to update the position information of the clutter patch before utilizing the prior information.

As shown in Figure 6, points $O^{\prime }$ and $O^{″}$ are the positions of the platform in the sensing stage and the detection stage respectively. Suppose that clutter patch C is the j-th clutter patch within the i-th range ring. In the detection stage, the updated formula for the slant range and the cone angle of clutter patch C can be written as

(37)$\begin{array}{cc}\hfill & {\overline{r}}_{ij}=\sqrt{r_{ij}^2+\Delta y^2-2sin\left({\phi }_{ij}\right)·\Delta y·r_{ij}},\hfill \\ \hfill & {\overline{\phi }}_{ij}=arccos\left({\displaystyle \frac{r_{ij}cos\left({\phi }_{ij}\right)}{{\overline{r}}_{ij}}}\right).\hfill \end{array}$

where ${\overline{r}}_{ij}$ and ${\overline{\phi }}_{ij}$ are the slant range and cone angle of clutter patch C in the detection stage respectively, and $\Delta y$ is the distance between $O^{\prime }$ and $O^{″}$.

In the detection stage, the radar transmits a linear frequency-modulated (LFM) signal $u\left(t\right)$. The prior CNCM of the CUT can be estimated according to the RCS information, updated slant range, and updated cone angle of the clutter patches, and the radar operating parameters. The prior CNCM can be constructed as follows:

(38)$\begin{array}{cc}\hfill {\mathbf{R}}_{\mathrm{prior},i}=& \sum _{j=1}^{N_c}\left\{{\left|{\xi }_{ij}\right|}^2\left[{\mathbf{b}}_d\left({\overline{\phi }}_{ij}\right){\mathbf{b}}_d^{\mathrm{H}}\left({\overline{\phi }}_{ij}\right)\right]\otimes \left[{\left|{\mathbf{g}}_r\left({\overline{\phi }}_{ij}\right)\right|}^2{\mathbf{a}}_R\left({\overline{\phi }}_{ij}\right){\mathbf{a}}_R^{\mathrm{H}}\left({\overline{\phi }}_{ij}\right)\right]\right\}+{\widehat{\sigma }}_n^2{\mathbf{I}}_{N_r{K}_d}\hfill \end{array},$

(39)${\xi }_{ij}={\displaystyle \frac{{\widehat{\gamma }}_{ij}{G}_t\left({\overline{\phi }}_{ij}\right)\sqrt{P_t}}{4\pi {{\overline{r}}_{ij}}^2}}{e}^{{\displaystyle \frac{-j4\pi f_0{r}_i}{c}}}u\left(t_i-{\tau }_i\right)\ast u^{\ast }\left(-t_i\right),$

where ${\xi }_{ij}$ represents the complex amplitude of the echo signal in the detection stage. ${\widehat{\gamma }}_{ij}$ can be obtained by (23). $t_i$ indicates the sampling time of the CUT. ${\widehat{\sigma }}_n^2$ is the noise power estimated in the detection stage. $G_t\left({\overline{\phi }}_{ij}\right)={\mathbf{a}}_t^{\mathrm{H}}\left({\phi }_0\right){\mathbf{a}}_t\left({\overline{\phi }}_{ij}\right)·{\mathbf{a}}_{\mathrm{T}}^{\mathrm{H}}\left({\phi }_0\right){\mathbf{a}}_{\mathrm{T}}\left({\overline{\phi }}_{ij}\right)$ represents the beam gain of the transmitting array. ${\mathbf{I}}_{N_r{K}_d}$ represents the $N_r{K}_d$-dimension unit matrix. ${\mathbf{b}}_d\left({\overline{\phi }}_{ij}\right)$ represents the temporal steering vector of the detection stage, which can be expressed as

(40)${\mathbf{b}}_d\left({\overline{\phi }}_{ij}\right)=\left[1,e^{j{\displaystyle \frac{4\pi v_{0}sin{\overline{\phi }}_{ij}}{\lambda f_r}}},\cdots ,e^{j{\displaystyle \frac{4\pi v_0\left(K_d-1\right)sin{\overline{\phi }}_{ij}}{\lambda f_r}}}\right],$

where $K_d$ represents the number of pulses contained in a CPI in the detection stage.

4.2. Design of Space-Time Adaptive Filter

The estimated clutter RCS vector ${\widehat{\gamma }}_i$ inevitably has some degree of error, and this error will affect the estimation accuracy of the clutter subspace. If the prior CNCM in (38) is directly used to calculate the space-time adaptive filter, this would lead to the degradation of clutter suppression performance. Therefore, the color-loading (CL) algorithm is introduced to correct clutter in the prior CNCM using the SCM.

Let ${\left\{{\mathbf{y}}_l\right\}}_{l=1}^L$ represent the sample data adjacent to the CUT; then, the SCM can be expressed as

(41)${\mathbf{R}}_{\mathrm{smi},i}={\displaystyle \frac{1}{L}}\sum _{l=1}^L{\mathbf{z}}_l{\mathbf{z}}_l^{\mathrm{H}},$

where ${\mathbf{z}}_l$ is the space-time snapshot of the detection stage, and L represents the number of training samples, usually taken as twice the system degree of freedom $2N_r{K}_d$.

According to [18], the CNCM of the CUT corrected by the CL algorithm can be expressed as

(42)${\mathbf{R}}_{\mathrm{CL},i}=\beta {\mathbf{R}}_{\mathrm{prior},i}+\left(1-\beta \right){\mathbf{R}}_{\mathrm{smi},i},0\le \beta \le 1,$

where $\beta$ represents the color-loading coefficient, which is used to adjust the weight of the prior CNCM ${\mathbf{R}}_{\mathrm{prior},i}$ in the color-loading CNCM ${\mathbf{R}}_{\mathrm{CL},i}$. Then, the STAP weight vector is given by

(43)${\mathbf{w}}_{\mathrm{opt},i}={\displaystyle \frac{{\mathbf{R}}_{\mathrm{CL},i}^{-1}{\mathbf{b}}_0}{{\mathbf{b}}_0^{\mathrm{H}}{\mathbf{R}}_{\mathrm{CL},i}^{-1}{\mathbf{b}}_0}},$

where ${\mathbf{b}}_0$ represents the desired spatial-temporal steering vector. The main steps of the proposed approach are summarized in Table 1.

4.3. Computational Complexity Analysis

Table 2 gives the computational complexity analysis of the entire algorithm process. $L_r$ represents the number of sampling points of a single echo pulse, $L_t$ represents the number of points of the reference signal during pulse compression, and $\epsilon >0$ represents the specified solution accuracy of the interior-point method. $L_r$ represents the number of sampling points of a single echo pulse, $L_t$ represents the number of points of the reference signal during pulse compression, and $\epsilon >0$ represents the specified solution accuracy of the interior-point method.

5. Simulation Experiments

In this section, the efficacy of the proposed method is validated through numerical experiments, and its performance is benchmarked against other algorithms. For the experiment, a synthetic aperture radar (SAR) image of a specified area, depicted in Figure 7, serves as a realistic clutter environment, with each pixel representing a scattering source. The platform parameters and radar system specifications are detailed in Table 3.

5.1. Sparsity of Actual Clutter Scene

In this experiment, the clutter scene is divided into different range rings based on range resolution, and each range ring is divided into different clutter patches at intervals of ${0.25}^{\circ }$. As shown in Figure 8, the sum of the intensity of SAR image pixels contained in the clutter patch is taken as the RCS of the clutter patch, and the RCS of the j-th clutter patch on the i-th range ring can be expressed as

(44)${\gamma }_{ij}=\sum _{q=1}^Q{\vartheta }_q,$

where Q represents the number of SAR image pixels contained in the clutter patch, and ${\vartheta }_q$ represents the intensity of the pixels.

Assume that the initial position of the platform is at x = 2500 and y = 0 in the clutter scene in Figure 6. According to (44), the real clutter scene, as shown in Figure 9, can be obtained. It can be seen that in the actual clutter scene, only some of the clutter patches have a large RCS, while most of the clutter patches have a small RCS. Figure 10 shows a box plot of the RCS amplitude of all clutter patches in the actual clutter scene. The maximum value of the upper edge of the box plot, 1.778, is taken as the threshold to distinguish the common clutter patch and the discrete and strong clutter patch. The quantity ratio and power ratio of the common clutter patch and the discrete and strong clutter patch are finally analyzed, as shown in Table 4. Discrete and strong clutter patches with a quantity ratio of 8.59% contribute 97.25% of the power, so it can be considered that the actual clutter scene has sparse characteristics, that is, the hypothesis in this paper is consistent with the actual situation.

5.2. Sub-Array Beam Optimization

In this experiment, the transmitting sub-array was evenly divided into eight non-overlapping sub-arrays, and each receiving sub-array contained eight elements. The mainlobe region was defined as ${\Theta }_{\mathrm{m}}=\left[-{70}^{\circ },{70}^{\circ }\right]$, while the sidelobe region was designated as ${\Theta }_{\mathrm{s}}=\left[-{90}^{\circ },-{70}^{\circ }\right]\cup \left[{70}^{\circ },{90}^{\circ }\right]$. The mainlobe’s upper and lower bounds were set at $U_{\mathrm{m}}=2$ dB and $L_{\mathrm{m}}=-2$ dB, respectively. The angle-sampling interval was established at $0.5^{\circ }$. As illustrated in Figure 11a, the beam optimized by the proposed method demonstrates a more even gain distribution across the designated mainlobe area when compared to the static transmitting sub-array beam pattern, which exhibits high gain solely in the beam’s direction. The normalized gain within the mainlobe is approximately −10 dB.

The receiving array was evenly divided into eight overlapping sub-arrays, and each receiving sub-array contained 57 elements. The mainlobe region and sidelobe region were also set to ${\Theta }_{\mathrm{m}}=\left[-{70}^{\circ },{70}^{\circ }\right]$ and ${\Theta }_{\mathrm{s}}=\left[-{90}^{\circ },-{70}^{\circ }\right]\cup \left[{70}^{\circ },{90}^{\circ }\right]$, and the angle sampling interval was set to ${0.5}^{\circ }$. Parameter ${\epsilon }_{\mathrm{m}}$ was set to 1. The simulation results are shown in Figure 11b. It can be seen from the simulation results that the beam optimized by the proposed method can also obtain relatively uniform gain within the expected mainlobe region, and the normalized beam gain is about −17.5 dB. In addition, we compared the optimization results of the proposed transmit sub-array beam-forming optimization method with those obtained using the alternating optimization (AO) algorithm [39] and the genetic algorithm (GA). The differences between the maximum and minimum gains within $\left[-{70}^{\circ },{70}^{\circ }\right]$ for the three methods were 4.02 dB, 6.62 dB, and 7.68 dB, respectively. The proposed method exhibited the smallest gain variation within the main lobe.

Similarly, we also compared the optimization results of the proposed receive sub-array beam-forming optimization method with those of the gradient descent (GD) algorithm [40] and GA. The differences between the maximum and minimum gains of the three algorithms in $\left[-{70}^{\circ },{70}^{\circ }\right]$ were 0.62 dB, 6.32 dB, and 14.26 dB, respectively. The proposed method exhibited the smallest gain fluctuation within the main lobe.

Figure 12 depicts the gain of the optimized sub-array joint beam across various angles. When contrasted with the static sub-array joint beam, there is a marked enhancement in the sidelobe region for the gain of the optimized sub-array joint beam. Upcoming simulation results will demonstrate that achieving uniform gain across the region of interest for both the transmit and receive sub-arrays is crucial for effective clutter sensing in a sub-array-level digital array radar system.

In actual radar systems, placement errors on the array are inevitable. The impact of placement errors can be reflected as phase errors superimposed on the spatial steering vector. As shown in Figure 13, we simulated the transmit and receive sub-array beam patterns under phase errors of $0^{\circ }$, $2^{\circ }$, $4^{\circ }$, $6^{\circ }$, $8^{\circ }$, and ${10}^{\circ }$. We can see that the placement errors will cause some disturbance to the sub-array beam pattern, but overall, the sub-array beam still maintains the flat-top characteristic. The experimental results show that when the placement error is within ${10}^{\circ }$, the sub-array beam will not change significantly.

5.3. Clutter-Sensing Performance

In this subsection, clutter sensing is performed using both a digital array at the sub-array level with a static sub-array joint beam and with an optimized sub-array joint beam, to validate the efficacy of the proposed methodology. As per Equation (20), the dictionary matrix is formulated with an angular interval of $0.{25}^{\circ }$, and the BPDN estimator is employed to address the SR problem outlined in Equation (23). Figure 14 shows the clutter-sensing results under different conditions and their mapping to the joint beam gain of the transmit–receive sub-array. Figure 14a shows the real clutter scene, while Figure 14b,c represent the clutter-sensing results and the joint transmit–receive sub-array beam pattern without and with sub-array beam optimization, respectively.

As shown in Figure 14b, when the sub-array beam is not optimized, the joint transmit–receive sub-array beam exhibits high gain only near the main lobe, while the low gain in the sidelobe region results in a low clutter-to-noise ratio (CNR) in the sidelobe area. This, in turn, leads to an inability to effectively reconstruct the clutter patches in the sidelobe region, allowing only the reconstruction of clutter patches near the main lobe. In Figure 14c, the optimized joint transmit–receive sub-array beam broadens the main lobe over a wide spatial area, ensuring high beam gain across a large spatial region. This improves the CNR in the sidelobe area. The comparison of Figure 14a,c further demonstrates that the proposed method, which optimizes the transmit–receive sub-array beam into a wide flat-topped beam, can effectively improve the reconstruction of clutter patches in the sidelobe region.

5.4. Clutter Suppression Performance

In this part, the 450th range cell is selected as the CUT, the color load factor is set to 0.8, and the CNR is set to 60 dB. The clutter spectrum distributions obtained by different algorithms are analyzed. Figure 15 shows the clutter spectrum estimated for the actual CNCM and loading sample matrix inversion (LSMI), and the clutter-sensing result with the static sub-array joint beam (CSSB) and the clutter-sensing result with the optimized sub-array beam (CSOB). It can be seen that the clutter spectra estimated by LSMI and the CSSB, as shown in Figure 15b,c have large differences with the actual clutter spectrum in the distribution region of sidelobe clutter, and then, the clutter spectrum estimated by the CSOB is close to the actual clutter spectrum.

The clutter suppression performance of different methods, which is defined as the ratio of the output signal-to-clutter-plus-noise ratio (SCNR) to the input SCNR, can be described by the improvement factor (IF):

(45)$\mathrm{IF}={\displaystyle \frac{{\left|{\mathbf{w}}^{\mathrm{H}}\mathbf{s}\right|}^2}{{\mathbf{w}}^{\mathrm{H}}{\mathbf{R}}_{\mathrm{actual}}\mathbf{w}}}{\displaystyle \frac{\mathrm{tr}\left({\mathbf{R}}_{\mathrm{actual}}\right)}{{\mathbf{s}}^{\mathrm{H}}\mathbf{s}}},$

where ${\mathbf{R}}_{\mathrm{actual}}$ is the actual CNCM of the CUT, $\mathbf{w}$ is the spatial–temporal adaptive weight vector, and $\mathbf{s}$ is the target signal.

When $\Delta y=0$, the IF curves obtained by the optimal processor, CSSB-STAP, CSOB-STAP, and LSMI-STAP are shown in Figure 16a. We can see that the clutter suppression performance of the CSOB-STAP method is close to that of the optimal processor. Due to the lack of valid samples, the performance of LSMI-STAP decreases by about 20 dB compared to the optimal processor. CSSB-STAP performs better than LSMI-STAP, but its performance loss is about 15dB greater than that of the optimal processor.

To ensure the effectiveness of the prior CNCM in the detection stage, the observation of the same clutter patch by the radar must be correlated in the sensing stage and the detection stage. However, due to the movement of the platform, the phase difference between the scattering sources inside the clutter patch will change, which will lead to a gradual decrease in the correlation between the observation of the same clutter patch by the radar in the sensing stage and the detection stage, that is, the prior knowledge described in this paper is time-sensitive.

Figure 16a–c shows the IF curves when the displacement of the platform is 0 m, 30 m, and 60 m, respectively. From the experimental results, it can be seen that the prior knowledge still matches the clutter environment even if there is some displacement of the carrier platform. Therefore, after the radar obtains the prior knowledge of the clutter environment by transmitting the sensing signal once, it can continue to use this prior knowledge to assist STAP in completing heterogeneous clutter suppression in the next few CPIs.

To further verify the target detection performance of the proposed algorithm, we added a target in the received echo. The target is positioned at the 385th range unit, with a radial velocity of 150 m/s. We simulated the clutter suppression results for the optimal processor, CSSB-STAP, CSOB-STAP, and LSMI-STAP. Figure 17 shows the filtering output for the Doppler channel where the target is located.

The experimental result shows the normalized peak sidelobe level (PSL) of the filtering outputs. The PSL values for the four filters are −13.75 dB, 1.9 dB, 23.65 dB, and 34.48 dB, respectively. It is evident that the normalized PSL of the proposed algorithm is significantly improved compared to LSMI-STAP and CSSB-STAP.

5.5. Algorithm Running Time Analysis

In addition, we calculated the running time of each main step of the proposed algorithm. The simulation software was MATLAB 2021b, the simulation platform was Windows 11 system (developed by Microsoft Corporation, Washington, DC, USA), the processor was a 12th Gen Intel(R) Core(TM) i5-12500H (developed by Intel Corporation, Santa Clara, CA, USA), and the number of Monte Carlo experiments was 200. As can be seen from Table 5, the construction of the prior CNCM takes the longest time in the main steps of the proposed algorithm. The main reason is that in order to accurately estimate the prior CNCM, we need to calculate the echo snapshots of each clutter patch. The computational time of the clutter-sensing stage is comparable to that of the space-time adaptive processing, and its main time consumption lies in solving the sparse recovery problem.

6. Conclusions

In this paper, we address the clutter-sensing STAP challenge for sub-array digital arrays within heterogeneous clutter environments. By examining the clutter signal model specific to sub-array digital array radars, we assess the impact of the conventional sub-array digital array radar’s low-sidelobe characteristics on clutter-sensing efficacy. Consequently, we formulate an optimization function targeting the enhancement of sub-array joint beams with flat-top characteristics. To simplify the complexity of the optimization variables, we divide the original optimization problem into two distinct subproblems, one for the transmitting sub-array beam optimization and another for the receiving sub-array beam optimization. The solutions to these sub-problems are then sought using the PDD algorithm and the interior-point method, respectively. It is important to note that in practical applications, the orthogonal waveforms emitted by each sub-array during the sensing phase may not be completely incoherent. As such, the transmitting waveforms can influence the radar beam to a certain degree. Recognizing this factor, it becomes necessary to consider an optimization model that jointly addresses the transmitting sub-array weight vector, the receiving sub-array weight vector, and the transmitting waveform. This is an avenue for future research.

Author Contributions

Conceptualization, B.J., W.P., K.L. and H.L.; Methodology, Y.W., W.P. and H.Z.; Software, Y.W.; Formal analysis, H.Z.; Data curation, B.J.; Writing—original draft, Y.W. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The data presented in this study are available on request from the corresponding author. Due to the privacy of the data, it is inconvenient to disclose.

Conflicts of Interest

The authors declare no conflicts of interest.

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. Geometric configuration of airborne radar.

Figure 2. Sub-array structure diagram. (a) Transmitting sub-array. (b) Receiving sub-array.

Figure 3. Signal processing flow of proposed method.

Figure 4. Sensing stage/detection stage timing diagram.

Figure 5. Sub-array joint beam pattern ([Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] m, [Forumla omitted. See PDF.] m, [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.]).

Figure 6. Geometrical configuration of airborne radar in sensing/detection stage.

Figure 7. SAR image of a certain area.

Figure 8. Relationship between clutter patch and SAR image pixels.

Figure 9. Actual clutter scene.

Figure 10. Box plot of RCS amplitude for all clutter patches.

Figure 11. Optimization result of sub-array beam. (a) Optimized transmitting sub-array beam pattern. (b) Optimized receiving sub-array beam pattern.

Figure 12. Optimized sub-array joint beam pattern.

Figure 13. Sub-array beam patterns under different phase errors. (a) Transmitting sub-array beam pattern under different phase errors. (b) Receiving sub-array beam pattern under different phase errors.

Figure 14. Comparison of clutter-sensing results. (a) Actual clutter scene. (b) Clutter-sensing result with static sub-array joint beam. (c) Clutter-sensing result with optimized sub-array joint beam.

Figure 15. Clutter spectrum of CUT. (a) Clutter spectrum estimated with actual CNCM. (b) Clutter spectrum estimated by LSMI. (c) Clutter spectrum estimated by CSSB. (d) Clutter spectrum estimated by CSOB.

Figure 16. IF curves at different platform displacements. (a) [Forumla omitted. See PDF.] m. (b) [Forumla omitted. See PDF.] m. (c) [Forumla omitted. See PDF.] m.

Figure 17. The filtering output for the Doppler channel where the target is located.

Table 1

Steps of proposed method.

Sensing Stage
step1:	Optimize the transmitting sub-array weight vector ${\mathbf{w}}_t$ and receiving sub-array weight vector ${\mathbf{w}}_r$ using (30) and (36);
step2:	The radar transmits orthogonal signals $\mathbf{p}\left(t\right)$;
step3:	Construct the dictionary matrix ${\mathbf{H}}_i$ for each range cell according to (20);
step4:	Estimate the RCS vector ${\widehat{\gamma }}_i$ of each range ring using (23);
step5:	Store the prior information of the clutter environment in the database.
Detection Stage
step1:	Update the slant range $r_{ij}$ and cone angle ${\phi }_{ij}$ of each clutter patch using (37);
step2:	The radar transmits LFM signals $u\left(t\right)$;
step3:	Construct the prior CNCM ${\mathbf{R}}_{\mathrm{prior},i}$ using (38);
step4:	Fuse the prior CNCM ${\mathbf{R}}_{\mathrm{prior},i}$ and SCM ${\mathbf{R}}_{\mathrm{smi},i}$ to obtain the color-loading CNCM ${\mathbf{R}}_{\mathrm{CL},i}$ using (42);
step5:	Calculate the optimal space-time adaptive processing weight vector ${\mathbf{w}}_{\mathrm{opt},i}$ using (43);
step6:	Calculate the clutter suppression output.

Table 2

The computational complexity of the main steps proposed.

Algorithm Steps			Computational Complexity
Sub-array beam optimization (offline processing part)	Transmitting sub-array beam optimization		$O\left(N_{sub1}^3+I_{in}·I_{out}max\left(N_{sub1}\left(L_m+L_s\right),N_{sub1}^2\right)\right)$
	Receiving beam sub-array optimization		$O\left(\left(N_{sub2}+1\right)\left(L_m+L_s\right)N_{sub2}^{4.5}log\left({\epsilon }^{-1}\right)\right)$
Clutter sensing and target detection (online processing part)	Sensing stage	Pulse compression	$O\left(N_1{N}_2{K}_s{L}_r{L}_t\right)$
		Clutter sensing	$O\left(\left(N_1{N}_2{K}_s{N}_c+N_c\right)N_c^{4.5}log\left({\epsilon }^{-1}\right)\right)$
	Detection stage	Pulse compression	$O\left(N_2{K}_d{L}_r{L}_t\right)$
		Construct prior CNCM (include location update)	$O\left(\left(2N_2{M}_d+L_t^2+3L_t+14\right)L_t{N}_c\right)$
		Space-time adaptive processing	$O\left(2N_2^3{K}_s^3+2N_2^2{K}_s^2\right)$

Table 3

Radar system parameters.

Symbols	Parameters	Value
$\lambda$	Wavelength	0.24 m
B	Bandwidth	2 MHz
N	Number of elements in whole array	64
$N_1$	Number of transmitting sub-arrays	8
$N_2$	Number of receiving sub-arrays	8
$N_{sub1}$	Number of elements in transmitting sub-array	8
$N_{sub2}$	Number of elements in receiving sub-array	57
d	Element spacing	0.12 m
$f_r$	Pulse repetition frequency	2500 Hz
$K_s$	Number of pulses in sensing stage	8
$K_d$	Number of pulses in detection stage	32
h	Platform height	5 km
$v_0$	Platform velocity	140 m/s

Table 4

Amplitude and number statistics of clutter patches.

	Quantity Ratio	Power Ratio
Common clutter patch	91.41%	2.75%
Discrete strong clutter patch	8.59%	97.25%

Table 5

Algorithm running time analysis.

Algorithm Steps			Running Time/s
Sub-array beam optimization (offline processing part)	Transmitting sub-array beam optimization		$5.87$
	Receiving beam sub-array optimization		$0.74$
Clutter sesning and target detection (online processing part)	Sensing stage	Pulse compression	$5.8\times {10}^{-3}$
		Clutter sensing	$5.95$
	Detection stage	Pulse compression	$2.8\times {10}^{-3}$
		Construct prior CNCM (include location update)	$13.58$
		Space-time adaptive processing	$5.12$

References

1. Ward, J. Space-time adaptive processing for airborne radar. IET Conf. Proc.; 1998; [DOI: https://dx.doi.org/10.1049/ic:19980240]

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

3. Melvin, W. A STAP overview. IEEE Aerosp. Electron. Syst. Mag.; 2004; 19, pp. 19-35. [DOI: https://dx.doi.org/10.1109/MAES.2004.1263229]

4. Reed, I.; Mallett, J.; Brennan, L. Rapid Convergence Rate in Adaptive Arrays. IEEE Trans. Aerosp. Electron. Syst.; 1974; AES-10, pp. 853-863. [DOI: https://dx.doi.org/10.1109/TAES.1974.307893]

5. Rabideau, D.; Steinhardt, A. Improved adaptive clutter cancellation through data-adaptive training. IEEE Trans. Aerosp. Electron. Syst.; 1999; 35, pp. 879-891. [DOI: https://dx.doi.org/10.1109/7.784058]

6. Kogon, S.; Zatman, M. STAP adaptive weight training using phase and power selection criteria. Proceedings of the Conference Record of Thirty-Fifth Asilomar Conference on Signals, Systems and Computers (Cat.No.01CH37256); Pacific Grove, CA, USA, 4–7 November 2001; Volume 1, pp. 98-102. [DOI: https://dx.doi.org/10.1109/ACSSC.2001.986887]

7. Melvin, W.; Wicks, M.; Brown, R. Assessment of multichannel airborne radar measurements for analysis and design of space-time processing architectures and algorithms. Proceedings of the 1996 IEEE National Radar Conference; Ann Arbor, MI, USA, 13–16 May 1996; pp. 130-135. [DOI: https://dx.doi.org/10.1109/NRC.1996.510669]

8. Guo, Q.; Liu, L.; Kaliuzhnyi, M.; Wang, Y.; Qi, L. STAP Training Samples Selection Based on GIP and Volume Cross Correlation. IEEE Geosci. Remote Sens. Lett.; 2022; 19, 4028205. [DOI: https://dx.doi.org/10.1109/LGRS.2022.3218670]

9. 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. Antennas Propag.; 2001; 49, pp. 91-103. [DOI: https://dx.doi.org/10.1109/8.910535]

10. Sun, K.; Zhang, H.; Li, G.; Meng, H.; Wang, X. A novel STAP algorithm using sparse recovery technique. Proceedings of the 2009 IEEE International Geoscience and Remote Sensing Symposium; Cape Town, South Africa, 12–17 July 2009; [DOI: https://dx.doi.org/10.1109/IGARSS.2009.5417664]

11. Yang, Z.; de Lamare, R.C.; Li, X. L₁-Regularized STAP Algorithms With a Generalized Sidelobe Canceler Architecture for Airborne Radar. IEEE Trans. Signal Process.; 2012; 60, pp. 674-686. [DOI: https://dx.doi.org/10.1109/TSP.2011.2172435]

12. Zhang, W.; An, R.; He, N.; He, Z.; Li, H. Reduced Dimension STAP Based on Sparse Recovery in Heterogeneous Clutter Environments. IEEE Trans. Aerosp. Electron. Syst.; 2020; 56, pp. 785-795. [DOI: https://dx.doi.org/10.1109/TAES.2019.2921141]

13. Cui, N.; Xing, K.; Yu, Z.; Duan, K. Tensor-Based Sparse Recovery Space-Time Adaptive Processing for Large Size Data Clutter Suppression in Airborne Radar. IEEE Trans. Aerosp. Electron. Syst.; 2023; 59, pp. 907-922. [DOI: https://dx.doi.org/10.1109/TAES.2022.3192223]

14. Wang, D.; Wang, T.; Cui, W.; Zhang, X. A Clutter Suppression Algorithm via Enhanced Sparse Bayesian Learning for Airborne Radar. IEEE Sens. J.; 2023; 23, pp. 10900-10911. [DOI: https://dx.doi.org/10.1109/JSEN.2023.3263919]

15. Liu, C.; Wang, T.; Zhang, S.; Ren, B. A Clutter Suppression Algorithm via Weighted ℓ₂-norm Penalty for Airborne Radar. IEEE Signal Process. Lett.; 2022; 29, pp. 1522-1525. [DOI: https://dx.doi.org/10.1109/LSP.2022.3187347]

16. Sun, K.; Meng, H.; Wang, Y.; Wang, X. Direct data domain STAP using sparse representation of clutter spectrum. Signal Process.; 2011; 91, pp. 2222-2236. [DOI: https://dx.doi.org/10.1016/j.sigpro.2011.04.006]

17. Song, D.; Chen, S.; Li, H.; Xi, F. Space-Time Adaptive Processing via Random Matrix Theory for Finite Training Samples. IEEE Sens. J.; 2023; 23, pp. 7334-7344. [DOI: https://dx.doi.org/10.1109/JSEN.2023.3245581]

18. Guerci, J.; Baranoski, E. Knowledge-aided adaptive radar at DARPA: An overview. IEEE Signal Process. Mag.; 2006; 23, pp. 41-50. [DOI: https://dx.doi.org/10.1109/MSP.2006.1593336]

19. Melvin, W.; Guerci, J. Knowledge-aided signal processing: A new paradigm for radar and other advanced sensors. IEEE Trans. Aerosp. Electron. Syst.; 2006; 42, pp. 983-996. [DOI: https://dx.doi.org/10.1109/TAES.2006.248215]

20. Capraro, C.; Capraro, G.; Bradaric, I.; Weiner, D.; Wicks, M.; Baldygo, W. Implementing digital terrain data in knowledge-aided space-time adaptive processing. IEEE Trans. Aerosp. Electron. Syst.; 2006; 42, pp. 1080-1099. [DOI: https://dx.doi.org/10.1109/TAES.2006.248199]

21. Xiong, Y.; Xie, W.; Li, H.; Gao, X. Colored-Loading Factor Optimization for Airborne KA-STAP Radar. IEEE Sens. J.; 2023; 23, pp. 23317-23326. [DOI: https://dx.doi.org/10.1109/JSEN.2023.3303264]

22. Tao, F.; Wang, T.; Wu, J.; Lin, X. A novel KA-STAP method based on Mahalanobis distance metric learning. Digit. Signal Process.; 2020; 97, 102613. [DOI: https://dx.doi.org/10.1016/j.dsp.2019.102613]

23. Kang, N.; Shang, Z.; Du, Q. Knowledge-Aided Structured Covariance Matrix Estimator Applied for Radar Sensor Signal Detection. Sensors; 2019; 19, 664. [DOI: https://dx.doi.org/10.3390/s19030664]

24. Gao, Z.; Tao, H. Robust STAP algorithm based on knowledge-aided SR for airborne radar. IET Radar Sonar Navig.; 2017; 11, pp. 321-329.

25. Adve, R.; Hale, T.; Wicks, M. Knowledge based adaptive processing for ground moving target indication. Digit. Signal Process.; 2007; 17, pp. 495-514. [DOI: https://dx.doi.org/10.1016/j.dsp.2005.06.005]

26. Xiong, Y.; Xie, W.; Wang, Y. Space time adaptive processing for airborne MIMO radar based on space time sampling matrix. Signal Process.; 2023; 211, 109119. [DOI: https://dx.doi.org/10.1016/j.sigpro.2023.109119]

27. Fang, M.; Liu, H.; Dai, F.; Wang, X. Space-time Adaptive Processing via Dynamic Environment Sensing. J. Electron. Inf. Technol.; 2015; 37, 1786. [DOI: https://dx.doi.org/10.11999/JEIT141505]

28. Wu, Y.; Jiu, B.; Guo, Z.; Liu, H. Space-time Adaptive Processing via Fast Environment Sensing. Proceedings of the 2022 IEEE International Conference on Signal Processing, Communications and Computing (ICSPCC); Xi’an, China, 25–27 October 2022; pp. 1-6. [DOI: https://dx.doi.org/10.1109/ICSPCC55723.2022.9984609]

29. Brookner, E. Phased-array and radar astounding breakthroughs—An update. Proceedings of the 2008 IEEE Radar Conference; Rome, Italy, 26–30 May 2008; pp. 1-6. [DOI: https://dx.doi.org/10.1109/RADAR.2008.4720771]

30. Lin, Z.; Hu, H.; Lei, S.; Li, R.; Tian, J.; Chen, B. Low-Sidelobe Shaped-Beam Pattern Synthesis With Amplitude Constraints. IEEE Trans. Antennas Propag.; 2022; 70, pp. 2717-2731. [DOI: https://dx.doi.org/10.1109/TAP.2021.3125319]

31. Ye, Y. Interior Point Algorithms: Theory and Analysis; John Wiley & Sons: Hoboken, NJ, USA, 2011.

32. Melvin, W. Space-time adaptive processing for radar. Academic Press Library in Signal Processing; Elsevier: Amsterdam, The Netherlands, 2014; Volume 2, pp. 595-665. [DOI: https://dx.doi.org/10.1016/B978-0-12-396500-4.00012-0]

33. Donoho, D. Compressed sensing. IEEE Trans. Inf. Theory; 2006; 52, pp. 1289-1306. [DOI: https://dx.doi.org/10.1109/TIT.2006.871582]

34. Tropp, J.; Gilbert, A. Signal Recovery From Random Measurements via Orthogonal Matching Pursuit. IEEE Trans. Inf. Theory; 2007; 53, pp. 4655-4666. [DOI: https://dx.doi.org/10.1109/TIT.2007.909108]

35. Chen, S.; Donoho, D.; Saunders, M. Atomic decomposition by basis pursuit. SIAM Rev.; 2001; 43, pp. 129-159. [DOI: https://dx.doi.org/10.1137/S003614450037906X]

36. Melman, A.; Polyak, R. The Newton modified barrier method for QP problems. Ann. Oper. Res.; 1996; 62, pp. 465-519. [DOI: https://dx.doi.org/10.1007/BF02206827]

37. Luo, Z.Q.; Ma, W.K.; So, A.M.C.; Ye, Y.; Zhang, S. Semidefinite Relaxation of Quadratic Optimization Problems. IEEE Signal Process. Mag.; 2010; 27, pp. 20-34. [DOI: https://dx.doi.org/10.1109/MSP.2010.936019]

38. Page, D.; Scarborough, S.; Crooks, S. Improving knowledge-aided STAP performance using past CPI data [radar signal processing]. Proceedings of the 2004 IEEE Radar Conference; Philadelphia, PA, USA, 29 April 2004; pp. 295-300. [DOI: https://dx.doi.org/10.1109/NRC.2004.1316438]

39. Cao, P.; Thompson, J.S.; Haas, H. Constant Modulus Shaped Beam Synthesis via Convex Relaxation. IEEE Antennas Wirel. Propag. Lett.; 2017; 16, pp. 617-620. [DOI: https://dx.doi.org/10.1109/LAWP.2016.2594141]

40. Ruder, S. An overview of gradient descent optimization algorithms. arXiv; 2016; arXiv: 1609.04747