Space-Time Cascaded Processing-Based Adaptive Transient Interference Mitigation for Compact HFSWR

Abstract

In high-frequency (HF) radar systems, transient interference is a common phenomenon that dramatically degrades the performance of target detection and remote sensing. Up until now, various suppression algorithms of transient interference have been proposed. They mainly concentrate on the skywave over-the-horizon radar on the basis of the assumption that the interference is sparse in a coherent processing interval (CPI). However, HF surface wave radar (HFSWR) often faces more complex transient interference due to various extreme types of weather, such as thunderstorm and typhoon, etc. The above algorithms usually suffer dramatic performance loss when transient interference contaminates the enormously continuous parts of a CPI. Especially for the compact HFSWR, which suffers from severe beam broadening and fewer array degrees of freedom. In order to solve the above problem, this study developed a two-dimensional interference suppression algorithm based on space-time cascaded processing. First, according to the spatial correlation of the compact array, the statistical samples of the main-lobe transient interference are estimated using a rotating spatial beam method. Next, an adaptive selection strategy is developed to obtain the optimal secondary samples based on information geometry distance. Finally, based on a quadratic constraint approximation, a precise estimation method of the optimal weight is developed when the interference covariance matrix is singular. The experimental results of simulation and measured data demonstrate that the proposed approach provides far superior suppression performance.

1. Introduction

High-frequency surface wave radar (HFSWR) plays an important role in maritime targets detection and wide-range ocean remote sensing over the horizon, working at 3–30 MHz. Based on the theoretical basis of Barrick’s first- and second-order cross-sectional equations, the ocean dynamics parameters can be effectively extracted using the high-frequency (HF) waves backscattered from the surface of the sea, such as the speed and direction of current and wind and wave [1,2]. However, in the HF band, the performance of the HFSWR system is seriously affected by external transient interferences, e.g., lightning, meteor echoes, man-made impulse bursts, and so on [3]. Transient interferences are uncorrelated with the transmitted signal, and even increase the noise level by as much as 20 dB in the Doppler spectrum [4]. Therefore, the weak sea echoes and targets would be submerged by the transient interference because it usually has very powerful energy and can occupy all the range units and most Doppler units.

Depending on its duration in the slow-time domain (STD), transient interference can be roughly classified into two types: short transient interference (STI), and long continuous interference (LCI). The duration of STI is no more than 3% that of echoes, while LCI has a relatively long duration [5]. Several previous studies have been investigated the suppression of transient interference. The most considerable interference excision methods for HF radar systems are based on spectrum reconstruction. According to the location result of transient interferences, the corrupted data are firstly set to zeros in the STD. Next, removed data samples are reconstructed using a linear prediction model [6], compressed sensing method [7], complex empirical mode decomposition (CEMD) [8], time-frequency analysis [9,10], and so on. The procedure recovers the lost data using the non-contaminated samples. However, the reconstruction error is not acceptable in situations with complex echoes. A robust principal component analysis (RPCA) technology is proposed by constructing the STD signals into a Hankel matrix in which sea echoes and transient interferences are considered as a low-rank and sparse matrix, respectively, and it separates transient interference components and the sea echoes directly [4,11]. This method has demonstrated good performance in STI excision when interference satisfies the sparsity in a coherent processing interval (CPI), but it is difficult to remove the LCI. In [12,13], a no-data interpolation method was proposed based on the time-domain filter technique. The basic principle of the time-domain filter method is to minimize the quadratic fitting errors and the energy orthogonal to the signal subspace, meanwhile, not requiring the interference detection step. In practice, we have found that the above types of methods only work well for transient interference when they occupy a small number of consecutive STD units. However, in cases where signals are complex and many consecutive STD units are contaminated by the interferences, they cannot effectively suppress the interference. Focusing on the problem, the adaptive beamforming method is usually used to simultaneously suppress the STI and LCI owing to the directional characteristics of transient interferences. Next, the transient interference components are subtracted from the original echoes [14,15,16]. Nevertheless, this kind of method cannot remove the interference received through the main-lobe beam and requires a large array aperture to ensure enough degrees of freedom (DOF) for the system. Especially for compact HFSWR, a smaller array aperture has few array DOFs and can severely widen the space beam, which will significantly increase the happening chance of main-lobe interferences.

The compact HFSWRs require less deployment space and lower manufacturing costs and are more appropriate for installation at the seashore due to covering less space. Up to now, various compact HFSWRs such as WERA using four antennas, OSMAR using eight antennas, and CODAR have represented the most commonly installed systems around the world [17,18]. Therefore, this study is concerned with two aspects of the performance improvement of the transient interference suppression approach: first, to effectively remove the STI and LCI simultaneously. Second, to solve the problem of few DOFs and main-lobe interference in compact HFSWR.

With regard to the aforementioned analysis, a novel approach for transient interference excision is proposed based on space-time two-dimensional cascaded processing technology, which is suitable for the STI and LCI cases in compact HFSWR. Based on the published literature, it is assumed that transient interferences have been effectively detected and located, based on prior information. First, according to the strong space correlation of a small aperture array, the rotating spatial beam method is proposed to improve system DOFs and obtains enough statistical samples of the main-lobe interferences. Next, based on information geometry distance theory, an adaptive selection strategy is developed to construct the optimal secondary samples of transient interference in the time domain. Finally, to solve the singularity of the interference covariance matrix problem, a precise estimation method of the optimal weight is developed according to the quadratic constraint approximation model.

2. Materials and Methods

2.1. Data Sets Received by the Compact HFSWR

The radar system, with an eight-element configuration, is located in Weihai, China. The radar frequency is 8.15 MHz. Figure 1a shows the amplitude of the range-STD LCI data shown in dB. It can be seen that the powers of interferences from the 656th to 1515th slow-time indexes are much larger than the sea echoes and background noise. This situation is usually caused by a large number of discrete lightning impulse interferences. Figure 1b presents the range-Doppler spectrum of the corresponding LCI, where it can be observed that all the range units and Doppler units are corrupted by transient interferences and the weak sea echoes are submerged. Figure 1c displays the range-STD spectrum of STI, where the transient interference becomes sparse in the STD. The range–Doppler spectrum in Figure 1d shows the same distribution of range and Doppler as Figure 1b, except the interference energy is lower.

2.2. Signal Model of Transient Interference and Sea Echoes

In the HFSWR system, the transient interference suppression is usually performed in the STD after pulse compression and array beamforming but prior to Doppler processing. For a given range unit and beam unit, the STD signal of HFSWR within a CPI is symbolized by $\mathit{x}(m)$. It can be generally modeled as an additive mixture of the sum of sea echoes $\mathit{c}(m)$, transient interference $\mathit{i}(t)$, and additive noise $\mathit{n}(m)$, which leads to the expression,

$${\displaystyle \mathit{x}}(m)=\mathit{c}(m)+\mathit{i}(m)+\mathit{n}(m)$$

(1)

where $m=1,2,\cdots ,M$ denotes the slow-time index, $M$ represents the number of pulses in a CPI. The receiving array of the HFSWR is a uniform linear array (ULA) with $N$ linear equispaced omnidirectional sensors. The spacing between adjacent sensors is $d$.

Generally, the dominant sea echoes received by the HFSWR system can be characterized by a pair of peaks in the Doppler domain, whose energy occupies numerous range bins and few Doppler bins symmetrically placed around 0 Hz. According to the Bragg scattering hypothesis, the Doppler frequencies of the first-order scattering Bragg waves are $\pm f_b=\pm \sqrt{g{f}_0/\pi \mathit{c}}$, where $f_0$ is radar working frequency, $g$ indicates the gravity acceleration, and $\mathit{c}$ is the speed of light. However, in practice, the surface current of the sea echoes may cause the Bragg lines to be shifted in the Doppler domain, which can be adequately modeled by,

$$\mathit{c}(m)=A_p{e}^{j\left(2\pi \left(f_s+\left|f_b\right|\right)mT+{\phi }_1\left(m\right)\right)}+A_n{e}^{j\left(2\pi \left(f_s-\left|f_b\right|\right)mT+{\phi }_2\left(m\right)\right)}$$

(2)

where $f_s=2v_s/\lambda$ denotes Doppler shift corresponding to the surface current of radial velocity $v_s$, $\lambda$ is radar wavelength, $A_p$ and $A_n$ denote the complex amplitudes of the advance and recede Bragg components, respectively. ${\phi }_1(m)$ and ${\phi }_2(m)$ indicate the slow-time-varying phase modulated by ocean turbulence and ionosphere, respectively, $T$ indicates the pulse repetition interval (PRI).

The once transient interference is usually short-lived in the STD, but its intensity is great and spread in the Doppler domain [19]. These physical characteristics of transient interference usually lead the weak sea echoes to be submerged, which fails the detection of the desired signals. A typical transient interference of array form can be modeled as,

$$\mathit{i}(m)={\displaystyle \sum _{k=1}^K{S}_k(m)\ast {\displaystyle \mathit{a}}({\theta }_k)}$$

(3)

where,

$$S_k(m)=A_k{e}^{-\frac{{(mT-t_k)}^2}{2{\sigma }_k^2}}{e}^{j(2\pi f_{k}mT+{\phi }_k(m))}$$

(4)

$$\mathit{a}({\theta }_k)=\left[1,e^{j2\pi d\sin {\theta }_k/\lambda },\cdots ,e^{j2\pi (N-1)d\sin {\theta }_k/\lambda }\right]$$

(5)

where $K$ indicates the number of transient interferences, ${\sigma }_k$ is the variance, $t_k$ is the location in the time domain, $f_k$ is the locations in the Doppler domain, $A_k$ is the complex envelope, ${\phi }_k(m)$ means the modulation phase, $\mathit{a}({\theta }_k)$ denotes the array steering vector in the direction of arrival (DOA) ${\theta }_k$ of the kth interference.

2.3. Signal Model of Transient Interference Excision

Based on the sliding sub-arrays technology, a HFSWR main-lobe interference canceller was developed in [20], which was called spread interference estimation canceller (SIEC) for space spread interference (SSI) mitigation. This is an effective interference mitigation method for the situation of the target submerged by the main-lobe SSI. The outputs of SIEC can be expressed as,

$$\mathit{y}(m)={({\mathit{w}}_q-\mathit{B}{\mathit{w}}_a)}^H\mathit{x}(m)$$

(6)

where ${\mathit{w}}_q$ is the $N\times 1$ complex weight, which indicates the static weight to keep the main-beam direction and side-lobe level. The block matrix $\mathit{B}$ is used for blocking the desired signal, which can be substituted by a single-notch space filter (SNSF). Let $\mathit{J}=E\left[{\left|\mathit{y}(m)\right|}^2\right]$ denote a cost function. The weight ${\mathit{w}}_{\mathit{a}}$ can be obtained by solving the following unconstrained optimization problem:

$$\underset{{\mathit{w}}_a}{\min }\mathit{J}=\underset{{\mathit{w}}_a}{\min }{({\mathit{w}}_q-\mathit{B}{\mathit{w}}_a)}^H{\mathit{R}}_{\mathit{x}}({\mathit{w}}_q-\mathit{B}{\mathit{w}}_{\mathit{a}})$$

(7)

Next, we can find the optimum ${\mathit{w}}_{\mathit{a}}$ as,

$${\mathit{w}}_{\mathit{a},opt}={({\mathit{B}}^H{\mathit{R}}_{\mathit{x}}\mathit{B})}^{-1}{\mathit{B}}^H{\mathit{R}}_{\mathit{x}}{\mathit{w}}_q$$

(8)

Let ${\mathit{w}}_{opt}={\mathit{w}}_q-\mathit{B}{\mathit{w}}_{\mathit{a},opt}$. The minimum mean-squared error (MMSE) for Equation (6), denoted as ${\mathit{J}}_{\min }$, is written as,

$${\mathit{J}}_{\min }={\mathit{w}}_{opt}^H{\mathit{R}}_{\mathit{x}}{\mathit{w}}_{opt}^{}={\mathit{w}}_q^H{\mathit{R}}_{\mathit{x}}{\mathit{w}}_{opt}^{}$$

(9)

It can be clearly seen that ${\mathit{J}}_{\min }$ is just the minimum output power of SSI, which is denoted as $P_{o,\min }$. Whole interferences can be effectively removed in the SIEC output when they meet the requirements of the following two cases simultaneously: first, the number of antennas is larger than the number of interferences, i.e., system DOFs are high enough; and second, the distortionless constraint is set to the DOA of the desired signal. The MMSE will then be dominated by the output desired signal power,

$${\mathit{J}}_{\min }\simeq {\sigma }_{s_0}^2{\left|{\mathit{w}}_q^H\mathit{a}({\theta }_0)\right|}^2\triangleq P_s$$

(10)

where ${\sigma }_{s_0}^2$ and $P_s$ denote the input and output of the desired signal power in the SIEC, respectively.

In the case of SIEC, the numerator of the signal-interference-noise-ratio (SINR) expression is just the power of the desired signal in the output, as shown in Equation (12), and the denominator of the SINR is the power of the interference-plus-noise in the output, which can be calculated as,

$$\begin{array}{l}E\left\{{\left|{\mathit{w}}_{opt}^H\mathit{x}(m)-{\mathit{w}}_{opt}^H\mathit{s}(m)\right|}^2\right\}=P_{o,\min }-P_s\\ ={\mathit{w}}_q^H{\mathit{R}}_{\mathit{x}}{\mathit{w}}_{opt}^{}-{\sigma }_{s_0}^2{\left|{\mathit{w}}_q^H\mathit{a}({\theta }_0)\right|}^2\end{array}$$

(11)

The general expression for the optimum SINR can then be written as follows:

$${\mathrm{SINR}}_{opt}=\frac{E\left\{{\left|{\mathit{w}}_{opt}^H\mathit{s}(m)\right|}^2\right\}}{E\left\{{\left|{\mathit{w}}_{opt}^H\mathit{x}(m)-{\mathit{w}}_{opt}^H\mathit{s}(m)\right|}^2\right\}}=\frac{P_s}{P_{o,\min }-P_s}$$

(12)

2.4. Proposed Interference Mitigation Method

For the compact HFSWR, the current adaptive beamforming methods make it very difficult to obtain enough statistical samples of main-lobe transient interferences, i.e., SIEC [21,22]. In this study, the statistical samples of main-lobe transient interferences are first estimated by spatial auxiliary beams, which have similar interference characteristics as the main beam and blocks the desired signals. Next, the optimal secondary samples of transient interference in the time domain are obtained using the information geometry distance method. At last, the quadratic constraint approximation model is developed for solving the problem of a singular covariance matrix of interferences and maintaining the main- and side-lobe pattern, simultaneously.

2.4.1. Adaptive Optimal Sample Selection

According to [23,24], the suppression performance of the transient interferences is associated with the covariance matrix of the secondary samples. Therefore, distances-based covariance matrix estimation algorithms have drawn a lot of attention, in the considered space. The information geometry method is a kind of multi-dimensional information detector and does not require the knowledge of the statistical characteristics of the secondary data, which has been widely applied in radar signal processing [25]. The geodesic distance calculated by the Riemannian metric is a generalization of the real scalars’ mean on the matrix manifold and an important metric for measuring distribution similarity between different matrixes. The new selection strategy for optimal secondary samples based on Riemannian distance is proposed in this section.

The symmetrized log-determinant divergence (SLD) is typically used to measure the difference between two probability distributions. The divergence can act as a measure of the distance between two matrixes, which can be expressed as,

$$\begin{array}{l}{\mathit{D}}_{\mathrm{SLD}}({\mathit{R}}_0,{\mathit{R}}_r)=\frac{1}{2}\left({\mathit{D}}_{\mathrm{LD}}({\mathit{R}}_0,{\mathit{R}}_r)+{\mathit{D}}_{\mathrm{LD}}({\mathit{R}}_r,{\mathit{R}}_0)\right)\\ =\frac{1}{2}\left({\displaystyle \sum _{i=1}^n({\lambda }_i-\ln {\lambda }_i-1)}+{\displaystyle \sum _{k=1}^n({\eta }_k-\ln {\eta }_k-1)}\right)\end{array}$$

(13)

where ${\lambda }_i$ is the eigenvalue of the matrix ${\mathit{R}}_r^{-1}{\mathit{R}}_0$, ${\eta }_k$ is the eigenvalue of the matrix ${\mathit{R}}_0^{-1}{\mathit{R}}_r$, $n$ is the number of the eigenvalue, ${\mathit{D}}_{\mathrm{LD}}({\mathit{R}}_0,{\mathit{R}}_r)$ denotes the log-determinant divergence of test matrix ${\mathit{R}}_0$ and ${\mathit{R}}_r$, which can be written as,

$$\begin{array}{l}{\mathit{D}}_{\mathrm{LD}}({\mathit{R}}_0,{\mathit{R}}_r)=\mathrm{tr}\left({\mathit{R}}_r^{-1}({\mathit{R}}_0-{\mathit{R}}_r)\right)-\ln \det ({\mathit{R}}_r^{-1}{\mathit{R}}_0)\\ ={\displaystyle \sum _{i=1}^n({\lambda }_i-\ln {\lambda }_i-1)}\end{array}$$

(14)

where ${\mathit{R}}_0$ and ${\mathit{R}}_r$ express the covariance matrixes of the reference sample at $r_0$ range unit and the tested samples at $r,r=1,2\cdots ,R$ range units, respectively. ${\mathit{R}}_0$ and ${\mathit{R}}_r$ can be obtained with a complete model of the structured Hermitian positive-definite matrix. Their detailed calculation method is presented in [26]. It can be defined as follows:

$${\mathit{R}}_r^{}=\left[\begin{array}{cccc}{\mathit{c}}_0& {\mathit{c}}_1^{\ast }& \cdots & {\mathit{c}}_{\delta -1}^{\ast }\\ {\mathit{c}}_1& {\mathit{c}}_0& \ddots & ⋮\\ ⋮& ⋮& \ddots & {\mathit{c}}_1^{\ast }\\ {\mathit{c}}_{\delta -1}& \cdots & {\mathit{c}}_1& {\mathit{c}}_0\end{array}\right],r=0,1,2,\cdots R$$

(15)

where ${\mathit{c}}_k$ is the correlation coefficient, which can be calculated by averaging in the time domain, as follows:

$${\mathit{c}}_k=\frac{1}{\delta }{\displaystyle \sum _{\delta =0}^{\delta -1-k}\mathit{y}(\delta ){\mathit{y}}^{\ast }(\delta +k)},k\le \delta -1,\delta ={\delta }_d{\delta }_a$$

(16)

where $\mathit{y}(\delta )$ indicates the test data constructed by vectoring an angle-Doppler local region (ADLR) consisting of three adjacent angles and Doppler units contaminated by the transient interference. ${\delta }_a$ and ${\delta }_d$ denote the number of the selected angle units and Doppler units in the ADLR. According to [23,27], ${\delta }_d$ and ${\delta }_a$ are both set to 3. The ADLR in the compact HFSWR is shown in Figure 2.

The test matrixes of the reference unit and all test units can be calculated by Equation (15) in the range domain. All the SLDs between the reference and test units can be obtained after implementing Equation (14) $R$ times. The smaller the SLD, the higher the similarity between the two matrixes. In order to realize the optimal mitigation performance of interference, the tested range units are regarded as the secondary samples, also called training samples when their SLDs are less than the predefined threshold. In order to avoid that the training samples contain sea echoes consistent with the interested one (i.e., the geophysical signal of interest), the protected units must be used around the interested unit. The optimal training sample structure is shown in Figure 3 where its mathematical expression is defined as,

$${\overline{\mathit{x}}}_{r_0}(m_0)={\mathit{X}}_{aux}(m_0)\ast \mathbf{\Gamma }\in {ℂ}^{P\times \kappa }$$

(17)

where ${\mathit{X}}_{aux}\in {ℂ}^{P\times R}$ denotes the signal received by the whole auxiliary beams and range units at $m_0$ slow-time index. $\mathbf{\Gamma }$ expresses a transform matrix employed to choose the optimal training samples, which can be shown as,

$$\mathbf{\Gamma }=[\mathit{\tau }\left(r_1\right),\mathit{\tau }\left(r_2\right),\cdots ,\mathit{\tau }\left(r_{\kappa -1}\right),\mathit{\tau }\left(r_{\kappa }\right)]\in {ℝ}^{R\times \kappa }$$

(18)

$$\mathit{\tau }\left(r\right)=\underset{1,\cdots ,(r-1)}{\underbrace{[0,\cdots ,0,}}1\underset{(r+1),\cdots ,R}{\underbrace{,0,\cdots ,0}}{]}^T\in {ℝ}^{R\times 1},r=1,2,\cdots ,R$$

(19)

2.4.2. Rotating Spatial Beam Method

It is assumed that the echo of the compact HFSWR is occupied by the transient interference in the slow-time index $m_0$ and reference range unit $r_0$. The training samples of the reference unit have been obtained based on Section 3.1. Let $r_1,r_2,\cdots r_{\kappa }$ denote the range units of selected training samples in the range domain. The main-lobe beam and auxiliary beams data in the $r_1,r_2,\cdots r_{\kappa }$ range cells are ${\tilde{\mathit{x}}}_{r_0}(m_0)$ and ${\overline{\mathit{x}}}_{r_0}(m_0)$, respectively. The auxiliary beam data of the reference range unit are ${\overline{\overline{\mathit{x}}}}_{r_0}(m_0)$. They can, respectively, be given as,

$${\tilde{\mathit{x}}}_{r_0}(m_0)=\left[{\mathit{x}}_{r_1}^{{\theta }_0}(m_0),{\mathit{x}}_{r_2}^{{\theta }_0}(m_0),\cdots ,{\mathit{x}}_{r_{\kappa }}^{{\theta }_0}(m_0)\right]\in {ℂ}^{1\times \kappa }$$

(20)

$${\overline{\mathit{x}}}_{r_0}(m_0)=\left[\begin{array}{l}{\mathit{x}}_{r_1}^{{\theta }_1}(m_0),{\mathit{x}}_{r_2}^{{\theta }_1}(m_0),\cdots ,{\mathit{x}}_{r_{\kappa }}^{{\theta }_1}(m_0)\hfill \\ {\mathit{x}}_{r_1}^{{\theta }_2}(m_0),{\mathit{x}}_{r_2}^{{\theta }_2}(m_0),\cdots ,{\mathit{x}}_{r_{\kappa }}^{{\theta }_2}(m_0)\hfill \\ \hspace{1em}⋮\hspace{1em}\hspace{1em}\hspace{1em}⋮\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}⋮\hspace{1em}\hspace{1em}⋮\hfill \\ {\mathit{x}}_{r_1}^{{\theta }_P}(m_0),{\mathit{x}}_{r_2}^{{\theta }_P}(m_0),\cdots ,{\mathit{x}}_{r_{\kappa }}^{{\theta }_P}(m_0)\hfill \end{array}\right]\in {ℂ}^{P\times \kappa }$$

(21)

$${\overline{\overline{\mathit{x}}}}_{r_0}(m_0)={\left[{\mathit{x}}_{r_0}^{{\theta }_1}(m_0),{\mathit{x}}_{r_0}^{{\theta }_2}(m_0),\cdots ,{\mathit{x}}_{r_0}^{{\theta }_P}(m_0)\right]}^T\in {ℂ}^{P\times 1}$$

(22)

where $\kappa$ denotes the number of training samples in the range domain, $P$ indicates the number of auxiliary beams, ${\theta }_0$ denotes the DOA of the main-lobe beam, and ${\theta }_1,{\theta }_2,\cdots ,{\theta }_P$ indicate the DOAs of auxiliary beams.

Now, the output of our proposed method with transient interference canceled can be expressed as,

$${\mathit{x}}_{out}(r_0,m_0,{\theta }_0)={\mathit{x}}_{r_0}^{{\theta }_0}(m_0)-{\mathit{w}}_a^H{\overline{\overline{\mathit{x}}}}_{r_0}(m_0)$$

(23)

where ${\mathit{w}}_a={\mathit{R}}_{\overline{\mathit{x}}\overline{\mathit{x}}}^{-1}{\tilde{\mathit{x}}}_{r_0}$ is the adaptive weight calculated for suppressing main-lobe interference in the form of a Wiener solution. ${\mathit{R}}_{\overline{\mathit{x}}\overline{\mathit{x}}}$ denotes the self-correlation matrix of the secondary samples received from auxiliary beams, it can be expressed as,

$${\mathit{R}}_{\overline{\mathit{x}}\overline{\mathit{x}}}=\frac{1}{\kappa }{\displaystyle \sum {\overline{\mathit{x}}}_{r_0}(m_0){\overline{\mathit{x}}}_{r_0}^H(m_0)}$$

(24)

2.4.3. Optimal Weight Estimation Method

In the class algorithms of the beamformer, the array pattern may distort if the interference component is contained in the main-lobe beam, because it leads to errors in estimating the target angle and Doppler frequency [28]. Furthermore, if the interference has a very strong correlation in the space domain, the self-correlation matrix may be a singular matrix, which would suffer from SINR attenuation. To avoid the above two problems, a quadratic constraint approximation model [29] was developed for adaptively calculating the optimal diagonal loading (DL) level, which is added to the ${\mathit{w}}_a$ for obtaining the optimal weight. The problem can be formulated as follows:

$$\begin{array}{l}{\mathit{w}}_{opt}=\arg \min ({\mathit{w}}^H{\mathit{R}}_{\overline{\mathit{x}}\overline{\mathit{x}}}\mathit{w}+\gamma {\mathit{w}}^H\mathit{w})\\ =\arg \min {\mathit{w}}^H({\mathit{R}}_{\overline{\mathit{x}}\overline{\mathit{x}}}+\gamma \mathit{I})\mathit{w} \mathit{s}.\mathit{t}. {\mathit{U}}_{\Theta }^H\mathit{w}={\mathit{U}}_{\Theta }^H{\mathit{w}}_q\end{array}$$

(25)

where ${\mathit{w}}_q$ denotes the static weight to keep the direction of the main-lobe beam and the side-lobe level, ${\mathit{U}}_{\Theta }^H$ indicates the main-lobe subspace. The quadratic constraint of Equation (25) can be appropriately relaxed as the following fast estimation of DL level γ [29], it can be written as,

$$\begin{array}{l}{\gamma }_{opt}=\left\{\gamma |{\Vert {\mathit{w}}_{opt}^{}(\gamma )\Vert }_2^2\le T_0^2{\Vert {\mathit{w}}_q^{}\Vert }_2^2=T_0^2\right\}\\ =\left\{\gamma |{\Vert {\mathit{w}}_{\mathit{a}}^{}(\gamma )\Vert }_2^2\le T_0^2-{\Vert {\mathit{w}}_0\Vert }_2^2\right\}\end{array}$$

(26)

where $T_0$ denotes the allowed antenna gain loss factor, which is always slightly larger than 1. For example, if the antenna gain loss is 0.423 dB, the factor $T_0$ should be set to 1.05. The square norm of the weight ${\mathit{w}}_a(\gamma )$ in Equation (26) can be written as Equation (27), and it decreases monotonically as $\gamma$ increases:

$${\Vert {\mathit{w}}_a(\gamma )\Vert }_2^2={\mathit{w}}_a^H(\gamma ){\mathit{w}}_a(\gamma )={\tilde{\mathit{x}}}_{r_0}^H{({\mathit{R}}_{\overline{\mathit{x}}\overline{\mathit{x}}}+\gamma \mathit{I})}^{-2}{\tilde{\mathit{x}}}_{r_0}$$

(27)

Taking the derivative with respect to $\gamma$, then,

$$\frac{d{\Vert {\mathit{w}}_a(\gamma )\Vert }_2^2}{d\gamma }=-2{\tilde{\mathit{x}}}_{r_0}^H{({\mathit{R}}_{\overline{\mathit{x}}\overline{\mathit{x}}}+\gamma \mathit{I})}^{-3}{\tilde{\mathit{x}}}_{r_0}$$

(28)

It can be seen that Equation (28) is negative if the DL factor $\gamma \ge 0$, and the weight norm is monotonically decreasing in $\gamma$. Defining the positive factor $\ell (\gamma )$ evaluates the difference between the square norm of the actual weight ${\mathit{w}}_{opt}(\gamma )$ and the allowed square norm value. The factor $\gamma$ can be calculated by the following iterative method:

$$\begin{array}{l}\ell (\gamma )=\frac{{\Vert {\mathit{w}}_{opt}(\gamma )\Vert }_2^2}{T_0^2}=\frac{\Vert {\mathit{w}}_0\Vert }{T_0^2}=\frac{{\Vert {\mathit{w}}_{\mathit{a}}(\gamma )\Vert }_2^2}{T_0^2}\\ ={\mu }_1+{\mu }_2{\Vert {\mathit{w}}_a(\gamma )\Vert }_2^2={\mu }_1+{\mu }_2{\tilde{\mathit{x}}}_{r_0}^H{({\mathit{R}}_{\overline{\mathit{x}}\overline{\mathit{x}}}+\gamma \mathit{I})}^{-2}{\tilde{\mathit{x}}}_{r_0}\end{array}$$

(29)

where ${\mu }_1$ and ${\mu }_2$ are constant real values. The optimal factor $\gamma$ can be obtained using the following iterative formula:

$${\gamma }_{i+1}=\max (0,\ell ({\gamma }_i)({\gamma }_i+1)-1)$$

(30)

when $\left|{\gamma }_{i+1}-{\gamma }_i\right|\le \epsilon$, the iteration is stopped. To reduce the computational complexity in Equation (29), the ${({\mathit{R}}_{\overline{\mathit{x}}\overline{\mathit{x}}}+\gamma \mathit{I})}^{-2}$ can be approximately solved by,

$${({\mathit{R}}_{\overline{\mathit{x}}\overline{\mathit{x}}}+\gamma \mathit{I})}^{-2}=\mathit{V}{(\mathit{D}+\gamma \mathit{I})}^{-2}{\mathit{V}}^H$$

(31)

where $\mathit{D}$ and $\mathit{V}$ denote eigenvalues and corresponding eigenvectors of Eigen-decomposing the self-correlation matrix ${\mathit{R}}_{\overline{\mathit{x}}\overline{\mathit{x}}}$, respectively. Now, the ${\mathit{w}}_a$ in Equation (23) can be replaced by the following final optimal weight:

$$\begin{array}{l}{\mathit{w}}_{\mathit{a},opt}({\gamma }_{opt})={({\mathit{R}}_{\overline{\mathit{x}}\overline{\mathit{x}}}+{\gamma }_{opt}I)}^{-1}{\mathit{r}}_{\overline{\mathit{x}}\tilde{\mathit{x}}}\\ =\mathit{V}{(\mathit{D}+{\gamma }_{opt}\mathit{I})}^{-1}{\mathit{V}}^H{\mathit{r}}_{\overline{\mathit{x}}\tilde{\mathit{x}}}\end{array}$$

(32)

After processing every occupied range unit, STD, and beam unit, the proposed cancellation algorithm can be used to obtain a corresponding three-dimensional output result of compact HFSWR. Based on the adaptive selection method of the optimal training samples in range unit (i.e., time-domain process) and the main-lobe interference mitigation with auxiliary beams (i.e., space-domain process), it guarantees that there are few sea echo components and the estimated transient interferences are as similar to those of main-lobe beam as possible in ${\overline{\mathit{x}}}_{r_0}(m_0)$, simultaneously. Therefore, the suppression result ${\mathit{x}}_{out}(r_0,m_0,{\theta }_0)$ will approximate the optimal value.

2.4.4. The Procedure of the Proposed Algorithm

The algorithm flow can be summarized as follows:

(1)

Input: compact HFSWR original echo data preprocessed by pulse compression, coherent integration and digital beamforming (DBF), location index of the transient interferences in the slow-time region;

(2)

Initialize the allowed error $\eta$, the number of auxiliary beams $P$, and the corresponding directions ${\theta }_1,{\theta }_2,\cdots ,{\theta }_P$;

(3)

Suppose the range unit for the cell under test is $r_0$, calculate the geometric distances ${\mathit{D}}_{\mathrm{SLD}}({\mathit{R}}_0,{\mathit{R}}_r)$ for all the training data in the range domain using Equation (13);

(4)

Set two guard units to prevent the target self-elimination and sort ${\mathit{D}}_{\mathrm{SLD}}$ in ascending order;

(5)

Select the $\kappa =2P$, which indicates as the number of training samples that correspond to $\kappa$ lowest values of ${\mathit{D}}_{\mathrm{SLD}L\mathit{D}}$;

(6)

Calculate the self-correlation matrix ${\mathit{R}}_{\overline{\mathit{x}}\overline{\mathit{x}}}$ and the crosscorrelation matrix ${\mathit{r}}_{\overline{\mathit{x}}\tilde{\mathit{x}}}$;

(7)

Calculate the DL level ${\gamma }_i$ using equations from Equation (29) to Equation (31) until the predefined convergence criteria is satisfied;

(8)

Calculate the final optimal weight using Equation (32);

(9)

Calculate the final output using Equation (23);

(10)

Repeat step 2 to step 9 for processing the whole occupied range units, STD units and beam units;

(11)

Output: the same three-dimensional matrix as input echo data, but interferences have been removed.

3. Results

In this section, simulated and experimental data are applied to validate the suppression feasibility of our proposed method for the compact HFSWR. The collected experimental data as shown in Figure 1 and the simulation data are conducted according to [4]. The robust principal component analysis (RPCA) method proposed in [4], the bicriterion optimization suppression (BOS) method proposed in [14], and the complex empirical mode decomposition (CEMD) method proposed in [5], main-lobe cancellation method (MLCM) proposed in [15], are applied for comparing the suppression performance of transient interferences. Distinct from the HF skywave over-the-horizon radar, the measured data of the compact HFSWR have a sufficiently long coherent accumulation time (hundreds of seconds). Considering the accuracy and complexity of the algorithm, the output STD data of the proposed method is processed using the windowed fast Fourier transform method to estimate the Doppler information of the submerged sea echoes. The windowed function is the Hanning window.

3.1. Simulation Data

As stated in [4], the simulation parameters are set as follows: the signal PRI is T = 25 ms, the CPI contains M = 512 pulses, the radar working frequency is 15 MHz, the number of the array elements is 8, the additive noise is assumed to be temporally and spatially white, the target’s Doppler frequency is −4.3 Hz, the sea-to-noise ratio is set to 20.23 dB, the target-to-noise ratio is set to 4.62 dB, and two simulated transient interferences are introduced at time instants of 1.2 s and 5 s. The STD amplitude of the selected range unit data with two transient interferences is shown in Figure 4a, where it can be found that the interference amplitudes change fiercely and are clearly larger than the target and sea echoes in the slow time index. The Doppler profile of the DBF result corrupted by two transient interferences is shown in Figure 4b. Clearly, the background noise level is significantly increased, and the sea-to-noise ratio and target-to-noise ratio suffer from serious attenuation. The transient-removed Doppler profile obtained by the whole compared methods is shown in Figure 4b. It can be seen that all the interference is suppressed, as expected, which confirms that our proposed method has at least 7 dB improvement in the sea-to-noise ratio and target-to-noise ratio compared with DBF, RPCA, BOS, CEMD, and MLCM. The MLCM suffers dramatic performance degradation owing to the loss of DOFs in compact array. This is because the MLCM needs enough DOFs to form deep nulling filter avoiding signal self-cancellation and obtain the sufficient secondary samples of the main-lobe interference. Thus, we can conclude that our method outperforms the whole methods on transient interference excision and is also effective in noise suppression.

3.2. Training Sample Selection for the Experimental Data

In a HF radar system, the distribution characteristics of interference and clutter are significantly heterogeneous [15,27]. Based on the above analysis, the mitigation of interference performance depends on the accuracy of the selected training samples, and the SLDs can reflect the homogeneity between the reference unit and test units. When the SLD values are smaller, the homogeneity is stronger. Consequently, the SLDs are first calculated between the reference unit and every test unit using Equation (13) for adaptively selecting the most homogeneous units as training samples. In our proposed method, the number of SLDs should be twice as large as the auxiliary beams for approaching the optimal performance of 3 dB loss. The number of auxiliary beams is set to 21, which was proved in our previous study [24]. Based on this, 42 training samples with the lowest SLDs need to be chosen. Figure 5a displays the normalized SLD results of LCI where its reference unit $r_0$ is set to the 12th unit and the STD index is the 881th. The reference unit of STI is also set to the 12th unit. The STD index is the 881th corresponding to the left transient interference in Figure 1c. The normalized SLD results are shown in Figure 5b. In Figure 5, the red point denotes the reference unit, the cyan points are the projected units, the blue points denote the selected training samples, the gray points are the heterogeneous units.

3.3. Doppler Profile of Transient Interference Suppression for the Experimental Data

In this section, the transient-removed Doppler profiles with the real-measured data are shown in Figure 6. Based on [5], the constant scale factor is set to 9 dB. The subarray number is set to 6 for forming the block filter in the MLCM. The numbers of snapshots and signal subspace in the BOS method are then set to 1000 and 6, respectively. Figure 6a,c show the LCI and STI mitigation results of the whole comparison methods at 2nd range unit, the 12th range unit data are shown in Figure 6b,d. As shown in Figure 6, the transient interference significantly raises the background noise, and even submerges the weak sea echo of the negative frequency. Comparing Figure 6a,b with Figure 6c,d, it is obvious that our proposed method can not only mitigate the interference but also better protect the sea echoes, including the weak and strong echoes, and it is the only one that well works for both the LCI and STI simultaneously. Therefore, our proposed method makes the effective inversion of ocean dynamics parameters using the submerged sea echoes possible. In addition, the proposed method is clearly superior to the MLCM in compact HFSWR. It can demonstrate that our method can more effectively solve the problem of few DOFs and main-lobe interference, simultaneously.

3.4. Range-Doppler Map of Transient Interference Suppression for the Experimental Data

Figure 7 and Figure 8 display the range-Doppler maps after transient interference excision with different methods. The methods are the RPCA, BOS, our proposed method and CEMD. The corresponding range-Doppler maps corrupted by the transient interferences are shown in Figure 1b,d. From Figure 1, it is clear that the energy of the interference is spread over the entire range of units and Doppler units, which completely submerges the weak sea echoes. Figure 7 and Figure 8 show range-Doppler maps of the LCI and STI, respectively. The spectrum of the LCI-removed signal obtained by our proposed method is shown in Figure 7c, where the sea echoes are more visible and the interference is suppressed more clearly than the other methods. As shown in Figure 8, it can be clearly observed that the proposed method can still achieve excellent performance of background level reduction and SINR increment. Furthermore, its performance is almost the same as the RPCA method.

4. Discussion

In Section 3, the configuration of the experiment was described in detail, and the mitigation performance of transient interference using simulation and experiment data was tested. Next, the SINR improvement experiment against different range units was carried out and the various methods were compared. The SINR comparison of sea echoes with varying range units as shown in Figure 9 and Figure 10, and the preset parameters of the different methods are the same as in Section 3. In this section, the SINR improvement of interference mitigation methods is evaluated using a quantized value (MSINR) obtained by calculating the mean value of the output SINR of all the range units. The experimental part is analyzed in detail below. As shown in Figure 1b,d, the sea echoes of the negative frequency are weak and only exist within the 22nd range unit. Thus, the number of range units was set to 22 for the sea echoes of the negative frequency.

4.1. SINR Improvement for the LCI Data

The output SINRs of the compared methods after suppressing the transient interference are shown in Figure 9. The result of sea echo on the LCI negative frequency is shown in Figure 9a. It can be seen that the MSINR of our proposed method has increased as much as 10.88 dB, 3.79 dB, 5.03 dB, and 3.18 dB compared with the DBF, RPCA, BOS and CEMD, respectively. The MSINR improvements of approximately 16.29 dB, 7.11 dB, 5.99 dB, and 18.96 dB are presented in Figure 9b, where it shows the result of the LCI positive frequency. Moreover, we can also observe from Figure 9, in general, that our proposed method almost has the optimal output SINRs in all range units. These advantages come from the accurate estimation and construction methods of statistical samples, secondary samples, and interference covariance matrix developed in this study, instead of requiring the transient interference to meet the sparsity and low-frequency property and so on.

4.2. SINR Improvement for the STI Data

In this way, the output SINRs of STI sea echoes after interference mitigation are shown in Figure 10a,b. Figure 10a displays the negative frequency result, the MSINRs are 11.25 dB, −2.04 dB, 4.78 dB, 2.35 dB corresponding to DBF, RPCA, BOS and CEMD methods, respectively. The range profile of positive frequency is shown in Figure 10b, the MSINRs are 13.46 dB, −2.41 dB, 5.31 dB, 18.28 dB. As shown in Figure 10, our proposed method is also significantly better than DBF, BOS, and CEMD methods, but slightly less than the RPCA method. This is because when transient interferences are sparse in the STD, RPCA can more effectively separate two components of the sparse interferences and low-rank sea echoes. Therefore, while the results are encouraging, it is clear that further work is needed. Potential sources leading to performance degradation are listed below.

Case 1: Further improving the estimation accuracy of training samples for more effectively suppressing transient interference.

Case 2: More accurately constructing the interference covariance matrix for avoiding the target self-cancellation.

Case 3: Developing a phase calibration method for adjusting the nonlinear phase of complex signals caused by the interference mitigation procedure.

5. Conclusions

In this paper, a novel space-time cascaded processing algorithm is proposed for LCI and STI suppression in the background of the compact HFSWR, based on the assumption that the transient interferences have been accurately located using published literature. The proposed method is implemented with the following steps: the rotating spatial beam method is developed for estimating the statistical samples of the main-lobe transient interference, and the optimal secondary samples are constructed by using an adaptive selection strategy to extract the homogeneous samples as the interested unit based on information geometry distance, and the estimation procedure of the optimal weight is considered as solving a quadratic constraint approximation problem when the interference covariance matrix is singular. The simulation and experimental results indicate that the proposed method can effectively improve the detection performance of the submerged sea echoes, and the sea echo of negative frequency has at least an average 3.01 dB output SINR increment than the RPCA, BOS, and CEDM methods for the LCI.