*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 ±*ωB*, defined in Equation (25). The contribution comes from two wave vectors, ±*k* - *<sup>B</sup>*, 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 *kB* includes a factor of cos *ϕ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 *ϕbi* → 90° are referred to as *forward scatter* and when this occurs, both *ω<sup>B</sup>* and *kB* tend to 0. Consequently, the wavelengths of the Bragg waves in this region become infinitely long and in addition, the cos<sup>4</sup> *ϕ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 *k* -<sup>1</sup> and *k* -2, and the nonlinear interaction between the same two ocean waves. The wave vector pair sum to give the Bragg wavevector *k* - *<sup>B</sup>* 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 *k* -<sup>1</sup> and *k* -2, the solution to the delta function in Equation (28), such that

$$
\omega - m\sqrt{gk\_1 \tanh(k\_1 d)} - m'\sqrt{gk\_2 \tanh(k\_2 d)} = 0 \tag{36}
$$

is sought. For each value of *ω*, the set of solutions of *k* -<sup>1</sup> and *k* -<sup>2</sup> defines a frequency contour in the *p*, *q* plane. As *m* and *m* can both take the values of 1 or −1, Equation (36) has four different forms.

**Figure 6.** Geometry of the second order scattering wave vectors, *k* - <sup>1</sup> and *k* -2, at angles *θ*<sup>1</sup> and *θ*2, respectively.

#### Case *m* = *m* :

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*<sup>1</sup> + *k*<sup>2</sup> > 2*k*<sup>0</sup> cos *ϕbi*, and therefore

$$
\omega^2 > 2\lg k\_0 \cos \varphi\_{bi} \tanh(2\lg k\_0 \cos \varphi\_{bi} d).
$$

Taking the square root gives

$$|\omega| > \sqrt{2\lg k\_0 \cos\varrho\_{bi} \tanh(2\lg k\_0 \cos\varrho\_{bi} d)}$$

which is equal to the bragg frequency *ω<sup>B</sup>* given in Equation (25). This leads to two conditions;

$$\begin{cases} \omega > \omega\_B & m = m' = 1 \\ \omega \le -\omega\_B & m = m' = -1 \end{cases}$$

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 *ω* increases, the contours become less circular, until |*ω*| = 2 *gk*<sup>0</sup> cos *<sup>ϕ</sup>bi* tanh(*k*0*<sup>d</sup>* cos *<sup>ϕ</sup>bi*), shown in white in Figure 7, where they separate.

**Figure 7.** The frequency contours of Equation (36) for two values of *ϕbi*, (**a**) monostatic *ϕbi* = 0 and (**b**) bistatic angle *ϕbi* = 25° when *m* = *m* = 1. The normalised frequency, *η* = *ω*/*ωB*, is shown by the colour, in the *p*, *q* plane.

When *ϕbi* = 0, the contours are symmetrical about the *p* and *q* axes, however, when *ϕbi* > 0, the contours rotate clockwise in the *p*, *q* plane, becoming symmetrical about some other axes, say *p* and *q* , shown by the additional black lines in Figure 7b.

Case *m* = *m* :

> When *m* = *m* , the square of Equation (36) is

$$
\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 \le \lg\left( (k\_2 - k\_1) \tanh((k\_2 - k\_1)d) \right) \dots
$$

In the right hand plane, when *k* -<sup>1</sup> and *k* -<sup>2</sup> lie in opposite directions along the *p* axis, meeting at a point past the bragg frequency, *k*<sup>2</sup> − *k*<sup>1</sup> reaches its maximum value of 2*k*<sup>0</sup> cos *ϕbi*. Therefore,

$$
\omega^2 \le \lg(2k\_0 \cos \,\upvarphi\_{bi} \tanh(2k\_0 \cos \,\upvarphi\_{bi} d))
$$

and this leads to the conditions

$$\begin{cases} 0 < \omega \le \omega\_B & m = -1, m' = 1 \\ -\omega\_B < \omega \le 0 & m = 1, m' = -1. \end{cases}$$

In the left hand plane, where *k*<sup>2</sup> < *k*1, the result is reversed giving

$$\begin{cases} 0 < \omega \le \omega\_B & m = 1, m' = -1 \\ -\omega\_B < \omega \le 0 & m = -1, m' = 1. \end{cases}$$

Figure 8 shows the normalised frequency contours when *m* = *m* for different values of *ϕbi*. Like when *m* = *m* , the contours are symmetric about the *p* and *q* axes, however, in this case, they do not cross the *q* 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.

**Figure 8.** The frequency contours of Equation (36) shown for *m* = *m* (where *m* = 1) with two different values for *ϕbi*: (a) monostatic *ϕbi* = 0 and (b) bistatic angle *ϕbi* = 25°. The colour shows the value of the normalised frequency, *η* = *ω*/*ωB*, in the *p*, *q* plane.

As the frequency contours for all four possible combinations of *m* and *m* are symmetrical about the *q* 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* plane, and double the result and hence

$$\sigma\_2(\omega) = 2^6 \pi k\_0^4 \cos^4 \varphi\_{bi} \sum\_{m, m' = \pm 1} \int\_{-\infty}^{\infty} \int\_{q'}^{\infty} |\Gamma\_E - i\Gamma\_H|^2 \mathcal{S}(m\vec{k\_1}) \mathcal{S}(m'\vec{k\_2}) \delta(\omega - m\omega\_1 - m'\omega\_2) \, dp \, dq. \tag{37}$$

In the right hand *p* plane, where *k*<sup>1</sup> ≤ *k*2, we can calculate the integral in polar coordinates *k*<sup>1</sup> and *θ*1, where *θ*<sup>1</sup> is the angle between *k* -<sup>1</sup> and the *p* axis, as shown Figure 6. Explicitly,

$$\sigma\_2(\omega) = 2^6 \pi k\_0^4 \cos^4 \eta\_{\rm li} \sum\_{m, n' = \pm 1} \int\_{-\theta\_L^{-}}^{\theta\_L^{+}} \int\_0^\infty |\Gamma\_E - i\Gamma\_H|^2 S(m\vec{k}\_1) S(m'\vec{k}\_2) \delta(\omega - m\omega\_1 - m'\omega\_2) k\_1 \, dk\_1 \, d\theta\_1,\tag{38}$$

where *θ*− *<sup>L</sup>* and *<sup>θ</sup>*<sup>+</sup> *<sup>L</sup>* are the integration limits of *θ*<sup>1</sup> and vary with frequency as well as *ϕbi*. In general, we find that the integration limits are

$$
\theta\_L^- = -(\pi + \varphi\_{bi}) \quad \text{and} \quad \theta\_L^+ = \pi - \varphi\_{bi}. \tag{39}
$$

however a particular set of *ω* values require different limits. This happens when *m* = *m* and |*ω*| > 2 *<sup>k</sup>*0*<sup>g</sup>* cos *<sup>ϕ</sup>bi* tanh(*k*0*<sup>d</sup>* cos *<sup>ϕ</sup>bi*), as the frequency contours cross the *<sup>q</sup>* axis and no longer complete full rotations in the right hand *p* plane. By symmetry, when the contours cross the *q* axis, *k*<sup>1</sup> and *k*<sup>2</sup> are the same length. Therefore, when *k*<sup>2</sup> = *k*<sup>1</sup> and *m* = *m* , the delta constraint of Equation (36) becomes

$$
\omega = \pm 2\sqrt{\text{g}k\_1 \tanh(k\_1 d)}\tag{40}
$$

and then squaring Equation (40) gives

$$\frac{\omega^2}{4g} = k\_1 \tanh(k\_1 d),\tag{41}$$

which can be solved numerically for *k*1.

Introducing a term, *β*, as the angle between *k* -<sup>1</sup> and *k* - *<sup>B</sup>* (see Figure 6), the integration limits are given by

$$
\theta\_L^+ = \pi - \varrho\_{bi} - \beta \quad \text{and} \quad \theta\_L^- = -\left(2\pi - \theta\_L^+ - 2\beta\right). \tag{42}
$$

Therefore, as *β* = cos−<sup>1</sup> *kB* 2*k*<sup>1</sup> , where the solution for *k*<sup>1</sup> from Equation (41) is used,

$$
\theta\_L^+ = \pi - \varphi\_{hi} - \cos^{-1}\left(\frac{k\_B}{2k\_1}\right) \quad \text{and} \quad \theta\_L^- = -\left(\pi + \varphi\_{hi} - \cos^{-1}\left(\frac{k\_B}{2k\_1}\right)\right).
\tag{43}
$$

In terms of *k*<sup>1</sup> and *θ*1,

$$k\_2 = \sqrt{k\_1^2 + k\_B^2 + 2k\_B k\_1 \cos(\theta\_1 + q\_{bi})} \quad \text{and} \quad \theta\_2 = \pi + \theta\_1 - \cos^{-1}\left(\frac{k\_1 + k\_B \cos(\theta\_1 + q\_{bi})}{k\_2}\right).$$

To calculate *σ*(2)(*ω*), we now reduce the double integral in Equation (38) to a single integral using the delta function. By defining

$$\begin{aligned} y\_s &= \sqrt{k\_1} \\ h(y\_s, \theta\_1) &= my\_s \sqrt{g \tanh(y\_s^2 d)} + m' \sqrt{gk\_2 \tanh(k\_2 d)} \\ I(y\_s, \theta\_1) &= 2^7 \pi |\Gamma\_T|^2 k\_0^4 \cos^4 \eta\_{\text{bi}} S(m\vec{k\_1}) S(m'\vec{k\_2}) y\_{s'}^3 \end{aligned}$$

Equation (38) can be written as

$$
\sigma\_2(\omega) = \int\_{-\theta\_L^-}^{\theta\_L^+} \int\_0^\infty I(y\_{s\prime}\theta\_1) \delta\left(\omega - h(y\_{s\prime}\theta\_1)\right) \, dy\_s \, d\theta\_1. \tag{44}
$$

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

$$
\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) \left| \frac{\partial y\_s}{\partial h} \right|\_\theta dh \, d\theta\_1,\tag{45}
$$

*Remote Sens.* **2020**, *12*, 313

where

$$\begin{split} \left| \frac{\partial h}{\partial y\_{s}} \right|\_{\theta\_{1}} &= \sqrt{\mathfrak{F}} \left\{ m \left( \sqrt{\tanh(y\_{s}^{2}d)} + \frac{y\_{s}^{2}d(\mathrm{sech}^{2}(y\_{s}^{2}d))}{\sqrt{\tanh(y\_{s}^{2}d)}} \right) \right. \\ &\left. + \frac{m'(y\_{s}^{3} + y\_{s}k\_{B}\cos(\theta\_{1} + \varphi\_{\mathrm{bi}}))}{k\_{2}^{3/2}} \left\{ \sqrt{\tanh(k\_{2}d)} + k\_{2}d \frac{\mathrm{sech}^{2}(k\_{2}d)}{\sqrt{\tanh(k\_{2}d)}} \right\} \right\}, \end{split}$$

whose reciprocal is *∂ys ∂h θ*1 . To integrate over the delta function, the solution, *y*∗, to

$$
\omega - h(y^\*, \theta\_1) = 0 \tag{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 *mm* = 1 and *θ*<sup>1</sup> = −*ϕbi*, the solution of *f*(*y*) is

$$y\_0^\* = \frac{\omega^2 - gk\_B}{2m\sqrt{\mathcal{g}}\omega}.\tag{47}$$

• When *mm* = −1 and *θ*<sup>1</sup> = *π* − *ϕbi*, the solution is

$$y\_0^\* = \frac{m\omega + \sqrt{2gk\_B - \omega^2}}{2\sqrt{\mathcal{g}}}.\tag{48}$$

Upon finding *y*∗, the second order cross section calculation in Equation (45) reduces to

$$
\sigma\_2(\omega) = \int\_{-\theta\_L^-}^{\theta\_L^+} I(y\_{s\prime}\theta\_1) \left| \frac{\partial y\_s}{\partial h} \right|\_{\theta\_1} \Big|\_{y\_s = y^\*} d\theta\_{1\prime} \tag{49}
$$

which can be calculated using a numerical integration method. For speed and convergence we update the value of *y*∗ <sup>0</sup> to the previously found solution for *y*∗, as *θ*<sup>1</sup> incrementally increases.
