# Spread Sea Clutter Suppression in HF Hybrid Sky-Surface Wave Radars Based on General Parameterized Time-Frequency Analysis

Authors: Zhuo-qun Wang; Ya-jun Li; Jun-nan Shi; Peng-fei Wang; Ai-hua Liu; Xin-fan Xia; De-hong Chen
Journal: International Journal of Antennas and Propagation
 (2020)

---

## Abstract

In high-frequency (HF) hybrid sky-surface wave radar, the first-order sea clutter broadening is severe under the action of ionospheric phase disturbance and bistatic angles. In this paper, a cascaded method is described to suppress the spread sea clutter. Firstly, the radar configuration and sea clutter broadening model are introduced based on the newly developed integrated HF sky-surface wave experimental system. In the cascaded processing method, a new ionospheric decontamination method based on general parameterized time-frequency (GPTF) analysis is proposed to estimate or correct the ionospheric phase distortion with large amplitude. Then, the forward-backward linear prediction (FBLP) algorithm is used to suppress the spread sea clutter caused by bistatic angle. Simulation results show that such ionospheric decontamination method based on GPTF is helpful for the large-amplitude ionospheric contamination when the target masking effect happens even after ionospheric phase decontamination. Finally, the proposed method is examined by the measured data. Experimental results indicate that the proposed method can well suppress the broadening sea clutter for HF hybrid sky-surface wave radars.


---

## Body

## 1. Introduction

Ship detection is an important mission of the HF over-the-horizon radar. Based on the propagation mode associated with sky-wave transmitting and surface wave receiving, HF hybrid sky-surface wave radar not only maintains the superiority of HF sky-wave radar which has a long detection range and wide coverage, the advantage of HF surface wave radar which has a stable propagation channel but also keeps a good invisibility and anti-interference ability [1–4]. However, the sea clutter under this combined propagation mode with sky-wave transmitting and shore-based receiving is influenced by composite factors such as the bistatic angle and ionospheric phase disturbance. Among the influencing factors of the spectrum spread of sea clutter, it is unclear whether the ionospheric disturbance or the bistatic angle dominates. In addition, sea clutter is also affected by marine dynamics factors. Therefore, under this hybrid propagation mode radar, the spectrum of first-order sea clutter is complicated [5, 6].Based on the hybrid operating mode of sky-wave transmitting and ground-wave receiving, HF hybrid sky-surface wave radar is expected to improve detection probability for ships by overcoming defects of existing sky-wave and ground-wave OTH radar and making complementary advantages of both. Ionosphere is a dispersive, hierarchical, and nonstationary medium. Unfortunately, the ionospheric disturbance often causes the sea clutter spectrum and target to spread in the frequency domain, rendering extended coherent integration pointless. On the other hand, HF hybrid sky-surface wave radar is actually a bistatic radar system. The receiving beam width is wide for the array aperture of our radar experimental system. Thus, the sea clutter spectrum will show different broadening characteristics in the different resolution cells. What is worse, ionospheric phase disturbance further contributes to the sea clutter spectrum broadening. Thus, the broadening of the first-order sea clutter spectrum is very severe under the influence of ionosphere and bistatic angles. In the process of cascade processing of broadened sea clutter in this paper, ionospheric phase decontamination was first processed. At this time, decontaminated broadening sea clutter is only caused by bistatic angles. Therefore, the broadening sea clutter suppression was processed after decontamination.For the ionospheric phase path disturbance, the method of extracting the contamination function is often used for compensation. The current ionospheric phase disturbance suppression method extracts the frequency modulation function by estimating the instantaneous frequency variation of the broadened echo spectrum, thereby obtaining the ionospheric phase disturbance correction function, and then correcting the echo signal by using the obtained correction function. It can sharpen the broadened echo spectrum and improve the target detection performance of the radar. For ionospheric contamination with small amplitude, the maximum entropy spectrum estimation method, the phase gradient algorithm (PGA) method [7], the minimum entropy searching method, the eigendecomposition (ED) method, and the piecewise polynomial phase modeling (PPPM) method [8, 9] are typical methods. These methods require that the broadened Bragg peak is extracted to estimate the instantaneous frequency. Time-frequency analysis is a more effective method for ionospheric decontamination, especially large-amplitude ionospheric contamination. Currently, nonparametric time-frequency analysis methods are mainly used, such as short-time Fourier transform (STFT), Winger–Ville distribution (WVD), pseudo-WVD distribution (PWVD) method, and smooth pseudo-WVD distribution (SPWVD) method. However, those methods are sensitive to the signal-to-noise ratio and are affected by cross-term interference and resolution limitations. Especially for the ionospheric contamination with large amplitude, these methods often fail. Currently, the main technical approach to detect the ship target is the sea clutter suppression method. Representative algorithms include sea clutter cancellation method, subspace class method, and singular value decomposition (SVD) method, but these algorithms are difficult to effectively suppress nonstationary broadened sea clutter [8–10].In 1997, Melyanovski [1] mentioned the HF hybrid sky-surface wave radar system and explored the feasibility. Riddolls [3, 4] discussed the ship target detection performance of this radar and analyzed the limits on the radar resolution capability in 2007 and 2008. Zhao et al. [5] presented an experimental study of HF passive bistatic radar via hybrid sky-surface wave mode based on Digital Radio Mondiale digital amplitude modulation broadcasting in 2013 and analyzed the target detection performance. In 2014, Li et al. [6] analyzed the characteristics of frequency shift and broadening sea clutter based on the newly developed integrated HF sky-surface wave experimental system. Experimental results show that some sea clutter broadening is serious and the ship target is difficult to effectively detect. In 2015, Li et al. [11] proposed a cascaded processing method for broadening sea clutter by S2-method and FBLP. The S2 transform method is a nonparametric method, which can suppress cross-term interference to a certain extent, but the resolution and the time-frequency characteristics of the signals cannot be accurately described. In this paper, a cascaded method for ionospheric decontamination and sea clutter suppression to enable detection of ship target is proposed. The method utilizes an important property of the GPTF and FBLP to suppress the broadening sea clutter. GPTF is a parameterized time-frequency analysis method. Compared with the S2-method, it has better resolution and no cross-term interference. And a novel estimator for the instantaneous frequency of the Bragg lines based on GPTF is proposed.The content of this paper is organized as follows. Firstly, the radar configuration and sea clutter broadening model are described based on the newly developed integrated HF sky-surface wave experimental system. Secondly, a cascaded processing method for ionospheric decontamination and sea clutter suppression is presented. In this method, the time-frequency analysis method based on GPTF is proposed to correct the ionospheric phase contamination with large amplitude. Then, the FBLP algorithm is used to suppress the broadening bistatic sea clutter caused by bistatic angle, which is based on the prior knowledge of distribution characteristics and the multidimensional feature of first-order sea clutter [12]. Finally, the proposed method is examined by the measured data. Compared with the traditional nonparametric time-frequency analysis method, the ionospheric pollution function extracted by the GPTF method is more accurate. Experimental results indicate that the proposed method can well suppress the broadening sea clutter for HF hybrid sky-surface wave radars. The research in this dissertation will provide the theory and technology foundation to detect the target over the horizon under the broadening sea clutter in HF hybrid sky-surface wave radar.


## 2. First-Order Sea Clutter Broadening Characteristics and Model

### 2.1. System Layout

Based on the propagation mode associated with sky and surface wave, the HF hybrid sky-surface wave radar not only maintains the superiority of HF sky-wave radar which has a long detection range and wide coverage, the advantage of HF surface radar which has a stable propagation channel but also keeps a good invisibility and anti-interference ability. The system layout of the HF hybrid sky-surface wave radar is shown in Figure1 [6]. The beams generated by the sky-wave radar are reflected by ionosphere to monitored sea area (path P1andP2). Then, the sea echo of the monitored sea area arrives at the surface wave radar station by ground-wave diffraction (path Rr1andRr2). In Figure 1, TR = L is the baseline length between the transmitting and receiving stations; Rri is the distance between the target and receiving station; θri is the angle between the Rri and L; h is the ionospheric height; D is the ground distance; and βHi is bistatic angle.Figure 1
Geometry and layout diagram of HF sky-surface wave radar.


### 2.2. First-Order Sea Clutter Bragg Frequency

Thus, the total Bragg frequency expression for HF hybrid sky-surface wave radar can be written as follows [6, 13]:(1)fBragg=±0.102×10−3f0cosβi2cosΔi+12+fiono,where fiono is the Doppler shift caused by ionospheric disturbance, βi is bistatic angle, and Δi is the grazing angle.Therefore, the maximum covering range of the first-order sea clutter Bragg frequency can be obtained by the following formula [6, 13]:(2)−0.102×10−3f0⋅cosΔimin+12⋅12+Rri−LcosθrimaxRri2+L2−2RriLcosθrimax+fiono_min,−0.102×10−3f0⋅cosΔimax+12⋅12+Rri−LcosθriminRri2+L2−2RriLcosθrimin+fiono_max+0.102×10−3f0⋅cosΔimax+12⋅12+Rri−LcosθriminRri2+L2−2RriLcosθrimin+fiono_min,+0.102×10−3f0⋅cosΔimin+12⋅12+Rri−LcosθrimaxRri2+L2−2RriLcosθrimax+fiono_max,where fiono_max and fiono_min represent the maximum and minimum Doppler shift value caused by ionosphere. Thus, the Bragg frequency of first-order sea clutter is a function of frequency, azimuth, range, and ionospheric state.


### 2.3. First-Order Sea Clutter Broadening Model

For a sea clutter scattering cell at a certain range and beam direction in HF hybrid sky-surface wave radar, the broadening sea clutter model can be written as follows [6, 13], considering the combined actions of bistatic angle and ionosphere:(3)st=∑i=1Kait⋅expj2πfBraggt+exp−j2πfBraggt⋅expjϕit,where ait is the amplitude of sea clutter signal, which meets Gaussian distribution; fBraggi=1,…,K is the K complex frequency component within the Bragg frequency band, and fBragg_min≤fBragg≤fBragg_max; and ϕit is the ionospheric phase contamination function.


## 2.1. System Layout

Based on the propagation mode associated with sky and surface wave, the HF hybrid sky-surface wave radar not only maintains the superiority of HF sky-wave radar which has a long detection range and wide coverage, the advantage of HF surface radar which has a stable propagation channel but also keeps a good invisibility and anti-interference ability. The system layout of the HF hybrid sky-surface wave radar is shown in Figure1 [6]. The beams generated by the sky-wave radar are reflected by ionosphere to monitored sea area (path P1andP2). Then, the sea echo of the monitored sea area arrives at the surface wave radar station by ground-wave diffraction (path Rr1andRr2). In Figure 1, TR = L is the baseline length between the transmitting and receiving stations; Rri is the distance between the target and receiving station; θri is the angle between the Rri and L; h is the ionospheric height; D is the ground distance; and βHi is bistatic angle.Figure 1
Geometry and layout diagram of HF sky-surface wave radar.


## 2.2. First-Order Sea Clutter Bragg Frequency

Thus, the total Bragg frequency expression for HF hybrid sky-surface wave radar can be written as follows [6, 13]:(1)fBragg=±0.102×10−3f0cosβi2cosΔi+12+fiono,where fiono is the Doppler shift caused by ionospheric disturbance, βi is bistatic angle, and Δi is the grazing angle.Therefore, the maximum covering range of the first-order sea clutter Bragg frequency can be obtained by the following formula [6, 13]:(2)−0.102×10−3f0⋅cosΔimin+12⋅12+Rri−LcosθrimaxRri2+L2−2RriLcosθrimax+fiono_min,−0.102×10−3f0⋅cosΔimax+12⋅12+Rri−LcosθriminRri2+L2−2RriLcosθrimin+fiono_max+0.102×10−3f0⋅cosΔimax+12⋅12+Rri−LcosθriminRri2+L2−2RriLcosθrimin+fiono_min,+0.102×10−3f0⋅cosΔimin+12⋅12+Rri−LcosθrimaxRri2+L2−2RriLcosθrimax+fiono_max,where fiono_max and fiono_min represent the maximum and minimum Doppler shift value caused by ionosphere. Thus, the Bragg frequency of first-order sea clutter is a function of frequency, azimuth, range, and ionospheric state.


## 2.3. First-Order Sea Clutter Broadening Model

For a sea clutter scattering cell at a certain range and beam direction in HF hybrid sky-surface wave radar, the broadening sea clutter model can be written as follows [6, 13], considering the combined actions of bistatic angle and ionosphere:(3)st=∑i=1Kait⋅expj2πfBraggt+exp−j2πfBraggt⋅expjϕit,where ait is the amplitude of sea clutter signal, which meets Gaussian distribution; fBraggi=1,…,K is the K complex frequency component within the Bragg frequency band, and fBragg_min≤fBragg≤fBragg_max; and ϕit is the ionospheric phase contamination function.


## 3. Ionosphere Decontamination Method Based on GPTF

### 3.1. Basic Principle of GPTF

Ionospheric disturbances have the characteristics of change rapidly over time. Compared with the nonparametric time-frequency analysis method, the parameterized time-frequency analysis method selects the appropriate kernel describing the nonstationary signal by introducing the model of the prior signal information. When the kernel form is consistent with the analyzed signal, the time-frequency resolution can be effectively improved. The typical methods have adaptive chirp wavelet decomposition, atomic decomposition, and polynomial Fourier transform, etc., but such methods mainly use polynomial kernels, which are not suitable for analyzing strong time-varying nonstationary signals that change faster with time.The method of short-time Fourier transform (STFT) is to window the signal based on the Fourier transform. The default window signal is approximately stationary. It can simultaneously describe the signal in the time domain and the frequency domain, thus reflecting the variation of the signal spectrum over time. Since the fixed-length window function cannot capture the change of the signal frequency in time, the time-frequency representation of the STFT is poorly concentrated, and the time-frequency characteristics of the signals cannot be accurately described. The generalized parameterized time-frequency analysis method uses the frequency rotation operator to rotate the time-frequency characteristics of the nonstationary signal by introducing the frequency rotation and translation operator, so that the signal tends to be stable, and then, the STFT is used for the rotated signal, and finally, the frequency is utilized. The translation operator shifts the signal time-frequency characteristics to the true ridge position. Since the signal analyzed by STFT is an approximate stationary signal, the generalized parameterized time-frequency analysis can effectively improve the time-frequency resolution, and there is no cross-term interference.Assuming that the instantaneous frequency of the signal is an arbitrary functionfmt, the complex signal form at this time is srt=ej2π∫fmtdt. The generalized parameterized time-frequency analysis of the signal can be expressed as follows [14]:(4)TFt0,f;P=∫−∞∞s¯τgσ∗τ−t0e−j2πfτdτs¯τ=srτ⋅ΦPRτ⋅Φt0,PSτΦPRτ=e−j2π∫κPτdτΦt0,PSτ=ej2⋅π⋅τ⋅κPt0,wherein the instantaneous frequency of srτ is an arbitrary function, κPτ is a transform kernel representing of the generalized parametric time-frequency analysis, P is a transform kernel parameter, ΦPRτ is a frequency rotation operator, Φt0,PSτ is a frequency shift operator, and gστ is a window function of the time-frequency transform.According to formula (4), the basic process of GPTF analysis can be described as follows. Firstly, the signal is rotated in the time-frequency plane, that is, the instantaneous frequency fmτ of the signal is subtracted from κPτ, and then, the signal is translated on the time-frequency plane by the frequency increment κPt0. Finally, the window function is used to make a STFT on the signal. When κPτ=fmτ, the signal forms after the signal srτ is rotated and translated as shown in the following equation:(5)s¯τ=srτ⋅ΦPRτ⋅Φt0,PSτ=ej2π∫fmτdτ⋅e−j2π∫κPτdτ⋅ej2π·τ·κPt0=ej2π⋅τ⋅κPt0.At this time, the instantaneous frequencyκPt0 of the signal s¯τ does not change with time, and the signal s¯τ is a stationary signal. At this time, the STFT is performed on the signal s¯τ, and the time-frequency concentration is significantly improved compared to the nonstationary signal. It can be seen from (4) and (5) that the frequency rotation operator and the translation operator are determined by selecting appropriate transformation kernel parameters, and the nonstationary signal is rotated and translated to stabilize the signal for the purpose of improving the time-frequency resolution by the determined frequency rotation operator and the frequency translation operator. It can be seen that the time-frequency resolution of the GPTF analysis is closely related to the accuracy of the transform kernel parameters. From the above analysis, the key to GPTF analysis is to estimate the parameters of the transform kernel.


### 3.2. Estimation Method for Transforming Kernel Parameters

According to the introduction of the principle of GPTF in Section3.1, the working mechanism is actually the STFT based on the frequency rotation operator and the translation operator. When the selected transform kernel matches the time-frequency characteristic of the analyzed signal, the rotated signal is close to the stationary signal, the time-frequency resolution is high after STFT. The description of the time-frequency characteristics of the stationary signal is accurate. It can be seen that the accuracy of the transform kernel parameter determines the performance of the GPTF method, which affects its performance in analyzing nonstationary signals. In [14], a method for estimating transform kernel parameters using the principle of time-frequency feature approximation is proposed. The time-frequency feature approximation principle indicates that the more accurate the transform kernel parameter estimation, the higher the time-frequency concentration. According to the time-frequency feature approximation principle, if the transform kernel parameter estimation is more accurate, the time-frequency concentration is higher. The principle of time-frequency approximation can be represented in Figure 2.Figure 2
Time-frequency approximation principle.It can be seen from Figure2 that the GPTF analysis method achieves the iterative optimization of the transform kernel parameters by combining the precision of the transform kernel parameters with the time-frequency concentration. The above principle is based on the generalized parametric time-frequency analysis to cyclically approximate and solve the signal time-frequency characteristics, so as to obtain the most suitable transform kernel parameters. According to the principle of time-frequency approximation, the detailed steps of transforming the kernel parameter estimation are as follows:Step1:determining the transform kernel form according to the signal priori form.Step 2: the number of iterationsi=0, initializing the transformation kernel parameters P0.Step4: performing parameterized time-frequency analysisTFt0,f;Pi on the signal according to equation (4).Step 4: detecting the peak ridge linef¯mit in the time-frequency diagram.Step 5: according to the ridge line form, select a suitable fitting method, and fit the ridge line and estimate transformation kernel parameterP˜i.Step 6: the number of iterations is increased by one, and the transformation kernel parameter is updatedPi=P˜i.Step7: calculating a termination conditionΛ=∫f˜mit−f˜mi−1tdt/∫f˜mitdt, where f˜mit is the peak ridge of the ith iteration and f˜mi−1t is the peak ridge of the i−1th iteration, which needs to be noted that f˜m0t=0.Step 8: setting the termination condition parameterζ. The size of this value depends on the situation. If it is too large, the calculation amount is too large, and if it is too small, the accuracy of the extracted time-frequency line is not high.Step 9: comparing the size of Λ andζ, when ζ > Λ, return to step 3. When ζ ≤ Λ, perform step 10.Step 10: outputting a transformation kernel parameter to obtain a time-frequency curve.According to the above steps of transforming the kernel parameter estimation, the flowchart of the design is shown in Figure3.Figure 3
The flowchart of transform kernel parameter estimation.


### 3.3. Ionospheric Decontamination Based on GPTF

Based on the above analysis process and steps, the derivation process of ionospheric decontamination based on GPTF is as follows:Step 1: assume that the signal after the ionospheric contamination is expressed as follows [7–9]:(6)sp=expjm2πsin2πf0texpj2πf1t,wherem is the amplitude coefficient, f0 is the modulation frequency, and f1 is the fixed frequency.Step 2: derive the signal and obtain the Doppler curve form, which can be expressed as follows:(7)ft=m2πf0cos2πf0t+f1.Since the cosine Doppler curve of the signal has periodicity, the ridge line is fitted using the Fourier series for the curve. The Doppler instantaneous frequency form is known. Since the Doppler curve is known to be a cosine form by prior knowledge, the transform kernel of the applied harmonics is as follows:(8)ft=∑m=1∞bmcosm2πfct+b0.There is a certain relationship between the cosine series coefficient of the curve and its Fourier transform coefficient. Consider the Fourier transform to find the coefficients of the cosine series of each order. The discrete form of the Doppler curve parameter estimation method is derived as follows.Radar sampling intervalΔt=1/PRF; equation (8) discretization is expressed as follows:(9)fn⋅Δt=fdn⋅Δt+b0,where fdn⋅Δt=∑m=1∞bmcosm⋅2πfcnΔt.After equation (9) is truncated by adding a rectangular window, the spectrum of each order cosine series becomes a sinc function form. In this case, in order to eliminate the influence of the sinc side lobes of the cosine series spectrum of each order on the spectrum of other order cosine series, the rectangular window length is set to TT=1/fc,2/fc,3/fc. In this way, the peak of the cosine series spectrum of other orders will be located at the zero point of the side lobe of the spectrum of each order cosine series, so as to eliminate the interference between the spectra. Therefore, in the latter experiment, the accumulation time is required to be an integer multiple of the period. When the accumulation time is not an integer multiple, the signal can be truncated, so this method requires at least one periodic signal. Without loss of generality, set TT=i/fc, at this time N=i/fcΔt. Similarly, the Fourier transform of the N point in the harmonic form fdn⋅Δt in (9) is obtained as follows:(10)Fdk=∑n=0N−1fdnΔte−j2π/Nkn=∑n=0N−1∑m=1∞bmcosm2πiNne−j2π/Nkn=∑m=1∞bm∑n=0N−1cos2πiNn⋅micos2πNn⋅k−j∑m=1∞bm∑n=0N−1cos2πNn⋅misin2πNn⋅k.According to the characteristics of the orthogonal trigonometric function set, formula (10) can be reduced to(11)Fdk=bmN2,k=mi,0,else.The transformation kernel parameter can be obtained as follows:(12)bm=2NrealFdkk=mi,b0=λ22realFd1NΔt.


## 3.1. Basic Principle of GPTF

Ionospheric disturbances have the characteristics of change rapidly over time. Compared with the nonparametric time-frequency analysis method, the parameterized time-frequency analysis method selects the appropriate kernel describing the nonstationary signal by introducing the model of the prior signal information. When the kernel form is consistent with the analyzed signal, the time-frequency resolution can be effectively improved. The typical methods have adaptive chirp wavelet decomposition, atomic decomposition, and polynomial Fourier transform, etc., but such methods mainly use polynomial kernels, which are not suitable for analyzing strong time-varying nonstationary signals that change faster with time.The method of short-time Fourier transform (STFT) is to window the signal based on the Fourier transform. The default window signal is approximately stationary. It can simultaneously describe the signal in the time domain and the frequency domain, thus reflecting the variation of the signal spectrum over time. Since the fixed-length window function cannot capture the change of the signal frequency in time, the time-frequency representation of the STFT is poorly concentrated, and the time-frequency characteristics of the signals cannot be accurately described. The generalized parameterized time-frequency analysis method uses the frequency rotation operator to rotate the time-frequency characteristics of the nonstationary signal by introducing the frequency rotation and translation operator, so that the signal tends to be stable, and then, the STFT is used for the rotated signal, and finally, the frequency is utilized. The translation operator shifts the signal time-frequency characteristics to the true ridge position. Since the signal analyzed by STFT is an approximate stationary signal, the generalized parameterized time-frequency analysis can effectively improve the time-frequency resolution, and there is no cross-term interference.Assuming that the instantaneous frequency of the signal is an arbitrary functionfmt, the complex signal form at this time is srt=ej2π∫fmtdt. The generalized parameterized time-frequency analysis of the signal can be expressed as follows [14]:(4)TFt0,f;P=∫−∞∞s¯τgσ∗τ−t0e−j2πfτdτs¯τ=srτ⋅ΦPRτ⋅Φt0,PSτΦPRτ=e−j2π∫κPτdτΦt0,PSτ=ej2⋅π⋅τ⋅κPt0,wherein the instantaneous frequency of srτ is an arbitrary function, κPτ is a transform kernel representing of the generalized parametric time-frequency analysis, P is a transform kernel parameter, ΦPRτ is a frequency rotation operator, Φt0,PSτ is a frequency shift operator, and gστ is a window function of the time-frequency transform.According to formula (4), the basic process of GPTF analysis can be described as follows. Firstly, the signal is rotated in the time-frequency plane, that is, the instantaneous frequency fmτ of the signal is subtracted from κPτ, and then, the signal is translated on the time-frequency plane by the frequency increment κPt0. Finally, the window function is used to make a STFT on the signal. When κPτ=fmτ, the signal forms after the signal srτ is rotated and translated as shown in the following equation:(5)s¯τ=srτ⋅ΦPRτ⋅Φt0,PSτ=ej2π∫fmτdτ⋅e−j2π∫κPτdτ⋅ej2π·τ·κPt0=ej2π⋅τ⋅κPt0.At this time, the instantaneous frequencyκPt0 of the signal s¯τ does not change with time, and the signal s¯τ is a stationary signal. At this time, the STFT is performed on the signal s¯τ, and the time-frequency concentration is significantly improved compared to the nonstationary signal. It can be seen from (4) and (5) that the frequency rotation operator and the translation operator are determined by selecting appropriate transformation kernel parameters, and the nonstationary signal is rotated and translated to stabilize the signal for the purpose of improving the time-frequency resolution by the determined frequency rotation operator and the frequency translation operator. It can be seen that the time-frequency resolution of the GPTF analysis is closely related to the accuracy of the transform kernel parameters. From the above analysis, the key to GPTF analysis is to estimate the parameters of the transform kernel.


## 3.2. Estimation Method for Transforming Kernel Parameters

According to the introduction of the principle of GPTF in Section3.1, the working mechanism is actually the STFT based on the frequency rotation operator and the translation operator. When the selected transform kernel matches the time-frequency characteristic of the analyzed signal, the rotated signal is close to the stationary signal, the time-frequency resolution is high after STFT. The description of the time-frequency characteristics of the stationary signal is accurate. It can be seen that the accuracy of the transform kernel parameter determines the performance of the GPTF method, which affects its performance in analyzing nonstationary signals. In [14], a method for estimating transform kernel parameters using the principle of time-frequency feature approximation is proposed. The time-frequency feature approximation principle indicates that the more accurate the transform kernel parameter estimation, the higher the time-frequency concentration. According to the time-frequency feature approximation principle, if the transform kernel parameter estimation is more accurate, the time-frequency concentration is higher. The principle of time-frequency approximation can be represented in Figure 2.Figure 2
Time-frequency approximation principle.It can be seen from Figure2 that the GPTF analysis method achieves the iterative optimization of the transform kernel parameters by combining the precision of the transform kernel parameters with the time-frequency concentration. The above principle is based on the generalized parametric time-frequency analysis to cyclically approximate and solve the signal time-frequency characteristics, so as to obtain the most suitable transform kernel parameters. According to the principle of time-frequency approximation, the detailed steps of transforming the kernel parameter estimation are as follows:Step1:determining the transform kernel form according to the signal priori form.Step 2: the number of iterationsi=0, initializing the transformation kernel parameters P0.Step4: performing parameterized time-frequency analysisTFt0,f;Pi on the signal according to equation (4).Step 4: detecting the peak ridge linef¯mit in the time-frequency diagram.Step 5: according to the ridge line form, select a suitable fitting method, and fit the ridge line and estimate transformation kernel parameterP˜i.Step 6: the number of iterations is increased by one, and the transformation kernel parameter is updatedPi=P˜i.Step7: calculating a termination conditionΛ=∫f˜mit−f˜mi−1tdt/∫f˜mitdt, where f˜mit is the peak ridge of the ith iteration and f˜mi−1t is the peak ridge of the i−1th iteration, which needs to be noted that f˜m0t=0.Step 8: setting the termination condition parameterζ. The size of this value depends on the situation. If it is too large, the calculation amount is too large, and if it is too small, the accuracy of the extracted time-frequency line is not high.Step 9: comparing the size of Λ andζ, when ζ > Λ, return to step 3. When ζ ≤ Λ, perform step 10.Step 10: outputting a transformation kernel parameter to obtain a time-frequency curve.According to the above steps of transforming the kernel parameter estimation, the flowchart of the design is shown in Figure3.Figure 3
The flowchart of transform kernel parameter estimation.


## 3.3. Ionospheric Decontamination Based on GPTF

Based on the above analysis process and steps, the derivation process of ionospheric decontamination based on GPTF is as follows:Step 1: assume that the signal after the ionospheric contamination is expressed as follows [7–9]:(6)sp=expjm2πsin2πf0texpj2πf1t,wherem is the amplitude coefficient, f0 is the modulation frequency, and f1 is the fixed frequency.Step 2: derive the signal and obtain the Doppler curve form, which can be expressed as follows:(7)ft=m2πf0cos2πf0t+f1.Since the cosine Doppler curve of the signal has periodicity, the ridge line is fitted using the Fourier series for the curve. The Doppler instantaneous frequency form is known. Since the Doppler curve is known to be a cosine form by prior knowledge, the transform kernel of the applied harmonics is as follows:(8)ft=∑m=1∞bmcosm2πfct+b0.There is a certain relationship between the cosine series coefficient of the curve and its Fourier transform coefficient. Consider the Fourier transform to find the coefficients of the cosine series of each order. The discrete form of the Doppler curve parameter estimation method is derived as follows.Radar sampling intervalΔt=1/PRF; equation (8) discretization is expressed as follows:(9)fn⋅Δt=fdn⋅Δt+b0,where fdn⋅Δt=∑m=1∞bmcosm⋅2πfcnΔt.After equation (9) is truncated by adding a rectangular window, the spectrum of each order cosine series becomes a sinc function form. In this case, in order to eliminate the influence of the sinc side lobes of the cosine series spectrum of each order on the spectrum of other order cosine series, the rectangular window length is set to TT=1/fc,2/fc,3/fc. In this way, the peak of the cosine series spectrum of other orders will be located at the zero point of the side lobe of the spectrum of each order cosine series, so as to eliminate the interference between the spectra. Therefore, in the latter experiment, the accumulation time is required to be an integer multiple of the period. When the accumulation time is not an integer multiple, the signal can be truncated, so this method requires at least one periodic signal. Without loss of generality, set TT=i/fc, at this time N=i/fcΔt. Similarly, the Fourier transform of the N point in the harmonic form fdn⋅Δt in (9) is obtained as follows:(10)Fdk=∑n=0N−1fdnΔte−j2π/Nkn=∑n=0N−1∑m=1∞bmcosm2πiNne−j2π/Nkn=∑m=1∞bm∑n=0N−1cos2πiNn⋅micos2πNn⋅k−j∑m=1∞bm∑n=0N−1cos2πNn⋅misin2πNn⋅k.According to the characteristics of the orthogonal trigonometric function set, formula (10) can be reduced to(11)Fdk=bmN2,k=mi,0,else.The transformation kernel parameter can be obtained as follows:(12)bm=2NrealFdkk=mi,b0=λ22realFd1NΔt.


## 4. Bistatic Sea Clutter Suppression Method Based on FBLP Algorithm

If the dominant factor of sea clutter broadening is the bistatic angle, or the sea clutter disturbed by the ionospheric phase path has been dedisturbed, then the suppression of bistatic sea clutter is required. At present, there are few studies on the problem of HF radar bistatic sea clutter suppression. The bistatic sea clutter shows different characteristics from monostatic sea clutter, and it makes the traditional Bragg line theoretical prediction or monostatic sea clutter characteristic suppression methods fail. At present, there are two main ideas for the bistatic sea clutter suppression method [15–17]: one is the sea clutter method based on singular value decomposition (SVD) [15], and the other is based on the sinusoidal signal parameter estimation method, such as the root cancellation method and the CLEAN method. The root cancellation method and the CLEAN algorithm are simple to implement, but the number of iterations is difficult to determine accurately. The SVD algorithm can suppress sea clutter better, but the amount of computation is large. When the target is located near the first-order Bragg peak, the three algorithms may lose part of the target information.Based on this, this paper uses the FBLP algorithm to suppress the sea clutter. The FBLP algorithm of this paper was previously a conference paper by authors, which was published in [12]. This algorithm uses the first-order sea clutter Bragg band distribution feature as a priori knowledge and combines the multidimensional characteristics of the measured sea clutter data. The method can effectively and accurately identify the first-order sea clutter and can better suppress the sea clutter while retaining the target. And the calculation amount is small. Firstly, the bistatic first-order sea clutter signal model is introduced. Secondly, the theoretical Bragg band range is determined in advance on the first-order sea clutter distance-Doppler spectrum. Then the FBLP algorithm is used to establish the forward and backward prediction matrix to estimate the sea clutter signal parameters. Then, based on the multidimensional characteristics of the measured sea clutter, the first-order sea clutter signal is better identified and extracted by the signal-to-noise ratio threshold method and the symmetry detection matrix. Finally, the amplitude of each sea clutter frequency component is zeroed to achieve the purpose of suppressing sea clutter.


### 4.1. Multidimensional Feature Analysis of Sea Clutter Spectrum

When the CLEAN or root cancellation algorithms are used to suppress sea clutter for traditional SVD, it is necessary to detect sea clutter firstly in order to suppress sea clutter. However, due to the influence of wind direction or ocean current in the actual marine environment, the characteristics of the sea clutter spectrum are complex, which makes the identification and extraction of the first-order sea clutter difficult, resulting in insufficient sea clutter suppression and the formation of false targets and tracks.Figure4(a) shows the range-Doppler spectrum of the first-order sea clutter measured by high-frequency hybrid sky-surface wave radar. From the measured graphs, we can also see that the first-order sea clutter exhibits multidimensional characteristics. It can be seen from Figure 4(a) that the first-order sea clutter spectrum of the close-range gate shows good continuity and symmetry in both the range and the Doppler dimension. Figure 4(a) shows the sea clutter spectrum of continuous distance gates. Due to the influence of wind direction, the amplitude of the first-order sea clutter positive and negative Bragg peaks shows asymmetry in Figure 4(b).Figure 4
The range-Doppler spectrum and asymmetry of the measured sea clutter in HF sky-surface wave radar. (a) Range-Doppler spectrum of actual sea clutter. (b) The asymmetry of actual sea clutter spectrum.
(a)
(b)


### 4.2. Parameter Estimation of Sea Clutter Signal Based on Linear Prediction Model

For a bistatic HF radar with the hybrid operating mode of sky-wave transmitting and ground-wave receiving, it is assumed that there areK complex sinusoidal components in the Bragg band in the scattering unit. Therefore, the bistatic first-order sea clutter model after ionospheric phase path disturbance suppression can also be written as follows [6, 12]:(13)yn=∑i=1Kaiej2πfBin−1Ts+jφi,n=1,2,…,N,where ai, fBi, and φi are, respectively, the amplitude of each signal component, first-order sea clutter Bragg frequency, and phase. Ts is the signal sampling period, and N is the signal sequence length. Here, fBmin≤fBi≤fBmax, fBmin,fBmax is the boundary of the frequency band caused by the bistatic angle in the scattering unit.The parameters of each sea clutter signal are estimated below based on a linear prediction method. When the signal-to-noise ratio (SNR) is sufficiently high, the sea clutter signal can be written as the sum of the firstL weights:(14)yn=−∑i=1Lyn−i⋅αi,where αi stands for the weight coefficients.Besides, so as to track the time-varying behavior of the sea clutter, the coefficients of the prediction error filter must be estimated over short data segments so that the filter coefficients could be updated adaptively. And the prediction equation matrix, which is defined in [18], could be expressed as follows:(15)yLyL−1…y1yL+1yL…y2⋮⋮⋱⋮yN−1yN−2…yN−1y∗2y∗3…y∗L+1y∗3y∗4…y∗L+2⋮⋮⋱⋮y∗N−Ly∗N−L+1…y∗Nα1α2⋮αL=−yL+1yL+2⋮yNy∗1y∗2⋮y∗N−1.Simply, we denote as(16)A⋅α=−y.And the weight coefficient matrixα could be estimated by the following equation:(17)α=−A−1y.As the linear equations above are overdetermined, the total least square method is used to solve this problem [18], namely:(18)α=−A−1y=−AHA−1AHy,where H indicates the conjugate transpose. Simultaneously, the order of prediction error filter polynomial L could be determined in [18], which satisfies the inequality as follows:(19)M≤L≤N−M2.That is to say, the order of prediction error filter should exceed the estimated signal number. Afterwards, we define frequency estimation matrix as follows:(20)f=f1,f2,…,fM,where(21)fk=1,e−sk,e−2sk,…,e−Lsk,k=1,2,…,M.It is easy to observe that each row inA is a linear combination of L linearly independent vectors in fk. That is to say, the rank of A is M as long as A has at least M rows. Thus, the dimension of null space in A is L+1−M dimension. In addition, as αk lies in the null space of A, we have(22)α0+α1e−5k+α2e−2sk+⋯+αLe−Lsk=0.The signal frequency could be estimated from the roots of (22). Besides, in order to obtain the amplitude and initial phase of each signal, we define the following matrix equation as follows:(23)β10β20…βM0β11β22…βM1⋮⋮⋮⋮β1N−1β2N−1…βMN−1h1h2h3h4=y1y2⋮yN.Simply, we denote as(24)R⋅h=Y,where βik=ej2πfik/fs, h=RHR−1RHY, in which fs represents the sampling frequency and N is the sampling number. In addition, the amplitude and initial phase could be obtained after taking the manipulation of the absolute and angular value of h, respectively:(25)ak=hk,k=1,…,M,φk=arctanImhkRehkk=1,…,M,where Im⋅ and Re⋅ represent taking the imaginary and real parts of the signal, respectively.


### 4.3. First-Order Sea Clutter Suppression

Through the previous analysis, in order to better achieve the purpose of suppressing sea clutter based on the FBLP algorithm, we use the first-order sea clutter obtained in the theoretical Bragg band of the range-Doppler spectrum, combined with the multidimensional features of the measured first-order sea clutter to better identify and suppress sea clutter. First, the range of the sea clutter Bragg band is determined in advance on the range-Doppler spectrum, and then, the FBLP algorithm is used to extract and estimate the parameters of the sea clutter signal. Based on the previous processing, we use the SNR threshold method defined below to determine the number of first-order sea clutter signals of interest. Specially, we supposeS as the number of signals, and the total power of sea clutter echo could be obtained by the following equation:(26)P=∑n=1Nyn2,where N is the number of signal sampling points. Then the power of each signal and noise is given by(27)Pk=1N∑n=1Nakβkn−12=ak21−βk2NN1−βk2k=1,2,…,S,Pnoise=1N∑n=1Nyn−∑k=1Sakβkn−12.The SNR of each signal is given by(28)SNRk=PkPnoise,k=1,2,…,S.According to the local dominant characteristics of sea clutter, the SNR of each signal component can be compared with the threshold to screen out the suspected sea clutter signal. Then, according to the frequency domain symmetry property of the sea clutter, combined with the frequency estimation matrixf, each estimated frequency is preprocessed as follows:(29)fi′=fi+fB,fi<0,fi−fB,fi>0,where fB is the theoretical Bragg scattering frequency, λ indicates the wave length, and g is the acceleration of gravity. Afterwards, we define the feature detection matrix (FDM) as follows [12]:(30)F=f1,2f1,3f1,4…f1,M0f2,3f2,4…f2,M00f3,4…f3,M000…fM−1,M,where fi,j=fi′−fj′,i=1,…,S−1;j=i+1,…,S.The processing diagram of bistatic sea clutter suppression based on FBLP algorithm is shown in Figure5. By performing a two-dimensional search on the triangular matrix F and finding the minimum value, the frequency component fx,fy satisfying the symmetry is obtained, that is, the component is the Bragg frequency corresponding to the first-order sea clutter. Finally, the corresponding amplitude in the amplitude matrix A corresponding to the component is zeroed Ax=0andAy=0. This completes the cancellation of the first-order sea clutter. In addition, considering the asymmetry of the first-order sea clutter positive and negative peaks caused by the wind direction, since the sea clutter maintains good continuity in the range dimension, it can be adopted as a reference by using adjacent distance gates. The above FDM detection step is then corrected.Figure 5
The processing diagram of bistatic sea clutter suppression based on the FBLP algorithm.


## 4.1. Multidimensional Feature Analysis of Sea Clutter Spectrum

When the CLEAN or root cancellation algorithms are used to suppress sea clutter for traditional SVD, it is necessary to detect sea clutter firstly in order to suppress sea clutter. However, due to the influence of wind direction or ocean current in the actual marine environment, the characteristics of the sea clutter spectrum are complex, which makes the identification and extraction of the first-order sea clutter difficult, resulting in insufficient sea clutter suppression and the formation of false targets and tracks.Figure4(a) shows the range-Doppler spectrum of the first-order sea clutter measured by high-frequency hybrid sky-surface wave radar. From the measured graphs, we can also see that the first-order sea clutter exhibits multidimensional characteristics. It can be seen from Figure 4(a) that the first-order sea clutter spectrum of the close-range gate shows good continuity and symmetry in both the range and the Doppler dimension. Figure 4(a) shows the sea clutter spectrum of continuous distance gates. Due to the influence of wind direction, the amplitude of the first-order sea clutter positive and negative Bragg peaks shows asymmetry in Figure 4(b).Figure 4
The range-Doppler spectrum and asymmetry of the measured sea clutter in HF sky-surface wave radar. (a) Range-Doppler spectrum of actual sea clutter. (b) The asymmetry of actual sea clutter spectrum.
(a)
(b)


## 4.2. Parameter Estimation of Sea Clutter Signal Based on Linear Prediction Model

For a bistatic HF radar with the hybrid operating mode of sky-wave transmitting and ground-wave receiving, it is assumed that there areK complex sinusoidal components in the Bragg band in the scattering unit. Therefore, the bistatic first-order sea clutter model after ionospheric phase path disturbance suppression can also be written as follows [6, 12]:(13)yn=∑i=1Kaiej2πfBin−1Ts+jφi,n=1,2,…,N,where ai, fBi, and φi are, respectively, the amplitude of each signal component, first-order sea clutter Bragg frequency, and phase. Ts is the signal sampling period, and N is the signal sequence length. Here, fBmin≤fBi≤fBmax, fBmin,fBmax is the boundary of the frequency band caused by the bistatic angle in the scattering unit.The parameters of each sea clutter signal are estimated below based on a linear prediction method. When the signal-to-noise ratio (SNR) is sufficiently high, the sea clutter signal can be written as the sum of the firstL weights:(14)yn=−∑i=1Lyn−i⋅αi,where αi stands for the weight coefficients.Besides, so as to track the time-varying behavior of the sea clutter, the coefficients of the prediction error filter must be estimated over short data segments so that the filter coefficients could be updated adaptively. And the prediction equation matrix, which is defined in [18], could be expressed as follows:(15)yLyL−1…y1yL+1yL…y2⋮⋮⋱⋮yN−1yN−2…yN−1y∗2y∗3…y∗L+1y∗3y∗4…y∗L+2⋮⋮⋱⋮y∗N−Ly∗N−L+1…y∗Nα1α2⋮αL=−yL+1yL+2⋮yNy∗1y∗2⋮y∗N−1.Simply, we denote as(16)A⋅α=−y.And the weight coefficient matrixα could be estimated by the following equation:(17)α=−A−1y.As the linear equations above are overdetermined, the total least square method is used to solve this problem [18], namely:(18)α=−A−1y=−AHA−1AHy,where H indicates the conjugate transpose. Simultaneously, the order of prediction error filter polynomial L could be determined in [18], which satisfies the inequality as follows:(19)M≤L≤N−M2.That is to say, the order of prediction error filter should exceed the estimated signal number. Afterwards, we define frequency estimation matrix as follows:(20)f=f1,f2,…,fM,where(21)fk=1,e−sk,e−2sk,…,e−Lsk,k=1,2,…,M.It is easy to observe that each row inA is a linear combination of L linearly independent vectors in fk. That is to say, the rank of A is M as long as A has at least M rows. Thus, the dimension of null space in A is L+1−M dimension. In addition, as αk lies in the null space of A, we have(22)α0+α1e−5k+α2e−2sk+⋯+αLe−Lsk=0.The signal frequency could be estimated from the roots of (22). Besides, in order to obtain the amplitude and initial phase of each signal, we define the following matrix equation as follows:(23)β10β20…βM0β11β22…βM1⋮⋮⋮⋮β1N−1β2N−1…βMN−1h1h2h3h4=y1y2⋮yN.Simply, we denote as(24)R⋅h=Y,where βik=ej2πfik/fs, h=RHR−1RHY, in which fs represents the sampling frequency and N is the sampling number. In addition, the amplitude and initial phase could be obtained after taking the manipulation of the absolute and angular value of h, respectively:(25)ak=hk,k=1,…,M,φk=arctanImhkRehkk=1,…,M,where Im⋅ and Re⋅ represent taking the imaginary and real parts of the signal, respectively.


## 4.3. First-Order Sea Clutter Suppression

Through the previous analysis, in order to better achieve the purpose of suppressing sea clutter based on the FBLP algorithm, we use the first-order sea clutter obtained in the theoretical Bragg band of the range-Doppler spectrum, combined with the multidimensional features of the measured first-order sea clutter to better identify and suppress sea clutter. First, the range of the sea clutter Bragg band is determined in advance on the range-Doppler spectrum, and then, the FBLP algorithm is used to extract and estimate the parameters of the sea clutter signal. Based on the previous processing, we use the SNR threshold method defined below to determine the number of first-order sea clutter signals of interest. Specially, we supposeS as the number of signals, and the total power of sea clutter echo could be obtained by the following equation:(26)P=∑n=1Nyn2,where N is the number of signal sampling points. Then the power of each signal and noise is given by(27)Pk=1N∑n=1Nakβkn−12=ak21−βk2NN1−βk2k=1,2,…,S,Pnoise=1N∑n=1Nyn−∑k=1Sakβkn−12.The SNR of each signal is given by(28)SNRk=PkPnoise,k=1,2,…,S.According to the local dominant characteristics of sea clutter, the SNR of each signal component can be compared with the threshold to screen out the suspected sea clutter signal. Then, according to the frequency domain symmetry property of the sea clutter, combined with the frequency estimation matrixf, each estimated frequency is preprocessed as follows:(29)fi′=fi+fB,fi<0,fi−fB,fi>0,where fB is the theoretical Bragg scattering frequency, λ indicates the wave length, and g is the acceleration of gravity. Afterwards, we define the feature detection matrix (FDM) as follows [12]:(30)F=f1,2f1,3f1,4…f1,M0f2,3f2,4…f2,M00f3,4…f3,M000…fM−1,M,where fi,j=fi′−fj′,i=1,…,S−1;j=i+1,…,S.The processing diagram of bistatic sea clutter suppression based on FBLP algorithm is shown in Figure5. By performing a two-dimensional search on the triangular matrix F and finding the minimum value, the frequency component fx,fy satisfying the symmetry is obtained, that is, the component is the Bragg frequency corresponding to the first-order sea clutter. Finally, the corresponding amplitude in the amplitude matrix A corresponding to the component is zeroed Ax=0andAy=0. This completes the cancellation of the first-order sea clutter. In addition, considering the asymmetry of the first-order sea clutter positive and negative peaks caused by the wind direction, since the sea clutter maintains good continuity in the range dimension, it can be adopted as a reference by using adjacent distance gates. The above FDM detection step is then corrected.Figure 5
The processing diagram of bistatic sea clutter suppression based on the FBLP algorithm.


## 5. Cascaded Decontamination and Suppression Algorithm Performance Analysis

### 5.1. Simulation Data Analysis

Simulation parameter setting is as follows: the signal type is a linear frequency modulated (LFM) signal, pulse repetition period is PRT = 40 ms, pulse accumulation time isT = 40 s, operating frequency is f0 = 13 MHz, sea clutter Bragg frequency is fB = 0.3750 Hz, signal-to-noise ratio (SNR) is SNR = 10 dB, and ionosphere contamination function is spt=0.5∗2∗pi∗sin2∗pi∗0.1∗t. Figure 6(a) shows the contaminated sea clutter and decontaminated sea clutter spectrum by S2-method, STFT, and GPTF, Figure 6(b) shows the estimated phase contamination by S2-method, STFT, and GPTF, Figure 6(c) shows decontaminated sea clutter spectrum by GPTF, Figure 6(d) shows the broadening sea clutter suppression based on the FBLP algorithm.Figure 6
Sea clutter suppression with ionospheric phase disturbance based on GPTF and FBLP. (a) Contaminated sea clutter and decontaminated sea clutter spectrum by S2-method, STFT, and GPTF. (b) Estimated phase contamination by PGA, STFT, and GPTF. (c) Decontaminated sea clutter spectrum by GPTF. (d) The broadening sea clutter suppression based on the FBLP algorithm.
(a)
(b)
(c)
(d)Figure6 shows the simulation results of sea clutter suppression with ionospheric phase disturbance based on GPTF and FBLP and the target Doppler frequency ftarget = 0.75 Hz. It can be seen from Figure 6 that the sea clutter can still be effectively suppressed by GPTF and FBLP and the target can be clearly seen. As can be seen from the above results, the approach based on GPTF can be used for correcting ionospheric contamination with large amplitude. And compared to the traditional time-frequency analysis methods, this method has high time-frequency resolution and no cross-term interference.


### 5.2. Experimental Data Analysis

The performance of ionosphere decontamination based on GPTF is verified by the measured data of HF sky-wave radar. The measured data have been processed by the traditional range-Doppler-azimuth processing. The echo accumulation time is 40 ms. The results of ionospheric decontamination based on GPTF at range gate = 60 and azimuth angle DBF = 3° and at range gate = 100 and azimuth angle DBF = 5° are shown in Figures7(a)-7(b). The green curve in the figure represents the contaminated echo spectrum. The red curve in the figure represents the decontaminated echo spectrum after the GPTF. The black curve in the figure represents the decontaminated echo spectrum after the S2-method. It can be seen from the above processing results that the GPTF method has better ionospheric decontamination performance than the S2-method.Figure 7
Decontamination results based on the GPTF method using the HF sky-wave radar. Decontamination result by the GPTF method (a) at range gate = 60 and azimuth angle DBF = 3° and (b) at range gate = 100 and azimuth angle DBF = 5°.
(a)
(b)In order to justify the effectiveness of the cascade suppression processing method based on GPTF and FBLP proposed in this paper, the following uses the measured data of HF sky-surface wave radar to verify. The measured data were acquired using the newly developed “high-frequency sky-surface wave radar experimental platform.” The measured data have been processed by the traditional range-Doppler-azimuth processing. The echo accumulation time is 50 s. The suppressed sea clutter spectrum by GPTF and FBLP algorithm at range gate = 35 and azimuth angle DBF = 10° is shown in Figures8(a)-8(b) among them, the added target Doppler velocity ftarget = −0.6 Hz.Figure 8
The broadening sea clutter suppression based on GPTF and FBLP algorithm. (a) Decontamination result by GPTF method. (b) Sea clutter suppression by FBLP.
(a)
(b)As shown in Figures8(a)-8(b), the contaminated sea clutter can be effectively suppressed by GPTF and the broadening sea clutter is obviously sharpened, as indicated by the red line. On this basis, the FBLP algorithm can effectively suppress the sea clutter, as indicated by the green line. From Figures 8(a)-8(b), we can see that the target can be effectively detected and SCR increased. This method can steadily and effectively sharpen severe sea clutter spectrum contaminated by ionospheric phase perturbation with large amplitude. Finally, it should be pointed out that we use different measured data to make many tests for the cascaded suppression method. And we found that, in most cases, this method can lead to good decontamination and suppression effect, which sufficiently shows its robustness.


## 5.1. Simulation Data Analysis

Simulation parameter setting is as follows: the signal type is a linear frequency modulated (LFM) signal, pulse repetition period is PRT = 40 ms, pulse accumulation time isT = 40 s, operating frequency is f0 = 13 MHz, sea clutter Bragg frequency is fB = 0.3750 Hz, signal-to-noise ratio (SNR) is SNR = 10 dB, and ionosphere contamination function is spt=0.5∗2∗pi∗sin2∗pi∗0.1∗t. Figure 6(a) shows the contaminated sea clutter and decontaminated sea clutter spectrum by S2-method, STFT, and GPTF, Figure 6(b) shows the estimated phase contamination by S2-method, STFT, and GPTF, Figure 6(c) shows decontaminated sea clutter spectrum by GPTF, Figure 6(d) shows the broadening sea clutter suppression based on the FBLP algorithm.Figure 6
Sea clutter suppression with ionospheric phase disturbance based on GPTF and FBLP. (a) Contaminated sea clutter and decontaminated sea clutter spectrum by S2-method, STFT, and GPTF. (b) Estimated phase contamination by PGA, STFT, and GPTF. (c) Decontaminated sea clutter spectrum by GPTF. (d) The broadening sea clutter suppression based on the FBLP algorithm.
(a)
(b)
(c)
(d)Figure6 shows the simulation results of sea clutter suppression with ionospheric phase disturbance based on GPTF and FBLP and the target Doppler frequency ftarget = 0.75 Hz. It can be seen from Figure 6 that the sea clutter can still be effectively suppressed by GPTF and FBLP and the target can be clearly seen. As can be seen from the above results, the approach based on GPTF can be used for correcting ionospheric contamination with large amplitude. And compared to the traditional time-frequency analysis methods, this method has high time-frequency resolution and no cross-term interference.


## 5.2. Experimental Data Analysis

The performance of ionosphere decontamination based on GPTF is verified by the measured data of HF sky-wave radar. The measured data have been processed by the traditional range-Doppler-azimuth processing. The echo accumulation time is 40 ms. The results of ionospheric decontamination based on GPTF at range gate = 60 and azimuth angle DBF = 3° and at range gate = 100 and azimuth angle DBF = 5° are shown in Figures7(a)-7(b). The green curve in the figure represents the contaminated echo spectrum. The red curve in the figure represents the decontaminated echo spectrum after the GPTF. The black curve in the figure represents the decontaminated echo spectrum after the S2-method. It can be seen from the above processing results that the GPTF method has better ionospheric decontamination performance than the S2-method.Figure 7
Decontamination results based on the GPTF method using the HF sky-wave radar. Decontamination result by the GPTF method (a) at range gate = 60 and azimuth angle DBF = 3° and (b) at range gate = 100 and azimuth angle DBF = 5°.
(a)
(b)In order to justify the effectiveness of the cascade suppression processing method based on GPTF and FBLP proposed in this paper, the following uses the measured data of HF sky-surface wave radar to verify. The measured data were acquired using the newly developed “high-frequency sky-surface wave radar experimental platform.” The measured data have been processed by the traditional range-Doppler-azimuth processing. The echo accumulation time is 50 s. The suppressed sea clutter spectrum by GPTF and FBLP algorithm at range gate = 35 and azimuth angle DBF = 10° is shown in Figures8(a)-8(b) among them, the added target Doppler velocity ftarget = −0.6 Hz.Figure 8
The broadening sea clutter suppression based on GPTF and FBLP algorithm. (a) Decontamination result by GPTF method. (b) Sea clutter suppression by FBLP.
(a)
(b)As shown in Figures8(a)-8(b), the contaminated sea clutter can be effectively suppressed by GPTF and the broadening sea clutter is obviously sharpened, as indicated by the red line. On this basis, the FBLP algorithm can effectively suppress the sea clutter, as indicated by the green line. From Figures 8(a)-8(b), we can see that the target can be effectively detected and SCR increased. This method can steadily and effectively sharpen severe sea clutter spectrum contaminated by ionospheric phase perturbation with large amplitude. Finally, it should be pointed out that we use different measured data to make many tests for the cascaded suppression method. And we found that, in most cases, this method can lead to good decontamination and suppression effect, which sufficiently shows its robustness.


## 6. Conclusions

Ship detection is an important mission of the HF hybrid sky-surface wave radar. The sea clutter spectrum is influenced by composite factors such as the bistatic angle and ionospheric phase disturbance. Thus, an ionospheric decontamination and sea clutter suppression method for HF hybrid sky-surface wave radars based on GPTF analysis is proposed. In this method, the time-frequency analysis method based on GPTF is proposed to correct the ionospheric phase contamination with large amplitude. Compared with the traditional nonparametric time-frequency analysis method, the phase contamination function extracted by this GPTF method has higher accuracy. At this time, the FBLP algorithm is used to suppress the sharped sea clutter spectrum. Finally, the proposed method is examined by measured data. Experimental results indicate that the proposed method can well suppress the broadening sea clutter for HF hybrid sky-surface wave radar. Although some positive results have been obtained, it must be pointed that there is still a need for more studies and improvements.



---

*Source: 7627521-2020-08-19.xml*