# Time-Varying Ocean-Like Surface Scattering at Grazing Incidence: Numerical Analysis of Doppler Spectrum at HF/VHF/UHF Bands

Authors: Yidong Hou; Biyang Wen; Caijun Wang; Yonghuai Yang
Journal: International Journal of Antennas and Propagation
 (2019)

---

## Abstract

This paper numerically analyzes the characteristics of the Doppler spectrum at HF/VHF/UHF bands from 1D time-varying ocean-like surfaces at grazing incidence in vertical polarization mode. The rough surface is transformed into a local perturbation plane which has its roughness flattened at the edges. The scattering waves include coherent reflected wave and incoherent scattering waves. The surface currents exciting the incoherent scattering waves are regarded as the unknowns which can be solved from the improved surface integral equation using the method of moments (MoM). The incident plane wave allows the incident angle to reach up to 90° (grazing incidence). Then the backscattering wave in the far field can be calculated, and the Doppler spectrum is obtained by coherent Monte-Carlo simulation. Firstly, the validity of the method is verified by comparing with the mature small perturbation method at the HF band. Then the incident wave frequency is asymptotically increased from HF to UHF, and the application range of the SPM is quantitatively evaluated in the Doppler spectrum domain. Finally, the paper focuses on analyzing the characteristics of Doppler spectrum in different bands and different sea states and comparing the influence of nonlinear ocean waves on the Doppler spectrum at different frequencies.


---

## Body

## 1. Introduction

In the issue of remote sensing of sea state using radar, the echo Doppler spectrum can provide more details than the radar cross-section (RCS). The Doppler spectrum from time-varying ocean-like rough surface has been studied both theoretically and numerically for the past several decades [1–4]. Benefiting from the small perturbation method (SPM) and Kirchhoff approximation (KA) theory, it has been possible to extract the sea state information from the echo Doppler spectrum at HF and X band [5–8]. However, these analytical methods will be invalid at VHF and UHF band [9, 10]. In fact, the sea echoes in the VHF and UHF bands also contain abundant sea state information. As early as 1965, [11] had used the UHF coherent Doppler radar to observe the echo Doppler spectrum of the Pacific region. In recent years, our team has used shore-based all-digital UHF radar to measure near-shore ocean echo Doppler spectrum [12], hoping to extract near-shore sea state information from the sea clutter of UHF radar. However, due to lack of theoretical basis, this work has been stagnant. Because accuracy is independent of frequency, the numerical method is a reasonable way to analyze the relationship between the Doppler spectrum and sea state at VHF and UHF band.When using numerical methods to solve rough surface scattering problems, it is inevitable to deal with the influence of the finite rough surface edges on the scattering wave. The classical rigorous approach is to constrain incident wave with tapered beam to make it has a Gaussian footprint on the rough surface [9]. As the grazing angle becomes smaller, the tapered wave width g needs to be increased exponentially to satisfy the Helmholtz equation [13]. Especially at the grazing incidence, g can reach thousands of times of radio wavelength [14], which greatly increases the calculating burden. Therefore, the problem of low grazing angle (LGA) has become a stumbling block for numerically analyzing the electromagnetic scattering and Doppler spectrum from rough ocean-like surface. In order to deal with this kind of problem, many scholars have made a lot of attempts and improvements. First, a number of acceleration algorithms, including the banded matrix iterative approach (BMIA) [15], the steepest descent fast multipole method (SDFMM) [16], and the method of ordered multiple interaction (MOMI) [17], were proposed to accelerate Monta-Carlo simulations of large-scale random rough surface scattering. Paper [14] used the MOMI to simulate the Doppler spectrum from time-varying ocean-like surface at a grazing angle of 5∘ and L band. Paper [18] studied 14-GHz backscattering from ocean-like surfaces at a grazing angle of 1∘ with the BMIA. Another kind of approach is to process the rough surface. Paper [19] assumed that the rough surface is periodic and then used periodic boundary condition and incident plane wave to analyze rough surface scattering at LGA. However, the actual rough sea surface is aperiodic, and the effect of this assumption on the Doppler spectrum has not been evaluated. Paper [20] added the resistive loading on rough surface edges to suppress the edge current, which effectively eliminated edge effects when the LGA was greater than 20∘. Paper [21] proposed the grazing method of moments (GMoM), a novel and rigorous approach to deal with the scattering from rough surfaces under grazing illumination. GMoM represents the rough surfaces as a bounded perturbation of a plane which enlightened by a plane wave and solves the tangential components of the scattered field with improved boundary integral equation method (IEM). The effectiveness of the GMoM at LGA has been verified by comparison with both classical numerical and analytical methods. More implementation details and performance evaluations can also refer to these works [22–24].In this paper, considering that most shore-based ocean state remote sensing radars work in the vertical polarization mode (TM case), we focus on TM case and analyze backscattering Doppler spectral characteristics from time-varying ocean-like surfaces at HF/VHF/UHF bands and grazing incidence. For simplicity and to reduce the computational burden, we restrict our study to 1D surface. Following the GMoM, ocean-like rough surface is considered as a local perturbation on an infinite plane and the roughness of the perturbation region close to zero at the edges. The surface currents producing the incoherent scattered waves, which we call scattering currents, are seen as the surface unknowns. Then the improved surface integral equation with impedance boundary condition is established to solve surface unknowns using the moment method. Using incident plane wave as excitation allows grazing incidence. The surface integral equation acts on unbounded local perturbation plane, so there is no edge effect. The scattering currents on the perturbation region also close to zero at the edges, which can be seen as windowing. The windowed scattering currents can effectively suppress the side lobes of the scattered waves. Finally, the Doppler spectrum is obtained through coherent Monte-Carlo simulations. This calculation process is very similar to the classical moment method [25]. The differences are treating the rough surface as a locally perturbed plane and replacing the total surface currents with the surface scattering currents as surface unknowns. To verify the effectiveness of these processes, the simulated Doppler spectra under different levels of sea state are compared with that calculated by the SPM at HF band. Then the simulation frequency is increased gradually from HF to UHF to quantitatively analyze the application scope of the SPM in the Doppler domain. This paper focuses on numerical analysis of the Doppler spectral characteristics of the linear and nonlinear sea surface (including Creamer model [26] and Lagrange model [27]) at different wavebands and the variation of the Doppler spectrum with sea state levels.The rest of the paper is organized as follows. Section2 describes the generation of linear and nonlinear time-varying ocean-like rough surface and the processing of rough surface edges. Section 3 gives the basic electromagnetic scattering formulas and the Monte-Carlo simulation of backscattering Doppler spectrum. In Section 4, the validity of the numerical Doppler spectrum at grazing incidence is verified by comparing with the SPM at HF band, and then the Doppler spectrum characteristics of linear and nonlinear ocean surface at HF/VHF/UHF bands are discussed in detail. Section 5 summarizes the conclusions.


## 2. Linear and Nonlinear Gravity Waves

From HF to UHF band, sea waves interacting with electromagnetic waves mainly belong to gravity waves [28]. Therefore, the contribution of the capillary waves to the scattered field is not considered in this paper. In this section, we will generate one-dimensional linear and nonlinear gravity waves for the follow-up study of the paper. The nonlinear waves contain the Creamer waves [26] and the Lagrange waves [27].


### 2.1. Linear Waves

The wave height of the linear sea surface satisfies the Gaussian distribution, which means that the linear waves can be regarded as a sum of a series of sine waves, and the phase and amplitude of each component obey the uniform distribution and Rayleigh distribution, respectively. The intensity of the sine wave with different wavenumbers is determined by the sea spectrum. Therefore, the linear gravity waves can be expressed as(1)hlx,t=Re∑nAkn,texp⁡jknx=Re∑nαnSkndkexp⁡-jwntexp⁡jknxwhere kn=ndk is the wavenumber, dk=2π/L depends on the rough surface length L, αn is the independent complex random Gauss variable with zero mean and unit variance, wn=gkn is the angular frequency of sea wave with a wavenumber of kn determined by the dispersion relation under the deep water condition, g is the gravitational acceleration, and S(kn) denotes the sea spectrum. The rough surface is sampled as N discrete points, and n is an integer in the interval of [-N/2,N/2-1]. Then the discrete sea surface profile can be calculated using the fast Fourier transform (FFT).In this paper, JONSWAP spectrum is employed to drive the sea waves, which is expressed by(2)Skw=0.038kw3U102gxw0.22exp⁡-1.25wpgkw4γrwhere kw is the water wavenumber, U10 represents the wind speed at 10 meters above the sea level, xw is a constant representing the range of wind area, γ is usually 3.3, and(3a)wp=22g2U10xw1/3(3b)r=exp⁡-gkw-wp22σ2wp2(3c)σ=0.07,gkw<wp0.09,gkw≥wpAs an inadequately developed spectrum, JONSWAP spectrum is more realistic than the PM spectrum. When xw increases, JONSWAP spectrum and PM spectrum will be closer.


### 2.2. Creamer Waves

The actual sea waves are nonlinear, and there are weak interactions between waves with different frequencies. These weak interactions can be considered as higher-order perturbation solutions to the hydraulic motion equation, which will significantly affect the characteristics of the Doppler spectrum, and the influence degree will also change with the sea state. Therefore, the nonlinearity of the sea waves cannot be neglected when studying the Doppler spectrum of the sea surface no matter which band the radar works in. The second-order nonlinear waves, generated by two wave interactions, have been described by [29, 30], which have been applied in the second-order RCS of ocean surface at HF band [31]. Although the perturbation expansion method is clear, it is not suitable for numerical implementation due to the high calculating burden especially for higher-order nonlinear waves. Another alternative method to generate nonlinear sea surfaces is the Hamiltonian formalism based on the weak wave-turbulence theory developed by [26], which has been widely used to study the electromagnetic scattering of 1D or 2D time-varying sea surfaces [3, 4].To begin with, the Creamer method expresses the nonlinear term of the sea surface as a Hilbert transform of the linear sea surface. In the 1D case, the Creamer nonlinear term is denoted as(4)htx,t=Re∑n-jknknAkn,texp⁡jknxwhere the complex amplitude A(kn,t) is as defined in (1). This transform can also be implemented using FFT to reduce the computation cost to Nlog⁡N. In the frequency domain, the complex amplitude of Creamer nonlinear waves can be expressed as(5)Ctkn,t=1N∑xexp⁡jknhtx,t-1knexp⁡-jknxSince (5) cannot be calculated using FFT directly, it needs approximate processing to improve computational efficiency. Expanding exp⁡jknhtx,t to Taylor series, themth expansion of (5) can be obtained as(6)Ctmkn,t=1N∑xjknhtx,tmm!knexp⁡-jknx=jknmm!NknFhtx,tmwhere F· denotes Fourier transform. Thenmth-order Creamer waves can be calculated using inverse Fourier transform, which is (7)hcmx,t=Re1NFi∑nCtmkn,tIt can be proved that the first-order Creamer waves just are linear waves.


### 2.3. Lagrange Waves

Both perturbation method and the Creamer method describe the vertical skewness of sea waves in the Euler coordinate system. There is horizontal skewness in the actual nonlinear waves, which directly affects the slope distribution of the sea surface and induces more remarkable influence on backscattering signals than vertical skewness [32, 33]. References [32, 33] have attempted to modify the perturbation model to add horizontal skewness on the nonlinear waves and discussed its effect on microwave scattering and emission. Recently [27, 34] proposed a more concise method derived in the Lagrange coordinate system to generate nonlinear waves with horizontal skewness also called Lagrange waves. A disturbance is added in the horizontal direction so that an additional horizontal offset is produced in the linear waves, which is(8)x,hlx,t→x+xLx,t,hlx,twhere hl(x,t) and xL(x,t) describe the vertical and horizontal movements of individual particles with time t, which are two correlated Gauss processes. The vertical process of hl(x,t) is given by (1), and horizontal process of xL(x,t) is completely determined by the vertical process, given by(9)xLx,t=Re∑nHknAkn,texp⁡jknxwhere H(kn) is the response function which determines the nonlinear characters of the Lagrange waves, in particular crest-trough and front-back asymmetry. In Lindgrens derivation [27], H(kn) consists of two parts:(10)Hkn=jtanh⁡hkn+sgnknβωn2where h represents the water depth, sgn· is the sign function, and β is obtained from the relation between the horizontal acceleration of water particles and vertical displacement [27]. Following the study of [27, 34], we set β=0.4 in this paper. The first part of H(kn) determines the vertical skewness and results in waves with more peaked crests and shallower troughs (crest-trough asymmetry) compared to the linear waves, and the second part of H(kn) determines the horizontal skewness and gives waves front-back asymmetry. When the water depth is assumed to be infinite, the response function can be simplified to Hkn=j+sgnknβ/ωn2.


### 2.4. Wave Height and Slope Distribution

Figure1 gives the comparison of linear and nonlinear surfaces at the wind speed of 10 m/s. The Creamer waves are accurate to the third order. In order to highlight the front-back asymmetry of Lagrange waves, the generated sea surfaces contain only waves spreading forward (A(kn,t)=0forkn<0). We can see that both Cramer and Lagrange waves have more flattened troughs and sharper crests, which belongs to vertical skewness. In addition, Lagrange waves also have front-back asymmetry compared to Creamer waves. The crests tilt in the direction of waves spreading.Figure 1
Comparison of linear and nonlinear wave profiles.In terms of wave profiles, the difference between linear and nonlinear waves is very small. This weak difference is mainly reflected in the wave height and slope distribution and can significantly affect the radar echo Doppler spectrum. Figure2 gives the comparison of linear and nonlinear wave height and slope distribution. The comparison shows that the probability density of the nonlinear wave troughs is slightly higher than that of the nonlinear wave crests in spatial distribution, and the slope of the nonlinear waves has a smaller variance than the linear waves, all of which directly indicate that the nonlinear wave crests are sharper and nonlinear wave troughs are more flattened. In addition, the probability density of the positive slope of the Lagrange waves is higher than that of the negative slope, which is just because the Lagrange waves tilt in the direction of waves spreading. These nonlinear characteristics make the wavenumber and wave frequency no longer accurately satisfy the first-order gravity wave dispersion relation, which is beneficial to the remote sensing of sea state using radar because the nonlinear characteristics make the radar Doppler spectrum contain more abundant sea state information. However, the nonlinearity of sea waves also makes the Doppler spectrum of the sea surface more complicated.Figure 2
Comparison of linear and nonlinear wave height and slope distribution.
(a)
Wave height distribution (b)
Slope distribution


## 2.1. Linear Waves

The wave height of the linear sea surface satisfies the Gaussian distribution, which means that the linear waves can be regarded as a sum of a series of sine waves, and the phase and amplitude of each component obey the uniform distribution and Rayleigh distribution, respectively. The intensity of the sine wave with different wavenumbers is determined by the sea spectrum. Therefore, the linear gravity waves can be expressed as(1)hlx,t=Re∑nAkn,texp⁡jknx=Re∑nαnSkndkexp⁡-jwntexp⁡jknxwhere kn=ndk is the wavenumber, dk=2π/L depends on the rough surface length L, αn is the independent complex random Gauss variable with zero mean and unit variance, wn=gkn is the angular frequency of sea wave with a wavenumber of kn determined by the dispersion relation under the deep water condition, g is the gravitational acceleration, and S(kn) denotes the sea spectrum. The rough surface is sampled as N discrete points, and n is an integer in the interval of [-N/2,N/2-1]. Then the discrete sea surface profile can be calculated using the fast Fourier transform (FFT).In this paper, JONSWAP spectrum is employed to drive the sea waves, which is expressed by(2)Skw=0.038kw3U102gxw0.22exp⁡-1.25wpgkw4γrwhere kw is the water wavenumber, U10 represents the wind speed at 10 meters above the sea level, xw is a constant representing the range of wind area, γ is usually 3.3, and(3a)wp=22g2U10xw1/3(3b)r=exp⁡-gkw-wp22σ2wp2(3c)σ=0.07,gkw<wp0.09,gkw≥wpAs an inadequately developed spectrum, JONSWAP spectrum is more realistic than the PM spectrum. When xw increases, JONSWAP spectrum and PM spectrum will be closer.


## 2.2. Creamer Waves

The actual sea waves are nonlinear, and there are weak interactions between waves with different frequencies. These weak interactions can be considered as higher-order perturbation solutions to the hydraulic motion equation, which will significantly affect the characteristics of the Doppler spectrum, and the influence degree will also change with the sea state. Therefore, the nonlinearity of the sea waves cannot be neglected when studying the Doppler spectrum of the sea surface no matter which band the radar works in. The second-order nonlinear waves, generated by two wave interactions, have been described by [29, 30], which have been applied in the second-order RCS of ocean surface at HF band [31]. Although the perturbation expansion method is clear, it is not suitable for numerical implementation due to the high calculating burden especially for higher-order nonlinear waves. Another alternative method to generate nonlinear sea surfaces is the Hamiltonian formalism based on the weak wave-turbulence theory developed by [26], which has been widely used to study the electromagnetic scattering of 1D or 2D time-varying sea surfaces [3, 4].To begin with, the Creamer method expresses the nonlinear term of the sea surface as a Hilbert transform of the linear sea surface. In the 1D case, the Creamer nonlinear term is denoted as(4)htx,t=Re∑n-jknknAkn,texp⁡jknxwhere the complex amplitude A(kn,t) is as defined in (1). This transform can also be implemented using FFT to reduce the computation cost to Nlog⁡N. In the frequency domain, the complex amplitude of Creamer nonlinear waves can be expressed as(5)Ctkn,t=1N∑xexp⁡jknhtx,t-1knexp⁡-jknxSince (5) cannot be calculated using FFT directly, it needs approximate processing to improve computational efficiency. Expanding exp⁡jknhtx,t to Taylor series, themth expansion of (5) can be obtained as(6)Ctmkn,t=1N∑xjknhtx,tmm!knexp⁡-jknx=jknmm!NknFhtx,tmwhere F· denotes Fourier transform. Thenmth-order Creamer waves can be calculated using inverse Fourier transform, which is (7)hcmx,t=Re1NFi∑nCtmkn,tIt can be proved that the first-order Creamer waves just are linear waves.


## 2.3. Lagrange Waves

Both perturbation method and the Creamer method describe the vertical skewness of sea waves in the Euler coordinate system. There is horizontal skewness in the actual nonlinear waves, which directly affects the slope distribution of the sea surface and induces more remarkable influence on backscattering signals than vertical skewness [32, 33]. References [32, 33] have attempted to modify the perturbation model to add horizontal skewness on the nonlinear waves and discussed its effect on microwave scattering and emission. Recently [27, 34] proposed a more concise method derived in the Lagrange coordinate system to generate nonlinear waves with horizontal skewness also called Lagrange waves. A disturbance is added in the horizontal direction so that an additional horizontal offset is produced in the linear waves, which is(8)x,hlx,t→x+xLx,t,hlx,twhere hl(x,t) and xL(x,t) describe the vertical and horizontal movements of individual particles with time t, which are two correlated Gauss processes. The vertical process of hl(x,t) is given by (1), and horizontal process of xL(x,t) is completely determined by the vertical process, given by(9)xLx,t=Re∑nHknAkn,texp⁡jknxwhere H(kn) is the response function which determines the nonlinear characters of the Lagrange waves, in particular crest-trough and front-back asymmetry. In Lindgrens derivation [27], H(kn) consists of two parts:(10)Hkn=jtanh⁡hkn+sgnknβωn2where h represents the water depth, sgn· is the sign function, and β is obtained from the relation between the horizontal acceleration of water particles and vertical displacement [27]. Following the study of [27, 34], we set β=0.4 in this paper. The first part of H(kn) determines the vertical skewness and results in waves with more peaked crests and shallower troughs (crest-trough asymmetry) compared to the linear waves, and the second part of H(kn) determines the horizontal skewness and gives waves front-back asymmetry. When the water depth is assumed to be infinite, the response function can be simplified to Hkn=j+sgnknβ/ωn2.


## 2.4. Wave Height and Slope Distribution

Figure1 gives the comparison of linear and nonlinear surfaces at the wind speed of 10 m/s. The Creamer waves are accurate to the third order. In order to highlight the front-back asymmetry of Lagrange waves, the generated sea surfaces contain only waves spreading forward (A(kn,t)=0forkn<0). We can see that both Cramer and Lagrange waves have more flattened troughs and sharper crests, which belongs to vertical skewness. In addition, Lagrange waves also have front-back asymmetry compared to Creamer waves. The crests tilt in the direction of waves spreading.Figure 1
Comparison of linear and nonlinear wave profiles.In terms of wave profiles, the difference between linear and nonlinear waves is very small. This weak difference is mainly reflected in the wave height and slope distribution and can significantly affect the radar echo Doppler spectrum. Figure2 gives the comparison of linear and nonlinear wave height and slope distribution. The comparison shows that the probability density of the nonlinear wave troughs is slightly higher than that of the nonlinear wave crests in spatial distribution, and the slope of the nonlinear waves has a smaller variance than the linear waves, all of which directly indicate that the nonlinear wave crests are sharper and nonlinear wave troughs are more flattened. In addition, the probability density of the positive slope of the Lagrange waves is higher than that of the negative slope, which is just because the Lagrange waves tilt in the direction of waves spreading. These nonlinear characteristics make the wavenumber and wave frequency no longer accurately satisfy the first-order gravity wave dispersion relation, which is beneficial to the remote sensing of sea state using radar because the nonlinear characteristics make the radar Doppler spectrum contain more abundant sea state information. However, the nonlinearity of sea waves also makes the Doppler spectrum of the sea surface more complicated.Figure 2
Comparison of linear and nonlinear wave height and slope distribution.
(a)
Wave height distribution (b)
Slope distribution


## 3. Formulation of the Electromagnetic Scattering

Limited by the edge effect of finite rough surface, the MoM cannot directly deal with the rough surface scattering problem under near-grazing incidence. Although adopting the tapered wave as incident wave can solve the large incident angle problem well, the computation load increases exponentially with the increase of the incident angle [13]. And in theory, the tapered wave cannot handle the complete grazing incidence problem. Other processing methods, such as using the impedance rough surfaces [20] and the periodic boundary condition [19], are not robust or have not been validated. Following the GMoM [21], in this section, we first treat the rough surface as a locally perturbed plane and then use the impedance boundary condition and the plane incident wave to simulate the grazing incident backscatter Doppler spectrum.


### 3.1. Construction of Locally Perturbed Surface

In this paper, we mainly study the characteristics of backscattering Doppler spectrum from one-dimensional time-varying ocean-like rough surface. Considering that most sea state monitoring radar is deployed on the coast and operates in vertical polarization mode, we focus on the vertical polarization incident wave at the grazing incidence. When using numerical methods to analyze rough surface scattering problems, the rough surface needs to be truncated to finite area. On the one hand, the discontinuous edge current of the truncated rough surface will affect the distribution of the scattered field. On the other hand, the scattering field can be approximated as the Fourier transform of the rough surface current, and the truncated surface current will generate side lobe at the undesired scattering angles, which is similar to the spectrum leakage.In order to deal with the edge effect, the rough surface is seen as a local perturbation on an infinite plane. We window the truncated rough surface with a Tukey-Hanning window function to construct the locally perturbed surface. As shown in Figure3, the width of the random and time-varying rough sea surface is L, and its roughness is determined by the wave spectrum and wind speed. The Tukey-Hanning window function consists of three parts: rectangular window part in [-L2/2,L2/2], Hanning window part in [-L1/2,-L2/2] and [L2/2,L1/2], and zero part in the both ends of the rough surface, which can be expressed as(11)wx=121+cos2πx+L2/2L1-L2,-L12≤x<-L221,-L22≤x<L22121+cos2πx-L2/2L1-L2,L22≤x<L120,elseThe windowed rough surface f(x,t) for the MoM is(12)fx,t=wxhl,c,Lx,twhere hl(x,t), hc(x,t), and hL(x,t) are linear wave, Creamer wave, and Lagrange wave, respectively. The windowed rough surface contains three areas: scattering area with a width of L2, asymptotical area with a width of L1-L2, and plane area with a width of L-L1. The edges of the windowed rough surface are asymptotically flat, and the roughness near the edges is zero. The scattering field is mainly excited by the surface currents in the scattering area, and the width of L2 is recommended to be greater than 32λ. Asymptotical area ensures the continuity of the surface currents and its width L1-L2 should be greater than 6λ. Although the roughness of the plane area is zero, there are still incoherent surface currents produced by the near-field interaction of current source in the asymptotical area. Considering the range of near-field interaction will not exceed λ [35], so its width L-L1 should be at least set to 2λ. The coherent currents always exist on the whole infinite surface, but they are known according to Fresnel formula. The incoherent current near the edges of the windowed rough surface is asymptotically reduced to zero, which implies that the current is also windowed, although this window is not strictly analytic and constant. The windowed currents will effectively suppress the side lobe of the scattering wave.Figure 3
Geometry of the scattering problem and the locally perturbed surface.


### 3.2. Doppler Spectrum Simulation

For vertically polarized wave, the magnetic field has only component in they^ direction. Therefore, the magnetic field can be represented as(13)H¯r¯,t=y^ψr¯,t=y^ψincr¯,t+ψrr¯,t+ψsr¯,twhere r¯=xx^+zz^, ψr¯,t is the total scalar magnetic field, ψincr¯,t=exp⁡ikixx-ikizz is the incident scalar magnetic field, ψrr¯,t=Rexp⁡ikixx+ikizz is the reflected scalar magnetic field, ψsr¯,t is the scattering scalar magnetic field, and the e-iωt time dependence is assumed and suppressed. The Fresnel reflection coefficient(14)R=k2/k1-kizkiz+k2/k1determines the intensity of reflected wave, where kix=ksin⁡θi, kiz=kcos⁡θi, k is the incident radio wavenumber in the vacuum, θi is the incident angle as shown in Figure 3, k1=ω2μ0ε1~ is the radio wavenumber in sea water, and ε1~=80ε0+i4/ω is the complex permittivity of sea water. Note that J=n^×H¯=ψ and M=E¯×n^=(n^·∇ψ)/iωε0 are the total scalar surface electric current and magnetic current, Jr=ψinc+ψr and Mr=(n^·∇ψinc+ψr)/iωε0 are the coherent currents, and Js=ψs and Ms=(n^·∇ψs)/iωε0 are the incoherent currents.In general, after the surface equivalent electric currentψr¯,t and magnetic current n^·∇ψr¯,t are calculated using the integral equation method (IEM) and the moment method [14], the scattering field can be obtained directly by using the Huygens principle, i.e.,(15)ψrr¯,t+ψsr¯,t=∫Cdl′Jr¯′,tn^·∇Gr¯,r¯′-iωε0Gr¯,r¯′Mr¯′,tfor r¯′ and r¯′ on the rough surface, where the integral region is C on the rough surface fx,t and n^=-fxx^+z^/1+fx2 is the unit normal vector of the rough surface. The two-dimensional Green function Gr¯,r¯′=iH01kr¯-r¯′/4, where H01kr¯-r¯′ is the Hankel function of the first kind. Note that the scattered field calculated by IEM directly includes both reflected and scattered waves, and it is affected by the truncation of finite rough surface, which must be eliminated by tapering the incident wave to make it have a Gaussian footprint on the rough surface [14], i.e.,(16)ψincr¯=exp⁡ikxsin⁡θi-zcos⁡θi1+wr¯×exp⁡-x+ztan⁡θi2q2where q is the tapering parameter, which must be greater than 6λ/cos1.5θi to make the tapering wave satisfy the Helmholtz equation [13], and also the rough surface size is usually set to 4q, which will increase rapidly with the increase of incident angle. One can see that the calculation amount of solving the IEM using the moment method is enormous when approaching grazing incidence, and it is theoretically impossible to deal with the case of full grazing incidence.Now we consider the directly truncated rough surface as a perturbation on an infinite plane and use the Tukey-Hanning window to restrain the rough surface so that its roughness fades to zero at the edge, as described in (11) and (12). After such treatment, the integral region of (15) can be extended to negative infinity to positive infinity. Next substituting J=Jr+Js and M=Mr+Ms into (15) and according to the linear characteristics of the integral equation, we have(17)ψrr¯,t+ψsr¯,t=∫-∞+∞dl′Jsr¯′,tn^·∇Gr¯,r¯′-iωε0Gr¯,r¯′Msr¯′,t+∫-∞+∞dl′Jrr¯′,tn^·∇Gr¯,r¯′-iωε0Gr¯,r¯′Mrr¯′,tObviously, coherent surface currents Jr and Mr are known, and ψrr¯,t=∫-∞∞dl′Jrn^·∇Gr¯,r¯′-iωε0Gr¯,r¯′Mr. Therefore, eliminating the reflected wave in the scattered field, we immediately have(18)ψsr¯,t=∫Cdl′ψsr¯′,tn^·∇Gr¯,r¯′-Gr¯,r¯′n^·∇ψsr¯′,twhere the integral region returns to the rough surface again because the incoherent scattering currents only exist on the rough surface. It is predictable that there are only the coherent currents associated with the reflected wave when the rough surface is an infinite plane. There is a singular value in ∇Gr¯,r¯′ and r¯=r¯′ which can be approximated using the principal value integral [36]. Then (18) can be expressed in an abbreviated form:(19)12-Dψsr¯′,t+Sn^·∇ψsr¯′,t=0where(20)D=⨏Cdl′n^·∇Gr¯,r¯′,S=∫Cdl′Gr¯,r¯′where ⨏ denotes the principal value integral with a vanishingly small piece subtracted out from the domain of the integration. The electric current ψsr¯,t and the magnetic current n^·∇ψsr¯,t are the surface unknowns, which cannot be solved directly from the single integral (19). Considering that seawater is a good conductor, the impedance boundary condition n^×E¯=-μ/ε1~H¯ can be employed to establish the relationship between ψsr¯,t and n^·∇ψsr¯,t, which is(21)n^·∇ψr¯,t=-ik2k1ψr¯,tSubstituting (13) into (21), we have(22)n^·∇ψs=-ik2k1ψs-ik2k1ψinc+ψr-n^·∇ψinc+ψrThe simultaneous (19) and (22) give(23)12-D-ik2k1Sψsr¯′,t=Sn^′·∇′+ik2k1×ψincr¯′,t+ψrr¯′,tBecause ψincr¯′,t and ψrr¯′,t are known, the surface incoherent electric current ψsr¯,t can be solved directly by moment method, and then the surface magnetic current n^·∇ψsr¯,t is determined from (22). The scattering far field can be obtained from the surface current, which is [36](24)ψsr¯,t=i42πkreikr-π/4ψsNθs,twhere r=r¯ and θs are the distance and scattering angle between the rough surface and the observed point, and(25)ψsNθs,t=∫Cdxe-iksin⁡θs+fx,tcos⁡θs×n^·∇ψsrs¯,t1+fx2x,t-ψsrs¯,tikfxx,tsin⁡θs-cos⁡θswhere rs¯=xx^+fx,tz^ is on the rough surface. The interaction of the radio waves with the time-varying ocean-like rough surface can produce the Doppler effect, which means that the scattering field ψsr¯,t changes over time. Because the phase velocity of the sea wave can be negligibly relative to that of the radio wave, the rough sea surface can be considered frozen when calculating the scattering far field. The continuous time t is discretized as tn=nΔT,n=0,1,…N-1. The discrete periodic ΔT=π/ωmax is adapted to the maximum simulation Doppler frequency of ωmax. For each rough surface frozen at tn, we can obtain a scattering wave ψsr¯,tn at the observed point of r¯. Finally, the Doppler spectrum can be calculated by the discrete Fourier transform (DFT):(26)Ψwk=∑n=0N-1ψsr¯,tne-i2π/Nnwk,wk=0,1,…,N-1Since the rough sea surface is random, one-time simulation of the Doppler spectrum is rough. We need to take the average of Monte-Carlo simulation of the Doppler spectrum of the independent random sea surfaces as the final result.As shown in Figure4, all the processing steps are summarized as follows:(1)
Generate time-varying linear waves, Creamer waves, and Lagrange waves(2)
Use the Tukey-Hanning window to restrain the rough surfaces to construct local perturbation planes(3)
Employ the moment method to solve the improved magnetic field integral equation to calculate incoherent backscattered waves(4)
The time-ordered sequence of backscattered waves is collected to produce the Doppler spectrum using Fourier transform(5)
Repeat the above steps to obtain multiple Doppler spectra for smoothingFigure 4
Sketch map of the Doppler spectrum simulation.


## 3.1. Construction of Locally Perturbed Surface

In this paper, we mainly study the characteristics of backscattering Doppler spectrum from one-dimensional time-varying ocean-like rough surface. Considering that most sea state monitoring radar is deployed on the coast and operates in vertical polarization mode, we focus on the vertical polarization incident wave at the grazing incidence. When using numerical methods to analyze rough surface scattering problems, the rough surface needs to be truncated to finite area. On the one hand, the discontinuous edge current of the truncated rough surface will affect the distribution of the scattered field. On the other hand, the scattering field can be approximated as the Fourier transform of the rough surface current, and the truncated surface current will generate side lobe at the undesired scattering angles, which is similar to the spectrum leakage.In order to deal with the edge effect, the rough surface is seen as a local perturbation on an infinite plane. We window the truncated rough surface with a Tukey-Hanning window function to construct the locally perturbed surface. As shown in Figure3, the width of the random and time-varying rough sea surface is L, and its roughness is determined by the wave spectrum and wind speed. The Tukey-Hanning window function consists of three parts: rectangular window part in [-L2/2,L2/2], Hanning window part in [-L1/2,-L2/2] and [L2/2,L1/2], and zero part in the both ends of the rough surface, which can be expressed as(11)wx=121+cos2πx+L2/2L1-L2,-L12≤x<-L221,-L22≤x<L22121+cos2πx-L2/2L1-L2,L22≤x<L120,elseThe windowed rough surface f(x,t) for the MoM is(12)fx,t=wxhl,c,Lx,twhere hl(x,t), hc(x,t), and hL(x,t) are linear wave, Creamer wave, and Lagrange wave, respectively. The windowed rough surface contains three areas: scattering area with a width of L2, asymptotical area with a width of L1-L2, and plane area with a width of L-L1. The edges of the windowed rough surface are asymptotically flat, and the roughness near the edges is zero. The scattering field is mainly excited by the surface currents in the scattering area, and the width of L2 is recommended to be greater than 32λ. Asymptotical area ensures the continuity of the surface currents and its width L1-L2 should be greater than 6λ. Although the roughness of the plane area is zero, there are still incoherent surface currents produced by the near-field interaction of current source in the asymptotical area. Considering the range of near-field interaction will not exceed λ [35], so its width L-L1 should be at least set to 2λ. The coherent currents always exist on the whole infinite surface, but they are known according to Fresnel formula. The incoherent current near the edges of the windowed rough surface is asymptotically reduced to zero, which implies that the current is also windowed, although this window is not strictly analytic and constant. The windowed currents will effectively suppress the side lobe of the scattering wave.Figure 3
Geometry of the scattering problem and the locally perturbed surface.


## 3.2. Doppler Spectrum Simulation

For vertically polarized wave, the magnetic field has only component in they^ direction. Therefore, the magnetic field can be represented as(13)H¯r¯,t=y^ψr¯,t=y^ψincr¯,t+ψrr¯,t+ψsr¯,twhere r¯=xx^+zz^, ψr¯,t is the total scalar magnetic field, ψincr¯,t=exp⁡ikixx-ikizz is the incident scalar magnetic field, ψrr¯,t=Rexp⁡ikixx+ikizz is the reflected scalar magnetic field, ψsr¯,t is the scattering scalar magnetic field, and the e-iωt time dependence is assumed and suppressed. The Fresnel reflection coefficient(14)R=k2/k1-kizkiz+k2/k1determines the intensity of reflected wave, where kix=ksin⁡θi, kiz=kcos⁡θi, k is the incident radio wavenumber in the vacuum, θi is the incident angle as shown in Figure 3, k1=ω2μ0ε1~ is the radio wavenumber in sea water, and ε1~=80ε0+i4/ω is the complex permittivity of sea water. Note that J=n^×H¯=ψ and M=E¯×n^=(n^·∇ψ)/iωε0 are the total scalar surface electric current and magnetic current, Jr=ψinc+ψr and Mr=(n^·∇ψinc+ψr)/iωε0 are the coherent currents, and Js=ψs and Ms=(n^·∇ψs)/iωε0 are the incoherent currents.In general, after the surface equivalent electric currentψr¯,t and magnetic current n^·∇ψr¯,t are calculated using the integral equation method (IEM) and the moment method [14], the scattering field can be obtained directly by using the Huygens principle, i.e.,(15)ψrr¯,t+ψsr¯,t=∫Cdl′Jr¯′,tn^·∇Gr¯,r¯′-iωε0Gr¯,r¯′Mr¯′,tfor r¯′ and r¯′ on the rough surface, where the integral region is C on the rough surface fx,t and n^=-fxx^+z^/1+fx2 is the unit normal vector of the rough surface. The two-dimensional Green function Gr¯,r¯′=iH01kr¯-r¯′/4, where H01kr¯-r¯′ is the Hankel function of the first kind. Note that the scattered field calculated by IEM directly includes both reflected and scattered waves, and it is affected by the truncation of finite rough surface, which must be eliminated by tapering the incident wave to make it have a Gaussian footprint on the rough surface [14], i.e.,(16)ψincr¯=exp⁡ikxsin⁡θi-zcos⁡θi1+wr¯×exp⁡-x+ztan⁡θi2q2where q is the tapering parameter, which must be greater than 6λ/cos1.5θi to make the tapering wave satisfy the Helmholtz equation [13], and also the rough surface size is usually set to 4q, which will increase rapidly with the increase of incident angle. One can see that the calculation amount of solving the IEM using the moment method is enormous when approaching grazing incidence, and it is theoretically impossible to deal with the case of full grazing incidence.Now we consider the directly truncated rough surface as a perturbation on an infinite plane and use the Tukey-Hanning window to restrain the rough surface so that its roughness fades to zero at the edge, as described in (11) and (12). After such treatment, the integral region of (15) can be extended to negative infinity to positive infinity. Next substituting J=Jr+Js and M=Mr+Ms into (15) and according to the linear characteristics of the integral equation, we have(17)ψrr¯,t+ψsr¯,t=∫-∞+∞dl′Jsr¯′,tn^·∇Gr¯,r¯′-iωε0Gr¯,r¯′Msr¯′,t+∫-∞+∞dl′Jrr¯′,tn^·∇Gr¯,r¯′-iωε0Gr¯,r¯′Mrr¯′,tObviously, coherent surface currents Jr and Mr are known, and ψrr¯,t=∫-∞∞dl′Jrn^·∇Gr¯,r¯′-iωε0Gr¯,r¯′Mr. Therefore, eliminating the reflected wave in the scattered field, we immediately have(18)ψsr¯,t=∫Cdl′ψsr¯′,tn^·∇Gr¯,r¯′-Gr¯,r¯′n^·∇ψsr¯′,twhere the integral region returns to the rough surface again because the incoherent scattering currents only exist on the rough surface. It is predictable that there are only the coherent currents associated with the reflected wave when the rough surface is an infinite plane. There is a singular value in ∇Gr¯,r¯′ and r¯=r¯′ which can be approximated using the principal value integral [36]. Then (18) can be expressed in an abbreviated form:(19)12-Dψsr¯′,t+Sn^·∇ψsr¯′,t=0where(20)D=⨏Cdl′n^·∇Gr¯,r¯′,S=∫Cdl′Gr¯,r¯′where ⨏ denotes the principal value integral with a vanishingly small piece subtracted out from the domain of the integration. The electric current ψsr¯,t and the magnetic current n^·∇ψsr¯,t are the surface unknowns, which cannot be solved directly from the single integral (19). Considering that seawater is a good conductor, the impedance boundary condition n^×E¯=-μ/ε1~H¯ can be employed to establish the relationship between ψsr¯,t and n^·∇ψsr¯,t, which is(21)n^·∇ψr¯,t=-ik2k1ψr¯,tSubstituting (13) into (21), we have(22)n^·∇ψs=-ik2k1ψs-ik2k1ψinc+ψr-n^·∇ψinc+ψrThe simultaneous (19) and (22) give(23)12-D-ik2k1Sψsr¯′,t=Sn^′·∇′+ik2k1×ψincr¯′,t+ψrr¯′,tBecause ψincr¯′,t and ψrr¯′,t are known, the surface incoherent electric current ψsr¯,t can be solved directly by moment method, and then the surface magnetic current n^·∇ψsr¯,t is determined from (22). The scattering far field can be obtained from the surface current, which is [36](24)ψsr¯,t=i42πkreikr-π/4ψsNθs,twhere r=r¯ and θs are the distance and scattering angle between the rough surface and the observed point, and(25)ψsNθs,t=∫Cdxe-iksin⁡θs+fx,tcos⁡θs×n^·∇ψsrs¯,t1+fx2x,t-ψsrs¯,tikfxx,tsin⁡θs-cos⁡θswhere rs¯=xx^+fx,tz^ is on the rough surface. The interaction of the radio waves with the time-varying ocean-like rough surface can produce the Doppler effect, which means that the scattering field ψsr¯,t changes over time. Because the phase velocity of the sea wave can be negligibly relative to that of the radio wave, the rough sea surface can be considered frozen when calculating the scattering far field. The continuous time t is discretized as tn=nΔT,n=0,1,…N-1. The discrete periodic ΔT=π/ωmax is adapted to the maximum simulation Doppler frequency of ωmax. For each rough surface frozen at tn, we can obtain a scattering wave ψsr¯,tn at the observed point of r¯. Finally, the Doppler spectrum can be calculated by the discrete Fourier transform (DFT):(26)Ψwk=∑n=0N-1ψsr¯,tne-i2π/Nnwk,wk=0,1,…,N-1Since the rough sea surface is random, one-time simulation of the Doppler spectrum is rough. We need to take the average of Monte-Carlo simulation of the Doppler spectrum of the independent random sea surfaces as the final result.As shown in Figure4, all the processing steps are summarized as follows:(1)
Generate time-varying linear waves, Creamer waves, and Lagrange waves(2)
Use the Tukey-Hanning window to restrain the rough surfaces to construct local perturbation planes(3)
Employ the moment method to solve the improved magnetic field integral equation to calculate incoherent backscattered waves(4)
The time-ordered sequence of backscattered waves is collected to produce the Doppler spectrum using Fourier transform(5)
Repeat the above steps to obtain multiple Doppler spectra for smoothingFigure 4
Sketch map of the Doppler spectrum simulation.


## 4. Results and Discussion

In this section, we will discuss in detail the characteristics of the Doppler spectrum backscattering from one-dimensional linear and nonlinear time-varying ocean-like surfaces at grazing incidence. The simulation frequencies range from HF to UHF bands. Table1 gives the detailed simulation conditions. Unless otherwise stated, all results in this section are simulated under these conditions. Before analysis and discussion, the Doppler spectra from linear wave are compared with the one derived by the SPM to verify the effectiveness of the proposed numerical method at grazing incidence. At the same time, we will analyze the applicable range of the SPM in the Doppler domain. Then we discuss the characteristics of the Doppler spectrum in HF/VHF/UHF bands and the effect of nonlinear waves on the Doppler spectrum.Table 1
Simulation conditions.
Parameters Value Simulation frequency HF to UHF bands Incident angle 90∘ Scattering angle - 9 0 ∘ (Backscattering) Δx λ / 24 (λ is the radio wavelength) Rough surface widthL 1024Δx Scattering area widthL2 816Δx Asymptotical area widthL1-L2 144Δx Plane area widthL-L1 64Δx Maximum Doppler frequencyωmax 2.8 ω b (ωb=2gksin⁡θi is the Bragg frequency) Simulation frameN 1024 Monte Carlo simulation times 32


### 4.1. Surface Incoherent Current Distribution

Figure5 shows the rough surface current distribution over the sea wave profile in HF band. As expected, the incoherent scattered currents exist only in the region of the local perturbation plane where the roughness is not zero. At the edges of the local perturbation surface, the current is asymptotically reduced to zero. The envelope of the current can be equivalent to being windowed, which will effectively suppress the side lobes of the scattering waves. In sharp contrast, the total surface current calculated using the conventional integral equation method contains both coherent and incoherent components, and the coherent current, which excites the reflected wave, occupies the vast majority. Moreover, the surface current near the edge of rough surface is inaccurate due to the edge effect. The surface integral equation in this paper is established on the unbounded local perturbation plane, and the incoherent scattered current is separated from the total current as the surface unknowns for calculating the scattering field without reflected waves. The effectiveness of this method at grazing incidence will be verified by comparison with the SPM.Figure 5
Surface current distribution over wave profile. (a) The incoherent surface currentψs calculated by the method in this paper. (b) The total surface current ψ calculated by the conventional integral equation method with incident plane wave.


### 4.2. Doppler Spectrum Comparison with the SPM

The SPM is widely recognized for analyzing sea surface scattering problems in HF band. The first- and second-order backscattering RSC of the rough sea surface have been derived by Barrick using the SPM [1, 2], which is the theoretical basis for remote sensing of the sea state using HF radar. The perturbation solution of the scattered field is not restricted by the incidence angle, even at grazing incidence. Therefore, it is reasonable to compare the Doppler spectrum with the SPM to verify the effectiveness of the simulation method. For the sake of simplicity, seawater is assumed to be perfect electric conductor and the linear waves are used. The first- and second-order RCS of the one-dimensional linear rough sea surface are given by [37](27a)σvv1θs=4πk3sin⁡θssin⁡θi-12Sksin⁡θs-ksin⁡θi(27b)σvv2θs=2πk3∫-∞+∞S-ksin⁡θi-qSq-ksin⁡θi×q2sin2θi-k2k2-q2dqThe third- and higher-order perturbation scattering waves are ignored. For backscattering at grazing incidence, θs=-θi=-90∘. In addition, according to the first-order gravitational wave dispersion relation, a more general spectrum Skw,ω including wave frequency can be expressed by the undirectional spectrum Skw[2]:(28)Skw,ω=Skwδω-sgnkwgkw+Skwδω+sgnkwgkwwhere δω is the Dirac impulse function of the variable ω and sgnkw=±1 depends on the ± of variable kw. The two items on the right of (28) represent the spatial/temporal spectrum of the waves moving into the +x space and the waves moving into the -x space, respectively. Putting (28) into (27a) and (27b), the backscattering Doppler spectrum at grazing incidence derived by the SPM is obtained:(29)ΨSPMω=∑m=±116πk3S2mkδω-mgk+2πk3∑m1,m2=±1∫-∞+∞Sm1k1Sm2k2k2-q2δω-m1gk1-m2gk2dqwhere k1=-k-q and k2=q-k. Note that ΨSPMω is calculated from the RCS, so it is the power spectrum.Figure6 gives the Doppler spectrum comparison between the simulation and the SPM in different sea states. We can see that the simulated Doppler spectra are highly consistent with that calculated by the SPM. As we all know, the sea surface backscattering Doppler spectrum is composed of first-order peak, also called Bragg peak, and high-order peaks. In the HF band, the third- and higher-order echo intensity is small enough to be ignored. Note that the simulation uses linear waves, which means the high-order peaks in the Figure 6 are all from the high-order electromagnetic (EM) scattering. In fact, the interaction between the gravity waves will produce nonlinear waves with high-order dispersion relation [30]. If the wavelengths of the produced waves satisfy the Bragg scattering condition, they will have a nonnegligible contribution to the backscattering echo intensity, which will be discussed in detail in the following sections. The high-order peaks from high-order EM scattering and nonlinear sea waves are attributed to the EM contribution and hydrodynamic contribution by [2]. The purpose of comparison with the SPM is to verify the effectiveness of the numerical method in this paper at grazing incidence, so at this point, we do not need to consider the hydrodynamic contribution. From Figure 6 we can also see that the high-order peaks are more sensitive to the wind speed in the low sea state, and when the wind speed exceeds 10 m/s, the electromagnetic contribution to the Doppler spectrum almost does not change. In addition, there is an obvious peak at 2ωb in low sea state. All of these characteristics are consistently reflected in the simulated and calculated Doppler spectra.Figure 6
Comparison of Doppler spectrum between the numerical simulation and the SPM in different sea states. The frequency is 10 MHz and the sea surfaces are linear. The Doppler frequency is normalized by the Bragg frequencyωb, and the spectral intensity is normalized by Bragg peak.
(a)
Wind speed 5 m/s (b)
Wind speed 10 m/s (c)
Wind speed 15 m/s (d)
Wind speed 20 m/sThe consistency of numerical simulation and the SPM in HF band shows that the method in this paper is effective for grazing incidence scattering. Next we will demonstrate the effective frequency range of the SPM by numerical simulation. At the same time, in order to verify the effectiveness of the proposed method in VHF and UHF bands, the classical integral equation method (IEM) is also employed for comparison, which calculates the surface equivalent currents by directly solving the EM integral equation using the moment method, and uses tapered incident wave, expressed by (16), to suppress the discontinuity of surface currents at the edges. Since the tapering parameter q must be greater than 6λ/cos1.5θi to make the tapered wave satisfy Helmholtz equation  [13],the width of the rough surface will increase rapidly with the incident angle, which results in a huge computation amount of classical method at low grazing angle. Here we set the grazing angle of 2∘, the tapering parameter of 1024λ, and the rough surface width of 4096λ for the classical method. The simulation conditions in Table 1 are still applicable for the proposed method. The comparison results of Doppler spectrum of linear waves among the proposed method, the classical integral equation method, and the SPM in different bands are shown in Figure 7. As we can see, the Doppler spectra simulated by the proposed method are almost coincident with that of the classical IEM from HF to UHF bands but with less computational amount, while with the increase of the operating frequency, the intensity of the Doppler spectrum calculated by the SPM is gradually smaller than that of the numerical simulation, which is because, on the one hand, the higher-order perturbation solution to scattering waves cannot be ignored as the frequency increases, and on the other hand, the gradually decreasing radio wavelength does not satisfy the applicable condition of the SPM. For the former, one might be able to improve the accuracy of the SPM by retaining the third- or higher-order perturbation solutions, but this is very troublesome because the form of the higher-order perturbation solutions is too complicated. For the latter, this is restricted by the condition of the SPM and cannot be solved. The mean absolute error of the Doppler spectral intensity over operating frequency between the SPM and the numerical method is shown in Figure 8, which provides a new perspective to measure the applicable frequency of the SPM in the Doppler domain. The Doppler spectral error of the SPM in HF band is very small and can be ignored, but it grows rapidly in VHF band.Figure 7
Comparison of Doppler spectrum of linear waves among the proposed method, the classical integral equation method, and the SPM in different frequency bands. The wind speed is 10 m/s, the Doppler frequency is normalized by the Bragg frequencyωb, and the spectral intensity is normalized by Bragg peak.
(a)
10 MHz (b)
30 MHz (c)
60 MHz (d)
100 MHz (e)
200 MHz (f)
400 MHzFigure 8
The mean absolute error of the Doppler spectral intensity over operating frequency between the SPM and the numerical method. The sea waves are linear and the wind speed is 10 m/s.From Figure7, we can roughly see that the Doppler spectra of linear waves are different at different operating frequencies. Especially at frequencies below 100 MHz, the Doppler spectrum changes more significantly with operating frequency. Figure 9 gives a more detailed description of the numerical simulated Doppler spectrum of linear waves in different frequencies. As the frequency increases, the Doppler frequency and intensity of the higher-order peaks gradually approach the Bragg peak, and the approaching rate in the high waves is much greater than that in the low waves. The approaching rate gradually slows down as the frequency increases, and the higher-order peaks are almost overlapped with the Bragg peaks when the frequency exceeds 200 MHz. There is no significant difference in the shape of the Doppler spectrum in UHF band, which means that the EM contribution tends to stabilize with the increase of frequency.Figure 9
Numerical simulated Doppler spectrum of linear waves over operating frequency. The horizontal axis is the Doppler angular frequency normalized by the Bragg frequency, and the vertical axis represents different simulation operating frequencies from 5 MHz to 400 MHz.
(a)
Wind speed 3 m/s (b)
Wind speed 7 m/s


### 4.3. Doppler Spectrum of Linear and Nonlinear Waves

The actual sea surface is nonlinear, and the nonlinear components generated by the interaction of gravity waves will introduce hydrodynamic contribution in Doppler spectrum. Figure10 gives the comparison of Doppler spectrum among linear waves, Creamer waves, and Lagrange waves in different frequencies. We can see that, in all bands from HF to UHF, the high-order peak intensity of nonlinear waves is enhanced relative to that of linear waves, but the hydrodynamic contribution does not affect the Bragg peak. In general, the high-order peaks of the Lagrange waves and the Creamer waves are very close in the Doppler frequency range of [-1.5ωb,1.5ωb], and the high-order peak intensity of the Lagrange waves is gradually higher than that of the Creamer waves with the increase of the Doppler frequency, which implies that the nonlinear characteristics of the Lagrange waves are stronger than the Creamer waves. Consistent with the linear wave, after entering the VHF band, the shape of the Doppler spectrum of the nonlinear wave hardly changes with the increase of the operating frequency, except that the high-order peak intensity of the Creamer waves gradually approaches that of the Lagrange waves, which on the one hand shows that both EM and hydrodynamic contributions are stable in the VHF and UHF bands, and on the other hand, the frequency of the water wave interacting with the electromagnetic wave gradually increases with the increase of the operating frequency, and at this point the nonlinear characteristics of the Lagrange waves and the Creamer waves will gradually approach.Figure 10
Doppler spectrum comparison between linear and nonlinear waves in HF/VHF/UHF bands. The wind speed is 10 m/s, and the Doppler frequency is normalized by the Bragg frequencyωb.
(a)
10 MHz (b)
60 MHz (c)
200 MHz (d)
600 MHzIn the issue of remote sensing of sea state using radar, one is more concerned about the characteristics of echo Doppler spectrum in different sea states. Figure11 shows the variation of the Doppler spectrum of nonlinear waves with the sea state in different operating frequencies. We can see that, in the HF band, high-order peaks are very sensitive to the sea state, especially in the low- and medium-sea state, and the Bragg peak intensity, however, hardly changes with the sea state, which are consistent with Barrick’s sea surface first- and second-order scattering theory [1, 2]. Typically, HF radar uses the intensity changes of higher-order peaks relative to the Bragg peak to invert the wave height [5]. When the wind speed exceeds 15 m/s, the Doppler spectrum of the Creamer wave is almost no longer changed, which is explained by the fact that the sea waves interacting with the electromagnetic waves reach saturation. It is obvious that the Lagrange waves reach saturation more slowly. As the operating frequency increases, the wavelength of the sea waves interacting with the electromagnetic waves becomes smaller, and thus more easily saturated. As shown in Figures 11(c) and 11(f), the Doppler spectra in the UHF band have reached stability at very low wind speeds, which is actually not conducive to inversing the wave height through the Doppler spectrum characteristics. Another obvious feature is that the frequencies of the high-order peaks are very close to that of Bragg peak in the VHF and UHF bands. Considering the frequency shifts caused by surface currents in the actual ocean, the higher-order peaks will overlap with Bragg peak, and therefore the method of extracting wave height in the HF band will be ineffective in the VHF and UHF bands.Figure 11
Doppler spectrum of nonlinear waves in different sea states and different operating frequencies. The horizontal axis represents the Doppler frequency normalized by the Bragg frequencyωb, the vertical axis is the wind speed at 10 meters above sea level, and the color depth represents Doppler spectrum intensity.
(a)
Creamer waves at 10 MHz (b)
Creamer waves at 60 MHz (c)
Creamer waves at 600 MHz (d)
Lagrange waves at 10 MHz (e)
Lagrange waves at 60 MHz (f)
Lagrange waves at 600 MHzA more intuitive comparison of the Doppler spectra in the HF and UHF bands at different wind speeds is shown in Figure12. It is obvious that the shape of the Doppler spectrum in the UHF band changes very slightly with the wind speed relative to that in the HF band. The high-order peaks in the UHF band are continuously widened as the wind speed increases, and the intensity is also slightly enhanced, which may provide an idea for remote sensing of the wave height using UHF radar.Figure 12
Doppler spectrum of Creamer waves in different sea states and different operating frequencies. The Doppler frequency is normalized by the Bragg frequencyωb.
(a)
10 MHz (b)
600 MHz


## 4.1. Surface Incoherent Current Distribution

Figure5 shows the rough surface current distribution over the sea wave profile in HF band. As expected, the incoherent scattered currents exist only in the region of the local perturbation plane where the roughness is not zero. At the edges of the local perturbation surface, the current is asymptotically reduced to zero. The envelope of the current can be equivalent to being windowed, which will effectively suppress the side lobes of the scattering waves. In sharp contrast, the total surface current calculated using the conventional integral equation method contains both coherent and incoherent components, and the coherent current, which excites the reflected wave, occupies the vast majority. Moreover, the surface current near the edge of rough surface is inaccurate due to the edge effect. The surface integral equation in this paper is established on the unbounded local perturbation plane, and the incoherent scattered current is separated from the total current as the surface unknowns for calculating the scattering field without reflected waves. The effectiveness of this method at grazing incidence will be verified by comparison with the SPM.Figure 5
Surface current distribution over wave profile. (a) The incoherent surface currentψs calculated by the method in this paper. (b) The total surface current ψ calculated by the conventional integral equation method with incident plane wave.


## 4.2. Doppler Spectrum Comparison with the SPM

The SPM is widely recognized for analyzing sea surface scattering problems in HF band. The first- and second-order backscattering RSC of the rough sea surface have been derived by Barrick using the SPM [1, 2], which is the theoretical basis for remote sensing of the sea state using HF radar. The perturbation solution of the scattered field is not restricted by the incidence angle, even at grazing incidence. Therefore, it is reasonable to compare the Doppler spectrum with the SPM to verify the effectiveness of the simulation method. For the sake of simplicity, seawater is assumed to be perfect electric conductor and the linear waves are used. The first- and second-order RCS of the one-dimensional linear rough sea surface are given by [37](27a)σvv1θs=4πk3sin⁡θssin⁡θi-12Sksin⁡θs-ksin⁡θi(27b)σvv2θs=2πk3∫-∞+∞S-ksin⁡θi-qSq-ksin⁡θi×q2sin2θi-k2k2-q2dqThe third- and higher-order perturbation scattering waves are ignored. For backscattering at grazing incidence, θs=-θi=-90∘. In addition, according to the first-order gravitational wave dispersion relation, a more general spectrum Skw,ω including wave frequency can be expressed by the undirectional spectrum Skw[2]:(28)Skw,ω=Skwδω-sgnkwgkw+Skwδω+sgnkwgkwwhere δω is the Dirac impulse function of the variable ω and sgnkw=±1 depends on the ± of variable kw. The two items on the right of (28) represent the spatial/temporal spectrum of the waves moving into the +x space and the waves moving into the -x space, respectively. Putting (28) into (27a) and (27b), the backscattering Doppler spectrum at grazing incidence derived by the SPM is obtained:(29)ΨSPMω=∑m=±116πk3S2mkδω-mgk+2πk3∑m1,m2=±1∫-∞+∞Sm1k1Sm2k2k2-q2δω-m1gk1-m2gk2dqwhere k1=-k-q and k2=q-k. Note that ΨSPMω is calculated from the RCS, so it is the power spectrum.Figure6 gives the Doppler spectrum comparison between the simulation and the SPM in different sea states. We can see that the simulated Doppler spectra are highly consistent with that calculated by the SPM. As we all know, the sea surface backscattering Doppler spectrum is composed of first-order peak, also called Bragg peak, and high-order peaks. In the HF band, the third- and higher-order echo intensity is small enough to be ignored. Note that the simulation uses linear waves, which means the high-order peaks in the Figure 6 are all from the high-order electromagnetic (EM) scattering. In fact, the interaction between the gravity waves will produce nonlinear waves with high-order dispersion relation [30]. If the wavelengths of the produced waves satisfy the Bragg scattering condition, they will have a nonnegligible contribution to the backscattering echo intensity, which will be discussed in detail in the following sections. The high-order peaks from high-order EM scattering and nonlinear sea waves are attributed to the EM contribution and hydrodynamic contribution by [2]. The purpose of comparison with the SPM is to verify the effectiveness of the numerical method in this paper at grazing incidence, so at this point, we do not need to consider the hydrodynamic contribution. From Figure 6 we can also see that the high-order peaks are more sensitive to the wind speed in the low sea state, and when the wind speed exceeds 10 m/s, the electromagnetic contribution to the Doppler spectrum almost does not change. In addition, there is an obvious peak at 2ωb in low sea state. All of these characteristics are consistently reflected in the simulated and calculated Doppler spectra.Figure 6
Comparison of Doppler spectrum between the numerical simulation and the SPM in different sea states. The frequency is 10 MHz and the sea surfaces are linear. The Doppler frequency is normalized by the Bragg frequencyωb, and the spectral intensity is normalized by Bragg peak.
(a)
Wind speed 5 m/s (b)
Wind speed 10 m/s (c)
Wind speed 15 m/s (d)
Wind speed 20 m/sThe consistency of numerical simulation and the SPM in HF band shows that the method in this paper is effective for grazing incidence scattering. Next we will demonstrate the effective frequency range of the SPM by numerical simulation. At the same time, in order to verify the effectiveness of the proposed method in VHF and UHF bands, the classical integral equation method (IEM) is also employed for comparison, which calculates the surface equivalent currents by directly solving the EM integral equation using the moment method, and uses tapered incident wave, expressed by (16), to suppress the discontinuity of surface currents at the edges. Since the tapering parameter q must be greater than 6λ/cos1.5θi to make the tapered wave satisfy Helmholtz equation  [13],the width of the rough surface will increase rapidly with the incident angle, which results in a huge computation amount of classical method at low grazing angle. Here we set the grazing angle of 2∘, the tapering parameter of 1024λ, and the rough surface width of 4096λ for the classical method. The simulation conditions in Table 1 are still applicable for the proposed method. The comparison results of Doppler spectrum of linear waves among the proposed method, the classical integral equation method, and the SPM in different bands are shown in Figure 7. As we can see, the Doppler spectra simulated by the proposed method are almost coincident with that of the classical IEM from HF to UHF bands but with less computational amount, while with the increase of the operating frequency, the intensity of the Doppler spectrum calculated by the SPM is gradually smaller than that of the numerical simulation, which is because, on the one hand, the higher-order perturbation solution to scattering waves cannot be ignored as the frequency increases, and on the other hand, the gradually decreasing radio wavelength does not satisfy the applicable condition of the SPM. For the former, one might be able to improve the accuracy of the SPM by retaining the third- or higher-order perturbation solutions, but this is very troublesome because the form of the higher-order perturbation solutions is too complicated. For the latter, this is restricted by the condition of the SPM and cannot be solved. The mean absolute error of the Doppler spectral intensity over operating frequency between the SPM and the numerical method is shown in Figure 8, which provides a new perspective to measure the applicable frequency of the SPM in the Doppler domain. The Doppler spectral error of the SPM in HF band is very small and can be ignored, but it grows rapidly in VHF band.Figure 7
Comparison of Doppler spectrum of linear waves among the proposed method, the classical integral equation method, and the SPM in different frequency bands. The wind speed is 10 m/s, the Doppler frequency is normalized by the Bragg frequencyωb, and the spectral intensity is normalized by Bragg peak.
(a)
10 MHz (b)
30 MHz (c)
60 MHz (d)
100 MHz (e)
200 MHz (f)
400 MHzFigure 8
The mean absolute error of the Doppler spectral intensity over operating frequency between the SPM and the numerical method. The sea waves are linear and the wind speed is 10 m/s.From Figure7, we can roughly see that the Doppler spectra of linear waves are different at different operating frequencies. Especially at frequencies below 100 MHz, the Doppler spectrum changes more significantly with operating frequency. Figure 9 gives a more detailed description of the numerical simulated Doppler spectrum of linear waves in different frequencies. As the frequency increases, the Doppler frequency and intensity of the higher-order peaks gradually approach the Bragg peak, and the approaching rate in the high waves is much greater than that in the low waves. The approaching rate gradually slows down as the frequency increases, and the higher-order peaks are almost overlapped with the Bragg peaks when the frequency exceeds 200 MHz. There is no significant difference in the shape of the Doppler spectrum in UHF band, which means that the EM contribution tends to stabilize with the increase of frequency.Figure 9
Numerical simulated Doppler spectrum of linear waves over operating frequency. The horizontal axis is the Doppler angular frequency normalized by the Bragg frequency, and the vertical axis represents different simulation operating frequencies from 5 MHz to 400 MHz.
(a)
Wind speed 3 m/s (b)
Wind speed 7 m/s


## 4.3. Doppler Spectrum of Linear and Nonlinear Waves

The actual sea surface is nonlinear, and the nonlinear components generated by the interaction of gravity waves will introduce hydrodynamic contribution in Doppler spectrum. Figure10 gives the comparison of Doppler spectrum among linear waves, Creamer waves, and Lagrange waves in different frequencies. We can see that, in all bands from HF to UHF, the high-order peak intensity of nonlinear waves is enhanced relative to that of linear waves, but the hydrodynamic contribution does not affect the Bragg peak. In general, the high-order peaks of the Lagrange waves and the Creamer waves are very close in the Doppler frequency range of [-1.5ωb,1.5ωb], and the high-order peak intensity of the Lagrange waves is gradually higher than that of the Creamer waves with the increase of the Doppler frequency, which implies that the nonlinear characteristics of the Lagrange waves are stronger than the Creamer waves. Consistent with the linear wave, after entering the VHF band, the shape of the Doppler spectrum of the nonlinear wave hardly changes with the increase of the operating frequency, except that the high-order peak intensity of the Creamer waves gradually approaches that of the Lagrange waves, which on the one hand shows that both EM and hydrodynamic contributions are stable in the VHF and UHF bands, and on the other hand, the frequency of the water wave interacting with the electromagnetic wave gradually increases with the increase of the operating frequency, and at this point the nonlinear characteristics of the Lagrange waves and the Creamer waves will gradually approach.Figure 10
Doppler spectrum comparison between linear and nonlinear waves in HF/VHF/UHF bands. The wind speed is 10 m/s, and the Doppler frequency is normalized by the Bragg frequencyωb.
(a)
10 MHz (b)
60 MHz (c)
200 MHz (d)
600 MHzIn the issue of remote sensing of sea state using radar, one is more concerned about the characteristics of echo Doppler spectrum in different sea states. Figure11 shows the variation of the Doppler spectrum of nonlinear waves with the sea state in different operating frequencies. We can see that, in the HF band, high-order peaks are very sensitive to the sea state, especially in the low- and medium-sea state, and the Bragg peak intensity, however, hardly changes with the sea state, which are consistent with Barrick’s sea surface first- and second-order scattering theory [1, 2]. Typically, HF radar uses the intensity changes of higher-order peaks relative to the Bragg peak to invert the wave height [5]. When the wind speed exceeds 15 m/s, the Doppler spectrum of the Creamer wave is almost no longer changed, which is explained by the fact that the sea waves interacting with the electromagnetic waves reach saturation. It is obvious that the Lagrange waves reach saturation more slowly. As the operating frequency increases, the wavelength of the sea waves interacting with the electromagnetic waves becomes smaller, and thus more easily saturated. As shown in Figures 11(c) and 11(f), the Doppler spectra in the UHF band have reached stability at very low wind speeds, which is actually not conducive to inversing the wave height through the Doppler spectrum characteristics. Another obvious feature is that the frequencies of the high-order peaks are very close to that of Bragg peak in the VHF and UHF bands. Considering the frequency shifts caused by surface currents in the actual ocean, the higher-order peaks will overlap with Bragg peak, and therefore the method of extracting wave height in the HF band will be ineffective in the VHF and UHF bands.Figure 11
Doppler spectrum of nonlinear waves in different sea states and different operating frequencies. The horizontal axis represents the Doppler frequency normalized by the Bragg frequencyωb, the vertical axis is the wind speed at 10 meters above sea level, and the color depth represents Doppler spectrum intensity.
(a)
Creamer waves at 10 MHz (b)
Creamer waves at 60 MHz (c)
Creamer waves at 600 MHz (d)
Lagrange waves at 10 MHz (e)
Lagrange waves at 60 MHz (f)
Lagrange waves at 600 MHzA more intuitive comparison of the Doppler spectra in the HF and UHF bands at different wind speeds is shown in Figure12. It is obvious that the shape of the Doppler spectrum in the UHF band changes very slightly with the wind speed relative to that in the HF band. The high-order peaks in the UHF band are continuously widened as the wind speed increases, and the intensity is also slightly enhanced, which may provide an idea for remote sensing of the wave height using UHF radar.Figure 12
Doppler spectrum of Creamer waves in different sea states and different operating frequencies. The Doppler frequency is normalized by the Bragg frequencyωb.
(a)
10 MHz (b)
600 MHz


## 5. Conclusions

In this paper, we use the improved integral equation method to study the characteristics of Doppler spectra backscattered from linear and nonlinear sea waves in HF/VHF/UHF bands at grazing incidence. In order to eliminate the edge effect and suppress the side lobes of scattered waves, the time-varying rough sea surface is transformed into a local perturbation plane by multiplying a Tukey-Hanning window function, which flattened asymptotically at the edges. Unlike the classical integral equation method, we divide the surface currents into a known coherent current and unknown incoherent currents. The coherent current excites the reflected wave, while the incoherent currents excite the scattered waves. The incoherent currents are considered as the surface unknowns and calculated by solving the surface integral equation established on the unbounded local perturbation plane using the moment method. The incident plane wave ensures the accuracy of the solution of the integral equation at grazing incidence. The calculation result shows that the intensity of the incoherent currents decreases with the decrease of roughness of the local perturbation plane and is equal to zero at the edge. The comparisons with the SPM on the one hand verify the effectiveness of the method in this paper at grazing incidence, and on the other hand, it also measures the applicable frequency range of the SPM in the Doppler spectrum domain. In the HF band, the Doppler spectrum of the numerical simulation is completely consistent with that derived by the SPM whether in the low sea state or the high sea state. The difference between the Doppler spectra derived from SPM and simulated by numerical method increases rapidly with the increase of operating frequency, which is on the one hand because the slope and wave height of the sea surface gradually fail to meet the conditions of the SPM, and on the other hand, the third- and higher-order perturbation solution to the scattered waves cannot be ignored.The numerically simulated Doppler spectrum is robust in VHF and UHF bands, which facilitates further analysis of the characteristics of the Doppler spectrum. We compared the Doppler spectra of linear and nonlinear sea waves at different operating frequencies and different sea states. The nonlinear waves include Creamer waves with only crest-trough asymmetry and Lagrange waves with both crest-trough and front-back asymmetry. In general, the high-order peak intensity of the nonlinear waves is about 10 dB higher than that of the linear waves, and the intensity of the Bragg peak is not affected by the nonlinearity of the sea surface. In the Doppler frequency range of[-1.5ωb,1.5ωb], the Doppler spectrum of the Lagrange waves is very close to that of the Creamer waves, and outside this range, the higher-order peaks of the Lagrange wave are slightly stronger than the Creamer waves, which indicates that the nonlinearity of the Lagrange waves is more abundant. Furthermore, with the increase of simulation frequency, the Doppler spectrum becomes less and less sensitive to the change of the sea state, and the high-order peaks gradually overlap with the Bragg peak, which reveals that the mature sea state inversion methods applied in HF radar are no longer effective in VHF and UHF bands. Numerical simulations show that the higher-order peaks are continuously broadened as the wind speed increases, which may provide an access to invert the wave height in the VHF and UHF bands. We will extend the numerical simulation to the 2D sea surface in subsequent study to research the relationship between the Doppler spectrum and the sea state.



---

*Source: 5363264-2019-07-15.xml*