Measuring the Directional Ocean Spectrum from Simulated Bistatic HF Radar Data

Abstract

HF radars are becoming important components of coastal operational monitoring systems particularly for currents and mostly using monostatic radar systems where the transmit and receive antennas are colocated. A bistatic configuration, where the transmit antenna is separated from the receive antennas, offers some advantages and has been used for current measurement. Currents are measured using the Doppler shift from ocean waves which are Bragg-matched to the radio signal. Obtaining a wave measurement is more complicated. In this paper, the theoretical basis for bistatic wave measurement with a phased-array HF radar is reviewed and clarified. Simulations of monostatic and bistatic radar data have been made using wave models and wave spectral data. The Seaview monostatic inversion method for waves, currents and winds has been modified to allow for a bistatic configuration and has been applied to the simulated data for two receive sites. Comparisons of current and wave parameters and of wave spectra are presented. The results are encouraging, although the monostatic results are more accurate. Large bistatic angles seem to reduce the accuracy of the derived oceanographic measurements, although directional spectra match well over most of the frequency range.

1. Introduction

Coastal high frequency or HF radar (3–30 MHz) is a tool that has enabled users to remotely measure ocean currents, winds and waves, in real-time, since the initial observation of Crombie [1] in 1955, who realised the relationship between ocean wave and HF radar Doppler spectra (such as that shown in Figure 1). The measurements are important in a number of coastal engineering topics, including testing of and assimilation in operational wave models, sea vessel navigation, land/beach erosion, designing offshore structures, and in supporting marine activities. Additionally, collecting ocean data over a long period of time can be useful in climate changes studies; the same data can also be used to assess the potential of a coastal region to become a wave/wind farm as shown by Wyatt [2]. With such important applications, the accuracy of the measurements is imperative.

The majority of the existing theory is for monostatic radar, where the transmitter and receiver are co-located. Ocean surface current measurements from such a radar are robust, and hundreds of radars provide real-time current measurements all over the world; see the work of Paduan and Washburn [3] for an introduction to the subject. Methods for measuring wind direction and variability are also reliable, such as the method of Wyatt et al. [4], where the maximum likelihood method is used fit the wind distribution to the available radar data. Ocean wave measurements are also possible, however, they can be less robust as they are more vulnerable to noise (for limitations see Wyatt et al. [5]). Notable methods of Lipa [6,7], Howell and Walsh [8], and Hisaki [9] are all significant in the history of measuring the wave spectrum. Another key method is the Seaview method, presented by Wyatt [10], and Green and Wyatt [11], which will be explained in Section 3.

To obtain the ocean measurements using radar data, many researchers use the radar cross section of the ocean surface, which models what the radar output will be for a given ocean state and radar frequency. For monostatic radars, the most commonly used radar cross section is that of Barrick [12,13], derived in 1972, based on the perturbation method of Rice [14]. The expression is split into its first and second order components; the first order radar cross section, ${\sigma }^{(1)}(\omega )$, is due to resonance between the emitted radio waves and ocean waves of a particular length and direction, and the second order radar cross section, ${\sigma }^{(2)}(\omega )$, is due to double scattering of the emitted radio waves from two ocean waves and the non linear combination of the same two ocean waves. In full, for radar wavenumber $k_0$ and ocean spectrum $S(\overrightarrow{k})$,

$$\begin{array}{c}\hfill {\sigma }^{(1)}(\omega )=2^6\pi k_0^4\sum _{m=\pm 1}S(m\overrightarrow{k_B})\delta (\omega -m{\omega }_B),\end{array}$$

(1)

where $k_B=2k_0$, is known as the Bragg wavenumber and

$$\begin{array}{c}\hfill {\omega }_B=\sqrt{2g{k}_{0}tanh(2k_{0}d)}\end{array}$$

(2)

is known as the Bragg frequency, for ocean depth d and gravity g. The second order term is given by

$$\begin{array}{cc}\hfill {\sigma }^{(2)}(\omega )=2^6\pi k_0^4& \sum _{m,m^{\prime }=\pm 1}{\int \int }_{-\infty }^{\infty }{\left|{\mathrm{\Gamma }}_T\right|}^{2}S(\overrightarrow{m{k}_1})S(m^{\prime }\overrightarrow{k_2})\delta (\omega -m{\omega }_1-m^{\prime }{\omega }_2)dpdq,\hfill \end{array}$$

(3)

where $\overrightarrow{k_1}$ and $\overrightarrow{k_2}$, with respective angular frequencies ${\omega }_1$ and ${\omega }_2$, are the two contributing wave vectors, defined by the relationship $\overrightarrow{k_1}+\overrightarrow{k_2}=\overrightarrow{k_B}$. The $|{\mathrm{\Gamma }}_T{|}^2$ term is known as the coupling coefficient and contains the mathematics of the nonlinear combinations of the waves and, as such, is a function of $\overrightarrow{k_1}$ and $\overrightarrow{k_2}$. More detail on the monostatic coupling coefficient is given by Lipa and Barrick [15].

Numerical methods like that of Holden and Wyatt [16] can be used to simulate monostatic Doppler spectra for given ocean conditions and radar settings, using Equations (24) and (28); a comparison of a radar measured Doppler spectrum and a simulated Doppler spectrum, using co-spatial and co-temporal wave buoy data as input to the simulation, is shown in Figure 2 where good agreement is shown.

Recently, bistatic radar—where the transmitter and receiver are separated by a notable distance—is on the ascendancy. Therefore, conversely to monostatic radar which receives backscatter, the detected radio waves in a bistatic radar have been scattered at a non-zero angle. A traditional coastal HF radar site consists of two monostatic radars, each providing data from ocean backscatter. However, a third dataset can be obtained at no additional cost if one of the receivers also receives bistatic scatter from the other transmitter. In this case, the radar site is called multistatic. Each different radar setup is shown in Figure 3. The advantages of employing a bistatic/multistatic radar setup are that (1) it can reduce the cost of setting up/maintaining a HF radar and (2) it can increase spatial coverage and data quality as shown by Whelan Hubbard [17].

In this paper, the aim is to obtain directional wave measurements from bistatic radar data. Previously, for bistatic radar, Zhang and Gill [18] developed an inversion algorithm to obtain the nondirectional wave spectrum from bistatic radar data and, when tested on simulated data, they obtained good results. Silva [19] has also presented results of simulated bistatic HF radar data inversion where the directional wave spectrum was estimated using Tikhonov regularization. They achieved good results for simulated data, however, the method is limited as a model is assumed for the direction of the spectrum, and this assumption may not always be appropriate.

In the existing numerical methods for extracting wave measurements from monostatic radar data, the aim is to invert Equation (3), to obtain the ocean spectrum $S(\overrightarrow{k})$ in terms of the measured $\sigma (\omega )$. Anderson [20] stated that the existing algorithms should work for a bistatic system if the inverted radar cross section is changed to the bistatic expression. In this work, we test this hypothesis and modify the Seaview method to measure the directional wave spectrum from bistatic HF radar data. Therefore, to do this, the bistatic radar cross section must be known.

In 1975, Johnstone [21] presented a bistatic radar cross section of the ocean surface and then, in 2001, Gill and Walsh [22] presented an alternative expression. However, under monostatic conditions, neither of the expressions reduce exactly to the monostatic term of Barrick [12,13] (which the Seaview method depends on). The derivation of Johnstone appears to have an error which causes the difference in the resulting expressions; Gill and Walsh followed a more complicated method, however, it has been shown that the monostatic form of their radar cross section is similar to Barrick’s and it is, therefore, unnecessary to change the existing operational inversion programs to use theirs instead. Another recent derivation is given in Chen et al. [23].

A bistatic radar cross section that reduces exactly to the monostatic term of Barrick [12,13] would be beneficial to systems based on Barrick’s expression (such as the Seaview inversion) as, in the radar coverage area, the bistatic angle can vary between 0° and 90°, so the discontinuity between the monostatic and bistatic radar cross section expressions would cause a discontinuity in the inversion program used and perhaps, then, the results. Therefore, in this work, we follow the method of Barrick, whilst retaining the bistatic angle, to derive the bistatic radar cross section of the ocean surface. A reviewer of this paper has drawn our attention to similar work by Hisaki and Tokuda [24] who allowed for a finite scattering area and showed that their equations reduced to those of Barrick for an infinite scattering area and a monostatic geometry.

We begin with an overview of the derivation of the bistatic radar cross section in Section 2.1, before presenting the numerical solution of the resulting expression in Section 2.2. Details of the Seaview inversion method (for which details of the cross-section equations and numerical simulations are a pre-requisite) are then given in Section 3. The results of the modified Seaview inversion, when tested on simulated bistatic data, are given in Section 4 and these are discussed in Section 5 which also includes some concluding remarks.

2. Materials and Methods

2.1. Bistatic Radar Cross Section of the Ocean Surface

To derive the bistatic radar cross section of the ocean surface, we follow the method of Barrick [12], where the equivalent monostatic radar cross section was derived. In his work, Barrick used the perturbation analysis of Rice [14] where, by assuming small waveheights and slopes, the electric field scattered from the ocean surface, $\overrightarrow{E_s}$, was calculated. They key points of the derivation follow, however, more details can be found in the work of Hardman [25].

The value of $\overrightarrow{E_s}$ depends on both the incident radio waves and the properties of the scattering surface. Firstly, the incident waves will propagate as vertically polarised ground waves. Secondly, as the ocean varies in both time and space, by assuming that these variations are periodic and that the surface is of infinite extent, we can define the surface, $f(x,y,t)$, as a Fourier series expansion, such that

$$\begin{array}{c}\hfill z=f(x,y,t)=\sum _{mnl=-\infty }^{\infty }P(m,n,l){\mathrm{e}}^{ia(mx+ny)-iwlt},\end{array}$$

(4)

for wavenumber $a=2\pi /L$, angular frequency $w=2\pi /T$ and Fourier coefficients $P(m,n,l)$ which are dependent on the integers m, n and l.

The electric field scattered from this surface, will also be periodic with the same fundamental spatial and temporal periods and hence Rice defined ${\overrightarrow{E}}_s$ as a Fourier series. For a perfectly conducting flat surface, the exact solution can be found. Therefore, we perturb the solution around the flat surface, with ordering parameter $k_{0}f$ (for radar wavenumber $k_0$), using Maxwell’s equations and the tangential boundary condition to obtain the first and second order Fourier coefficients. Details of the calculations for vertically polarised waves can be found in the work of Hardman [25]; for horizontally polarised waves, the details are provided by Rice [14]. The resulting scattered electric field has components

$$\begin{array}{cc}\hfill E_x=& 2\sum _{mnl}E(m,n,z,l){\mathrm{e}}^{-i{\omega }_{0}t}\left[i(k_0-am)P(m-\nu ,n,l)\phantom{\sum _{kl}}\right.\hfill \\ \hfill & +\sum _{qrs}\left\{a^2(m-q)(\nu -q)k_0+(k_0-am)b^2(q,r)\right\}\left.Q(m,n,l,q,r,s)\phantom{\sum _{kl}}\right],\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill E_y=& 2a\sum _{mnl}E(m,n,z,l){\mathrm{e}}^{-i{\omega }_{0}t}\left[-inP(m-\nu ,n,l)\phantom{\sum _{qrs}}\right.\hfill \\ \hfill & \left.+\sum _{qrs}\left\{a(n-r)(\nu -q)k_0-n{b}^2(q,r)\right\}Q(m,n,l,q,r,s)\right]\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill E_z=& 2{\mathrm{e}}^{i{k}_{0}x}{\mathrm{e}}^{-i{\omega }_{0}t}+2\sum _{mnl}{\displaystyle \frac{E(m,n,z,l)}{b(m,n)}}{\mathrm{e}}^{-i{\omega }_{0}t}\left[\phantom{\sum _{klJ}}\left(\phantom{b^2}-i(a(m-\nu )k_0+b^2(m,n)\right)P(m-\nu ,n,l)\right.\hfill \\ \hfill & +\sum _{qrs}\left\{(a^3(q-\nu )(m^2+n^2-qm-rn)k_0)\right.\hfill \\ \hfill & +\left.\left.a(a(m^2+n^2)-m{k}_0)b^2(q,r)\right\}Q(m,n,l,q,r,s)\phantom{\sum _{kl}}\right],\hfill \end{array}$$

(7)

where ${\omega }_0$ is the angular frequency of the emitted radio waves,

$$E(m,n,z,l)={\mathrm{e}}^{i(a(mx+ny)+b(m,n)z)}{\mathrm{e}}^{-iwlt},$$

in which

$$\begin{array}{c}\hfill b(m,n)=\left\{\begin{array}{cc}{(k_0^2-a^2{m}^2-a^2{n}^2)}^{1/2}\hfill & \mathrm{if}{m}^2+n^2{k}_0^2/a^2\hfill \\ i{(a^2{m}^2+a^2{n}^2-k_0^2)}^{1/2}\hfill & \mathrm{if}{m}^2+n^2{k}_0^2/a^2\hfill \end{array}\right.,\end{array}$$

and

$$Q(m,n,l,q,r,s)={\displaystyle \frac{P(q-\nu ,r,s)P(m-q,n-r,l-s)}{b(q,r)}}.$$

By definition, the electric field in Equations (5)–(7) corresponds to the scattering of infinite plane waves from a surface of infinite extent. To transform the fields to finitely scattered fields, which is necessary as only a portion of the whole ocean will be illuminated by the radar, we can use the equation presented by Johnstone [21] who followed the work of Stratton [26]. He showed that if the electric field is known on a finite section of an infinitely large volume such as the surface $S_1$ on the hemisphere shown in Figure 4, then the scattered electric field at a point $(x^{\prime },y^{\prime },z^{\prime })$ inside of the volume can be calculated.

The equation is given by

$$\begin{array}{cc}\hfill \overrightarrow{E}(x^{\prime },y^{\prime },z^{\prime })=& {\displaystyle \frac{{\mathrm{e}}^{i{k}_{0}R}}{4\pi R}}{\int }_{S_1}\left\{{\left[\left({\displaystyle \frac{\partial E_x}{\partial z}}-{\displaystyle \frac{\partial E_z}{\partial x}}\right){\overrightarrow{\widehat{a}}}_x+\left({\displaystyle \frac{\partial E_y}{\partial z}}-{\displaystyle \frac{\partial E_z}{\partial y}}\right){\overrightarrow{\widehat{a}}}_y\right]}_{z=0}\right.\hfill \\ \hfill & +i{k}_0{\left[\phantom{{\displaystyle \frac{\partial E_z}{\partial x}}}{E}_{x}cos\theta {\overrightarrow{\widehat{a}}}_x+E_{y}cos\theta {\overrightarrow{\widehat{a}}}_y-(E_{x}sin\theta cos\phi +E_{y}sin\theta sin\phi ){\overrightarrow{\widehat{a}}}_z\right]}_{z=0}\hfill \\ \hfill & +{\left[\left(\phantom{{\displaystyle \frac{\partial E_z}{\partial x}}}{\displaystyle \frac{\partial E_z}{\partial x}}-{\displaystyle \frac{\partial E_x}{\partial z}}\right)sin\theta cos\phi +\left({\displaystyle \frac{\partial E_z}{\partial y}}-{\displaystyle \frac{\partial E_y}{\partial z}}\right)sin\theta sin\phi \right]}_{z=0}\hfill \\ \hfill & ·\left.\left[\phantom{{\displaystyle \frac{\partial E_z}{\partial x}}}sin\theta cos\phi {\overrightarrow{\widehat{a}}}_x+sin\theta sin\phi {\overrightarrow{\widehat{a}}}_y+cos\theta {\overrightarrow{\widehat{a}}}_z\right]\right\}{\mathrm{e}}^{-i\overrightarrow{k_0}·\mathbf{\rho }}d{S}_1,\hfill \end{array}$$

(8)

where $\mathbf{\rho }$ is the vector $(x,y,z)$ on $S_1$, $(R,\theta ,\phi )$ are the radius, polar angle and azimuthal angle measured from the origin to $(x^{\prime },y^{\prime },z^{\prime })$ and $\overrightarrow{k_0}$ is the radar wavevector given by

$$\overrightarrow{k_0}=k_0(sin\theta cos\phi ,sin\theta sin\phi ,cos\theta ).$$

Note that the integral can be evaluated on any plane and $z=0$ is used for convenience. However, as pointed out by Hisaki and Tokuda [24], this choice is in fact the infinite scattering surface limit. Evaluating the integral on z = f as they did changes the resulting power spectrum of the scatter but differences are very small at the Doppler frequencies used for inversion so z = 0 is sufficient for this work. To find the scattered field at a point $(x^{\prime },y^{\prime },z^{\prime })$, the values for $E_x$, $E_y$ and $E_z$ from Equations (5)–(7) are substituted into Equation (8) and the integral is calculated. The vertically polarised component, $E_{\theta }(t)$, is then identified in the resulting expression (as these are the radio waves that the receiver will detect) such that

$$\begin{array}{cc}\hfill E_{\theta }(t)={\displaystyle \frac{i{\mathrm{e}}^{i{k}_{0}R}}{2\pi R}}{L}^2\sum _{mnl}& \left\{\phantom{\sum _q}\right.B(t)[-i{x}_{1}P(m-\nu ,n,l)+\sum _{qrs}{x}_{2}Q]+C(t)[-i{y}_{1}P(m-\nu ,n,l)+\sum _{qrs}{y}_{2}Q]\hfill \\ \hfill & +D(t)[-i{z}_{1}P(m-\nu ,n,l)+\sum _{qrs}{z}_{2}Q]\left.\phantom{\sum _q}\right\},\hfill \end{array}$$

(9)

where,

$$\begin{array}{ccc}\hfill x_1& =am-k_0;\hfill & x_2=a^2(m-q)(\nu -q)k_0+(k_0-am)b^2(q,r)\hfill \\ \hfill y_1& =an;\hfill & y_2=a^2(n-r)(\nu -q)k_0-an{b}^2(q,r)\hfill \\ \hfill z_1& =a(m-\nu )k_0+b^2(m,n);\hfill & z_2=\left[a^3(q-\nu )(m^2+n^2-qm-rn)k_0+a\left(a(m^2+n^2)-m{k}_0\right)b^2(q,r)\right].\hfill \end{array}$$

and,

$$\begin{array}{cc}\hfill B(t)& =\sum _{mnl}cos\phi (k_0+b(m,n)cos\theta )sinc(XR)sinc(YR){\mathrm{e}}^{-i(wl+{\omega }_0)t}\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill C(t)& =\sum _{mnl}sin\phi (k_0+b(m,n)cos\theta )sinc(XR)sinc(YR){\mathrm{e}}^{-i(wl+{\omega }_0)t}\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill D(t)& =\sum _{mnl}-cos\theta (amcos\phi +ansin\phi ){\displaystyle \frac{sinc(XR)sinc(YR)}{b(m,n)}}{\mathrm{e}}^{-i(wl+{\omega }_0)t},\hfill \end{array}$$

(12)

for $XR={\displaystyle \frac{L}{2}}\left(am-k_{0}sin\theta cos\phi \right)$ and $YR={\displaystyle \frac{L}{2}}\left(an-k_{0}sin\theta sin\phi \right)$.

In the scattered electric field in Equation (9), the first order components (namely, the terms including a single P term) represent the single scattering of one electromagnetic wave, to the receiver, from one ocean wave. The second order components (which include a factor of Q) represent doubly scattered electromagnetic waves, to the receiver, from two single ocean waves. The order of the ocean wave, currently denoted by $P(m-\nu ,n,l)$, has not yet been considered and is assumed to be first order. However, in making such an assumption, a second order contribution from first order scattering from second order oceans waves is missed, where a second order ocean wave is the result of the nonlinear interaction between two first order ocean waves.

To allow for the second order hydrodynamic effects in shallow water, Barrick and Lipa [27] used a perturbation method to relate the second order coefficients $P^{(2)}(\overrightarrow{k},\omega )$, of a surface defined by $z={\sum }_{\overrightarrow{k},\omega }P(\overrightarrow{k},\omega ){\mathrm{e}}^{i\overrightarrow{k}·\overrightarrow{r}-i\omega t}$, to the first order coefficients $P^{(1)}(\overrightarrow{k},\omega )$. Their method involved expanding the surface height Fourier coefficients around the flat surface, i.e.,

$$P(\overrightarrow{k},\omega )=P^{(1)}(\overrightarrow{k},\omega )+P^{(2)}(\overrightarrow{k},\omega )+\dots ,$$

alongside boundary conditions from the equations of motion, also expanded to second order. They showed that for ocean waves with wavevectors $\overrightarrow{k_1}$ and $\overrightarrow{k_2}$, with corresponding angular frequencies ${\omega }_1$ and ${\omega }_2$ (related by the dispersion relation of ocean waves given by $\omega =\sqrt{gktanh(kd)}$),

$$\begin{array}{c}\hfill P^{(2)}(\overrightarrow{k^{″}},{\omega }^{″})=\sum _{\overrightarrow{k_1}\overrightarrow{k_2}}\sum _{{\omega }_1{\omega }_2}{\mathrm{\Gamma }}_H(\overrightarrow{k_1},{\omega }_1,\overrightarrow{k_2},{\omega }_2)P^{(1)}(\overrightarrow{k_1},{\omega }_1)P^{(1)}(\overrightarrow{k_2},{\omega }_2),\end{array}$$

(13)

where $\overrightarrow{k^{″}}=\overrightarrow{k_1}+\overrightarrow{k_2}$, ${\omega }^{″}={\omega }_1+{\omega }_2$, and

$$\begin{array}{cc}\hfill {\mathrm{\Gamma }}_H={\displaystyle \frac{1}{2}}& \left\{\phantom{{\displaystyle \frac{(\omega )}{(B)}}}{k}_{1}tanh(k_{1}d)+k_{2}tanh(k_{2}d)\right.+\frac{{\omega }^{″}}{g}\frac{({\omega }_1^3{csch}^2(k_{1}d)+{\omega }_2^3{csch}^2(k_{2}d))}{({\omega }^{″2}-g{k}^{″}tanh(k^{″}d))}\hfill \\ \hfill & +\frac{(k_1{k}_{2}tanh(k_{1}d)tanh(k_{2}d)-\overrightarrow{k_1}·\overrightarrow{k_2})}{\sqrt{k_1{k}_{2}tanh(k_{1}d)tanh(k_{2}d)}}\left(\frac{g{k}^{″}tanh(k^{″}d)+{\omega }^{″2}}{g{k}^{″}tanh(k^{″}d)-{\omega }^{″2}}\right)\left.\phantom{\sum _{kl}}\right\},\hfill \end{array}$$

(14)

is called the hydrodynamic coupling coefficient.

To include the second order hydrodynamic effects in the scattered electric field, we expand $P(m-\nu ,n,l)$ into $P^{(1)}(m-\nu ,n,l)+P^{(2)}(m-\nu ,n,l)$ and then substitute in the value of $P^{(2)}(m-\nu ,n,l)$ using Equation (13) (by letting $\overrightarrow{k^{″}}=(m-\nu ,n)$ and $\omega =l$), and so Equation (9) becomes

$$E_{\theta }(t)={\displaystyle \frac{i{\mathrm{e}}^{i{k}_{0}R}}{2\pi R}}{L}^2\sum _{mnl}\left\{\phantom{\sum _{mnl}}-i\zeta (t)P^{(1)}(m-\nu ,n,l)\right.+\sum _{qrs}\left[-i\zeta (t){\mathrm{\Gamma }}_{H}b(p,q)+\xi (t)\right]\left.\phantom{\sum _{mnl}}Q(m,n,l,q,r,s)\right\},$$

(15)

where for brevity

$$\zeta (t)=B(t)x_1+C(t)y_1+D(t)z_1$$

and

$$\xi (t)=B(t)x_2+C(t)y_2+D(t)z_2.$$

Following Johnstone [21], we finally calculate the radar cross section by substituting Equation (15) into

$$\begin{array}{c}\hfill \sigma (\omega )=\underset{R\to \infty }{lim}4\pi R^2{\displaystyle \frac{\mathcal{F}\left[\langle E_{\theta }(t_1)E_{\theta }^{*}(t_2)\rangle \right]}{L^2}},\end{array}$$

(16)

where R is the distance from the scatter patch to the receiver and $\mathcal{F}$ denotes a Fourier transform, with definition

$$\begin{array}{c}\hfill \mathcal{F}\left[f(t)\right]={\displaystyle \frac{1}{2\pi }}{\int }_{\infty }^{\infty }f(t){\mathrm{e}}^{-i\omega t}dt.\end{array}$$

(17)

To calculate Equation (16), the following properties of the surface height Fourier coefficients are used:

• The Fourier coefficients are normally distributed about zero; hence

  $$\begin{array}{c}\hfill \langle P(m,n,l)\rangle =0.\end{array}$$

  (18)
• As the surface is real, $f(x,y,t)$ is equal to $f^{*}(x,y,t)$, which is true when

  $$\begin{array}{c}\hfill P(-m,-n,-l)=P^{*}(m,n,l).\end{array}$$

  (19)
• From Thomas [28],

  $$\begin{array}{c}\hfill \langle P_1{P}_2{P}_3\rangle =0\end{array}$$

  (20)

  and

  $$\langle P_1{P}_2{P}_3{P}_4\rangle =\langle P_1{P}_2\rangle \langle P_3{P}_4\rangle +\langle P_1{P}_3\rangle \langle P_2{P}_4\rangle +\langle P_1{P}_4\rangle \langle P_2{P}_3\rangle .$$

  (21)
• The surface roughness spectrum $S(p,q,wl)$, found by utilising the Wiener–Khinchin theorem, is related to the surface height Fourier coefficients by

  $$\langle P(m,n,l)P(q,r,s)\rangle =\left\{\begin{array}{cc}{\displaystyle \frac{{(2\pi )}^{3}S(p,q,wl)}{L^{2}T}}\hfill & \mathrm{if}q,r,s=-m,-n,-l\hfill \\ 0\hfill & \mathrm{if}\mathrm{else},\hfill \end{array}\right.$$

  (22)

  where $p=am$ and $q=an$.

Then, substituting $E_{\theta }(t)$ from Equation (15) into Equation (16) and using Equation (20) leads to

$$\begin{array}{cc}\hfill \sigma (\omega )& ={\displaystyle \frac{1}{\pi }}\mathcal{F}\left[L^2\sum _{\begin{array}{c}mnl\\ m^{\prime }{n}^{\prime }{l}^{\prime }\end{array}}\left\{\zeta (t_1){\zeta }^{\prime *}(t_2)〈P^{(1)}(m-\nu ,n,l)P^{{(1)}^{\prime }*}(m^{\prime }-\nu ,n^{\prime },l^{\prime })〉\right.\right.\hfill \\ \hfill & +\sum _{\begin{array}{c}qrs\\ q^{\prime }{r}^{\prime }{s}^{\prime }\end{array}}\left\{\phantom{\sum }\left[-i\zeta (t_1)b(q,r){\mathrm{\Gamma }}_H+\xi (t_1)\right]\left[i{\zeta }^{\prime *}(t_2)b^{*}(q^{\prime },r^{\prime }){\mathrm{\Gamma }}_H^{\prime }+{\xi }^{\prime *}(t_2)\right]〈Q{Q}^{\prime *}〉\phantom{\sum }\right\}\left.\phantom{〈P^{(1)}{P}^{{(1)}^{\prime }*}〉}\right\}\left.\phantom{\sum _{\begin{array}{c}mnl\\ m^{\prime }{n}^{\prime }{l}^{\prime }\end{array}}{\displaystyle \frac{1}{L^2}}}\right],\hfill \end{array}$$

(23)

where the arguments of Q are implied. The calculation of Equation (23) can be separated into its first and second order terms, such that

$$\sigma (\omega )={\sigma }^{(1)}(\omega )+{\sigma }^{(2)}(\omega ),$$

where the first order radar cross section, ${\sigma }^{(1)}(\omega )$, is defined by the term including the average $〈P^{(1)}(m-\nu ,n,l)P^{{(1)}^{\prime }*}(m^{\prime }-\nu ,n^{\prime },l^{\prime })〉$, and the second order, ${\sigma }^{(2)}(\omega )$, including $〈Q{Q}^{\prime *}〉$. To calculate each of ${\sigma }^{(1)}(\omega )$ and ${\sigma }^{(2)}(\omega )$, the properties in Equations (18)–(22) are used to enforce restrictions on the Fourier coefficients and to introduce the roughness spectrum $S(\overrightarrow{k})$. The mathematical details are spared here, but can be found in the work of Hardman [25].

2.1.1. First Order

The first order radar cross section is given by

$$\begin{array}{c}\hfill {\sigma }^{(1)}(\omega )=2^5\pi k_0^4{cos}^4{\phi }_{bi}\sum _{m=\pm 1}S(m\overrightarrow{k_B})\delta (\omega -m{\omega }_B),\end{array}$$

(24)

defined at the Bragg frequencies, $\pm {\omega }_B$, where

$$\begin{array}{c}\hfill {\omega }_B=\sqrt{2g{k}_{0}cos{\phi }_{bi}tanh(2k_{0}dcos{\phi }_{bi})},\end{array}$$

(25)

for the bistatic angle, ${\phi }_{bi}$, which is shown in Figure 3 and Figure 5, and is related to the azimuthal scatter angle $\phi$ by

$$\begin{array}{c}\hfill {\phi }_{bi}={\displaystyle \frac{1}{2}}\left(\pi -\phi \right).\end{array}$$

(26)

The value of ${\sigma }^{(1)}(\omega )$ depends on the ocean spectrum contribution for the Bragg wavevector, $\overrightarrow{k_B}$, which travels in the elliptical normal direction from the scatter point (as shown in Figure 3 and Figure 5), and is defined by

$$\begin{array}{c}\hfill \overrightarrow{k_B}=-2k_{0}cos{\phi }_{bi}(cos{\phi }_{bi},-sin{\phi }_{bi}).\end{array}$$

(27)

2.1.2. Second Order

The second order bistatic radar cross section is

$$\begin{array}{cc}\hfill {\sigma }^{(2)}(\omega )=& 2^5\pi k_0^4{cos}^4{\phi }_{bi}\sum _{m,m^{\prime }=\pm 1}\int \int {\left|{\mathrm{\Gamma }}_E-i{\mathrm{\Gamma }}_H\right|}^{2}S(\overrightarrow{m{k}_1})S(\overrightarrow{m^{\prime }{k}_2})\delta (\omega -m{\omega }_1-m^{\prime }{\omega }_2)dpdq,\hfill \end{array}$$

(28)

for wavevector pairs $\overrightarrow{k_1}$ and $\overrightarrow{k_2}$ (with respective angular frequencies ${\omega }_1$ and ${\omega }_2$) such that

$$\overrightarrow{k_1}+\overrightarrow{k_2}=\overrightarrow{k_B}.$$

Explicitly,

$$\begin{array}{c}\hfill \overrightarrow{k_1}=(p-k_0,q)\mathrm{and}\overrightarrow{k_2}=(-k_{0}cos(2{\phi }_{bi})-p,k_{0}sin(2{\phi }_{bi})-q),\end{array}$$

(29)

and ${\mathrm{\Gamma }}_E$ is the electromagnetic coupling coefficient given by

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_E={\displaystyle \frac{1}{2^2{cos}^2{\phi }_{bi}}}\left({\displaystyle \frac{a_1}{b_1-k_0\Delta }}+{\displaystyle \frac{a_2}{b_2-k_0\Delta }}\right),\end{array}$$

(30)

where

$$\begin{array}{c}\hfill \Delta =0.011-0.012i,\end{array}$$

(31)

is the normalized surface impedance derived by Barrick [29] and

$$\begin{array}{c}\hfill a_1=-k_{1x}(\overrightarrow{k_2}·\overrightarrow{\widehat{a}})-2{cos}^2{\phi }_{bi}\left(-k_2^2+2k_0(\overrightarrow{k_2}·\overrightarrow{\widehat{a}})\right),\end{array}$$

(32)

$$\begin{array}{c}\hfill a_2=-k_{2x}(\overrightarrow{k_1}·\overrightarrow{\widehat{a}})-2{cos}^2{\phi }_{bi}(-k_1^2+2k_0(\overrightarrow{k_1}·\overrightarrow{\widehat{a}})),\end{array}$$

(33)

$$\begin{array}{c}\hfill b_1=\sqrt{-k_2^2+2k_0(\overrightarrow{k_2}·\overrightarrow{\widehat{a}})}\end{array}$$

(34)

and

$$\begin{array}{c}\hfill b_2=\sqrt{-k_1^2+2k_0(\overrightarrow{k_1}·\overrightarrow{\widehat{a}})}\end{array}$$

(35)

(noting that both $b_1$ and $b_2$ can be real or imaginary depending on the argument), where $\overrightarrow{\widehat{a}}$ is a unit vector in the direction of the receiver from the scattering patch (see Figure 5), namely,

$$\overrightarrow{\widehat{a}}=(-cos(2{\phi }_{bi}),sin(2{\phi }_{bi})).$$

2.1.3. Monostatic Conditions

When ${\phi }_{bi}=0$, i.e., under monostatic conditions, the first and second order radar cross section expressions given in Equations (24) and (28), respectively, are equivalent to the commonly used first order monostatic radar cross section of Lipa and Barrick [15], except for a factor of 2. The difference is due to differing Fourier transform definitions, however, the factor is ultimately not important as when the inverse Fourier transform (which will also be different by a factor of 2) is taken to find the power in the spectrum, the two terms will be equal. Furthermore, when inverting the expression, the whole spectrum is normalised to removed the effects of propagation over the ocean and so the factor is again inconsequential.

2.2. Numerical Solution

The method for finding the numerical solution of the bistatic radar cross section given in Equations (24) and (28), is analogous to the method of Holden and Wyatt [16].

2.2.1. First Order

Due to the delta function in Equation (24), the first order contribution to the Doppler spectrum will appear as two peaks at $\pm {\omega }_B$, defined in Equation (25). The contribution comes from two wave vectors, $\pm \overrightarrow{k_B}$, travelling in the direction of the Bragg bearing, both toward and away from the radar set up. As Crombie [1] hypothesised (for monostatic radar), these particular signals are amplified by resonance. As $k_B$ includes a factor of $cos{\phi }_{bi}$, for one radar using a single carrier frequency, there will be a number of different Bragg waves, dependent on the radar beam range and angle. Scattering from locations where ${\phi }_{bi}\to 90°$ are referred to as forward scatter and when this occurs, both ${\omega }_B$ and $k_B$ tend to 0. Consequently, the wavelengths of the Bragg waves in this region become infinitely long and in addition, the ${cos}^4{\phi }_{bi}$ factors in Equations (24) and (28) will tend to zero, leading to a low SNR.

2.2.2. Second Order

The second order radar cross section contribution, given in Equation (28), is due to double electromagnetic scattering from two first order ocean waves, with wavevectors $\overrightarrow{k_1}$ and $\overrightarrow{k_2}$, and the nonlinear interaction between the same two ocean waves. The wave vector pair sum to give the Bragg wavevector $\overrightarrow{k_B}$ and, theoretically, large numbers of pairs exist in the p, q wavenumber plane; an example pair can be seen in Figure 6. To find the values of $\overrightarrow{k_1}$ and $\overrightarrow{k_2}$, the solution to the delta function in Equation (28), such that

$$\begin{array}{c}\hfill \omega -m\sqrt{g{k}_{1}tanh(k_{1}d)}-m^{\prime }\sqrt{g{k}_{2}tanh(k_{2}d)}=0\end{array}$$

(36)

is sought. For each value of $\omega$, the set of solutions of $\overrightarrow{k_1}$ and $\overrightarrow{k_2}$ defines a frequency contour in the p, q plane. As m and $m^{\prime }$ can both take the values of 1 or $-1$, Equation (36) has four different forms.

Case $m=m^{\prime }$:

Squaring Equation (36) gives

$${\omega }^2=g\left(k_{1}tanh(k_{1}d)+k_{2}tanh(k_{2}d)+2\sqrt{k_{1}tanh(k_{1}d)k_{2}tanh(k_{2}d)}\right)$$

and it can be shown that

$${\omega }^2\ge (k_1+k_2)tanh((k_1+k_2)d).$$

Now, as the sum of two sides of a triangle is greater than the third, $k_1+k_{2}2k_{0}cos{\phi }_{bi}$, and therefore

$${\omega }^{2}2g{k}_{0}cos{\phi }_{bi}tanh(2g{k}_{0}cos{\phi }_{bi}d).$$

Taking the square root gives

$$\left|\omega \right|\sqrt{2g{k}_{0}cos{\phi }_{bi}tanh(2g{k}_{0}cos{\phi }_{bi}d)},$$

which is equal to the bragg frequency ${\omega }_B$ given in Equation (25). This leads to two conditions;

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}\omega {\omega }_B\hfill & m=m^{\prime }=1\hfill \\ \omega \le -{\omega }_B\hfill & m=m^{\prime }=-1.\hfill \end{array}\right.\end{array}$$

Figure 7 shows the frequency contours, defined by Equation (36). At frequencies close to the Bragg frequencies the contours are circular in shape, centred around the Bragg frequency. As $\omega$ increases, the contours become less circular, until $\left|\omega \right|=2\sqrt{g{k}_{0}cos{\phi }_{bi}tanh(k_{0}dcos{\phi }_{bi})}$, shown in white in Figure 7, where they separate.

When ${\phi }_{bi}=0$, the contours are symmetrical about the p and q axes, however, when ${\phi }_{bi}0$, the contours rotate clockwise in the p, q plane, becoming symmetrical about some other axes, say $p^{\prime }$ and $q^{\prime }$, shown by the additional black lines in Figure 7b.

Case $m\ne m^{\prime }$:

When $m\ne m^{\prime }$, the square of Equation (36) is

$$\begin{array}{c}\hfill {\omega }^2=g\left(k_{1}tanh(k_{1}d)+k_{2}tanh(k_{2}d)-2\sqrt{k_{1}tanh(k_{1}d)k_{2}tanh(k_{2}d)}\right),\end{array}$$

and it can be shown that

$${\omega }^2\le g\left((k_2-k_1)tanh((k_2-k_1)d)\right).$$

In the right hand plane, when $\overrightarrow{k_1}$ and $\overrightarrow{k_2}$ lie in opposite directions along the $p^{\prime }$ axis, meeting at a point past the bragg frequency, $k_2-k_1$ reaches its maximum value of $2k_{0}cos{\phi }_{bi}$. Therefore,

$$\begin{array}{c}\hfill {\omega }^2\le g(2k_{0}cos{\phi }_{bi}tanh(2k_{0}cos{\phi }_{bi}d))\end{array}$$

and this leads to the conditions

$$\left\{\begin{array}{cc}0\omega \le {\omega }_B\hfill & m=-1,m^{\prime }=1\hfill \\ -{\omega }_B\omega \le 0\hfill & m=1,m^{\prime }=-1.\hfill \end{array}\right.$$

In the left hand plane, where $k_2{k}_1$, the result is reversed giving

$$\left\{\begin{array}{cc}0\omega \le {\omega }_B\hfill & m=1,m^{\prime }=-1\hfill \\ -{\omega }_B\omega \le 0\hfill & m=-1,m^{\prime }=1.\hfill \end{array}\right.$$

Figure 8 shows the normalised frequency contours when $m\ne m^{\prime }$ for different values of ${\phi }_{bi}$. Like when $m=m^{\prime }$, the contours are symmetric about the $p^{\prime }$ and $q^{\prime }$ axes, however, in this case, they do not cross the $q^{\prime }$ axis for any frequency. The contours depict small circles around the Bragg frequencies, growing in size and becoming less circular as the frequency approaches zero.

As the frequency contours for all four possible combinations of m and $m^{\prime }$ are symmetrical about the $q^{\prime }$ axis, the integration in Equation (28) can be taken over one half of the symmetric plane and doubled. Therefore, we integrate Equation (28) over the right hand $p^{\prime }$ plane, and double the result and hence

$${\sigma }_2(\omega )=2^6\pi k_0^4{cos}^4{\phi }_{bi}\sum _{m,m^{\prime }=\pm 1}{\int }_{-\infty }^{\infty }{\int }_{q^{\prime }}^{\infty }|{\mathrm{\Gamma }}_E-i{\mathrm{\Gamma }}_H{|}^{2}S(m\overrightarrow{k_1})S(m^{\prime }\overrightarrow{k_2})\delta (\omega -m{\omega }_1-m^{\prime }{\omega }_2)dpdq.$$

(37)

In the right hand $p^{\prime }$ plane, where $k_1\le k_2$, we can calculate the integral in polar coordinates $k_1$ and ${\theta }_1$, where ${\theta }_1$ is the angle between $\overrightarrow{k_1}$ and the p axis, as shown Figure 6. Explicitly,

$${\sigma }_2(\omega )=2^6\pi k_0^4{cos}^4{\phi }_{bi}\sum _{m,m^{\prime }=\pm 1}{\int }_{-{\theta }_L^{-}}^{{\theta }_L^{+}}{\int }_0^{\infty }{|{\mathrm{\Gamma }}_E-i{\mathrm{\Gamma }}_H|}^{2}S(m\overrightarrow{k_1})S(m^{\prime }\overrightarrow{k_2})\delta (\omega -m{\omega }_1-m^{\prime }{\omega }_2)k_{1}d{k}_{1}d{\theta }_1,$$

(38)

where ${\theta }_L^{-}$ and ${\theta }_L^{+}$ are the integration limits of ${\theta }_1$ and vary with frequency as well as ${\phi }_{bi}$. In general, we find that the integration limits are

$${\theta }_L^{-}=-(\pi +{\phi }_{bi})\mathrm{and}{\theta }_L^{+}=\pi -{\phi }_{bi},$$

(39)

however a particular set of $\omega$ values require different limits. This happens when $m=m^{\prime }$ and $\left|\omega \right|2\sqrt{k_{0}gcos{\phi }_{bi}tanh(k_{0}dcos{\phi }_{bi})}$, as the frequency contours cross the $q^{\prime }$ axis and no longer complete full rotations in the right hand $p^{\prime }$ plane. By symmetry, when the contours cross the $q^{\prime }$ axis, $k_1$ and $k_2$ are the same length. Therefore, when $k_2=k_1$ and $m=m^{\prime }$, the delta constraint of Equation (36) becomes

$$\begin{array}{c}\hfill \omega =\pm 2\sqrt{g{k}_{1}tanh(k_{1}d)}\end{array}$$

(40)

and then squaring Equation (40) gives

$$\begin{array}{c}\hfill {\displaystyle \frac{{\omega }^2}{4g}}=k_{1}tanh(k_{1}d),\end{array}$$

(41)

which can be solved numerically for $k_1$.

Introducing a term, $\beta$, as the angle between $\overrightarrow{k_1}$ and $\overrightarrow{k_B}$ (see Figure 6), the integration limits are given by

$$\begin{array}{c}\hfill {\theta }_L^{+}=\pi -{\phi }_{bi}-\beta \mathrm{and}{\theta }_L^{-}=-\left(2\pi -{\theta }_L^{+}-2\beta \right).\end{array}$$

(42)

Therefore, as $\beta ={cos}^{-1}\left({\displaystyle \frac{k_B}{2k_1}}\right)$, where the solution for $k_1$ from Equation (41) is used,

$$\begin{array}{c}\hfill {\theta }_L^{+}=\pi -{\phi }_{bi}-{cos}^{-1}\left({\displaystyle \frac{k_B}{2k_1}}\right)\mathrm{and}{\theta }_L^{-}=-\left(\pi +{\phi }_{bi}-{cos}^{-1}\left({\displaystyle \frac{k_B}{2k_1}}\right)\right).\end{array}$$

(43)

In terms of $k_1$ and ${\theta }_1$,

$$\begin{array}{c}\hfill k_2=\sqrt{k_1^2+k_B^2+2k_B{k}_{1}cos({\theta }_1+{\phi }_{bi})}\mathrm{and}{\theta }_2=\pi +{\theta }_1-{cos}^{-1}\left({\displaystyle \frac{k_1+k_{B}cos({\theta }_1+{\phi }_{bi})}{k_2}}\right).\end{array}$$

To calculate ${\sigma }^{(2)}(\omega )$, we now reduce the double integral in Equation (38) to a single integral using the delta function. By defining

$$\begin{array}{cc}\hfill y_s& =\sqrt{k_1}\hfill \\ \hfill h(y_s,{\theta }_1)& =m{y}_s\sqrt{gtanh(y_s^{2}d)}+m^{\prime }\sqrt{g{k}_{2}tanh(k_{2}d)}\hfill \\ \hfill I(y_s,{\theta }_1)& =2^7\pi {\left|{\mathrm{\Gamma }}_T\right|}^2{k}_0^4{cos}^4{\phi }_{bi}S(m\overrightarrow{k_1})S(m^{\prime }\overrightarrow{k_2})y_s^3,\hfill \end{array}$$

Equation (38) can be written as

$$\begin{array}{c}\hfill {\sigma }_2(\omega )={\int }_{-{\theta }_L^{-}}^{{\theta }_L^{+}}{\int }_0^{\infty }I(y_s,{\theta }_1)\delta \left(\omega -h(y_s,{\theta }_1)\right)d{y}_{s}d{\theta }_1.\end{array}$$

(44)

Now, in order to integrate over the delta function, Equation (44) should have an integration variable of h. Therefore, we calculate

$$\begin{array}{c}\hfill {\sigma }_2(\omega )={\int }_{-{\theta }_L^{-}}^{{\theta }_L^{+}}{\int }_0^{\infty }I(y_s,{\theta }_1)\delta \left(\eta -h(y_s,{\theta }_1)\right)|{\displaystyle \frac{\partial y_s}{\partial h}}{|}_{\theta }dhd{\theta }_1,\end{array}$$

(45)

where

$$\begin{array}{cc}\hfill {\left|\frac{\partial h}{\partial y_s}\right|}_{{\theta }_1}=\sqrt{g}& \left\{m\left(\sqrt{tanh(y_s^{2}d)}+{\displaystyle \frac{y_s^{2}d({sech}^2(y_s^{2}d))}{\sqrt{tanh(y_s^{2}d)}}}\right)\right.\hfill \\ & \left.+{\displaystyle \frac{m^{\prime }(y_s^3+y_s{k}_{B}cos({\theta }_1+{\phi }_{bi}))}{k_2^{3/2}}}\left\{\sqrt{tanh(k_{2}d)}+k_{2}d{\displaystyle \frac{sec{h}^2(k_{2}d)}{\sqrt{tanh(k_{2}d)}}}\right\}\right\},\hfill \end{array}$$

whose reciprocal is $|{\displaystyle \frac{\partial y_s}{\partial h}}{|}_{{\theta }_1}$. To integrate over the delta function, the solution, $y^{*}$, to

$$\begin{array}{c}\hfill \omega -h(y^{*},{\theta }_1)=0\end{array}$$

(46)

is required, which can be found using a numerical method. For timely convergence in the numerical method, a good initial guess for $y^{*}$ is important. As the solution for shallow water should not be too different to that for deep water, we find an initial solution for the deep water case and use that, as a starting point, for the shallow water case. The deep water equation can be solved exactly in two cases:

• When $m{m}^{\prime }=1$ and ${\theta }_1=-{\phi }_{bi}$, the solution of $f(y)$ is

  $$\begin{array}{c}\hfill y_0^{*}={\displaystyle \frac{{\omega }^2-g{k}_B}{2m\sqrt{g}\omega }}.\end{array}$$

  (47)
• When $m{m}^{\prime }=-1$ and ${\theta }_1=\pi -{\phi }_{bi}$, the solution is

  $$\begin{array}{c}\hfill y_0^{*}={\displaystyle \frac{m\omega +\sqrt{2g{k}_B-{\omega }^2}}{2\sqrt{g}}}.\end{array}$$

  (48)

Upon finding $y^{*}$, the second order cross section calculation in Equation (45) reduces to

$$\begin{array}{c}\hfill {\sigma }_2(\omega )={\int }_{-{\theta }_L^{-}}^{{\theta }_L^{+}}I(y_s,{\theta }_1)|{\displaystyle \frac{\partial y_s}{\partial h}}{|}_{{\theta }_1}{|}_{y_s=y^{*}}d{\theta }_1,\end{array}$$

(49)

which can be calculated using a numerical integration method. For speed and convergence we update the value of $y_0^{*}$ to the previously found solution for $y^{*}$, as ${\theta }_1$ incrementally increases.

2.2.3. Electromagnetic Singularities

The electromagnetic coupling coefficient given in Equation (30), contains two singularities; either when $b_1$ or $b_2$ is equal to zero. The singularities lie on two circles in the p, q plane, shown in Figure 9; explicitly,

$$\begin{array}{c}\hfill p^2+q^2=k_0^2\end{array}$$

(50)

and

$$\begin{array}{c}\hfill {(k_{0}cos\phi -p+k_0)}^2+{(k_{0}sin\phi -q)}^2=k_0^2.\end{array}$$

(51)

Each singularity will be most prominent when a frequency contour is tangential to the singular circle. In order to find the frequencies that this is true for, the solutions for p and q such that the gradient of the frequency contour expression of Equation (36) is equal to the gradient of the circle functions of Equations (50) and (51) are sought. Knowledge of the geometry of the contours and the radii of the circles is exploited to find the solutions for p and q and then the solutions are substituted into Equation (36) to give the tangential frequencies. The four solutions for p and q are

$$\begin{array}{c}\hfill p=k_{0}sin{\phi }_{bi}\mathrm{and}q=k_{0}cos{\phi }_{bi},\end{array}$$

(52)

$$\begin{array}{c}\hfill p=-k_{0}sin{\phi }_{bi}\mathrm{and}q=-k_{0}cos{\phi }_{bi},\end{array}$$

(53)

$$\begin{array}{c}\hfill p=k_{0}sin{\phi }_{bi}+2k_0{sin}^2{\phi }_{bi}\mathrm{and}q=k_{0}cos{\phi }_{bi}+2k_{0}sin{\phi }_{bi}cos{\phi }_{bi}\end{array}$$

(54)

and

$$\begin{array}{c}\hfill p=-k_{0}sin{\phi }_{bi}+2k_0{sin}^2{\phi }_{bi}\mathrm{and}q=-k_{0}cos{\phi }_{bi}+2k_{0}sin{\phi }_{bi}cos{\phi }_{bi}.\end{array}$$

(55)

Substituting the solutions for p and q, from Equations (52)–(55) into Equation (36) gives two distinct tangential frequencies:

$$\begin{array}{c}\hfill \omega =2^{3/4}{\omega }_B\sqrt{{\displaystyle \frac{\sqrt{1\pm sin{\phi }_{bi}}}{cos{\phi }_{bi}}}}{\displaystyle \frac{\sqrt{tanh(d\sqrt{2k_0^2(1\pm sin{\phi }_{bi})})}}{\sqrt{tanh(2k_{0}cos{\phi }_{bi}d)}}},\end{array}$$

(56)

where the solution with the + signs is for the p and q in Equations (52) and (53), and the − signs for the solutions of p and q in Equations (54) and (55).

These values for $\omega$ are highlighted in Figure 7 by the white contours and are shown to be tangential to the circles expressed in Equations (50) and (51). Both singularities are highlighted in a simulated Doppler spectrum in Figure 10 by dashed vertical lines. A low amplitude Gaussian noise spectrum has been added as can be seen at the extremities of the plot. For deep water, or when $d\to \infty$, the values for $\omega$ become

$$\begin{array}{c}\hfill \omega =2^{3/4}{\omega }_B\sqrt{{\displaystyle \frac{\sqrt{1\pm sin{\phi }_{bi}}}{cos{\phi }_{bi}}}},\end{array}$$

(57)

which agree with the results of Gill Walsh [22].

2.2.4. Currents

An ocean current affects how an ocean wave propagates, both in direction and speed. The change in speed means that, because of the Doppler effect, the entire spectrum is subject to an additional shift, $\Delta \omega$. The additional shift is

$$\begin{array}{c}\hfill \Delta \omega =2k_0{\nu }_E(\phi )cos{\phi }_{bi},\end{array}$$

(58)

where ${\nu }_E(\phi )$ is the component of the current velocity in the elliptical normal direction for beam angle $\phi$.

3. Inversion of $\mathbf{\sigma }(\mathbf{\omega })$ to Measure the Directional Wave Spectrum

The aim of inversion is to obtain the ocean wave directional spectrum, $S(\overrightarrow{k})$, from the power spectrum of the measured radar cross section, $\sigma (\omega )$ using Equations (24) and (28). The power spectrum of the backscattered radar signal is proportional to $\sigma (\omega )$ and therefore in principle needs to be calibrated to account for antenna gains, propagation losses and other factors. To avoid this the problem is usually framed in terms of the ratio of the second and first order backscatter power spectra. A number of inversion methods have been published e.g., [6,8,9,15,16,30]. Here we use the method of Wyatt [10,11] which is referred to in this paper as the Seaview inversion method since Seaview Sensing Ltd has an exclusive license from the University of Sheffield to commercialise the software package and continue its development.

The Seaview Inversion Method

The method used makes the assumption that the first order Bragg wave, $\overrightarrow{k_B}$, is generated by the local wind and that, by limiting the Doppler frequency range used in the inversion, the waves contributing to the second order scatter can be separated into long waves $\overrightarrow{k_1}$ and short waves $\overrightarrow{k_2}$ the latter also being locally wind-driven. These short wind waves are then modelled with a Pierson–Moskowitz spectrum [31], using an initial waveheight estimate obtained using an empirical method [32], and a sech${}^2$ directional distribution [33] the parameters of which are determined using a wind direction estimation model [4].

The method is iterative and is initialised assuming the wind–wave model applies at all wavenumbers. The ratio of the Equations (1) and (3), $\frac{{\sigma }_2(\omega )}{{\sigma }_1(\omega )}$, is then integrated to provide a simulated Doppler spectrum ratio which is compared with the measured Doppler spectrum ratio to obtain the difference between them at each Doppler frequency. The wave spectrum is then updated taking into account the calculated Doppler spectrum ratio difference and the value of the coupling coefficient at the long wave vector wavenumber, $\overrightarrow{k_1}$ relative to the maximum along the Doppler frequency contour. The iteration continues until convergence is achieved or a specified number of iterations, usually 100, has been reached. The quality of the convergence is measured by a quantity that reflects the difference between the measured and simulated ratio. The solution can be very different from the initialising spectrum and often shows bimodality in frequency or direction or both due to the presence of swell or changing wind conditions. More details of the method can be found in Wyatt [10], and Green and Wyatt [11]. The maximum frequency (or equivalently wavenumber) of the long waves that can be measured is dependant on the radio frequency [5,34]. The minimum frequency depends on the quality of the radar data which impacts on the ability to clearly separate first and second order parts of the measured Doppler spectrum. An independent validation of the method applied to monostatic data is presented in [35].

The inversion software, providing surface current, wind direction and wave information, was written for monostatic radar configurations only. This has been extended to bistatic configurations using the analysis in Section 2.1 and Section 2.2 above with some modifications to account for a difference in the coordinate system used in the Seaview software. The bistatic extension is not yet part of the commercial package. It has been tested using simulated data, using the methods given in Section 2.2, of two types: (a) two monostatic radars (in order to check that the bistatic extension to the Seaview package provides the same results for zero bistatic angle as the original monostatic package); (b) one monostatic and one bistatic radar sharing a transmitter at the monostatic site. Tests using two bistatic radars with a transmitter between the two sites are in progress but the results are not ready to be reported on at this time. The two receiver sites are those of the University of Plymouth wavehub WERA radar [36,37] and a limited number of cells from the coverage grid for that radar have been used to provide different bistatic angles for the simulations. Using an existing configuration made it easier to provide data for the inversions in standard Seaview formats. The wave parameters used in each simulation are the same at all cells and propagations losses and any antenna effects are not included so the only differences at each cell are the bistatic angle and the Bragg wave bearing relative to wave, wind and current directions.

The wave parameters for the different simulations are described in

Table 1. Wave, wind and current parameters used for the Doppler spectra simulations. The buoy data are not separated into wind–waves and swell but their peak period and direction are included in the swell columns.

Case	Type	Wind–Wave Hs m	${\mathit{\theta }}_{\mathit{w}}$	Spread	Wind Speed m/s	Current Speed m/s	${\mathit{\theta }}_{\mathit{c}}$	Swell Hs m	${\mathit{\theta }}_{\mathit{s}}$	Tp s	Spread
1	Model	3.07	0.0	3.0	12.0	1.4	70.0
2	Model	3.07	30.0	2.0	12.0	1.8	90.0	3.0	140.0	13.2	10.0
3	Buoy	1.72				1.8	90.0		68.24	12.8
4	Buoy	1.72				1.8	90.0		148.16	6.74
5	Buoy	2.70				1.8	90.0		109.01	14.22
6	Buoy	5.87				1.0	245.0		162.36	9.85

. They include modelled and buoy-measured wave spectra.

4. Results

Figure 11 and Figure 12 show inverted surface current speed and direction, wind direction and significant waveheight and spectral peak direction maps for the model cases in Table 1 for the monostatic and bistatic configurations. For the bimodal case 2, Figure 13 shows the long (swell, here defined as waves in 0.05-0.1Hz band) and short (wind–wave, 0.1–0.2Hz) contributions separately to confirm that the latter are aligned with the wind. Figure 14, Figure 15, Figure 16 and Figure 17 show current, wind and wave maps for the buoy cases. Figure 18, Figure 19, Figure 20 and Figure 21 show sample directional spectra compared with those measured by the buoy and used in the simulations to provide a qualitative validation of the radar measured spectra. They are from 4 selected locations to cover the key parameter ranges expected to be important in the accuracy of the inversion. The key parameters are the bistatic angles and the difference in angle between the two Bragg directions (a minimum of $30°$ is required for monostatic processing [38]) and are presented in

Table 2. Configuration parameters at selected cells.

Cell Number	Configuration	Bistatic Angle	Angle between Braggs
1664	monostatic	0	40.6
	1 mono, 1 bistatic	20.4	20.3
3116	monostatic	0	60.0
	1 mono, 1 bistatic	30.1	29.9
3128	monostatic	0	77.8
	1 mono, 1 bistatic	39.1	38.8
3140	monostatic	0	99.6
	1 mono, 1 bistatic	50.1	49.5

. Note that for some of the locations the Bragg angle difference is below the suggested monostatic threshold. Cell 1664 is the one on the left of the top row in the maps, cells 3116, 3128 and 3140 are the lower three going south along the column to the east of $-5°{36}^{\prime }$.

There is generally good agreement both for both the standard monostatic and the bistatic case. Differences will be discussed further in the next section.

Scatter plots and statistics of the comparisons for currents are presented in Figure 22 and for waves in Figure 23. The data are colour-coded according to the bistatic angle with red being the largest bistatic angle. There is some dependence of the accuracy of the wave measurements on this parameter as can be seen in the right hand column of Figure 23. Most of the larger differences in peak period and peak direction are associated with the bimodal model case where the swell and wind–waves peaks were of similar magnitude and small differences in these magnitudes can lead to differences in peak identification. This is also evident when comparing Figure 12 with Figure 13.

5. Discussion

Monostatic and bistatic radar data have been simulated using the backscatter cross-section formulations developed in Section 2. The Seaview software package has been modified to include bistatic configurations and inverted to provide current, wind direction and wave measurements.

Currents have been obtained with good accuracy and consistency over many different bistatic and Bragg angles as evidenced in the scatter plots, Figure 22, and maps.

Wind directions are consistent with modelling except for the case shown in Figure 17. Note that the colour-coding in the wind plots is the derived directional spreading of the short waves which we haven’t attempted to validate at this point. In Figure 17 this goes beyond the expected maximum ($80°$) for this parameter. The reason is that the simulation used a wind direction that was roughly aligned with the Bragg direction and the smaller Bragg peak was mostly lost in the simulated noise level. The Seaview algorithm has difficulty estimating a wind direction accurately in these circumstances which, in our experience, rarely occur in measured data. It is interesting to note that waves are still measurable in these conditions confirming that the inversion result is independent of the initial guess which uses the wind direction.

The wave inversions are not as uniform as those for currents and winds. The significant waveheights are in reasonable agreement, Figure 23, although there is some evidence of overestimation at the highest simulated waveheight. That figure also shows that waveheight is underestimated for the largest bistatic angles of $64°$. More work is needed to determine a bistatic angle threshold for accurate wave measurement. Case 5 shows particularly noisy peak wave directions, Figure 16, including at the selected cells for which directional spectra are shown in Figure 20. While the frequency spectra (top left in each case) show good agreement, the low frequency part of the mean directions as a function of frequency (bottom left in each case) are not good at the low frequency peaks for three of these cases particuarly for the bistatic case. The inversion seems to be oversensitive to noise at these frequencies, which correspond with Doppler frequencies near the first order peak, and this needs further work as has been noted in many other applications of this method ([34,39]). Although some locations, including cell 1664, have Bragg angle differences in the bistatic case that are below the suggested monostatic angle threshold, there is no evidence that the results are worse there. This aspect also needs further work.

The directional spectra in Figure 18, Figure 19, Figure 20 and Figure 21 use log scales for both the frequency spectra (top left) and the directional spectra (right column) so that differences at both high and low amplitudes can be identified. Apart from the low frequency issue referred to above, the shape of the spectra are in good agreement in all cases. As mentioned above, the maximum frequency for the radar measurements is variable and depends on geometrical factors such as Bragg direction relative to wind direction. At the frequency (12.355 MHz) of the examples presented her, the monostatic cases have maximum frequencies in the range of 0.22–0.283 Hz and the bistatic in the range of 0.187–0.277 Hz.

6. Conclusions

In this paper, the theory for the interpretation of bistatic radar Doppler spectra in terms of currents, winds and waves has been reviewed and methods to simulate and then invert such data have been developed. As far as the authors are aware this is the first time that ocean wave directional spectra, to a maximum frequency that depends on the geometrical parameters and without any prior assumptions about the shape of those spectra, have been obtained from bistatic, albeit only simulated, data. The next step will be to apply the method to measured radar data.

Current, wind and wave measurements from bistatic radar data have been obtained with reasonable accuracy. The statistics for the bistatic cases are not quite as good as the monostatic cases, although they are biased by the large bistatic angle cases. The exact limits on bistatic angle and angle between Braggs still need to be determined but are expected to be about $60°$ and $\approx 30°$ respectively.

The results in this paper compare a monostatic configuration with a combined monostatic and bistatic configuration with one transmitter at one of the receive sites. Work on a configuration involving two bistatic radars with one transmitter located between the two sites is in progress. This will help to determine suitable configurations for bistatic radar installations for oceanographic measurements.