Low-Altitude Windshear Wind Speed Estimation Method Based on KASPICE-STAP

Abstract

Aiming at the problem of low-altitude windshear wind speed estimation for airborne weather radar without independent identically distributed (IID) training samples, this paper proposes a low-altitude windshear wind speed estimation method based on knowledge-aided sparse iterative covariance-based estimation STAP (KASPICE-STAP). Firstly, a clutter dictionary composed of clutter space–time steering vectors is constructed using prior knowledge of the distribution position of ground clutter echo signals in the space–time spectrum. Secondly, the SPICE algorithm is used to obtain the clutter covariance matrix iteratively. Finally, the STAP processor is designed to eliminate the ground clutter echo signal, and the wind speed is estimated after eliminating the ground clutter echo signal. The simulation results show that the proposed method can accurately realize a low-altitude windshear wind speed estimation without IID training samples.

1. Introduction

Low-altitude windshear is an atmospheric phenomenon in which the air flow below 600 m suddenly changes its direction or speed within a small area. It has the characteristics of a short duration, a small area of action, a high instantaneous intensity, and a strong potential to cause harm [1]. When an aircraft suddenly encounters low-altitude windshear in the takeoff or landing stage, the pilot usually lacks sufficient time and space to control the attitude of the aircraft, which leads to the occurrence of flight accidents. Therefore, the detection and warning of low-altitude windshear has become an important topic in the current civil aviation field [2]. The estimation of low-altitude windshear wind speed determines the accuracy of low-altitude windshear detection, so the accurate estimation of low-altitude windshear wind speed is of much concern. As airborne weather radar operates in look-down mode, the clutter has a wide distribution range and a high intensity, and the aircraft’s movement will cause the clutter spectrum to expand seriously, which will lead to the low-altitude windshear echo signal in the airborne weather radar echo signal often being submerged in the ground clutter echo signal, so the suppression effect of the ground clutter echo signal will directly affect the accuracy of the wind speed estimation [3].

Space-–time adaptive processing (STAP) technology [4] is a signal-processing technology that is applied to airborne weather radar. The application of STAP technology makes it possible to accurately estimate the low-altitude windshear wind speed in strong clutter environments. The optimal weight vector of STAP relies on the accurate estimation of the clutter covariance matrix (CCM) [5]. According to the RMB (Reed–Mallet–Brennan) criterion proposed by Reed et al. [6], the used training samples must be independent identically distributed (IID) samples, and only when the number of IID training samples reaches two times the freedom degree of the system or more can the performance of the traditional STAP method be reduced by less than 3 dB. However, in practice, the lack of IID training samples makes the traditional STAP technology unable to accurately obtain the clutter covariance matrix, which leads to the degradation of the ground clutter echo signal suppression performance, and this limits the application of the traditional STAP technology in low-altitude windshear wind-speed estimation. Although the rank-reduction STAP [7] and dimensionality-reduction STAP [8], which were proposed to estimate the low-altitude wind speed of windshear when the number of IID training samples is low, can effectively suppress the ground clutter echo signal and estimate the wind speed, they still need a certain number of IID training samples.

In recent years, scholars have proposed a variety of STAP methods that do not require IID training samples. In 2009, Sarkar et al. [9] proposed the direct data domain STAP (DDD-STAP) method. This method directly determines the window for the detection range units without IID training samples, which can effectively suppress clutter and achieve the purpose of target detection. However, the sliding window involve the loss of the freedom degree of the system, and the sliding window has high requirements on the degrees of freedom in the space and time domains. The development of knowledge-aided STAP (KA-STAP) technology brings a new round of innovation to STAP technology. Although KA-STAP technology can improve the clutter suppression performance to a certain extent, it requires a higher precision of prior knowledge. On this basis, in 2011, Sun Ke et al. from Tsinghua University introduced a direct data domain and prior knowledge into the sparse STAP method and proposed the DDD sparse recovery STAP (DDD-SR-STAP) method [10]. In this method, the detection range units were directly sparse recovered without a sliding window, but this requires prior knowledge of the specific position distribution of the echo signal in the target area in the two-dimensional space–time spectrum. Since the low-altitude windshear is a distributed target and the low-altitude windshear has the characteristics of a sudden occurrence and a short existence time, it is very difficult to know the specific position distribution of low-altitude windshear echo signals in a two-dimensional space–time spectrum in advance.

STAP methods that require few or no IID samples have been applied in clutter suppression and parameter estimation in many ways [11,12,13,14,15,16,17,18], but they are all aimed at point targets and cannot be directly applied to distributed targets such as low-altitude windshear. Therefore, it is of great significance to study how to accurately realize low-altitude windshear wind speed estimation without IID training samples.

Aiming at solving the above problems, this paper proposes a low-altitude windshear wind speed estimation method based on knowledge-aided sparse iterative covariance-based estimation STAP (KASPICE-STAP). Firstly, a clutter dictionary composed of clutter space–time steering vectors is constructed by using prior knowledge of the distribution position of the space–time spectrum of the ground clutter echo signals. Secondly, the SPICE algorithm is used to obtain the clutter covariance matrix iteratively. Finally, the STAP processor is designed to calculate the corresponding optimal weight vector for clutter suppression, the normalized Doppler frequency is obtained from the signal after clutter suppression, and then, the wind speed is estimated. Simulation results show that the proposed method can effectively and accurately estimate wind speed without IID training samples.

2. Echo Signal Model

The model of low-altitude windshear detected by an airborne weather radar forward-looking array is shown in Figure 1. The reference coordinate system is established with H meters directly below the particle point of the aircraft platform as the origin. It is assumed that the antenna array is composed of an N-element linear array, and it is placed evenly along the Y-axis with a spacing of $\mathrm{d}=0.5 \mathsf{\lambda }$, where $\mathsf{\lambda }$ is the wavelength. We assume that the pulse repetition rate is ${\mathrm{f}}_{\mathrm{r}}$ and the number of pulses within a coherent processing interval is K. In Figure 1, ${\mathsf{\theta }}_0$ is the horizontal azimuth angle of low-altitude windshear wind field echo signal, ${\mathsf{\phi }}_{\mathrm{l}}$ is its pitch angle, ${\mathsf{\psi }}_0$ is its spatial cone angle, and $\cos {\mathsf{\psi }}_0=\cos {\mathsf{\phi }}_{\mathrm{l}}\cos {\mathsf{\theta }}_0$ is valid.

This paper assumes that there are L range units in total; ${\mathrm{x}}_{\mathrm{l}}$ represents the $\mathrm{N}\mathrm{K}\times 1$-dimensional radar echo signal of the lth range unit to be detected ($\mathrm{l}=1,2,\dots \dots ,\mathrm{L}$), and its expression is:

$${\mathrm{x}}_{\mathrm{l}}={\mathrm{s}}_{\mathrm{l}}+{\mathrm{c}}_{\mathrm{l}}+{\mathrm{N}}_{\mathrm{l}}$$

(1)

where, ${\mathrm{s}}_{\mathrm{l}}$ is the low-altitude windshear wind field echo signal in the lth range unit to be detected; ${\mathrm{c}}_{\mathrm{l}}$ is the ground clutter echo signal in the lth range unit to be detected; ${\mathrm{n}}_{\mathrm{l}}$ is additive white Gaussian noise.

2.1. Low-Altitude Windshear Echo Signal Model

The echo signal ${\mathrm{s}}_{\mathrm{l}}$ of low-altitude windshear wind field in the lth range unit to be detected is:

$${\mathrm{s}}_{\mathrm{l}}=\mathsf{\sigma }.{\mathrm{v}}_{\mathrm{l}}{\left({\mathrm{f}}_{\mathrm{d}},{\mathsf{\psi }}_0\right)}_{\mathrm{N}\mathrm{K}\times 1}=\mathsf{\sigma }.{\mathrm{v}}_{\mathrm{d}}{\left({\mathrm{f}}_{\mathrm{d}}\right)}_{\mathrm{K}\times 1}\otimes {\mathrm{v}}_{\mathrm{s}}{\left({\mathsf{\psi }}_0\right)}_{\mathrm{N}\times 1}$$

(2)

where $\otimes$ is the Kronecker product; $\mathsf{\sigma }$ is the complex amplitude of the low-altitude windshear wind field echo signal of the lth range unit to be detected; ${\mathsf{\psi }}_0$ is the spatial cone angle of low-altitude windshear wind field echo signal. In this paper, it is considered that the spatial cone angle of low-altitude windshear wind field echo signal is prior information. ${\mathrm{f}}_{\mathrm{d}}$ is the normalized Doppler frequency of the wind field echo signal in the lth range unit to be detected, and the value range is $[-1,1]$; ${\mathrm{v}}_{\mathrm{l}}{\left({\mathrm{f}}_{\mathrm{d}},{\mathsf{\psi }}_0\right)}_{\mathrm{N}\mathrm{K}\times 1}$ is the space–time steering vector of wind field echo signal in the lth range unit to be detected; ${\mathrm{v}}_{\mathrm{d}}{\left({\mathrm{f}}_{\mathrm{d}}\right)}_{\mathrm{K}\times 1}$ and ${\mathrm{v}}_{\mathrm{s}}{\left({\mathsf{\psi }}_0\right)}_{\mathrm{N}\times 1}$ are, respectively, the time steering vector and space steering vector, constituting the space–time steering vector of the wind field echo signal in the range unit to be detected, and the expressions are as follows:

$${\mathrm{v}}_{\mathrm{d}}{\left({\mathrm{f}}_{\mathrm{d}}\right)}_{\mathrm{K}\times 1}={[1,{\mathrm{e}}^{\mathrm{J}\mathsf{\pi }{\mathrm{f}}_{\mathrm{d}}},\dots \dots ,{\mathrm{e}}^{\mathrm{J}\mathsf{\pi }\left(\mathrm{K}-1\right){\mathrm{f}}_{\mathrm{d}}}]}_{\mathrm{K}\times 1}^{\mathrm{T}}\odot \mathrm{g}({\mathsf{\sigma }}_{{\mathrm{f}}_{\mathrm{d}}})$$

(3)

$$\begin{array}{c}{\mathrm{v}}_{\mathrm{s}}{\left({\mathsf{\psi }}_0\right)}_{\mathrm{N}\times 1}={[1,{\mathrm{e}}^{\mathrm{J}2\mathsf{\pi }\frac{\mathrm{d}}{\mathsf{\lambda }}\cos {\mathsf{\theta }}_0\cos {\mathsf{\phi }}_{\mathrm{l}}},\dots \dots ,{\mathrm{e}}^{\mathrm{J}2\mathsf{\pi }\left(\mathrm{N}-1\right)\frac{\mathrm{d}}{\mathsf{\lambda }}\cos {\mathsf{\theta }}_0\cos {\mathsf{\phi }}_{\mathrm{l}}}]}_{\mathrm{N}\times 1}^{\mathrm{T}}\odot \mathrm{g}({\mathsf{\sigma }}_{{\mathsf{\psi }}_0})=\\ {[1,{\mathrm{e}}^{\mathrm{J}2\mathsf{\pi }\frac{\mathrm{d}}{\mathsf{\lambda }}\cos {\mathsf{\psi }}_0},\dots \dots ,{\mathrm{e}}^{\mathrm{J}2\mathsf{\pi }\left(\mathrm{N}-1\right)\frac{\mathrm{d}}{\mathsf{\lambda }}\cos {\mathsf{\psi }}_0}]}_{\mathrm{N}\times 1}^{\mathrm{T}}\odot \mathrm{g}({\mathsf{\sigma }}_{{\mathsf{\psi }}_0})\end{array}$$

(4)

In order to accurately describe the distributed model such as low-altitude windshear, Meng et al. [19] gradually established and improved the spatial distribution source model. By introducing an extension function to describe the angle and frequency expansion of the signal, the model extends the space–time steering vector of the point target to obtain the space–time steering vector of the distributed target. In Equations (3) and (4), $\odot$ is the Hadamard product; $\mathrm{g}({\mathsf{\sigma }}_{{\mathrm{f}}_{\mathrm{d}}})={[1,{\mathrm{e}}^{-2{\mathsf{\pi }}^2{\mathsf{\sigma }}_{{\mathrm{f}}_{\mathrm{d}}}},\dots \dots ,{\mathrm{e}}^{-2{\mathsf{\pi }}^2(\mathrm{K}-1){\mathsf{\sigma }}_{{\mathrm{f}}_{\mathrm{d}}}}]}^{\mathrm{T}}$ is the frequency expansion function of the low-altitude windshear wind field echo signal in the range unit to be detected; $\mathrm{g}({\mathsf{\sigma }}_{{\mathsf{\psi }}_0})={[1,{\mathrm{e}}^{-1/2{(2\mathsf{\pi }\frac{\mathrm{d}}{\mathsf{\lambda }})}^2{\mathsf{\sigma }}_{{\mathsf{\psi }}_0}},\dots \dots ,{\mathrm{e}}^{-1/2{(2\mathsf{\pi }(\mathrm{N}-1)\frac{\mathrm{d}}{\mathsf{\lambda }})}^2{\mathsf{\sigma }}_{{\mathsf{\psi }}_0}}]}^{\mathrm{T}}$ is the angle expansion function of the low-altitude windshear wind field echo signal in the range unit to be detected;${\mathsf{\sigma }}_{{\mathsf{\psi }}_0}={\mathsf{\sigma }}_{{\mathsf{\theta }}_0}^2{\cos }^2{\mathsf{\phi }}_{\mathrm{l}}{\sin }^2{\mathsf{\theta }}_0+{\mathsf{\sigma }}_{{\mathsf{\phi }}_{\mathrm{l}}}^2{\sin }^2{\mathsf{\phi }}_{\mathrm{l}}{\cos }^2{\mathsf{\theta }}_0$, ${\mathsf{\sigma }}_{{\mathsf{\theta }}_0}$ represents the extension of ${\mathsf{\theta }}_0$ in the horizontal azimuth angle of the low-altitude windshear wind field echo signal; ${\mathsf{\sigma }}_{{\mathsf{\phi }}_{\mathrm{l}}}$ represents the extension of ${\mathsf{\phi }}_{\mathrm{l}}$ in the direction of its pitch angle [20].

2.2. Ground Clutter Echo Signal Model

The ground clutter echo signal model of the airborne weather radar in this paper is established based on the Ward model [21] which was developed by Lincoln Laboratory of Massachusetts Institute of Technology. The ground clutter in the lth range unit to be detected is divided into M ground clutter patches, the mth ground clutter patch is shown in Figure 1. It is assumed that the ground clutter echo signal has no fluctuation and no fuzzy phenomenon [22], and the ground clutter echo signal ${\mathrm{c}}_{\mathrm{l}}$ in the lth range unit to be detected is:

$${\mathrm{c}}_{\mathrm{l}}={\displaystyle \sum _{\mathrm{m}=1}^{\mathrm{m}}\frac{{\mathrm{B}}_{\mathrm{m}}}{{\mathrm{R}}_{\mathrm{l}}^2}\mathrm{f}\left({\mathsf{\theta }}_{\mathrm{l},\mathrm{m}}\right)\mathrm{a}({\mathrm{f}}_{\mathrm{d}\mathrm{m}},{\mathsf{\psi }}_{\mathrm{s}\mathrm{m}})}$$

(5)

where ${\mathrm{B}}_{\mathrm{m}}$ is the reflection coefficient of the mth ground clutter patch; $\mathrm{F}\left({\mathsf{\theta }}_{\mathrm{l},\mathrm{m}}\right)$ is the radar antenna pattern; ${\mathrm{r}}_{\mathrm{l}}$ is the distance between the ground clutter patches of the lth range unit to be detected and the aircraft platform; $\mathrm{a}({\mathrm{f}}_{\mathrm{d},\mathrm{m}},{\mathsf{\psi }}_{\mathrm{s},\mathrm{m}})$ is the $\mathrm{N}\mathrm{K}\times 1$-dimensional space–time steering vector of the mth ground clutter patch. The expression is as follows:

$$\mathrm{a}({\mathrm{f}}_{\mathrm{dm}},{\mathsf{\psi }}_{\mathrm{sm}})={[1,{\mathrm{e}}^{\mathrm{J}\mathsf{\pi }{\mathrm{f}}_{\mathrm{d}\mathrm{m}}},\dots \dots ,{\mathrm{e}}^{\mathrm{J}\mathsf{\pi }(\mathrm{K}-1){\mathrm{f}}_{\mathrm{d}\mathrm{m}}}]}_{\mathrm{K}\times 1}^{\mathrm{T}}\otimes {[1,{\mathrm{e}}^{\mathrm{J}2\mathsf{\pi }\frac{\mathrm{d}}{\mathsf{\lambda }}\cos {\mathsf{\psi }}_{\mathrm{s}\mathrm{m}}},\dots \dots ,{\mathrm{e}}^{\mathrm{J}2\mathsf{\pi }(\mathrm{N}-1)\frac{\mathrm{d}}{\mathsf{\lambda }}\cos {\mathsf{\psi }}_{\mathrm{s}\mathrm{m}}}]}_{\mathrm{N}\times 1}^{\mathrm{T}}$$

(6)

where ${\mathsf{\psi }}_{\mathrm{sm}}$ is the corresponding spatial cone angle; ${\mathrm{f}}_{\mathrm{dm}}$ is the normalized Doppler frequency; when the pitch angle and horizontal azimuth angle of the clutter patch in this area are ${\mathsf{\phi }}_{\mathrm{l},\mathrm{m}}$ and ${\mathsf{\theta }}_{\mathrm{l},\mathrm{m}}$, respectively, then the spatial cone angle and the normalized Doppler frequency can be expressed as:

$${\mathsf{\psi }}_{\mathrm{s},\mathrm{m}}=\arccos (\cos {\mathsf{\theta }}_{\mathrm{l},\mathrm{m}}\cos {\mathsf{\phi }}_{\mathrm{l},\mathrm{c}})$$

(7)

$${\mathrm{f}}_{\mathrm{d},\mathrm{m}}=\frac{4\mathrm{V}}{\mathsf{\lambda }{\mathrm{f}}_{\mathrm{r}}}\cos (\mathsf{\alpha }-{\mathsf{\theta }}_{\mathrm{l},\mathrm{m}})\cos {\mathsf{\phi }}_{\mathrm{l},\mathrm{c}}$$

(8)

where $\mathsf{\alpha }$ is the included angle between the flight direction of the aircraft and the antenna array, and it is 90°.

3. Low-Altitude Windshear Wind Speed Estimation Using KASPICE-STAP Method

The estimation method of low-altitude windshear wind speed based on KASPICE-STAP is mainly divided into two parts: a clutter covariance matrix estimation based on SPICE algorithm and a low-altitude windshear wind speed estimation. The following are discussed, respectively.

3.1. Estimation of Clutter Covariance Matrix Based on SPICE Algorithm

To estimate the clutter covariance matrix by the SPICE algorithm, first of all, the clutter dictionary is constructed by using prior knowledge of the distribution position of the ground clutter echo signals in space–time two-dimensional spectrum, and then, the clutter covariance matrix is estimated by the SPICE algorithm. In the following two aspects, the construction of the clutter dictionary and the estimation of the clutter covariance matrix by the SPICE algorithm are described.

3.1.1. Construction of Clutter Dictionary

The position distribution of the ground clutter in space–time two-dimensional spectrum is shown in Figure 2. In the figure, the gray cell grids are the distribution position of the ground clutter echo signals, and the white cell grids are the distribution position of noise and low-altitude windshear echo signals. Since it is impossible to divide the clutter infinitely as in practical applications, a discrete approach is adopted to approximate the division, and each clutter cell grid is evenly divided into ${\mathrm{N}}_{\mathrm{s}}$ parts. Then, the $\mathrm{N}\mathrm{K}\times 1$-dimensional space–time steering vector corresponding to each clutter cell grid is $\mathrm{a}({\mathrm{f}}_{\mathrm{dm}},{\mathsf{\psi }}_{\mathrm{sm}})={\mathrm{s}}_{\mathrm{dm}}({\mathrm{f}}_{\mathrm{dm}})\otimes {\mathrm{s}}_{\mathrm{sm}}({\mathsf{\psi }}_{\mathrm{sm}}),\mathrm{m}=1,2,\dots ,{\mathrm{N}}_{\mathrm{s}}$, where ${\mathrm{s}}_{\mathrm{dm}}({\mathrm{f}}_{\mathrm{dm}})=[1,{\mathrm{e}}^{\mathrm{J}\mathsf{\pi }{\mathrm{f}}_{\mathrm{dm}}},\dots ,{\mathrm{e}}^{\mathrm{J}\mathsf{\pi }(\mathrm{K}-1){\mathrm{f}}_{\mathrm{dm}}}]$ is the time steering vector and ${\mathrm{s}}_{\mathrm{sm}}({\mathsf{\psi }}_{\mathrm{sm}})=[1,{\mathrm{e}}^{\mathrm{J}2\mathsf{\pi }\frac{\mathrm{d}}{\mathsf{\lambda }}{\mathsf{\psi }}_{\mathrm{dm}}},\dots ,{\mathrm{e}}^{\mathrm{J}2\mathsf{\pi }(\mathrm{N}-1)\frac{\mathrm{d}}{\mathsf{\lambda }}{\mathsf{\psi }}_{\mathrm{dm}}}]$ is the space steering vector. Then, the space–time steering vector of each clutter cell grid is combined into an $\mathrm{N}\mathrm{K}\times {\mathrm{N}}_{\mathrm{s}}$ dimensional clutter dictionary ${\mathrm{A}}_1$, which is expressed as:

$${\mathrm{A}}_1=[\mathrm{a}({\mathrm{f}}_{\mathrm{d}1},{\mathsf{\psi }}_{\mathrm{s}1}),\dots ,\mathrm{a}({\mathrm{f}}_{\mathrm{dm}},{\mathsf{\psi }}_{\mathrm{sm}}),\dots ,\mathrm{a}({\mathrm{f}}_{\mathrm{d}{\mathrm{N}}_{\mathrm{s}}},{\mathsf{\psi }}_{\mathrm{s}{\mathrm{N}}_{\mathrm{s}}})]$$

(9)

3.1.2. Estimation of Clutter Covariance Matrix Using SPICE Algorithm

If we let ${\mathrm{R}}_{\mathrm{l}}$ be the clutter covariance matrix of the lth range unit that is to be detected, its expression is [23]:

$${\mathrm{R}}_{\mathrm{l}}=\mathrm{E}\left[\left({\mathrm{c}}_{\mathrm{l}}+{\mathrm{n}}_{\mathrm{l}}\right){\left({\mathrm{c}}_{\mathrm{l}}+{\mathrm{n}}_{\mathrm{l}}\right)}^{\mathrm{H}}\right]$$

(10)

As it can be seen from Figure 2, since a single clutter cell grid is evenly divided into ${\mathrm{N}}_{\mathrm{s}}$ parts, $\mathrm{M}$ in Equation (5) is ${\mathrm{N}}_{\mathrm{s}}$, and ${\mathsf{\alpha }}_{\mathrm{l},\mathrm{m}}=\frac{{\mathrm{B}}_{\mathrm{m}}}{{\mathrm{R}}_{\mathrm{l}}^2}\mathrm{F}\left({\mathsf{\theta }}_{\mathrm{l},\mathrm{m}}\right)$ is the amplitude value of the ground clutter echo signal of the lth range unit that is to be detected and the mth clutter cell grid. In this case, ${\mathrm{c}}_{\mathrm{l}}$ in Equation (5) can be expressed as:

$${\mathrm{c}}_{\mathrm{l}}={\displaystyle \sum _{\mathrm{m}=1}^{{\mathrm{N}}_{\mathrm{s}}}{\mathsf{\alpha }}_{\mathrm{l},\mathrm{m}}.\mathrm{a}({\mathrm{f}}_{\mathrm{dm}},{\mathsf{\psi }}_{\mathrm{sm}})}$$

(11)

By substituting Equation (11) into Equation (10), we can obtain:

$$\begin{array}{l}{\mathrm{R}}_{\mathrm{l}}=\mathrm{E}\left[\left({\mathrm{c}}_{\mathrm{l}}+{\mathrm{n}}_{\mathrm{l}}\right){\left({\mathrm{c}}_{\mathrm{l}}+{\mathrm{n}}_{\mathrm{l}}\right)}^{\mathrm{H}}\right]\\ =\mathrm{E}[\left({\displaystyle \sum _{\mathrm{m}=1}^{{\mathrm{N}}_{\mathrm{s}}}{\mathsf{\alpha }}_{\mathrm{l}.\mathrm{m}}.\mathrm{a}({\mathrm{f}}_{\mathrm{dm}},{\mathsf{\psi }}_{\mathrm{sm}})+{\mathrm{n}}_{\mathrm{l}}}\right){\left({\displaystyle \sum _{\mathrm{m}=1}^{{\mathrm{N}}_{\mathrm{s}}}{\mathsf{\alpha }}_{\mathrm{l},\mathrm{m}}.\mathrm{a}({\mathrm{f}}_{\mathrm{d}\mathrm{m}},{\mathsf{\psi }}_{\mathrm{sm}})+{\mathrm{n}}_{\mathrm{l}}}\right)}^{\mathrm{H}})]\\ =\mathrm{E}\left[{\displaystyle \sum _{\mathrm{m}=1}^{{\mathrm{N}}_{\mathrm{s}}}{\left|{\mathsf{\alpha }}_{\mathrm{l},\mathrm{m}}\right|}^2}\mathrm{a}({\mathrm{f}}_{\mathrm{dm}},{\mathsf{\psi }}_{\mathrm{sm}}){(\mathrm{a}({\mathrm{f}}_{\mathrm{dm}},{\mathsf{\psi }}_{\mathrm{sm}}))}^{\mathrm{H}}\right]\\ +\mathrm{E}\left[{\mathrm{n}}_{\mathrm{l}}^{\mathrm{H}}{\displaystyle \sum _{\mathrm{m}=1}^{{\mathrm{N}}_{\mathrm{s}}}{\mathsf{\alpha }}_{\mathrm{l},\mathrm{m}}.\mathrm{a}({\mathrm{f}}_{\mathrm{dm}},{\mathsf{\psi }}_{\mathrm{sm}})}\right]+\mathrm{E}\left[{\mathrm{n}}_{\mathrm{l}}{\displaystyle \sum _{\mathrm{m}=1}^{{\mathrm{N}}_{\mathrm{s}}}{\mathsf{\alpha }}_{\mathrm{l},\mathrm{m}}.{(\mathrm{a}({\mathrm{f}}_{\mathrm{dm}},{\mathsf{\psi }}_{\mathrm{sm}}))}^{\mathrm{H}}}\right]+\mathrm{E}\left[{\mathrm{n}}_{\mathrm{l}}{\mathrm{n}}_{\mathrm{l}}^{\mathrm{H}}\right]\end{array}$$

(12)

It is assumed that since clutter and noise are independent of each other [24], $\mathrm{E}\left[{\mathrm{n}}_{\mathrm{l}}^{\mathrm{H}}{\displaystyle \sum _{\mathrm{m}=1}^{{\mathrm{N}}_{\mathrm{s}}}{\mathsf{\alpha }}_{\mathrm{l},\mathrm{m}}.\mathrm{a}({\mathrm{f}}_{\mathrm{dm}},{\mathsf{\psi }}_{\mathrm{sm}})}\right]=0$, $\mathrm{E}\left[{\mathrm{n}}_{\mathrm{l}}{\displaystyle \sum _{\mathrm{m}=1}^{{\mathrm{N}}_{\mathrm{s}}}{\mathsf{\alpha }}_{\mathrm{l},\mathrm{m}}.{(\mathrm{a}({\mathrm{f}}_{\mathrm{dm}},{\mathsf{\psi }}_{\mathrm{sm}}))}^{\mathrm{H}}}\right]=0$ in Equation (12) can be obtained by substituting it into Equation (12):

$${\mathrm{R}}_{\mathrm{l}}=\mathrm{E}\left[{\displaystyle \sum _{\mathrm{m}=1}^{{\mathrm{N}}_{\mathrm{s}}}{\left|{\mathsf{\alpha }}_{\mathrm{l},\mathrm{m}}\right|}^2}\mathrm{a}({\mathrm{f}}_{\mathrm{dm}},{\mathsf{\psi }}_{\mathrm{sm}}){(\mathrm{a}({\mathrm{f}}_{\mathrm{dm}},{\mathsf{\psi }}_{\mathrm{sm}}))}^{\mathrm{H}}\right]+\mathrm{E}\left[{\mathrm{n}}_{\mathrm{l}}{\mathrm{n}}_{\mathrm{l}}^{\mathrm{H}}\right]$$

(13)

In Equation (13), $\mathrm{E}\left[{\mathrm{n}}_{\mathrm{l}}{\mathrm{n}}_{\mathrm{l}}^{\mathrm{H}}\right]=\left[\begin{array}{cccc}{\mathsf{\sigma }}_1& 0& \cdots & 0\\ 0& {\mathsf{\sigma }}_2& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& \cdots & \cdots & {\mathsf{\sigma }}_{\mathrm{n}}\end{array}\right]={\mathsf{\sigma }}_{\mathrm{N}}{\mathrm{I}}_{\mathrm{N}\mathrm{K}\times \mathrm{N}\mathrm{K}}$. Therefore, Equation (13) can be expressed as:

$${\mathrm{R}}_{\mathrm{l}}={\displaystyle \sum _{\mathrm{m}=1}^{{\mathrm{N}}_{\mathrm{s}}}{\left|{\mathsf{\alpha }}_{\mathrm{l},\mathrm{m}}\right|}^2}\mathrm{a}({\mathrm{f}}_{\mathrm{dm}},{\mathsf{\psi }}_{\mathrm{sm}}){(\mathrm{a}({\mathrm{f}}_{\mathrm{dm}},{\mathsf{\psi }}_{\mathrm{sm}}))}^{\mathrm{H}}+{\mathsf{\sigma }}_{\mathrm{n}}{\mathrm{I}}_{\mathrm{N}\mathrm{K}\times \mathrm{N}\mathrm{K}}$$

(14)

where ${\left|{\mathsf{\alpha }}_{\mathrm{l},\mathrm{m}}\right|}^2$ is called the sparse coefficient vector, which represents the power value of the ground clutter echo signal of the lth range unit that is to be detected and the mth clutter cell grid. ${\mathsf{\sigma }}_{\mathrm{N}}$ is the receiver noise power, which is regarded as a known quantity.

In this case, the clutter covariance matrix can be estimated through Equation (14). In Equation (14), ${\left|{\mathsf{\alpha }}_{\mathrm{l},\mathrm{m}}\right|}^2$ is an unknown quantity, while in this paper, to solve the value of ${\left|{\mathsf{\alpha }}_{\mathrm{l},\mathrm{m}}\right|}^2$, the SPICE algorithm is used, and its objective function and constraint condition are as follows:

$$\{\begin{array}{c}\min \hspace{1em}\mathrm{J}={\Vert ({\displaystyle \sum _{\mathrm{m}=1}^{{\mathrm{N}}_{\mathrm{s}}}{\left|{\mathsf{\alpha }}_{\mathrm{l},\mathrm{m}}\right|}^2}\mathrm{a}({\mathrm{f}}_{\mathrm{dm}},{\mathsf{\psi }}_{\mathrm{s}\mathrm{m}}){(\mathrm{a}({\mathrm{f}}_{\mathrm{dm}},{\mathsf{\psi }}_{\mathrm{s}\mathrm{m}}))}^{\mathrm{H}}+{\mathsf{\sigma }}_{\mathrm{n}}{\mathrm{I}}_{\mathrm{N}\mathrm{K}\times \mathrm{N}\mathrm{K}})({\stackrel{⌢}{\mathrm{R}}}_{\mathrm{l}}-({\displaystyle \sum _{\mathrm{m}=1}^{{\mathrm{N}}_{\mathrm{s}}}{\left|{\mathsf{\alpha }}_{\mathrm{l},\mathrm{m}}\right|}^2}\mathrm{a}({\mathrm{f}}_{\mathrm{dm}},{\mathsf{\psi }}_{\mathrm{s}\mathrm{m}}){(\mathrm{a}({\mathrm{f}}_{\mathrm{dm}},{\mathsf{\psi }}_{\mathrm{sm}}))}^{\mathrm{H}}+{\mathsf{\sigma }}_{\mathrm{n}}{\mathrm{I}}_{\mathrm{N}\mathrm{K}\times \mathrm{N}\mathrm{K}}))\Vert }_{\mathrm{F}}^2\\ \mathrm{s}.\mathrm{t}.\frac{{\Vert {(\mathrm{diag}({\left|{\mathsf{\alpha }}_{\mathrm{l},\mathrm{m}}\right|}^2))}^{\mathrm{t}+1}-{(\mathrm{diag}({\left|{\mathsf{\alpha }}_{\mathrm{l},\mathrm{m}}\right|}^2))}^{\mathrm{t}}\Vert }_2^2}{{\Vert {(\mathrm{diag}({\left|{\mathsf{\alpha }}_{\mathrm{l},\mathrm{m}}\right|}^2))}^{\mathrm{t}}\Vert }_2^2}\le \mathsf{\kappa }\end{array}$$

(15)

In Equation (15), $\mathrm{J}$ is the objective function of the SPICE algorithm, which is to obtain the value of ${\left|{\mathsf{\alpha }}_{\mathrm{l},\mathrm{m}}\right|}^2$ by solving the minimum value problem of $\mathrm{J}$ and judging whether the constraint condition is valid. If the answer is yes, then ${\left|{\mathsf{\alpha }}_{\mathrm{l},\mathrm{m}}\right|}^2$ is the desired value; otherwise, we keep iterating until the value of ${\left|{\mathsf{\alpha }}_{\mathrm{l},\mathrm{m}}\right|}^2$ meets the constraint condition. ${\stackrel{⌢}{\mathrm{R}}}_{\mathrm{l}}$ is the maximum likelihood estimate of the covariance matrix; $\mathrm{diag}({\left|{\mathsf{\alpha }}_{\mathrm{l},\mathrm{m}}\right|}^2)$ is the power diagonal matrix of the clutter space–time two-dimensional spectrum formed by the sparse coefficient vector ${\left|{\mathsf{\alpha }}_{\mathrm{l},\mathrm{m}}\right|}^2$; $\mathsf{\kappa }$ represents the threshold of SPICE algorithm accuracy level, generally, the empirical value is 0.13; $\mathrm{t}+1$ is the current number of times that we need to use the SPICE algorithm, $\mathrm{t}=0,1,\dots ,\mathrm{N}$. For convenient representation, we substitute Equation (14) into the expression of J in Equation (15) and mark $\mathrm{diag}({\left|{\mathsf{\alpha }}_{\mathrm{l},\mathrm{m}}\right|}^2)$ as $\mathrm{B}$. Equation (15) can be expressed as:

$$\left\{\begin{array}{c}\mathrm{J}={‖{\mathrm{R}}_{\mathrm{l}}^{-\frac{1}{2}}({\stackrel{⌢}{\mathrm{R}}}_{\mathrm{l}}-{\mathrm{R}}_{\mathrm{l}})‖}_{\mathrm{F}}^2\\ \mathrm{s}.\mathrm{t}.\frac{{‖{\mathrm{B}}^{\mathrm{t}+1}-{\mathrm{B}}^{\mathrm{t}}‖}_2^2}{{‖{\mathrm{B}}^{\mathrm{t}}‖}_2^2}\le \mathsf{\kappa }\end{array}\right.$$

(16)

Since the proposed method does not use IID training samples, as it only uses the range unit ${\mathrm{x}}_{\mathrm{l}}$ to be detected, the ${\stackrel{⌢}{\mathrm{r}}}_{\mathrm{l}}$ expression is as follows:

$${\stackrel{⌢}{\mathrm{r}}}_{\mathrm{l}}={\mathrm{x}}_{\mathrm{l}}{\mathrm{x}}_{\mathrm{l}}^{\mathrm{H}}$$

(17)

We substitute Equation (17) into the objective function J expression of Equation (16) and expand and simplify the F norm to obtain:

$$\begin{array}{cc}\mathrm{J}\hfill & ={‖{\mathrm{R}}_{\mathrm{l}}^{-\frac{1}{2}}({\mathrm{x}}_{\mathrm{l}}{\mathrm{x}}_{\mathrm{l}}^{\mathrm{H}}-{\mathrm{R}}_{\mathrm{l}})‖}_{\mathrm{F}}^2\hfill \\ & =\mathrm{T}\mathrm{r}\left[{({\mathrm{x}}_{\mathrm{l}}{\mathrm{x}}_{\mathrm{l}}^{\mathrm{H}}-{\mathrm{R}}_{\mathrm{l}})}^{\mathrm{H}}{({\mathrm{R}}_{\mathrm{l}}^{-\frac{1}{2}})}^{\mathrm{H}}{\mathrm{R}}_{\mathrm{l}}^{-\frac{1}{2}}({\mathrm{x}}_{\mathrm{l}}{\mathrm{x}}_{\mathrm{l}}^{\mathrm{H}}-{\mathrm{R}}_{\mathrm{l}})\right]\hfill \\ & =\mathrm{T}\mathrm{r}\left[{\mathrm{R}}_{\mathrm{l}}^{-1}{({\mathrm{x}}_{\mathrm{l}}{\mathrm{x}}_{\mathrm{l}}^{\mathrm{H}}-{\mathrm{R}}_{\mathrm{l}})}^2\right]\hfill \\ & =\mathrm{T}\mathrm{r}({\mathrm{x}}_{\mathrm{l}}{\mathrm{x}}_{\mathrm{l}}^{\mathrm{H}}){\mathrm{x}}_{\mathrm{l}}^{\mathrm{H}}{\mathrm{R}}_{\mathrm{l}}^{-1}{\mathrm{x}}_{\mathrm{l}}+\mathrm{T}\mathrm{r}({\mathrm{R}}_{\mathrm{l}})-2{‖{\mathrm{x}}_{\mathrm{l}}‖}_2^2\hfill \\ & ={‖{\mathrm{x}}_{\mathrm{l}}‖}_2^2{\mathrm{x}}_{\mathrm{l}}^{\mathrm{H}}{\mathrm{R}}_{\mathrm{l}}^{-1}{\mathrm{x}}_{\mathrm{l}}+\mathrm{E}\left[{‖{\mathrm{x}}_{\mathrm{l}}‖}_2^2\right]-2{‖{\mathrm{x}}_{\mathrm{l}}‖}_2^2\hfill \end{array}$$

(18)

Strictly speaking, the mathematical expectation $\mathrm{E}\left[{‖{\mathrm{x}}_{\mathrm{l}}‖}_2^2\right]$ of L2-norm square of ${\mathrm{x}}_{\mathrm{l}}$ cannot be solved exactly [25]. However, since the radar echo signal ${\mathrm{x}}_{\mathrm{l}}$ of the lth range unit to be detected is superimposed by the ground clutter echo signal, noise, and low-altitude windshear echo signal, it can be written as follows:

$$\mathrm{E}\left[{‖{\mathrm{x}}_{\mathrm{l}}‖}_2^2\right]=\mathrm{E}\left[{‖{\mathrm{c}}_{\mathrm{l}}+{\mathrm{n}}_{\mathrm{l}}+{\mathrm{s}}_{\mathrm{l}}‖}_2^2\right]$$

(19)

Therefore, Equation (11) is substituted into Equation (19), and according to literature [24], we can obtain:

$$\begin{array}{ll}\mathrm{E}\left[{‖{\mathrm{x}}_{\mathrm{l}}‖}_2^2\right]\hfill & =\mathrm{E}\left[{‖{\mathrm{c}}_{\mathrm{l}}+{\mathrm{n}}_{\mathrm{l}}+{\mathrm{s}}_{\mathrm{l}}‖}_2^2\right]=\mathrm{E}\left[{‖{\displaystyle \sum _{\mathrm{m}=1}^{{\mathrm{N}}_{\mathrm{s}}}{\mathsf{\alpha }}_{\mathrm{l},\mathrm{m}}.\mathrm{a}({\mathrm{f}}_{\mathrm{dm}},{\mathsf{\psi }}_{\mathrm{sm}})}+{\mathrm{n}}_{\mathrm{l}}+{\mathrm{s}}_{\mathrm{l}}‖}_2^2\right]\hfill \\ & ={\displaystyle \sum _{\mathrm{m}=1}^{{\mathrm{N}}_{\mathrm{s}}}{\left|{\mathsf{\alpha }}_{\mathrm{l},\mathrm{m}}\right|}^2{‖\mathrm{a}({\mathrm{f}}_{\mathrm{dm}},{\mathsf{\psi }}_{\mathrm{sm}})‖}_2^2}+\mathrm{N}\mathrm{K}{\mathsf{\sigma }}_{\mathrm{n}}+{‖{\mathrm{s}}_{\mathrm{l}}‖}_2^2\hfill \end{array}$$

(20)

Substituting Equation (20) into Equation (18) to obtain:

$$\mathrm{J}={‖{\mathrm{x}}_{\mathrm{l}}‖}_2^2{\mathrm{x}}_{\mathrm{l}}^{\mathrm{H}}{\mathrm{R}}_{\mathrm{l}}^{-1}{\mathrm{x}}_{\mathrm{l}}+{\displaystyle \sum _{\mathrm{m}=1}^{{\mathrm{N}}_{\mathrm{s}}}{\left|{\mathsf{\alpha }}_{\mathrm{l},\mathrm{m}}\right|}^2{‖\mathrm{a}({\mathrm{f}}_{\mathrm{dm}},{\mathsf{\psi }}_{\mathrm{sm}})‖}_2^2}+\mathrm{N}\mathrm{K}{\mathsf{\sigma }}_{\mathrm{n}}+{‖{\mathrm{s}}_{\mathrm{l}}‖}_2^2-2{‖{\mathrm{x}}_{\mathrm{l}}‖}_2^2$$

(21)

The last two terms in Equation (21) are constants, which do not affect the result of solving its minimum value problem, so after ignoring the ${‖{\mathrm{s}}_{\mathrm{l}}‖}_2^2-2{‖{\mathrm{x}}_{\mathrm{l}}‖}_2^2$, the expression for ${\mathrm{J}}_1$ is:

$${\mathrm{J}}_1={‖{\mathrm{x}}_{\mathrm{l}}‖}_2^2{\mathrm{x}}_{\mathrm{l}}^{\mathrm{H}}{\mathrm{R}}_{\mathrm{l}}^{-1}{\mathrm{x}}_{\mathrm{l}}+{\displaystyle \sum _{\mathrm{m}=1}^{{\mathrm{N}}_{\mathrm{s}}}{\left|{\mathsf{\alpha }}_{\mathrm{l},\mathrm{m}}\right|}^2{‖\mathrm{a}({\mathrm{f}}_{\mathrm{dm}},{\mathsf{\psi }}_{\mathrm{sm}})‖}_2^2}+\mathrm{N}\mathrm{K}{\mathsf{\sigma }}_{\mathrm{n}}$$

(22)

Then, the ${\left|{\mathsf{\alpha }}_{\mathrm{l},\mathrm{m}}\right|}^2$ solution obtained by solving the minimum problem of $\mathrm{J}$ is the same as that obtained by solving the minimum problem of ${\mathrm{J}}_1$. Since ${‖{\mathrm{x}}_{\mathrm{l}}‖}_2^2$ is a constant, to simplify the calculation, we divide both of the sides of Equation (22) by ${‖{\mathrm{x}}_{\mathrm{l}}‖}_2^2$ at the same time to obtain:

$$\frac{{\mathrm{J}}_1}{{‖{\mathrm{x}}_{\mathrm{l}}‖}_2^2}={\mathrm{x}}_{\mathrm{l}}^{\mathrm{H}}{\mathrm{R}}_{\mathrm{l}}^{-1}{\mathrm{x}}_{\mathrm{l}}+{\displaystyle \sum _{\mathrm{m}=1}^{{\mathrm{N}}_{\mathrm{s}}}{\left|{\mathsf{\alpha }}_{\mathrm{l},\mathrm{m}}\right|}^2\frac{{‖\mathrm{a}({\mathrm{f}}_{\mathrm{dm}},{\mathsf{\psi }}_{\mathrm{sm}})‖}_2^2}{{‖{\mathrm{x}}_{\mathrm{l}}‖}_2^2}}+\mathrm{N}\mathrm{K}\frac{{\mathsf{\sigma }}_{\mathrm{n}}}{{‖{\mathrm{x}}_{\mathrm{l}}‖}_2^2}$$

(23)

Let the expression ${\mathrm{J}}_2$ be:

$${\mathrm{J}}_2={\mathrm{x}}_{\mathrm{l}}^{\mathrm{H}}{\mathrm{R}}_{\mathrm{l}}^{-1}{\mathrm{x}}_{\mathrm{l}}+{\displaystyle \sum _{\mathrm{m}=1}^{{\mathrm{N}}_{\mathrm{s}}}{\left|{\mathsf{\alpha }}_{\mathrm{l},\mathrm{m}}\right|}^2\frac{{‖\mathrm{a}({\mathrm{f}}_{\mathrm{dm}},{\mathsf{\psi }}_{\mathrm{sm}})‖}_2^2}{{‖{\mathrm{x}}_{\mathrm{l}}‖}_2^2}}+\mathrm{N}\mathrm{K}\frac{{\mathsf{\sigma }}_{\mathrm{n}}}{{‖{\mathrm{x}}_{\mathrm{l}}‖}_2^2}$$

(24)

Then, solving the minimum problem of ${\mathrm{J}}_1$ has the same solution as solving the minimum problem of ${\mathrm{J}}_2$. According to the study of Rojas et al. [26], the minimum problem of Equation (24) can be written as:

$$\left\{\begin{array}{l}\min {\mathrm{x}}_{\mathrm{l}}^{\mathrm{H}}{\mathrm{R}}_{\mathrm{l}}^{-1}{\mathrm{x}}_{\mathrm{l}}\\ \mathrm{s}.\mathrm{t}.{\displaystyle \sum _{\mathrm{m}=1}^{{\mathrm{N}}_{\mathrm{s}}}{\left|{\mathsf{\alpha }}_{\mathrm{l},\mathrm{m}}\right|}^2\frac{{‖\mathrm{a}({\mathrm{f}}_{\mathrm{dm}},{\mathsf{\psi }}_{\mathrm{sm}})‖}_2^2}{{‖{\mathrm{x}}_{\mathrm{l}}‖}_2^2}}+\mathrm{N}\mathrm{K}\frac{{\mathsf{\sigma }}_{\mathrm{n}}}{{‖{\mathrm{x}}_{\mathrm{l}}‖}_2^2}\le 1\end{array}\right.$$

(25)

According to the literature [23], to solve Equation (25), the clutter covariance matrix ${\mathrm{R}}_{\mathrm{l}}$ in Equation (25) is expressed as:

$${\mathrm{R}}_{\mathrm{l}}^{}={\Theta }_{}^{\mathrm{H}}\mathrm{B}\Theta$$

(26)

where $\Theta =[{\mathrm{A}}_1,{\mathrm{I}}_{\mathrm{N}\mathrm{K}}]$ is a $\mathrm{N}\mathrm{K}\times ({\mathrm{N}}_{\mathrm{s}}+\mathrm{N}\mathrm{K})$-dimensional sparse dictionary composed of clutter and noise. By substituting Equation (26) into Equation (25), we can obtain:

$$\left\{\begin{array}{l}\min {\mathrm{x}}_{\mathrm{l}}^{\mathrm{H}}{({\Theta }_{}^{\mathrm{H}}\mathrm{B}\Theta )}^{-1}{\mathrm{x}}_{\mathrm{l}}\\ \mathrm{s}.\mathrm{t}.{\displaystyle \sum _{\mathrm{m}=1}^{{\mathrm{N}}_{\mathrm{s}}}{\left|{\mathsf{\alpha }}_{\mathrm{l},\mathrm{m}}\right|}^2\frac{{‖\mathrm{a}({\mathrm{f}}_{\mathrm{dm}},{\mathsf{\psi }}_{\mathrm{sm}})‖}_2^2}{{‖{\mathrm{x}}_{\mathrm{l}}‖}_2^2}}+\mathrm{N}\mathrm{K}\frac{{\mathsf{\sigma }}_{\mathrm{n}}}{{‖{\mathrm{x}}_{\mathrm{l}}‖}_2^2}\le 1\end{array}\right.$$

(27)

In this case, in order to solve the minimum problem of Equation (27), a new matrix $\Psi$ is introduced, which satisfies:

$$\Psi ={({\Theta }^{-1})}^{\mathrm{H}}$$

(28)

Namely:

$$\left\{\begin{array}{l}\Theta ={\left({\Psi }^{\mathrm{H}}\right)}^{-1}\\ {\Psi }^{\mathrm{H}}\Theta ={\mathrm{I}}_{\mathrm{N}\mathrm{K}}\end{array}\right.$$

(29)

By substituting Equation (29) into Equation (27), we can obtain:

$$\left\{\begin{array}{l}\min {\mathrm{x}}_{\mathrm{l}}^{\mathrm{H}}{({\Theta }_{}^{\mathrm{H}}\mathrm{B}\Theta )}^{-1}{\mathrm{x}}_{\mathrm{l}}=\min {\mathrm{x}}_{\mathrm{l}}^{\mathrm{H}}{\Theta }^{-1}{\mathrm{B}}^{-1}{({\Theta }^{\mathrm{H}})}^{-1}{\mathrm{x}}_{\mathrm{l}}=\min {\mathrm{x}}_{\mathrm{l}}^{\mathrm{H}}{\Psi }^{\mathrm{H}}{\mathrm{B}}^{-1}\Psi {\mathrm{x}}_{\mathrm{l}}\\ \mathrm{s}.\mathrm{t}.{\Psi }^{\mathrm{H}}\Theta ={\mathrm{I}}_{\mathrm{N}\mathrm{K}}\\ {\displaystyle \sum _{\mathrm{m}=1}^{{\mathrm{N}}_{\mathrm{s}}}{\left|{\mathsf{\alpha }}_{\mathrm{l},\mathrm{m}}\right|}^2\frac{{‖\mathrm{a}({\mathrm{f}}_{\mathrm{dm}},{\mathsf{\psi }}_{\mathrm{sm}})‖}_2^2}{{‖{\mathrm{x}}_{\mathrm{l}}‖}_2^2}}+\mathrm{N}\mathrm{K}\frac{{\mathsf{\sigma }}_{\mathrm{n}}}{{‖{\mathrm{x}}_{\mathrm{l}}‖}_2^2}\le 1\end{array}\right.$$

(30)

According to the literature [27,28,29], the solution to the minimum value problem is:

$$\left\{\begin{array}{l}{\left|{\mathsf{\alpha }}_{\mathrm{l},\mathrm{m}}\right|}^2=\frac{\left|\Psi {\mathrm{x}}_{\mathrm{l}}\right|}{w_{\mathrm{l},\mathrm{m}}^{1/2}\mathsf{\rho }}\\ \Psi ={\mathrm{B}}_{(\mathrm{t})}\Theta {({\mathrm{r}}_{\mathrm{l}}{}_{(\mathrm{t})})}^{-1}\hspace{1em}(\mathrm{t}=0,1,\dots ,\mathrm{N})\end{array}\right.$$

(31)

where $\mathsf{\rho }={\displaystyle \sum _{\mathrm{m}=1}^{{\mathrm{N}}_{\mathrm{s}}}{w}_{\mathrm{l},\mathrm{m}}^{1/2}\left|\Psi {\mathrm{x}}_{\mathrm{l}}\right|}+{\displaystyle \sum _{\mathrm{N}={\mathrm{N}}_{\mathrm{s}}+1}^{{\mathrm{N}}_{\mathrm{s}}+\mathrm{N}\mathrm{K}}\frac{\left|\Psi {\mathrm{x}}_{\mathrm{l}}\right|}{{‖{\mathrm{x}}_{\mathrm{l}}‖}_2^{}}}$;$w_{\mathrm{l},\mathrm{m}}^{}=\frac{{‖\mathrm{a}({\mathrm{f}}_{\mathrm{dm}},{\mathsf{\psi }}_{\mathrm{sm}})‖}_2^2}{{‖{\mathrm{x}}_{\mathrm{l}}‖}_2^2}$. The solution of ${\mathrm{R}}_{\mathrm{l}}{}_{(\mathrm{t})}$ and ${\mathrm{B}}_{(\mathrm{t})}$ is as follows [24]:

$$\left\{\begin{array}{l}{\mathrm{R}}_{\mathrm{l}}{}_{(\mathrm{t})}={\displaystyle \sum _{\mathrm{m}=1}^{{\mathrm{N}}_{\mathrm{s}}}\frac{\left|\mathrm{a}{({\mathrm{f}}_{\mathrm{dm}},{\mathsf{\psi }}_{\mathrm{sm}})}^{\mathrm{H}}{\mathrm{x}}_{\mathrm{l}}\right|}{{‖\mathrm{a}({\mathrm{f}}_{\mathrm{dm}},{\mathsf{\psi }}_{\mathrm{sm}})‖}_2^2}}\mathrm{a}({\mathrm{f}}_{\mathrm{dm}},{\mathsf{\psi }}_{\mathrm{sm}}){(\mathrm{a}({\mathrm{f}}_{\mathrm{dm}},{\mathsf{\psi }}_{\mathrm{sm}}))}^{\mathrm{H}}+{\mathsf{\sigma }}_{\mathrm{n}}{\mathrm{I}}_{\mathrm{N}\mathrm{K}\times \mathrm{N}\mathrm{K}}\hspace{1em}(\mathrm{t}=0)\\ {\mathrm{B}}_{(\mathrm{t})}=\mathrm{diag}(\frac{\left|\mathrm{a}{({\mathrm{f}}_{\mathrm{dm}},{\mathsf{\psi }}_{\mathrm{sm}})}^{\mathrm{H}}{\mathrm{x}}_{\mathrm{l}}\right|}{{‖\mathrm{a}({\mathrm{f}}_{\mathrm{dm}},{\mathsf{\psi }}_{\mathrm{sm}})‖}_2^2})\hspace{1em}(\mathrm{t}=0)\end{array}\right.$$

(32)

$$\left\{\begin{array}{l}{\mathrm{R}}_{\mathrm{l}}{}_{(\mathrm{t})}={\displaystyle \sum _{\mathrm{m}=1}^{{\mathrm{N}}_{\mathrm{s}}}{\left|{\mathsf{\alpha }}_{\mathrm{l},\mathrm{m}}\right|}_{(\mathrm{t}-1)}^2}\mathrm{a}({\mathrm{f}}_{\mathrm{dm}},{\mathsf{\psi }}_{\mathrm{sm}}){(\mathrm{a}({\mathrm{f}}_{\mathrm{dm}},{\mathsf{\psi }}_{\mathrm{sm}}))}^{\mathrm{H}}+{\mathsf{\sigma }}_{\mathrm{n}}{\mathrm{I}}_{\mathrm{N}\mathrm{K}\times \mathrm{N}\mathrm{K}}\hspace{1em}(\mathrm{t}=1,2\dots ,\mathrm{N})\\ {\mathrm{B}}_{(\mathrm{t})}=\mathrm{diag}({\left|{\mathsf{\alpha }}_{\mathrm{l},\mathrm{m}}\right|}_{(\mathrm{t}-1)}^2)\hspace{1em}(\mathrm{t}=1,2,\dots ,\mathrm{N})\end{array}\right.$$

(33)

where Equation (32) is used to calculate ${\mathrm{R}}_{\mathrm{l}}{}_{(\mathrm{t})}$ and ${\mathrm{B}}_{(\mathrm{t})}$ when $\mathrm{t}=0$, and this is the first time the SPICE algorithm is used. Equation (33) is used when $\mathrm{t}=1,2,\dots ,\mathrm{N}$ to calculate ${\mathrm{R}}_{\mathrm{l}}{}_{(\mathrm{t})}$ and ${\mathrm{B}}_{(\mathrm{t})}$ in the subsequent iteration cycle and in the Equation, ${\left|{\mathsf{\alpha }}_{\mathrm{l},\mathrm{m}}\right|}_{(\mathrm{t}-1)}^2(\mathrm{t}=1,2,\dots ,\mathrm{N})$ is the ${\left|{\mathsf{\alpha }}_{\mathrm{l},\mathrm{m}}\right|}^2$ obtained from the previous cycle ${\mathrm{R}}_{\mathrm{l}}{}_{(\mathrm{t}-1)}(\mathrm{t}=1,2,\dots ,\mathrm{N})$ and ${\mathrm{B}}_{(\mathrm{t}-1)}(\mathrm{t}=1,2,\dots ,\mathrm{N})$ are obtained from the cycle $\mathrm{t}$.

Firstly, the constructed clutter dictionary is used to obtain ${\mathrm{R}}_{\mathrm{l}}{}_{(0)}$ and ${\mathrm{B}}_{(0)}$ when the SPICE algorithm is first used according to Equation (32), and ${\left|{\mathsf{\alpha }}_{\mathrm{l},\mathrm{m}}\right|}_{(0)}^2$ is obtained according to Equation (31). Secondly, it is judged whether ${\left|{\mathsf{\alpha }}_{\mathrm{l},\mathrm{m}}\right|}_{(0)}^2$ meets the constraint condition of Equation (15). If the constraint condition is satisfied, ${\left|{\mathsf{\alpha }}_{\mathrm{l},\mathrm{m}}\right|}_{(0)}^2$ is the required ${\left|{\mathsf{\alpha }}_{\mathrm{l},\mathrm{m}}\right|}^2$. Otherwise, we put ${\left|{\mathsf{\alpha }}_{\mathrm{l},\mathrm{m}}\right|}_{(0)}^2$ into Equation (33) and iterate the cycle continuously until the constraint condition of Equation (15) is met at a certain iteration, then, ${\left|{\mathsf{\alpha }}_{\mathrm{l},\mathrm{m}}\right|}_{(\mathrm{t}-1)}^2(\mathrm{t}=1,\dots ,\mathrm{N})$ that is obtained in this iteration is regarded as ${\left|{\mathsf{\alpha }}_{\mathrm{l},\mathrm{m}}\right|}^2$. Finally, we put ${\left|{\mathsf{\alpha }}_{\mathrm{l},\mathrm{m}}\right|}^2$ into Equation (14) to obtain the clutter covariance matrix ${\mathrm{R}}_{\mathrm{l}}^{}$.

The flow diagram of clutter covariance matrix estimation based on SPICE algorithm is shown in Figure 3 below.

3.2. Low-Altitude Windshear Wind Speed Estimation

After the clutter covariance matrix is obtained iteratively by the SPICE algorithm, the STAP processor is designed to suppress the ground clutter echo signal and estimate the wind speed. The optimal weight vector is solved by the following equation:

$$\left\{\begin{array}{l}\min w_{\mathrm{l}}^{\mathrm{H}}{\mathrm{R}}_{\mathrm{l}}{w}_{\mathrm{l}}^{}\hfill \\ \mathrm{s}.\mathrm{t}.w_{\mathrm{l}}^{\mathrm{H}}{\mathrm{G}}_{\mathrm{l}}=1\hfill \end{array}\right.$$

(34)

where ${\mathrm{G}}_{\mathrm{l}}={\mathrm{v}}_{\mathrm{l}}{\left({\mathrm{f}}_{\mathrm{d}},{\mathsf{\psi }}_0\right)}_{\mathrm{N}\mathrm{K}\times 1}$ is the extended low-altitude windshear space–time steering vector. According to Equation (34), the optimal weight vector is:

$$w_{\mathrm{l}}^{}=\mathsf{\mu }{\mathrm{R}}_{\mathrm{l}}^{-1}{\mathrm{G}}_{\mathrm{l}}^{}$$

(35)

After the optimal weight vector $w_{\mathrm{l}}^{}$ is used to process the radar echo signal, the ground clutter in the range unit ${\mathrm{x}}_{\mathrm{l}}$ that is to be detected is suppressed, and the normalized Doppler frequency is solved by the following equation:

$${\widehat{\mathrm{f}}}_{\mathrm{l}}=\underset{}{\mathrm{argmax}}{w}_{\mathrm{l}}^{\mathrm{H}}{\mathrm{x}}_{\mathrm{l}}$$

(36)

Then, the wind field target velocity estimation result of the lth range unit to be detected is:

$${\widehat{\mathrm{v}}}_{\mathrm{l}}=\frac{\mathsf{\lambda }{\widehat{\mathrm{f}}}_{\mathrm{l}}{\mathrm{f}}_{\mathrm{r}}}{4}$$

(37)

4. Process of the Method

Figure 4 shows the flow diagram of the estimation method of low-altitude windshear wind speed based on KASPICE-STAP.

The method proposed in this paper is able to realize a low-altitude windshear wind speed estimation without IID training samples. The specific steps are as follows:

Step 1: Construct the corresponding clutter dictionary by using prior knowledge of the position distribution of ground clutter echo signals in space–time two-dimensional spectrum.

Step 2: Obtain the clutter covariance matrix iteratively using the SPICE algorithm.

Step 3: Design the STAP processor, calculate the optimal weight vector, and then suppress the ground clutter.

Step 4: After suppressing the ground clutter, the normalized Doppler frequency is obtained, and the wind speed for the detection range unit is estimated.

5. Results

In this paper, it is assumed that the low-altitude windshear wind field is located 8.5 km~16.5 km in front of the aircraft platform, the center of the low-altitude windshear wind field is located 12.5 away from the 60° direction of the aircraft platform, the main lobe of the antenna is aligned with the low-altitude windshear wind field, and the width of the low-altitude windshear wind field is 8 km. System simulation parameters are shown in

Table 1. System simulation parameters.

Parameter	Value	Parameter	Value
Aircraft platform height (m)	600	Number of array elements	8
Aircraft platform speed (m/s)	75	Number of coherent pulse	32
Radar wavelength (m)	0.032	Angle of main lobe (°)	(60,0)
Pulse repetition frequency (Hz)	7000	Range resolution (m)	150
Signal-to-noise ratio (dB)	5	Clutter-to-noise ratio (dB)	40

.

Figure 5 shows the space–time two-dimensional spectrum of the radar echo signal of an airborne weather radar forward-looking array, wherein the ground clutter of the airborne forward-looking array has the space–time coupling characteristic, the ground clutter echo signal presents the semi-elliptical feature in the space–time two-dimensional spectrum, and the power spectrum of low-altitude windshear signal has a certain broadening in the normalized Doppler frequency domain and the spatial cone Angle cosine domain. It can be concluded from the analysis that the power of the ground clutter echo signal is much stronger than that of the low-altitude windshear signal, which makes the ground clutter echo signal completely cover up the normalized Doppler frequency information of the low-altitude windshear echo signal, which is an important factor causing the inaccurate performance of the low-altitude windshear estimation.

Figure 6a,b shows the space–time two-dimensional spectrum of the ground clutter obtained using the KASPICE-STAP method and DDD-SR-STAP method for the sparse recovery of the detection range units, respectively. As can be seen from Figure 6, although the DDD-SR-STAP method recovers the space–time two-dimensional spectrum of the ground clutter echo signal, it also presents semi-elliptical distribution characteristics, the space–time spectrum of ground clutter echo signal has a serious discontinuous distribution, and the recovery effect of the space–time two-dimensional spectrum of ground clutter echo signal is very poor. The space–time two-dimensional spectrum of the ground clutter echo signal obtained by the sparse recovery method in this paper is more accurate, and the semi elliptical distribution characteristics and the main lobe of the ground clutter echo signals can be clearly seen, which is more consistent with the space–time two-dimensional spectrum of the ground clutter echo signal in Figure 5 radar echo signal, showing the ultra-high resolution performance advantages of the sparse recovery technology.

Figure 7 shows the clutter characteristic spectra of the clutter covariance matrix estimated by the KASPICE-STAP method, the DDD-STAP method, and the DDD-SR-STAP method. It can be seen from Figure 7 that the number of large eigenvalues of the clutter covariance matrix estimated by the KASPICE-STAP method is smaller than that of the other two methods, that is, the clutter degree of freedom is low, so the clutter covariance matrix estimated by the KASPICE-STAP method is more accurate.

Figure 8 shows the factor of improvement curves of the KASPICE-STAP method, the DDD-STAP method, and the DDD-SR-STAP method. Taking the No. 60 range unit as an example, it can be seen from Figure 8 that compared with the other two methods, the KASPICE-STAP method has a deeper and narrower notch in the main clutter area, which is able to suppress the clutter echo signals more effectively.

Figure 9 shows the comparison of the low-altitude windshear wind speed estimations by the KASPICE-STAP method, the DDD-STAP method, and the DDD-SR-STAP method. As can be seen from Figure 9, the DDD-STAP method cannot accurately estimate the low-altitude windshear wind speed due to the loss of degrees of freedom in the sliding window and the limitation of the degrees of freedom in the space and time domain. The DDD-SR-STAP method is affected by only one single range unit data that are to be detected, and the accuracy of low-altitude windshear wind speed estimation is not high. The KASPICE-STAP method can obtain better wind speed estimation results with higher accuracy. At the same time, within the range of the 8.5 km–16.5 km low-altitude windshear wind field, the wind speed of the low-altitude windshear wind field shows an inverse “S” shape with distance.

Table 2. Comparison of root mean square error of low-altitude windshear wind speed estimation.

Wind Speed Estimation Method	Root Mean Square Error/(m/s)
DDD-STAP	13.2187
DDD-SR-STAP	4.7335
KASPICE-STAP	1.8692

shows the comparison of root mean square error of low-altitude windshear wind speed estimation by the KASPICE-STAP method, the DDD-STAP method, and the DDD-SR-STAP method. It can be seen from Table 2 that the root mean square error of low-altitude windshear wind speed estimation by the DDD-STAP method and the DDD-SR-STAP method is larger than that which was achieved by the KASPICE-STAP method. Therefore, the KASPICE-STAP method can achieve accurate estimation results of low-altitude windshear wind speed.

6. Conclusions

When one is estimating low-altitude windshear wind speed, being faced with no IID training samples is a common situation. The traditional STAP method is unable to accurately estimate the clutter covariance matrix and suppress the clutter, which severely degrades the performance of low-altitude windshear wind speed estimation. In this paper, a low-altitude wind-speed estimation method based on KASPICE-STAP was designed. This method was applied to airborne weather radar, and it can be used for the clutter suppression and parameter estimation of distributed targets without any IID training samples, and it shows a good performance when there is only a single range unit that is to be detected. It uses prior knowledge of the ground clutter echo signals, fully excavated the sparsity of the space–time two-dimensional spectrum of the range unit that is to be detected to construct the clutter dictionary, and then, it uses the SPICE algorithm to iteratively obtain the clutter covariance matrix. Compared with the traditional DDD-STAP method, this method obtains a more accurate clutter covariance matrix due to its excellent performance. Finally, the STAP processor was designed, and the corresponding optimal weight vector was calculated, the normalized Doppler frequency was proposed from the signal after eliminating the ground clutter afterwards, and ultimately, the wind speed was estimated. The simulation results showed that the KASPICE-STAP method can still effectively and accurately estimate the low-altitude windshear speed without IID training samples. At present, the limitation of the method is that the clutter dictionary is traversed in each iteration of the SPICE algorithm, which has a high computational complexity and leads to a long operation time. This will be our future research direction.