Introduction of vertical nonlinear features

The linear surface model is the prerequisite for nonlinear surface models. Common methods for modeling sea surfaces include fractal, bilinear superposition, and linear filtering [Reference Rino, Crystal, Koide, Ngo and Guthart15]. The linear filtering method treats the surface as a superposition of different harmonics, with the amplitude of the harmonics being independent Gaussian random variables that are specifically proportional to the wave spectrum. The linear sea surface can be realized by Fourier transform of the spectrum. The spatial Fourier components of the linear surface can be expressed as follows:

(1)\begin{align}{A_L}\left( {{k_x},{k_y}} \right) & = Y \left( {{k}} \right)\sqrt {W\left( {{k_x},{k_y}} \right)\Delta {k_x}\Delta {k_y}/2} \nonumber \\ & \qquad + {Y ^*}\left( { - {{\bf {k}}}} \right)\sqrt {W\left( {{k_x},{k_y}} \right)\Delta {k_x}\Delta {k_y}/2} \end{align}

(2)\begin{align}{h_L}\left( {x,y} \right) & = 4M \cdot N \cdot {F^{ - 1}}\left[ {{A_L}\left( {{k_x},{k_y}} \right)} \right] \nonumber \\ & \qquad = \sum\limits_{m = - M}^M {\sum\limits_{n = - N}^N {{A_L}\left( {{k_{xm}},{k_{yn}}} \right)} } \exp [\,j({k_{xm}}x + {k_{yn}}y)]\end{align}

where k represents the wave vector of the sea surface; ${k_x},{k_y}$ represent the wave vectors along the x-, y-direction, respectively; Y(k) is the independent Gaussian random variables; the * denotes the Hermitian form, which insures that ${h_L}(x,y)$is real; and ${h_L}(x,y)$ is the sea height at different points. The number of grid points along the surface in x- and y-directions are defined as 2 M and 2 N, respectively. ${F^{ - 1}}[ \cdot ]$ represents the inverse two-dimensional Fourier transform. $W\left( {{k_x},{k_y}} \right)$ represents the two-dimensional wave spectrum, composed of an omnidirectional and a directional spectrum. The wave spectrum describes the distribution of energy on the ocean surface with respect to wavelength and wave direction.

The sea surface height and velocity potential can be formulated as a pair of typical Hamiltonian variables. The nonlinear term of the second-order Creamer model can be written as the Hilbert transform of the linear surface [Reference Creamer, Henyey, Schult and Wright16]. In the two-dimensional case, the Hilbert transform vector of the linear term is

(3)\begin{equation}{{\bf \it {h}}_c}(x,y) = {\mathop{\textrm Re}\nolimits} \sum\limits_{\bf \it {k}} {( - j{{\bf \it {k}} \over k})} {A_L}({k_x},{k_y})\exp (\,j({k_x}x + {k_y}y))\end{equation}

where $k = \sqrt {k_x^2 + k_y^2} $, and the components of its Hilbert vector in the x-, y-directions are as follows:

(4)\begin{equation}{h_{cx}}(x,y) = {\mathop{\textrm {Re}}\nolimits} \sum\limits_{\bf \it {k}} {( - j{{{k_x}} \over k})} {A_L}({k_x},{k_y})\exp (\,j({k_x}x + {k_y}y))\end{equation}

(5)\begin{equation}{h_{cy}}(x,y) = {\mathop{\textrm {Re}}\nolimits} \sum\limits_{\bf \it {k}} {( - j{{{k_y}} \over k})} {A_L}({k_x},{k_y})\exp (\,j({k_x}x + {k_y}y))\end{equation}

The Creamer nonlinear Fourier components can be calculated by Equation (6) [Reference Fung, Li and Chen12].

(6)\begin{equation}{A_{NL}}({k_x},{k_y}) = {1 \over {MN}}\sum\limits_{x,y} {{{\exp \left[ {\,j{\bf \it {k}} \cdot {{\bf \it {h}}_c}(x,y)} \right] - 1} \over k}} \exp (\,j({k_x}x + {k_y}y))\end{equation}

where the exponential term in Equation (6) is dependent not only on k but also on the parameter ${{\boldsymbol{h}}_c}(x,y)$, which leads to complex model calculations. Soriano et al. retained the first two approximations by expanding the exponential term using Taylor series, which effectively reduces the computational complexity while preserving the nonlinear interaction between waves [Reference Soriano, Joelson and Saillard17]. The first term of the Taylor expansion is equivalent to the Fourier transform of the linear ocean surface, designated as ${A_L}\left( {{k_x},{k_y}} \right)$. The second term characterizes weakly nonlinear effects:

(7)\begin{equation}{C_2}{\textrm{(}}{k_x},{k_y}{\textrm{) = }} - {{k_x^2} \over {2k}}F\left[ {h_{cx}^2} \right] - {{{k_x}{k_y}} \over k}F\left[ {{h_{cx}}{h_{cy}}} \right] - {{k_y^2} \over {2k}}F\left[ {h_{cy}^2} \right]\end{equation}

where $F[ \cdot ]$ represents the Fourier transform, and ${C_2}\left( {{k_x},{k_y}} \right)$ defines the second term of Creamer model expansion. Then the Fourier transform of the nonlinear surface can be expressed as follows:

(8)\begin{equation}{A_{NL}}({k_x},{k_y}) = {A_L}({k_x},{k_y}) + {C_2}{\textrm{(}}{k_x},{k_y}{\textrm{)}}\end{equation}

The nonlinear sea surface of the second-order Creamer model can be obtained by replacing the ${A_L}\left( {{k_x},{k_y}} \right)$ of Equation (2) with the ${A_{NL}}({k_x},{k_y})$ and performing a Fourier inverse transform.

The shape of the sea surface will be different under different wind conditions. The wind serves as a significant environmental parameter in simulating the sea surface. It can be described by the wind speed at 10 m above sea level and the wind direction ($\phi $). Figure 1 shows the two-dimensional linear and nonlinear sea surfaces established using the Elfouhaily spectrum with a wind speed of 5 m/s and wind direction along the x-direction [Reference Elfouhaily, Chapron, Katsaros and Vandemark21].

Figure 1. Two-dimensional sea surface of (a) linear (b) nonlinear.

To facilitate a more detailed examination of the disparity between linear and nonlinear surfaces, a one-dimensional surface can be derived by sampling from a two-dimensional surface. In this study, Fig. 2(a) displays the one-dimensional linear and nonlinear surfaces obtained by sampling from Fig. 1 at y = 250 m along the x-direction. Additionally, Fig. 2(b) illustrates the cumulative distribution function of the slopes computed from the linear and nonlinear surfaces depicted in Fig. 2(a).

Figure 2. One-dimensional sea surface: (a) Height distribution; (b) Slope distribution.

Figure 2(a) shows that the nonlinear sea surface is sharper, with steeper crests and gentler troughs than the linear sea surface. The sea surface is commonly regarded as a superposition of waves with different wavelengths, which can be classified into shorter wavelength short waves and longer wavelength long waves. Wavelengths of short waves are usually less than 1 m, which is close to radar wavelengths and has a strong correlation with remote-sensing imaging results. On the other hand, wavelengths of long waves range from a few meters to several hundred meters. Due to the interaction between the waves in the sea surface, short waves are suppressed at the troughs of long waves and enhanced at the crests. By analyzing the cumulative distribution characteristics of the slope for linear and nonlinear sea surfaces in Fig. 2(b), it can be observed that the slope distribution of the nonlinear sea surface is steeper compared to the linear surface.

Introduction of horizontal nonlinear features

The sea surface is a randomly rough surface affected by wind and waves, with different height at different positions. The sea surface height $h(x,y)$ is a random variable that can be described by the probability density function of the sea surface height. A random variable can be characterized by its statistical distribution by using moments. Let the surface function $h(x,y)$ be a real stationary random process with zero mean. The first-order moments of the random variable represent its mathematical expectation $\left\langle h \right\rangle = 0$, the second-order moments represent the variance $\left\langle {{h^2}} \right\rangle = {\delta ^2}$, and the third-order moments represent the skewness $\left\langle {{h^3}} \right\rangle = {\mu _3}$ [Reference Su, Zhang, Dang and Tan18]. The second-order moments of the sea surface height are related to the autocorrelation function of the sea surface.

(9)\begin{equation}\left\langle {h\left( {x,y} \right)h\left( {x + {\xi _x},y + {\xi _y}} \right)} \right\rangle = {\delta ^2}\rho \left( {{\xi _x},{\xi _y}} \right)\end{equation}

(10)\begin{equation}W\left( {{k_x},{k_y}} \right) = {1 \over {2\pi }}\int {{\delta ^2}\rho \left( {{\xi _x},{\xi _y}} \right)\exp \left( { - j{k_x}{\xi _x} - j{k_y}{\xi _y}} \right)} d{\xi _x}d{\xi _y}\end{equation}

where $\delta $ is the root mean square height of the surface, ${\xi _x}$ represents the displacement in the x-direction, and ${\xi _y}$ represents the displacement in the y-direction, $\rho ({\xi _x},{\xi _y})$ is the autocorrelation function of surface.

Equation (10) shows that the wave spectrum $W({k_x},{k_y})$ is the Fourier transform of the autocorrelation function $\rho ({\xi _x},{\xi _y})$. The autocorrelation function is the second-order statistic that characterizes the statistical distribution of sea surface height. The third-order moments of the sea surface height are related to the bicorrelation function of the sea surface. The bispectrum is defined as the Fourier transform of the bicorrelation function.

(11)\begin{equation}\left\langle {h\left( {x,y} \right)h\left( {x + {\xi _x},y + {\xi _y}} \right)h(x + {\zeta _x},y + {\zeta _y})} \right\rangle = {\delta ^3}S\left( {{\xi _x},{\xi _y};{\zeta _x},{\zeta _y}} \right)\end{equation}

where ${\xi _x},{\zeta _x}$ denotes the displacement in the x-direction and ${\xi _y},{\zeta _y}$ denotes the displacement in the y-direction, and $S\left( {{\xi _x},{\xi _y};{\zeta _x},{\zeta _y}} \right)$ is the bicorrelation function, which can characterize both the linear and nonlinear features of the sea surface [Reference Chen, Fung and Amar13].

when ${\zeta _x} = {\xi _x}$, ${\zeta _y} = {\xi _y}$, we have

(12)\begin{equation}\left\langle {h\left( {x,y} \right){h^2}\left( {x + {\xi _x},y + {\xi _y}} \right)} \right\rangle = {\delta ^3}S\left( {{\xi _x},{\xi _y}} \right)\end{equation}

and when ${\zeta _x} = 0$, ${\zeta _y} = 0$, we have

(13)\begin{align}\left\langle {{h^2}\left( {x,y} \right)h\left( {x + {\xi _x},y + {\xi _y}} \right)} \right\rangle & = \left\langle {{h^2}\left( {x' - {\xi _x},y' - {\xi _y}} \right)h\left( {x',y'} \right)} \right\rangle \nonumber \\ & \quad = {\delta ^3}S\left( { - {\xi _x}, - {\xi _y}} \right)\end{align}

The variables have reduced from four to two. The bicorrelation function $S\left( {{\xi _x},{\xi _y}} \right)$ can be decomposed into a symmetric part and an asymmetric part.

(14)\begin{equation}S\left( {{\xi _x},{\xi _y}} \right) = {S_s}\left( {{\xi _x},{\xi _y}} \right) + {S_a}\left( {{\xi _x},{\xi _y}} \right)\end{equation}

(15)\begin{equation}{S_s}\left( {{\xi _x},{\xi _y}} \right) = {{S\left( {{\xi _x},{\xi _y}} \right) + S\left( { - {\xi _x}, - {\xi _y}} \right)} \over 2}\end{equation}

(16)\begin{equation}{S_a}\left( {{\xi _x},{\xi _y}} \right) = {{S\left( {{\xi _x},{\xi _y}} \right) - S\left( { - {\xi _x}, - {\xi _y}} \right)} \over 2}\end{equation}

According to the properties of Fourier transforms, the Fourier transform of an even function is a real-valued even function, and the Fourier transform of an odd function is a purely imaginary odd function. The bispectrum $B({k_x},{k_y})$ can be divided into real and imaginary parts.

(17)\begin{align} B\left( {{k_x},{k_y}} \right) & = {B_s}\left( {{k_x},{k_y}} \right) + j{B_a}\left( {{k_x},{k_y}} \right) \nonumber \\ & = {1 \over {{{\left( {2\pi } \right)}^2}}}\int {{\delta ^3}} S({\xi _x},{\xi _y})\exp \left( { - j{k_x}{\xi _x} - j{k_y}{\xi _y}} \right)d{\xi _x}d{\xi _y} \end{align}

where ${B_{\textrm{s}}}({k_x},{k_y})$ is the real part of the bispectrum, which is the Fourier transform of ${S_s}\left( {{\xi _x},{\xi _y}} \right)$, and ${B_{\textrm{a}}}({k_x},{k_y})$ is the imaginary part of the bispectrum, which is the Fourier transform of ${S_a}\left( {{\xi _x},{\xi _y}} \right)$.

In the IEM model, ${B_{\textrm{s}}}({k_x},{k_y})$ does not contribute to backscattering, while ${B_{\textrm{a}}}({k_x},{k_y})$ can be derived from Equation (18) [Reference Fung, Li and Chen12–Reference Chen, Wu, Tsang, Li, Shi and Fung14].

(18)\begin{equation}{B_a}(k,\phi ) = - {1 \over {16}}k{s_0}\left( {6 - {k^2}s_0^2{{\,\cos }^2}\phi } \right)\cos \phi \exp \left( { - {{k{s_0}} \over 4}} \right)\end{equation}

where ${s_0}$ is the correlation distance of the bispectrum function [Reference Creamer, Henyey, Schult and Wright16].

Establishment of nonlinear electromagnetic scattering model

In the electromagnetic scattering model of the random rough surface, radar waves are incident to the sea surface and received by the radar after undergoing scattering from the sea surface. In this paper, the IEM-C model is developed based on the IEM model, combining the Creamer nonlinear sea-surface and bispectrum, to calculate the backscattering coefficients. The backscattering coefficients of the IEM-C model can be divided into the contributions of the wave spectrum and bispectrum, as shown in Equation (19).

(19)\begin{equation}\sigma _{pp}^0 = \sigma _{pp}^0(N) + \sigma _{pp}^0(S)\end{equation}

The $\sigma _{pp}^0(N)$ represents the contribution from the wave spectrum, while the $\sigma _{pp}^0(S)$ represents the contribution from the bispectrum [Reference Chen, Fung and Amar13]. The backscattering coefficient generated by the wave spectrum$\sigma _{pp}^0(N)$ is given by Equation (20).

(20)\begin{equation}{\sigma ^0}_{pp}(N) = {{k_i^2} \over 2}\exp \left[ { - 2{\delta ^2}{k_i}^2{{\cos }^2}{\theta _i}} \right]\sum\limits_{n = 1}^\infty {{{\left| {I_{pp}^n} \right|}^2}{{{W^{\left( n \right)}}\left( {2{k_i}\sin {\theta _i},0} \right)} \over {n!}}} \end{equation}

(21)\begin{equation}I_{pp}^n = {\left( {2{k_i}\cos {\theta _i}\delta } \right)^n}{f_{pp}}\exp ( - k_i^2{\delta ^2}{\cos ^2}{\theta _i}) + {({k_i}\delta \cos {\theta _i})^n}{F_{pp}}\end{equation}

where ${k_i}$ is the wave number of incident radar, ${\theta _i}$ is the incidence angle of radar wave, ${f_{pp}}$ is Kirchhoff field coefficient, and ${F_{pp}}$ is the compensation field coefficient [Reference Chen, Fung and Weissman10]. The subscript pp represents the polarization state of radar waves. $W_{}^{(n)}(2{k_i}\sin {\theta _i},0)$ is the nth-order wave spectrum which is the Fourier transform of the nth autocorrelation function ${\rho ^{(n)}}(x,y)$.

$\sigma _{pp}^0(S)$ consists of three components representing the contribution of the bispectrum to the Kirchhoff field, the cross-field, and in the compensating field, as shown in Equation (22). The detailed derivation process can be found in paper [Reference Fung, Li and Chen12], specifically equations (31)–(53).

(22)\begin{align}& \sigma _{pp}^0(S) = {{k_i^2} \over 2}{\left| {{\,f_{pp}}} \right|^2}\exp \left( { - 4k_i^2{\delta ^2}{{\,\cos }^2}{\,\theta _i}} \right) \nonumber \\ & \qquad \qquad \quad \sum\limits_{n = 1}^\infty {{{{{\left( { - 8k_i^3{{\,\cos }^3}{\,\theta _i}} \right)}^n}} \over {n!}}} B_a^{(n)}(2{k_i}\sin {\theta _i},0) \nonumber \\ & \qquad \qquad + {{k_i^2} \over 2}{\mathop{\textrm Re}\nolimits} \left( {{f_{pp}}^*{F_{pp}}} \right)\exp \left( { - 3k_i^2{\delta ^2}{{\cos }^2}{\theta _i}} \right)\nonumber \\ & \qquad \qquad \quad \sum\limits_{n = 1}^\infty {{{{{\left( { - 3k_i^3{{\,\cos }^3}{\,\theta _i}} \right)}^n}} \over {n!}}B_a^{(n)}(2{k_i}\sin {\theta _i},0)} \nonumber \\ & \qquad\qquad + {{k_i^2} \over 8}{\left| {{F_{pp}}} \right|^2}\exp \left( { - 2k_i^2{\delta ^2}{{\,\cos }^2}{\,\theta _i}} \right)\nonumber \\ & \qquad \qquad \quad \sum\limits_{n = 1}^\infty {{{{{\left( { - k_i^3{{\,\cos }^3}{\,\theta _i}} \right)}^n}} \over {n!}}} B_a^{(n)}(2{k_i}\sin {\theta _i},0)\end{align}

where $B_a^{(n)}(2{k_i}\sin {\theta _i},0)$ represents the nth-order bispectrum, which is the Fourier transform of the nth-order bicorrelation function.

Different wind angles

First, the variation of the backscattering coefficient of the IEM-C model with the wind direction is analyzed. The wind direction angle ϕ is defined as the angle between the wind direction and the radar line of sight direction. With a fixed radar wave frequency of 5.3 GHz, an incident angle of 30°, the wind speed of 4 m/s, and the wind direction of 0° (north), the relationship between the backscattering coefficient of the sea surface and the wind direction angle is illustrated in Fig. 3.

Figure 3. Backscattering coefficients for different wind directions: (a) HH polarization; (b) VV polarization.

From Fig. 3, it can be observed that the backscattering coefficients of HH polarization and VV polarization exhibit cosine-like variations with respect to ϕ. The backscattering coefficients of the IEM model are the same in upwind (ϕ = 0°) and downwind (ϕ = 180°) conditions. The backscattering coefficients of the IEM-C model exhibit differences in upwind (ϕ = 0°) and downwind (ϕ = 180°) conditions, with a 0.6 dB increase in upwind compared to downwind. The IEM-C model demonstrates differences in the backscattering coefficients between upwind and downwind conditions, which are not present in the IEM model. For 0 < ϕ < 90°, the backscattering coefficient of the IEM-C model is greater than that of the IEM model. At ϕ = 90°, the backscattering coefficient of the IEM-C model is the same as the IEM model. For 90 < ϕ < 270°, the backscattering coefficient of the IEM-C model is smaller than that of the IEM model. The IEM-C model adds the bispectrum of the surface on the basis of the second-order Creamer nonlinear sea surface, considering the nonlinear skewness characteristics in the horizontal direction caused by wind direction. In upwind cases, the roughness increases, leading to an enhancement the backscattering coefficient. In downwind cases, the roughness decreases, reducing the backscattering coefficient. The model simulation results show that the IEM-C model exhibits wind-dependent characteristics that are closer to the real sea surface scattering conditions.

Based on Fig. 3(a), it is evident that the disparity between the backscattering coefficients of the IEM-C model and the IEM model decreases as ϕ increases from 0° to 90°. Conversely, the disparity between the backscattering coefficients of the IEM-C model and the IEM model gradually increases as ϕ increases from 90° to 180°.

The closer the wind direction, the stronger the nonlinear scattering due to the wind direction. The nonlinear scattering properties are most pronounced in the case of upwind and downwind, while they have no effect in the case of crosswind. The variation of nonlinear scattering contributions observed in VV polarization shows the same trend as in HH polarization. Comparing Fig. 3(a), (b), the backscattering coefficient of HH polarization (0.31 dB) is larger than that of VV polarization (0.15 dB) for the same wind angle. In the IEM-C model, the difference of backscattering coefficients due to wind direction is more pronounced in HH polarization.

Different radar wave incidence angles

Under the premise of wind speed at 5 m/s and radar frequency at 5.3 GHz, the backscattering coefficient of the IEM-C model and the IEM model was calculated for upwind and downwind conditions with different incidence angles. The results are shown in Fig. 4.

Figure 4. Backscattering coefficients at different incidence angles (a) HH – upwind; (b) VV – upwind; (c) HH – downwind; (d) VV – downwind.

Figure 4 illustrates a gradual decrease in the backscattering coefficients of both the IEM model and the IEM-C model as the incident angle increases. At small incidence angles, the dominant mechanism responsible for the backscattering is specular reflection, contributing to a higher backscatter coefficient. As the incidence angle increases, the contribution of specular reflection decreases, and Bragg scattering becomes dominant. As the angle of incidence increases further, the Bragg scattering decreases.

Based on the data presented in Fig. 4(a), (b), it can be observed that the backscattering coefficient of the IEM-C model is higher than that of the IEM model under upwind conditions. Similarly, from the observations in Fig. 4(c), (d), it is evident that the backscattering coefficient of the IEM-C model is lower compared to the IEM model under downwind conditions. As the incident angle increases, the difference of backscattering coefficients between the IEM-C model and the IEM model initially increases and then decreases, with the most significant difference (2.3 dB) observed at an incident angle of 30°.

Different wind speeds

Figure 5 illustrates the relationship between the wind speed and backscattering coefficient (upwind and downwind) for the radar wave frequency of 5.3 GHz and incident angle of 30°.

Figure 5. Backscattering coefficients at different wind speeds (a) HH – upwind; (b) VV – upwind; (c) HH – downwind; (d) VV – downwind.

A comprehensive analysis of Fig. 5 reveals that, irrespective of the polarization mode or the specific model employed, there is a discernible upward trend in the backscattering coefficient of the sea surface with increasing wind speed. This observed behavior can be attributed to the heightened instability of the sea surface as wind speed escalates. This leads to the increase in the roughness of the sea surface and the consequent increase in the backscattering coefficient.

By analyzing the backscattering results of the IEM and IEM-C models under the same polarization and wind direction conditions depicted in Fig. 5, it becomes evident that the IEM model exhibits a linear increase in the backscattering coefficient with increasing wind speed. Conversely, the backscattering coefficient of the IEM-C model demonstrates a nonlinear variation with wind speed. As the wind speed intensifies, the disparity between the IEM-C and IEM models initially enlarges and then diminishes, reaching its most pronounced level at a wind speed of 5 m/s.

Observing the backscattering coefficients for different polarizations and wind directions in Fig. 5(a), (c) and Fig. 5(b), (d), significant differences are observed. The backscattering coefficients of the IEM-C model for the downwind case are smaller than those of the IEM model for the same wind speeds. The scattering coefficients of the IEM model are smaller than those of the IEM-C model for the upwind case. Under the same wind speed conditions, the difference in backscattering coefficients between the IEM-C model and the IEM model is larger in the HH polarization compared to the VV polarization.