# Bistatic SAR Mapping of Ocean-Wave Spectra

Marcel Kleinherenbrink, Paco Lopez-Dekker, Frederic Nouguier, Bertrand Chapron

► To cite this version:

Marcel Kleinherenbrink, Paco Lopez-Dekker, Frederic Nouguier, Bertrand Chapron. Bistatic SAR Mapping of Ocean-Wave Spectra. IEEE Transactions on Geoscience and Remote Sensing, 2024, 62, 4205812 (12p.). 10.1109/TGRS.2024.3394245. hal-04669256

HAL Id: hal-04669256

https://hal.science/hal-04669256v1

Submitted on 25 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Bistatic SAR Mapping of Ocean-Wave Spectra

Kleinherenbrink Marcel <fnref n="1" />,*, Lopez-Dekker Paco <fnref n="1" />, Nouguier Frederic <fnref n="2" />, Chapron Bertrand <fnref n="2" />

**Abstract :**

Earth Explorer 10 mission Harmony will consist of two satellites that fly in formation with Sentinel-1. It will operate as a multistatic radar in which Sentinel-1 transmits signals and all three satellites receive signals from different lines-of-sight. To prepare for Harmony and other possible future bistatic missions, transforms are derived to map the ocean-wave spectrum into bistatic synthetic aperture radar (SAR) spectra. The SAR mapping follows the standard derivation using the multidimensional characteristic function, but with adjustments for the modulation transfer functions compared to the monostatic case. This article focuses on the SAR modulations caused by velocity bunching as it is the dominant distortion mechanism. We argue that a multistatic system, such as Harmony, leads to an inversion that constrains the real aperture radar (RAR) response on a scene-by-scene basis. A benefit of having additional receivers for wave spectra estimation is that the three lines-of-sight enable to capture a larger fraction of the wave spectrum. Improvements are especially expected in high wind speed conditions such as tropical cyclones, where large energetic surface motions strongly deteriorate the (azimuth) resolution of the SAR data. Enhanced directional wave spectral characteristics will further help to improve the interpretation of the new bistatic Harmony high-resolution scatter and Doppler combined directional measurements.

**Keywords** : Synthetic aperture radar, Transfer functions, Receivers, Modulation, Transmitters, Spaceborne radar, Sea surface, Bistatic synthetic aperture radar (SAR), cross-spectral analysis, Harmony, ocean-wave spectra, SAR spectra

<footnote n="1">¹ Department of Geoscience and Remote Sensing, Delft University of Technology, Delft, The Netherlands</footnote>

<footnote n="2">² Laboratoire d'Océanographie Physique et Spatiale, CNRS, Institut Francais de Recherche pour l'Exploitation de la Mer, IRD, Université de Brest, Plouzané, France</footnote>

<footnote>* Corresponding author : Marcel Kleinherenbrink, email address : m.kleinherenbrink@Tudelft.nl</footnote>
# I. INTRODUCTION

Synthetic aperture radar (SAR) processing assumes the received echoes result from reflections of the transmitted signals from static scatterers. Non-zero motions in the range direction of scatterers lead to shifts in the azimuth direction and defocusing [1]. During the seventies, airborne SAR imagery over the ocean revealed the complexity of radar signals reflecting from waves of various scales. Around the same time, the concept of velocity bunching was introduced, which is the azimuthal variation of scatter density in a SAR image as a consequence of the varying orbital motions within long waves [2]. With the launch of Seasat in 1978, the interest in SAR signals over the ocean increased, as it was the first opportunity to study ocean waves on a large scale. Alpers & Rufenach [3] then theoretically described the velocity bunching phenomenon and the loss of resolution due to surface accelerations.

Further developments allowed the interpretation SAR image modulations via transfer functions directly with the purpose of inferring ocean-wave properties. These transfer functions incorporated harmonic modulations caused by tilt, hydrodynamic effects, and velocity bunching [4], [5]. Hasselmann & Hasselmann [6] then developed a non-linear mapping function,

which directly relates a directional ocean wave spectrum to a SAR image intensity spectrum. Generalizations of the theory were made by others [7]–[10], leading to the definition of practical algorithms [11] to retrieve estimates of the underlying directional wave spectra from the observed intensity spectra.

In recent years, there has been an increased interest in bistatic spaceborne SAR missions, in which the transmitters and receivers fly on-board different platforms. Recently, Earth Explorer 10 mission Harmony got approved, consisting of two passive bistatic companions, which will fly in formation with Sentinel-1D [12]–[14]. Flying multiple platforms with receivers in formation with a monostatic system is an efficient way to provide simultaneous observations with multiple lines of sight, enhancing the ability to retrieve directional information. The theoretical mapping of wave modulations in bistatic SAR images has only been studied to a limited extent. Ying et al. [15] derived a velocity bunching transfer function and elaborated the theory of shifts and resolution degradation comparable to [3]. Transfer functions for the bistatic hydrodynamic and tilt modulations were not considered. A full non-linear bistatic SAR mapping function has therefore never been derived. Since the sensitivity of a SAR system to the different spectral components in an ocean surface wave field is highly dependent on the direction, we can expect that a multi-directional system will help recover larger parts of the directional surface wave spectrum.

In the present study, the bistatic SAR spectral mapping transformation is presented for systems consisting of one transmitter and multiple receivers with along-track separations up to $O(500 \text{ km})$. Starting from the use of the characteristic function [7], [8], [10], [11], an extension for the bistatic case is derived. Using simple parametric wind-sea and swell wave spectra, properties of the bistatic SAR spectrum are analysed. Having multiple bistatic companions, like Harmony, is demonstrated to help recover a larger fraction of the spectrum compared to a monostatic system. Lastly, we compare our closed-form solution to end-to-end numerical simulations.

# II. THEORY

This section addresses the SAR mapping of an ocean-surface spectrum for a bistatic radar system like Earth Explorer 10 Harmony. First, the Doppler shifts, consequences of motions in the range direction, are discussed. Then, transfer functions are derived to infer cross-track shifts, azimuth shifts and scatter variations. Finally, we reintroduce the full non-linear mapping between ocean-wave and SAR spectra.

## A. Doppler shift vectors

Under a static-surface assumption, the received signal phase history is fully determined by the evolving distance between
the receiving radar and the target. Surface motions during a SAR acquisition causes this phase, or range, history to deviate from its expected one. Typically described by a second-order function, multiple deviations are quadratic to cause smearing effects and a loss of resolution, while linear terms will cause shifts in the SAR image.

Let us considering the range-history of a target, i.e. the time-varying distance between a point on the surface and the moving radar. For a monostatic radar with the radar beam oriented in the zero-Doppler direction, this range history is well approximated by

$$R(t) \approx R_t + \frac{U^2}{2R_t} t^2, \qquad (1)$$

with $R_t$ the shortest range which happens at the zero-Doppler time, $t = 0$, and $U$ the effective velocity of the radar (which can account for a curved trajectory of a satellite). A radial velocity $v_t$ (in the direction of the transmitter) of the target will lead to a modified range-history

$$R'(t) = R(t) - v_t \cdot t. \qquad (2)$$

Solving for the time of shortest range, i.e. the zero-Doppler time ($\frac{1}{dt}R'(\Delta t) = 0$), gives

$$\Delta t = \frac{v_t R_t}{U^2}, \qquad (3)$$

which is the time-offset at which the target is imaged with respect to its true position. This time offset translates to a spatial offset in the flight direction given by

$$d_y \approx \frac{v_t R_t v_g}{U^2}, \qquad (4)$$

with $v_g$ the ground projection of the velocity of the radar. In a flat-Earth geometry we have $v_g = U$ and (4) simplifies to [1], [3]

$$d_y \approx \frac{v_t R_t}{U}. \qquad (5)$$

A more complex situation arises for bistatic observations. Fig. 1 shows the observation geometry for one of Harmony's receivers. The non-squinted transmitter (Sentinel-1) transmits a signal which hits the surface at incident angle $\theta_t$. The transmitter ground-range direction is aligned with the cross-track direction $x$. Part of the signal gets reflected towards the receiver at an angle $\theta_r$ with respect to the normal of the surface and bistatic angle $\alpha_s$. In the equations below we use the ground-projected bistatic angle $\alpha$, which is the azimuth angle between the transmitter and receiver ground-projected range directions. The non-symmetric bistatic geometry around the x-axis leads to a rotation of the ground range, which is proportional to the bistatic angle. For the same reason, a bistatic geometry also exhibits a small rotation of the Doppler gradient. As the geometry changes within the Sentinel-1 - Harmony swath, the range and Doppler isolines are slightly curving. However, radial surface velocities are rather small which result in displacements of O(100 m). We therefore assume the Doppler shifts are linear. A more general solution is provided in the appendix.

A shift in a SAR image displaces the scatterer along the iso-range line so that the resulting Doppler shift compensates

<figure><img src="image_2.png" /><figcaption>Fig. 1. Observing geometry of the considered bistatic system. The transmitted and reflected signals are indicated with orange and red lines. The blue lines are related to the bistatic ground-range direction and the purple line indicates the direction of maximum Doppler gradient.</figcaption></figure>

the Doppler due to the velocity of the target. Therefore we approximate

$$\frac{\delta}{\delta x} \dot{R}_b d_x + \frac{\delta}{\delta y} \dot{R}_b d_y = v_b, \qquad (6)$$

where the right side represents the spatial changes in the bistatic Doppler velocity $\dot{R}_b$ as a consequence of a bistatic line-of-sight velocity $v_b$. The velocity components split between transmitter and receiver components, to yield

$$\frac{\delta}{\delta x}(\dot{R}_t + \dot{R}_r)d_x + \frac{\delta}{\delta y}(\dot{R}_t + \dot{R}_r)d_y = v_b. \qquad (7)$$

Under the assumption that the platform velocity vectors $\vec{U}_t$ and $\vec{U}_r$ are constant over the integration time, this is rewritten as

$$(\vec{U}_t \cdot \frac{\delta \hat{r}_t}{\delta x} + \vec{U}_r \cdot \frac{\delta \hat{r}_r}{\delta x})d_x + (\vec{U}_t \cdot \frac{\delta \hat{r}_t}{\delta y} + \vec{U}_r \cdot \frac{\delta \hat{r}_r}{\delta y})d_y = v_b. \qquad (8)$$

The equation shows that the Doppler difference is related to a change in the line-of-sight unit vectors $\hat{r}_t$ and $\hat{r}_r$ from the transmitter and receiver, respectively. The partial derivatives are given as

$$\frac{\delta \hat{r}_i}{\delta x} = \frac{\hat{x} - (\hat{r}_i \cdot \hat{x}) \hat{r}_i}{r_i} \qquad (9) \\ \frac{\delta \hat{r}_i}{\delta y} = \frac{\hat{y} - (\hat{r}_i \cdot \hat{y}) \hat{r}_i}{r_i}.$$

By inserting the partial derivative in the previous equation, the four terms within brackets represent angular velocity changes of the two platforms in both directions, leading to

$$(\frac{U_{xt}}{R_t} + \frac{U_{xr}}{R_r})d_x + (\frac{U_{yt}}{R_t} + \frac{U_{yr}}{R_r})d_y = v_b, \qquad (10)$$

a function of projected platform velocities ($U_{xt}, U_{xr}, U_{yt}, U_{yr}$), resembling to the monostatic configuration.

The function is related to displacements in both along- and across-track directions. Note that a moving scatterer should not only carry the same range history, but should also have
approximately the same reference range as the target location. It therefore approximately experiences a shift over a range isoline. The spatial derivatives of the bistatic range are given as

$$ 
\begin{aligned}
\frac{\delta r_b}{\delta x} &= (\hat{r}_t + \hat{r}_r) \cdot \hat{x} \\
\frac{\delta r_b}{\delta y} &= (\hat{r}_t + \hat{r}_r) \cdot \hat{y}.
\end{aligned}
\qquad (11) $$

The combination of both functions leads to a dependency for both coordinates

$$ \frac{\delta y}{\delta x} = \frac{(\hat{r}_t + \hat{r}_r) \cdot \hat{x}}{(\hat{r}_t + \hat{r}_r) \cdot \hat{y}}, \qquad (12) $$

which reduces to

$$ \frac{\delta y}{\delta x} = \frac{\sin \theta_t + \cos \alpha \sin \theta_r}{\sin \alpha \sin \theta_r} \qquad (13) $$

for a bistatic system with a non-squinted transmitter and a squinted receiver (Fig. 1). By inserting the assumption of $d_y = \frac{\delta y}{\delta x} d_x$ into Eq. 10 two equations for the shifts

$$ 
\begin{aligned}
d_x &= \frac{v_b}{(\vec{U}_t \cdot \frac{\delta \hat{r}_t}{\delta x} + \vec{U}_r \cdot \frac{\delta \hat{r}_r}{\delta x}) + (\vec{U}_t \cdot \frac{\delta \hat{r}_t}{\delta y} + \vec{U}_r \cdot \frac{\delta \hat{r}_r}{\delta y}) \frac{\delta y}{\delta x}} \\
d_y &= \frac{v_b}{(\vec{U}_t \cdot \frac{\delta \hat{r}_t}{\delta x} + \vec{U}_r \cdot \frac{\delta \hat{r}_r}{\delta x}) \frac{\delta x}{\delta y} + (\vec{U}_t \cdot \frac{\delta \hat{r}_t}{\delta y} + \vec{U}_r \cdot \frac{\delta \hat{r}_r}{\delta y})}
\end{aligned}
\qquad (14) $$

are obtained. For the monostatic case $\vec{U}_t = \vec{U}_r$, $\vec{U}_t \frac{\delta \hat{r}_t}{\delta x} = \vec{U}_r \frac{\delta \hat{r}_r}{\delta x} = 0$ and $\frac{\delta y}{\delta x} = \infty$, which solely leads to a shift in the along-track direction. The velocity is perpendicular to the range direction, and

$$ d_y = \frac{v_b}{\vec{U}_t \cdot \frac{\delta \hat{r}_t}{\delta y} + \vec{U}_r \cdot \frac{\delta \hat{r}_r}{\delta y}} = \frac{v_b R_t}{2U_t}, \qquad (15) $$

is the monostatic azimuth shift as consequence of velocity $v_b = 2v_t$. This expression is similar to Eq. 5.

## B. Transfer functions

We consider transfer functions of the non-linear SAR modulation mapping based on the Doppler-shift vectors, and separately, a linear RAR modulation based on the derivatives of bistatic Normalized Radar Cross Section (NRCS) with respect to the incident angle for a rough surface.

1) **SAR modulation:** An elegant way of computing the SAR spectral transform is to consider the use of cross-correlation functions. Transfer functions, used to compute the shift and scatter modulations, are expressed in the spectral domain. For shifts, it simply requires the computation of the bistatic range rate $v_b$ for each wave number $(k_x, k_y)$ with amplitude $\zeta_k$, to enter in Eq. 14. The bistatic range rate is the sum of the velocities in the transmitter and receiver line-of-sights, but a cleaner solution is obtained by projecting the surface velocity $\vec{v}_s$ in the bistatic line-of-sight

$$ v_b = v_t + v_r = \vec{v}_s \cdot (\hat{r}_t + \hat{r}_r), \qquad (16) $$

which can be rewritten as a function of the angles of the bistatic system

$$ v_b = \begin{pmatrix} v_{s,x} \\ v_{s,y} \\ v_{s,z} \end{pmatrix} \begin{pmatrix} \sin \theta_t + \sin \theta_r \cos \alpha \\ \sin \theta_r \sin \alpha \\ \cos \theta_t + \cos \theta_r \end{pmatrix}. \qquad (17) $$

Assume the wave height at time $t = 0$

$$ \zeta(x, y) = \sum_k \zeta_k e^{i(k_x x + k_y y)}, \qquad (18) $$

with $\omega = \sqrt{gk}$ for deep ocean waves, then the velocity for a particular $(k_x, k_y)$ are readily computed as

$$ 
\begin{aligned}
v_b &= -\omega \left( \frac{k_x}{|k|} (\sin \theta_t + \sin \theta_r \cos \alpha) - \frac{k_y}{|k|} \sin \theta_r \sin \alpha + i(\cos \theta_t + \cos \theta_r) \right) \zeta_k e^{i(k_x x + k_y y)} \\
&= (-\omega'_x + \omega'_y - j\omega'_z) \zeta_k e^{i(k_x x + k_y y)},
\end{aligned}
\qquad (19) $$

where $(\omega'_x, \omega'_y, \omega'_z)$ are line-of-sight-projected components of the wave angular velocity. Entering these results into Eq. 14, leads to the transfer functions

$$ 
\begin{aligned}
T_x &= \frac{-\omega'_x + \omega'_y - j\omega'_z}{\Omega_x + \Omega_y \frac{\delta y}{\delta x}} \\
T_y &= \frac{-\omega'_x + \omega'_y - j\omega'_z}{\Omega_x \frac{\delta x}{\delta y} + \Omega_y}
\end{aligned}
\qquad (20) $$

with

$$ 
\begin{aligned}
\Omega_x &= \vec{U}_t \cdot \frac{\delta \hat{r}_t}{\delta x} + \vec{U}_r \cdot \frac{\delta \hat{r}_r}{\delta x} \\
\Omega_y &= \vec{U}_t \cdot \frac{\delta \hat{r}_t}{\delta y} + \vec{U}_r \cdot \frac{\delta \hat{r}_r}{\delta y}
\end{aligned}
\qquad (21) $$

so that the transfer functions depend on the platform velocity vectors in the local reference frame, the incident angles and the bistatic angle. Let $\alpha = 0$ and $\theta_r = \theta_t$ and using the same simplification in the denominator as in Eq. 15, we obtain

$$ 
\begin{aligned}
T_y &= -\frac{R_t \omega}{2U_t} \left( \frac{k_x}{|k|} 2 \sin \theta_t + i2 \cos \theta_t \right) \\
&= -\frac{R_t \omega}{U_t} \left( \frac{k_x}{|k|} \sin \theta_t + i \cos \theta_t \right).
\end{aligned}
\qquad (22) $$

the monostatic transfer functions [4], [6], [9] for azimuth shifts caused by surface velocity.

2) **RAR modulation:** Although small, the RAR is dominant in the ground-range direction, where velocity bunching is absent. A mapping of the SAR spectrum therefore requires a description of the RAR modulations. Note that the RAR for moderate incident angles (20° - 50°) is not always well modeled and constrained. Most inversion algorithms consider a RAR variant of two branches. The first branch inherently assumes a two-scale model, which considers tilting and hydrodynamic interaction of Bragg waves by larger waves in an analytical expression [10], [16]. The implementation of this model is relatively simple, but the applicability is very limited as the assumptions break when Bragg scattering is not dominant. The second branch is based on derivatives of geophysical model functions used for wind-speed estimation from scatterometry or SAR data [11], [17]. This second type of model is also used in the retrieval algorithm of Sentinel-1 ocean-wave spectra [18]. An advantage over the first branch is that the Sentinel-1 algorithm also accounts for polarization.

Since no long-baseline bistatic radar mission is currently flying and the RAR modulations for monostatic radar missions is poorly constraint, it is difficult to define a RAR transfer
function which accurately includes polarimetric rotations for a bistatic system [19]. We will therefore adopt a RAR modulation model of the first branch based on the incident-angle derivatives of the scattering functions described in Elfouhaily et al. (1999/2001) [20], [21]. Let us consider the description for the surface elevation in Eq. 18. For such a surface, by treating the system as a monostatic equivalent, a slope modulation transfer function is constructed as

$$T_I = -ik_b \frac{1}{\sigma_{pq}} \frac{\delta\sigma_{pq}}{\delta\theta_b}, \quad (23)$$

which depends on the NRCS $\sigma_{pq}$ in polarization pq, the derivative of the scatter with respect to the bistatic incident angle $\frac{\delta\sigma_{pq}}{\delta\theta_b}$ and the range wavenumber $k_b$. The range wavenumber is computed from the ground-range gradient in Eq. 13, such that

$$k_b = k_x \cos(\phi_r) - k_y \sin(\phi_r) \quad (24)$$

with

$$\phi_r = \arctan\left(\frac{\delta y}{\delta x}\right), \quad (25)$$

where $\frac{\delta y}{\delta x}$ as in Eq. 13. At the Small Perturbation Model (SPM) limit, the scatter is computed with radar wavenumber $k_r$, such that [21]

$$\sigma_{pq} = 16\pi k_r^2 g_{pq} (S(k_{Br}, \phi_r) + S(k_{Br}, \phi_r + \pi)), \quad (26)$$

where $S(k_{Br}, \phi_r)$ is the Elfouhaily spectrum evaluated at resonant Bragg wavenumber $k_{Br} = 2k_r \frac{|\vec{r}_i + \vec{r}_r|}{2} \sin(\theta_b)$ in the ground-range direction $\phi_r$. Further details on the spectrum are given in the appendix. The scattering parameter $g_{pq}$ is a function of the vertical components of the radar wavenumber $q_z = k_r(\cos(\theta_t) + \cos(\theta_r))$ and is computed as

$$g_{pq} = q_z^2 |\vec{P}_i \cdot \vec{P}_s|^2, \quad (27)$$

where $\vec{P}_i$ and $\vec{P}_s$ are the polarization vectors of the incoming and scattered electromagnetic fields as defined in Eq. 39 in [20]. Note that the derivative $\frac{\delta\sigma_{pq}}{\delta\theta_b}$ changes the Bragg wavenumber and requires for the bistatic case derivatives of the polarization vectors. A tilt, or local incident angle change in the line-of-sight $\delta\theta_b$ introduces $\delta\theta_t, \delta\theta_r$ and $\delta\alpha$, which modulates the polarization of the scattered signal. As this paper focuses on the bistatic SAR transform, an in-depth analysis of the RAR is omitted. Therefore we simplify the polarization vector product $|\vec{P}_i \cdot \vec{P}_s|$ to an equivalent monostatic case, which only depends on $\theta_b$. This has the consequence that we effectively assume VV polarized modulations, which is likely an underestimation for the bistatic case.

The RAR modulation is an ongoing topic of research. Under homogeneous wind conditions, with wind speeds below 15 m s⁻¹ the models perform reasonably well. They, however, soon break in situations with substantial wave breaking and irregular or asymmetric wave spectra, like in tropical storms. For these cases, it is difficult to generalize the RAR modulation into a model. Fortunately, in most of these scenarios, velocity bunching is the dominant modulating mechanism [17], but still a proper RAR modulation is required to accurately constrain the wave spectrum in the range direction. A multistatic system, like Harmony, will help to constrain the RAR model.

By incoherent processing of the bursts of pulses, we get a waveform that is only sensitive to range-direction waves. With two bistatic companions, the wave spectrum can be 'sliced' in three directions, which provides us with a RAR response in these three directions. Under the assumption of an initial wave spectrum that fits the (a-)symmetry a first-guess RAR modulation function can be estimated. An iterative procedure is then used to solve both the RAR transfer function and the wave spectrum, a comparable approach to Engen et al. [22] or Jacobsen & Høgda [23]. Further enhancements are expected by using Harmony's polarimetry.

### C. SAR spectral transform

Krogstad et al. [8] elegantly expressed the nonlinear SAR-spectral mapping in terms of the partial derivatives of the joint-characteristic function of radar modulations and shifts, i.e. the G-function. Ignoring a delta function at the origin that expresses the mean NRCS, the SAR image intensity spectrum is then described by

$$P(k_x, k_y) = \iint G(x, y, k_x, k_y) \exp(-i(k_x x + k_y y)) dxdy, \quad (28)$$

which is slightly modified compared to Eq. 6 in [10], where the integral should be computed for each $(k_x, k_y)$ separately. Assuming Gaussian statistics, the G-function is described with a set of cross-correlations as

$$G(k_x, k_y, x, y) = \exp(k_x^2 \mu_{xx} + k_y^2 \mu_{yy} + k_x k_y (\mu_{xy} + \mu_{yx})) \cdot \\ (1 + \rho_{II} + \\ ik_x(\rho_{Ix} - \rho_{xI}) + ik_y(\rho_{Iy} - \rho_{yI}) + \\ ik_x^2 \mu_{Ix} \mu_{xI} + ik_y^2 \mu_{Iy} \mu_{yI} + \\ ik_x k_y (\mu_{Ix} \mu_{yI} + \mu_{Iy} \mu_{xI})), \quad (29)$$

with

$$\mu_{ab}(x, y) = \rho_{ab}(x, y) - \rho_{ab}(0, 0) \quad (30)$$

and the cross-correlations computed as

$$\rho_{ab}(x, y) = \frac{1}{(2\pi)^2} \int \int \frac{1}{2} T_a(k_x, k_y) T_b^*(k_x, k_y) e^{-i\omega\Delta t} S(k_x, k_y) + \\ \frac{1}{2} (T_a(-k_x, -k_y) T_b^*(-k_x, -k_y) e^{-i\omega\Delta t})^* S(-k_x, -k_y) \\ \exp(-i(k_x x + k_y y)) dk_x dk_y. \quad (31)$$

The cross-correlations depend on the wave spectrum $S$ and transfer functions $T_{a,b}$ defined in the previous sections. In case of a co-spectrum the time difference $\Delta t = 0$, but in practice cross-spectra are used to suppress the speckle-noise bias [24] and to resolve the 180-degree ambiguity [10].

The computation of the full SAR spectrum involves calculating Eq. 31 for the different items and all wavenumber and then applying Eq. 28, which is computationally expensive. Therefore, we follow [6] and derive a series expansion. The cross-terms in $G$ (Eq. 29) outside of the exponential have in most cases a small effect on the SAR spectra [9]. Furthermore,
the exponential terms related to the cross-correlation at lag zero are taken outside the integral, which yields

$$ P(k_x, k_y) \approx \exp(-k_x^2 \rho_{xx}(0) - k_y^2 \rho_{yy}(0) - 2k_x k_y \rho_{xy}(0)) \cdot \\ \iint e^{k_x^2 \rho_{xx} + k_y^2 \rho_{yy} + k_x k_y (\rho_{xy} + \rho_{xy})} (1 + \rho_{II}) \exp(-i(k_x x + k_y y)) dxdy. \quad (32) $$

The exponential term outside the integral can be interpreted as a Gaussian-shaped low-pass filter in the two-dimensional wavenumber domain, where the width of the filter decreases if the velocity variance of the surface increases. The exponential terms within the integral represent the non-linear shifting as a consequence of a moving surface. Each of the exponential terms inside of the integral can be written as a Taylor series, resulting in

$$ P(k_x, k_y) = \exp(-k_x^2 \rho_{xx}(0) - k_y^2 \rho_{yy}(0) - 2k_x k_y \rho_{xy}(0)) \cdot \\ \iint \left( \frac{(k_x^2 \rho_{xx})^0}{0!} + \frac{(k_x^2 \rho_{xx})^1}{1!} + \dots \right) \cdot \\ \left( \frac{(k_y^2 \rho_{yy})^0}{0!} + \frac{(k_y^2 \rho_{yy})^1}{1!} + \dots \right) \cdot \\ \left( \frac{(k_x k_y (\rho_{xy} + \rho_{yx}))^0}{0!} + \frac{(k_x k_y (\rho_{xy} + \rho_{yx}))^1}{1!} + \dots \right) \cdot \\ (1 + \rho_{II}) \exp(-i(k_x x + k_y y)) dxdy, \quad (33) $$

or

$$ P(k_x, k_y) = \exp(-k_x^2 \rho_{xx}(0) - k_y^2 \rho_{yy}(0) - 2k_x k_y \rho_{xy}(0)) \cdot \\ \iint \sum_a \frac{(k_x^2 \rho_{xx})^a}{a!} \cdot \sum_b \frac{(k_y^2 \rho_{yy})^b}{b!} \cdot \\ \sum_c \frac{(k_x k_y (\rho_{xy} + \rho_{yx}))^c}{c!} (1 + \rho_{II}) \cdot \\ \exp(-i(k_x x + k_y y)) dxdy, \quad (34) $$

which eventually can be rearranged to leave a function in the form of a set of Fourier transforms, such that

$$ P(k_x, k_y) = \exp(-k_x^2 \rho_{xx}(0) - k_y^2 \rho_{yy}(0) - 2k_x k_y \rho_{xy}(0)) \cdot \\ \sum_a \frac{k_x^{2a}}{a!} \sum_b \frac{k_y^{2b}}{b!} \sum_c \frac{(k_x k_y)^c}{c!} \cdot \\ \iint \rho_{xx}^a \rho_{yy}^b (\rho_{xy} + \rho_{yx})^c (1 + \rho_{II}) \cdot \\ \exp(-i(k_x x + k_y y)) dxdy, \quad (35) $$

which should be evaluated up to a certain order $a+b+c \le O$. Typically, a fifth-order expansion provides sufficiently accurate results.

### III. MULTISTATIC SAR SPECTRA

Harmony's multistatic SAR system consisting of one transmitter (Sentinel-1) and three receivers (Harmony-A, Sentinel-1 and Harmony-B) are considered in this section as a reference. The data of multistatic SAR systems will be used to estimate stress-equivalent wind vectors, surface-current vectors and

ocean-wave spectra. In this section, Harmony's multistatic SAR spectra are discussed, which are the main input for for ocean-wave retrievals. To analyse the properties of the multistatic spectra, we separately elaborate on the frozen-surface response (RAR), the cut-off, the effect of wind waves and the directional response to swell. We will argue that having three lines-of-sight not only enhances the retrieval of the ocean-wave spectrum, but also enhances the retrieval of surface-current vectors by means of constraining the wave-Doppler.

#### A. The frozen-surface response

A frozen surface does not cause any Doppler shifts, and its accompanied velocity bunching would therefore be absent. The response from such a surface would yield an RAR response at the resolution of a SAR system. In that sense, a 'two-dimensional RAR' observation is comparable to an optical instantaneous snapshot (e.g. Sentinel-2 as in Kudryavtsev et al. (2017) [25]), which is not affected by any misregistrations. Such a two-dimensional RAR response will never be observed due to the fast-moving nature of the ocean, but will contribute to the full SAR response. Presently demonstrated, although one-dimensional, with the SWIM instrument on-board the CFOSAT mission [26], [27], RAR spectra only account for the tilt-hydro modulation term $1+\rho_{II}$. Since the transfer functions to define $\rho_{II}$ are linear, the RAR spectra solely scale the input spectra.

An example of two-dimensional multistatic RAR spectra is shown Fig. 2. For this particular scenario, the peak of the Gaussian swell wave spectrum (see appendix for the implementation) propagates at an angle of -30° with respect to the cross-track direction. The swell system is superimposed by a wind-wave system [28] (see appendix for the implementation) propagating at a mean angle is 45°. The maximum intensity modulation is reached for waves propagating in the ground-range directions of the transmitter-receiver systems. Therefore, the swell-signal spectral density is the highest for Harmony-A (the heading receiver). The spectral density of the wind-wave system is highest for Harmony-B (the trailing receiver) as its ground range is more closely aligned with the peak direction. Due to the large wavenumber-directional spread of the wind-wave system, the peak RAR spectral density of a swell system is often higher than that of a wind-wave system. However, in case of equal spreading, a short-wavelength system would have a larger spectral density due to steeper wave slopes than a long-wavelength system.

In the ground-range direction ($k_y = 0$), Harmony's three lines-of-sight provide three slices of the wave spectrum not affected by velocity, but solely by RAR modulations. This enables the computation of three MeAn Cross Spectrum (MACS) [29], defined as the spectral integral between wavelengths of 15 m and 20 m near the range direction. The integration yields a complex number of which the phase is related to the motion of the modulation wind waves. Its amplitude depends on the resolved wind-wave amplitudes and slopes, and on the Bragg waves amplitude and direction. As such, the imaginary part of the MACS (iMACS) is strongly correlated with the wind
and the wave-Doppler. With the considered multistatic system it is possible to compute iMACS in three lines-of-sight for two polarizations. It is expected to be especially useful in constraining the wave-Doppler as it respresents comparable wavelengths and associated relaxation scales.

## B. Cut-off wavelength

Kerboal et al. [30] and Schulz-Stellenfleth et al. [16] clearly showed that the monostatic zero-lag correlation $\rho_{yy}(0) \approx (\frac{R_t}{U})^2 f_{v,m}(0)$ is a function of the ocean's velocity variance $f_{v,m}(0)$ projected in the line-of-sight. Velocity variance is the dominant term for azimuth resolution loss for sea surface coherence times larger than 0.05 s [16]. The azimuth resolution, described as the cut-off wavelength, is to the first order approximated as

$$ \lambda_{cut,m} \approx \pi \sqrt{(\frac{R_t}{U})^2 f_{v,m}(0)} \quad (36) $$

for a monostatic system. In the spectral domain, the cut-off is described by the exponential terms outside of the integral of Eq. 32, known as the fall-off [31]. Note that the fall-off function for the side-looking case at incident angles larger than 10° is one dimensional and only filters signals in the $k_y$-direction. The ocean velocity variance causes the SAR processing to filter signals from waves shorter than the cut-off in the along-track direction.

For a bistatic system the $k_x$-terms outside of the integral in Eq. 32 cannot be ingored. In fact, the fall-off function is rotated and symmetric around the ground-range direction. The equation for the bistatic cut-off is more complex, but note that the denominators of the transfer functions in Eq. 20 are purely geometric. As the geometry is known, the bistatic cut-off wavelength in the major polarizations (in which Bragg scattering is dominant) is still a function of the line-of-sight projected velocity variance $f_{v,b}(0)$ and has the form

$$ \lambda_{cut,b} \approx \pi \sqrt{f(\theta_t, \theta_r, \alpha, \vec{U}_t, \vec{U}_r) f_{v,b}(0)}, \quad (37) $$

where $f(\theta_t, \theta_r, \alpha, \vec{U}_t, \vec{U}_r)$ is a function depending on the satellite velocity vectors and the geometry of the bistatic observation system. The cut-off is sensitive to two geometrical parameters: 1) the gradient of the Doppler and 2) the angle between the range isolines and the Doppler gradient (Fig. 1). As the first is close to that of the monostatic, the difference between monostatic and bistatic cut-off size primarily results from the non-perpendicularity between Doppler and range isolines. This leads to a narrower fall-off function in the spectral domain (Fig. 3) for any slanted observation. Line-of-sight diversity therefore comes at the cost of a relatively lower spectral coverage for a single observation. A practical way of computing the cut-off is by a rotation of the spectrum such that the wavenumber $k_x$ is aligned with the ground-range direction and then follow the traditional procedure in the new basis.

The cut-off is often used to constrain the unresolved wind-wave spectrum. Using the NRCS, the cut-off and an external estimate of wind direction, it is possible to constrain the wind speed (and direction) and the inverse wave age of a

parametrized wave spectrum [28]. Note that accurate estimation of the inverse wave age is essential for constraining the wave-Doppler and therefore important for surface-current estimation. Constraining the wind-wave spectrum if often done in a serial approach, where the wind speed is estimated first, because the NRCS is assumed to be only weakly dependent on the wave age. Although small, the cut-off is theoretically modulated by the wind direction due to the horizontal components of the velocity variance (as shown by the transfer functions of Eq. 20). It is foreseen that retrieval algorithms for Harmony use the NRCS and cutoffs from three lines-of-sight in a joint inversion to estimate stress-equivalent wind-speed vectors and inverse wave age. This inversion might be aided by the MACS parameter, particularly in cases where short and long wind waves are propagating in different directions (current, wind variability), and observations of both polarizations, particularly when wave breaking has an important role (currents, shallow bathymetry, high sea states).

The three lines-of-sight also provide a direct benefit for the ocean-wave inversion by increasing the spectral coverage (Fig. 3). The along-track resolution of a monostatic system in the orbit of Sentinel-1 is in the order of 150 m. High sea states in polar lows and tropical cyclones can further deteriorate the along-track resolution to 300 m or more. This prevents monostatic SAR systems to map waves propagating off the ground-range direction. The directional diversity of the multistatic fall-offs in the spectral domain show that a relatively larger part of the spectrum is captured. The relative improvement in spectral coverage depends on the sea state and varies with wavelength as shown in Fig. 4. In most cases, Harmony's bistatic companions double the spectral coverage with respect to Sentinel-1 alone. For swell waves travelling nearly parallel to the flight direction, however not much benefit of having bistatic companions, as the cut-offs overlap near $k_x = 0$ in the spectral domain. However, the part of the wave spectrum at intermediate angles are resolved better. The coverage of the wind-wave systems, which often have wavelengths shorter than 100 m, remains low 10%, but note that wind-wave systems typically have a broad spectrum. The 'slicing' of such a spectrum in three direction will still provide valuable information for constraints.

## C. Wind waves

In the current ocean-processing algorithms of Sentinel-1 [32], a non-linear wind-wave component is removed before a linear swell-wave estimation. To illustrate the response of wind waves in the SAR spectrum, Fig. 5 shows the modeled monostatic and bistatic SAR spectra for three different wind speeds. The input wave spectrum is based on [28] and assumes nearly fully developed waves (see appendix for a summarized description). Note that with increasing wind speed, the total energy contained in the spectrum is higher, the peak energy becomes higher and the energy is more spread in terms of wavelength and direction. With increasing wind speed also the wavelength of the peak becomes larger.

Wind waves tend to reveal the non-linear nature of the SAR mapping of ocean-wave spectra. At wavenumbers where
<figure><img src="image_3.png" /><figcaption>Fig. 2. Example of two-dimensional multistatic RAR spectra from an ocean spectrum consisting of swell waves and wind waves. The baseline between Sentinel-1 and Harmony-A (bistatic heading) and Harmony-B (bistatic trailing) is 350 km. On the right the input wave spectrum is shown, consisting of a wind-wave and a swell system.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. Width and orientation of the multistatic falloff function for two bistatic receivers at a 350 km along-track baseline and a monostatic system. The last panel shows the combined coverage in the wavenumber domain.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Fraction of the spectrum captured within 3 dB fall-off for a monostatic (Sentinel-1 only) and a multistatic (Harmony) system.</figcaption></figure>

the ocean-wave spectrum does not have any signal, the SAR spectrum exhibits signal. Except for the range direction, where no velocity bunching occurs, ocean-wave signals get smeared within the fall-off function. Is is therefore that wind waves play a critical role in the determination of the cut-off [33]. With increasing wind speed, the wind-wave energy becomes larger, thereby inducing a stronger SAR spectral response at the costs of a narrower spectral coverage. As the mean wave direction

is better aligned with the ground-range direction of the bistatic trailing satellite, its corresponding spectral response is stronger than that of the heading satellite.

Except for very low wind speeds, wind-sea systems have directional spreading larger wider than a single quadrant. With the intended line-of-sight diversity of Harmony, we get three slices, with a width of the fall-off function, through the wind-wave system. With respect to considering the wind-wave system as 'unresolved', like is currently done in ocean processing algorithm, Harmony would at least resolve part of the wind-wave system and at least can provide additional wind-wave-system properties, like peak wavelength and spreading. This would for example be beneficial in cases where long wind waves are not aligned with the local wind, which occurs, amongst others, in tropical cyclones [34].

## D. Directional response to swell

For the discussion of the SAR spectral response to swell, we consider a Gaussian swell-wave system with a peak wavelength of ~ 250 m, superimposed with a nearly fully-developed wind-wave system at 10 m s⁻¹. By varying the mean direction of the swell system, it allows to compute the directional response at its peak wavelength for the multistatic system of Harmony (Fig. 6). In the monostatic and bistatic
<figure><img src="image_6.png" /><figcaption>Fig. 5. Modeled monostatic (Sentinel-1) and bistatic (Harmony) SAR spectra for three wind speeds: 6 m s⁻¹ (top), 10 m s⁻¹ (middle) and 14 m s⁻¹ (bottom). The mean direction of the wind-wave system is kept at 45° with respect to cross track.</figcaption></figure>

ground-range directions, the swell response is weak due to the absence of velocity bunching and the relative small wave slopes of swell, which make the RAR response small. Fortunately, the line-of-sight diversity of Harmony ensures non-overlapping 'dead bands', where the spectral response is more than an order of magnitude weaker, as shown by the solid lines.

Once leaving the ground-range direction, the power quickly increases due to the non-linear velocity bunching. Due to larger Doppler shifts, as a consequence of the non- perpendicular iso-range and -Doppler lines, a slightly quicker increase is observed for the bistatic companions. This also results in a stronger response as shown by the difference in maxima between the yellow solid line (monostatic) and the green and purple lines (bistatic). At $\gtrsim 30°$ from the ground-range direction, the spectral response decays because of the fall-off function. As the cut-off wavelength is larger in the bistatic case, the fall-off is more rapid.

IV. VERIFICATION OF BISTATIC SAR SPECTRA

To verify the bistatic transformation, the multistatic closed-form spectra are compared to the results of a numerical simulation. As input we use WaveWatch3 ocean-wave spectra with a parametric (Elfouhaily) spectrum for the short wavelengths based on the wind speed. For the 'resolved' wave numbers we use a spectral grid of ranging from $-\frac{2\pi}{2.5}$ rad m⁻¹ to $\frac{2\pi}{2.5}$ rad m⁻¹ with a spectral resolution of $\frac{2\pi}{2000}$ rad m⁻¹ in both directions, which ensures that we capture most of the velocity variance. The RAR is based on the spectral energy at the 'unresolved' Bragg wavelengths. Note that there are small differences in RAR as the facet orientation is modelled in the numerical model, but not properly captured by the closed-form. The latter will affect primarily the ground-range direction, where the velocity bunching is absent and the higher wavenumbers near the ground-range direction, where the slopes are steep. Additionally, be aware that we use the fast (incomplete) transform of Eq. 35, which can introduce small discrepancies.
<figure><img src="image_7.png" /><figcaption>Fig. 6. Intensity in the SAR spectrum at the maximum for a swell system (~250 m peak wavelength) as a function of direction in the presence of a wind-wave system. An angle of zero represents the cross-track direction. The dashed lines show the RAR response only. The dotted lines are the ground-range directions of the observing systems.</figcaption></figure>

The comparison between the top and middle rows shows that the closed-form properly captures the geometrical transformation of the bistatic spectra. The peak signals match to within 10%, independent of the looks angle. Similar results are obtained for three other cases (supplementary material), with different surface properties. Compared to the numerical model, the spread of the signals in the closed-form model appears to be widened. By inflating the cut-off (bottom row) the match between the numerical and closed-form models improves. A cause for the increased cut-off, or quicker fall-off, is the irratic decorrelating behavior of short waves captured only by the numerical model. In reality, the cut-off always deviates from the velocity-variance driven model, with the discrepancy increasing for higher sea states. The discrepancy is primarily attributed to wave breaking, which increases signal decorrelation [17].

## V. CONCLUSIONS

This article presents developments to derive the mapping of an ocean-wave spectrum into a SAR spectrum for the bistatic configurations. We have provided an intuitive solution of the mapping for a bistatic geometry. The general solutions are also applicable to spectra from high-resolution scatterometers, such as SEASTAR. A simple approximation of the bistatic RAR was assumed for completeness. Monostatic RAR for moderate incident angles is still a topic of active research and the formulation of a bistatic solution is therefore beyond the scope of the present developments. We have carefully described the properties of the bistatic SAR spectra in terms of the Doppler shifts, the cut-off, the iMACS and its associated directional responses. The bistatic transformation was verified using a comparison with a numerical model.

Three benefits are identified for having two bistatic companions next to an active monostatic system, e.g. Earth Explorer 10 Harmony. Firstly, two additional lines-of-sight bring added spectral coverage, particularly for waves travelling at intermediate angles with respect to the azimuth direction. This helps to constrain the swell-wave spectrum and resolve a part of the wind-wave system. Secondly, an inversion of the geophysical parameters benefits from six cut-offs with two polarizations for each satellite. The cut-offs will provide additional information on the wind-wave system and potentially help to separate

the contributions from velocity variance and decorrelation (by wave breaking). Thirdly, six MACS parameters are available. The MACS primarily serve to constrain the wave-Doppler and are sensitive to the wind-driven part of the wave spectrum, with different relaxation scales than the observed NRCS.

Future investigations will concentrate on further testings and numerical developments to improve the determination of the RAR modulations under monostatic and bistatic configurations. Efforts will further be conducted to fully benefit from all measurable parameters to advance improved retrieval strategies to consistently derive surface wind, wave and current directional properties from combined mono- and bistatic SAR and RAR measurements.

## APPENDIX A: GENERAL DOPPLER SHIFT VECTORS

The position of a scatterer in a SAR image is given by the azimuth time at which the scatterer is seen with a reference Doppler, which is time at which $\dot{r}_b$ has as given value, and its two-way range at that azimuth time. Let us now consider a scatterer located at a position $x = x_s$ at $y = t_a = 0$. Let us also assume that it moves with a velocity $v_s$, adding a term to the bistatic range derivative

$$v_b = v_s \cdot (\hat{r}_t + \hat{r}_r), \qquad (38)$$

which will be considered to remain constant during the synthetic aperture. The azimuth position of the scatterer is now given by

$$\dot{R}_b(x_p, t_a) + v_b = \dot{R}_0(x), \qquad (39)$$

where the right hand expresses the fact that the reference velocity is, in general, range dependent. Using the time and spatial derivatives of the bistatic range this can be expanded as

$$\dot{R}_b(x_p) + \ddot{R}_b(x_p) \cdot t_a + v_b = \dot{R}_b(x_p) + \frac{\partial}{\partial x} \dot{R}_b(x_p) \cdot d_x, \quad (40)$$

which simplifies to

$$\ddot{R}_b(x_p) \cdot t_a + v_b = \frac{\partial}{\partial x} \dot{R}_b(x_p) \cdot d_x, \qquad (41)$$

and where

$$\ddot{R}_b \approx -\mathbf{U}_t \cdot \frac{\partial}{\partial t} \hat{r}_t - \mathbf{U}_r \cdot \frac{\partial}{\partial t} \hat{r}_r. \qquad (42)$$
<figure><img src="image_8.png" /><figcaption>Fig. 7. SAR spectra for multistatic system with a single transmitter and two bistatic receivers separated by 350 km along-track. The results are normalized to match the spectral density. The top panels show the numerical model, whereas the middle panels show the results of the closed form. As the part of the velocity variance or signal decorrelation is not captured by the closed-form model the cut-off is inflated to match that of the numerical model in the bottom panels.</figcaption></figure>

Likewise the apparent ground range position of the target is given by

$$ (\dot{R}_b(x_p) + v_b)t_a + \frac{1}{2}\ddot{R}_b(x_p)t_a^2 = \frac{\partial}{\partial x}R_b(x_p) \cdot d_x. \quad (43) $$

The non-linear system of equations formed by (41) and (43) can be readily solved by using (41) to express $d_x$ as a function of $t_a$:

$$ d_x = \frac{\ddot{R}_b t_a + v_b}{\partial \dot{R}_b / \partial x} \quad (44) $$

This can be substituted in (41) resulting in a second-order equation. The azimuth time is given by

$$ t_a = \frac{(\ddot{R}_b C - \dot{R}_b - v_b) - D \cdot \sqrt{(\ddot{R}_b C - \dot{R}_b - v_b)^2 + 2\ddot{R}_b C v_b}}{\ddot{R}_b} \quad (45) $$

with

$$ D = \operatorname{sign}(\ddot{R}_b C - \dot{R}_b - v_b), \quad (46) $$

which can be reorganized as

$$ t_a = \frac{(\ddot{R}_b C - \dot{R}_b - v_b) - D \cdot \sqrt{(\ddot{R}_b C - \dot{R}_b)^2 + v_b \cdot (v_b + 2\dot{R}_b)}}{\ddot{R}_b} \quad (47) $$

with the C,

$$ C = \frac{\partial R_b / \partial x}{\partial \dot{R}_b / \partial x} \quad (48) $$

The azimuth shift of the scatterer will be given by

$$ d_y = t_a \cdot v_g, \quad (49) $$

with $v_g$ the effective ground velocity of the observing system, which in a flat Earth approximation is equal to the linear velocity of the radars, $|U|$.

## APPENDIX B: GAUSSIAN WAVE SPECTRUM

Let the frequency $f$ be described as

$$ f = \frac{1}{2\pi} \sqrt{gk}, \quad (50) $$

then the one dimensional frequency spectrum is given by

$$ S_f = A e^{-(\frac{f-f_p}{2\sigma_f})^2} 0.5f, \quad (51) $$
with $f_p$ the peak frequency, $\sigma_f$ the frequency standard deviation and the amplitude as a function of significant wave height $H_s$, such that

$$A = \left(\frac{H_s}{4}\right)^2 \frac{1}{\sqrt{2\pi}\sigma_f}. \qquad (52)$$

The directional spreading function is given as

$$D_{\phi} = \frac{1}{\sqrt{2\pi}\sigma_{\phi}} e^{-\left(\frac{\phi-\phi_{w}}{2\sigma_{\phi}}\right)^2}, \qquad (53)$$

where $\phi$, $\phi_w$ and $\sigma_\phi$ are the direction, mean wave direction and its standard deviation. In Cartesian coordinates the full two-dimensional wave spectrum has to be scaled by $k^{-1}$ and therefore described by

$$S(k_x, k_y) = S_f(k_x, k_y) D_{\phi}(k_x, k_y)/k, \qquad (54)$$

where the wave direction and wave numbers can simply be computed from $k_x$ and $k_y$.

## APPENDIX C: ELFOUHAILY WAVE SPECTRUM

Wind waves in homogeneous wind conditions are more realistically modeled using the spectrum described by [28], which is briefly reviewed here. The Elfouhaily spectrum is composed of a short-wave $B_h$ and a long-wave curvature spectrum $B_l$, such that

$$S_k = k^{-3}(B_l + B_h). \qquad (55)$$

The equations and parameters to compute $B_l$ and $B_h$ can all be found in [28] and will not be further elaborated here. Directional spreading is taken into account by Eq. 57 in [28] given as

$$\Delta_k = \tanh \left( \alpha_0 + \alpha_p (c/c_p)^{2.5} + \alpha_m (c_m/c)^{2.5} \right), \qquad (56)$$

where the constants $\alpha_0 = \frac{\ln(2)}{4}$, $\alpha_p = 4$ and the minimum phase velocity $c_m = \sqrt{g/k_m + \gamma/\rho k_m}$, with $k_m$ the wave number at the minimum phase velocity. The terms $c_p$ and $c$ are the phase velocities at the peak $k_p$ and wave number $k$, and are computed using similar functions as for the minimum phase velocity. The equations for phase velocity depend on the surface tension $\gamma = 0.072 \text{ N m}^{-1}$ and density $\rho = 1000 \text{ kg m}^{-3}$ of water. The spreading is given as a function of $\Delta_k$, such that

$$D_{\phi,E} = \frac{1}{2\pi} (1 + \Delta_k \cos(2\phi)). \qquad (57)$$

In analogue to the Gaussian wave spectrum, the two dimensional Cartesian Elfouhaily spectrum is computed with

$$S_E(k_x, k_y) = S_k(k_x, k_y) D_{\phi,E}(k_x, k_y)/k. \qquad (58)$$

## REFERENCES

[1] R. Raney, “Synthetic Aperture Imaging Radar and Moving Targets,” *IEEE Transactions on Aerospace and Electronic Systems*, vol. AES-7, no. 3, pp. 499–505, May 1971. [Online]. Available: http://ieeexplore.ieee.org/document/4103740/

[2] T. Larson, L. Moskowitz, and J. Wright, “A note on SAR imagery of the ocean,” *IEEE Transactions on Antennas and Propagation*, vol. 24, no. 3, pp. 393–394, May 1976. [Online]. Available: http://ieeexplore.ieee.org/document/1141351/

[3] W. Alpers and C. Rufenach, “The effect of orbital motions on synthetic aperture radar imagery of ocean waves,” *IEEE Transactions on Antennas and Propagation*, vol. 27, no. 5, pp. 685–690, Sep. 1979. [Online]. Available: http://ieeexplore.ieee.org/document/1142163/

[4] W. R. Alpers, D. B. Ross, and C. L. Rufenach, “On the detectability of ocean surface waves by real and synthetic aperture radar,” *Journal of Geophysical Research*, vol. 86, no. C7, p. 6481, 1981. [Online]. Available: http://doi.wiley.com/10.1029/JC086iC07p06481

[5] K. Hasselmann, R. K. Raney, W. J. Plant, W. Alpers, R. A. Shuchman, D. R. Lyzenga, C. L. Rufenach, and M. J. Tucker, “Theory of synthetic aperture radar ocean imaging: A MARSEN view,” *Journal of Geophysical Research*, vol. 90, no. C3, p. 4659, 1985. [Online]. Available: http://doi.wiley.com/10.1029/JC090iC03p04659

[6] K. Hasselmann and S. Hasselmann, “On the nonlinear mapping of an ocean wave spectrum into a synthetic aperture radar image spectrum and its inversion,” *Journal of Geophysical Research*, vol. 96, no. C6, p. 10713, 1991. [Online]. Available: http://doi.wiley.com/10.1029/91JC00302

[7] H. Krogstad and H. Schyberg, “On Hasselmann’s Nonlinear Ocean-SAR Transformation,” in [*Proceedings*] IGARSS’91 Remote Sensing: Global Monitoring for Earth Management, vol. 2. Espoo, Finland: IEEE, 1991, pp. 841–846. [Online]. Available: http://ieeexplore.ieee.org/document/580017/

[8] H. E. Krogstad, “A simple derivation of Hasselmann’s nonlinear ocean-synthetic aperture radar transform,” *Journal of Geophysical Research*, vol. 97, no. C2, p. 2421, 1992. [Online]. Available: http://doi.wiley.com/10.1029/91JC03010

[9] H. E. Krogstad, O. Samset, and P. W. Vachon, “Generalizations of the non-linear ocean-SAR transform and a simplified SAR inversion algorithm,” *Atmosphere-Ocean*, vol. 32, no. 1, pp. 61–82, Mar. 1994. [Online]. Available: http://www.tandfonline.com/doi/abs/10.1080/07055900.1994.9649490

[10] G. Engen and H. Johnsen, “SAR-ocean wave inversion using image cross spectra,” *IEEE Transactions on Geoscience and Remote Sensing*, vol. 33, no. 4, pp. 1047–1056, Jul. 1995. [Online]. Available: http://ieeexplore.ieee.org/document/406690/

[11] B. Chapron, H. Johnsen, and R. Garello, “Wave and wind retrieval from sar images of the ocean,” *Annales Des Télécommunications*, vol. 56, no. 11-12, pp. 682–699, Nov. 2001. [Online]. Available: https://link.springer.com/10.1007/BF02995562

[12] P. Lopez-Dekker, H. Rott, P. Prats-Iraola, B. Chapron, K. Scipal, and E. D. Witte, “Harmony: an Earth Explorer 10 Mission Candidate to Observe Land, Ice, and Ocean Surface Dynamics,” in *IGARSS 2019 - 2019 IEEE International Geoscience and Remote Sensing Symposium*. Yokohama, Japan: IEEE, Jul. 2019, pp. 8381–8384. [Online]. Available: https://ieeexplore.ieee.org/document/8897983/

[13] M. Kleinherenbrink, A. Korosov, T. Newman, A. Theodosiou, A. S. Komarov, Y. Li, G. Mulder, P. Rampal, J. Stroeve, and P. Lopez-Dekker, “Estimating instantaneous sea-ice dynamics from space using the bi-static radar measurements of Earth Explorer 10 candidate Harmony,” *The Cryosphere*, vol. 15, no. 7, pp. 3101–3118, Jul. 2021. [Online]. Available: https://tc.copernicus.org/articles/15/3101/2021/

[14] A. Theodosiou, M. Kleinherenbrink, and P. López-Dekker, “Wide-Swath Ocean Altimetry Using Multisatellite Single-Pass Interferometry,” *IEEE Transactions on Geoscience and Remote Sensing*, vol. 61, pp. 1–21, 2023. [Online]. Available: https://ieeexplore.ieee.org/document/10155423/

[15] Y. Ying, W. Xiaoqing, Z. Minhui, and C. Jinsong, “Study on bistatic SAR ocean wave imaging mechanism,” in *2009 IEEE International Geoscience and Remote Sensing Symposium*. Cape Town, South Africa: IEEE, 2009, pp. II-41–II-44. [Online]. Available: http://ieeexplore.ieee.org/document/5417994/

[16] J. Schulz-Stellenfleth, “Spaceborne synthetic aperture radar observations of ocean waves traveling into sea ice,” *Journal of Geophysical Research*, vol. 107, no. C8, 2002. [Online]. Available: http://doi.wiley.com/10.1029/2001JC000837

[17] H. Li, A. Mouche, H. Wang, J. E. Stopa, and B. Chapron, “Polarization Dependence of Azimuth Cutoff From Quad-Pol SAR Images,” *IEEE Transactions on Geoscience and Remote Sensing*, vol. 57, no. 12, pp. 9878–9887, Dec. 2019. [Online]. Available: https://ieeexplore.ieee.org/document/8798973/

[18] R. Husson and P. Vincent, “Sentinel-1 Ocean Swell Wave Spectra (OSW) Algorithm Definition,” p. 56, 2020.

[19] L. Iannini, D. Comite, N. Pierdicca, and P. Lopez-Dekker, “Rough-Surface Polarimetry in Companion SAR Missions,” *IEEE Transactions on Geoscience and Remote Sensing*, vol. 60, pp. 1–15, 2022. [Online]. Available: https://ieeexplore.ieee.org/document/9751715/
[20] T. Elfouhaily, D. R. Thompson, D. Vandemark, and B. Chapron, "A new bistatic model for electromagnetic scattering from perfectly conducting random surfaces," *Waves in Random Media*, vol. 9, no. 3, pp. 281-294, Jul. 1999. [Online]. Available: http://www.tandfonline.com/doi/abs/10.1088/0959-7174/9/3/301

[21] T. Elfouhaily, D. R. Thompson, D. E. Freund, D. Vandemark, and B. Chapron, "A new bistatic model for electromagnetic scattering from perfectly conducting random surfaces: numerical evaluation and comparison with SPM," *Waves in Random Media*, vol. 11, no. 1, pp. 33-43, Jan. 2001. [Online]. Available: http://www.tandfonline.com/doi/abs/10.1088/0959-7174/11/1/303

[22] G. Engen, P. Vachon, H. Johnsen, and F. Dobson, "Retrieval of ocean wave spectra and RAR MTF's from dual-polarization SAR data," *IEEE Transactions on Geoscience and Remote Sensing*, vol. 38, no. 1, pp. 391-403, Jan. 2000. [Online]. Available: http://ieeexplore.ieee.org/document/823935/

[23] S. Jacobsen and K. A. Høgda, "Estimation of the real aperture radar modulation transfer function directly from synthetic aperture radar ocean wave image spectra without a priori knowledge of the ocean wave height spectrum," *Journal of Geophysical Research: Oceans*, vol. 99, no. C7, pp. 14291-14302, Jul. 1994. [Online]. Available: https://agupubs.onlinelibrary.wiley.com/doi/10.1029/94JC00633

[24] A. Goldfinger, "Estimation of Spectra from Speckled Images," *IEEE Transactions on Aerospace and Electronic Systems*, vol. AES-18, no. 5, pp. 675-681, Sep. 1982. [Online]. Available: http://ieeexplore.ieee.org/document/4102711/

[25] V. Kudryavtsev, M. Yurovskaya, B. Chapron, F. Collard, and C. Donlon, "Sun glitter imagery of ocean surface waves. Part 1: Directional spectrum retrieval and validation: SUN GLITTER IMAGERY OF SURFACE WAVES," *Journal of Geophysical Research: Oceans*, vol. 122, no. 2, pp. 1369-1383, Feb. 2017. [Online]. Available: http://doi.wiley.com/10.1002/2016JC012425

[26] D. Hauser, C. Tourain, L. Hermozo, D. Alraddawi, L. Aouf, B. Chapron, A. Dalphinet, L. Delaye, M. Dalila, E. Dormy, F. Gouillon, V. Gressani, A. Grouazel, G. Guitton, R. Husson, A. Mironov, A. Mouche, A. Ollivier, L. Oruba, F. Piras, R. Rodriguez Suquet, P. Schippers, C. Tison, and N. Tran, "New Observations From the SWIM Radar On-Board CFOSAT: Instrument Validation and Ocean Wave Measurement Assessment," *IEEE Transactions on Geoscience and Remote Sensing*, vol. 59, no. 1, pp. 5-26, Jan. 2021. [Online]. Available: https://ieeexplore.ieee.org/document/9104702/

[27] L. Aouf, D. Hauser, B. Chapron, A. Toffoli, C. Tourain, and C. Peureux, "New Directional Wave Satellite Observations: Towards Improved Wave Forecasts and Climate Description in Southern Ocean," *Geophysical Research Letters*, vol. 48, no. 5, Mar. 2021. [Online]. Available: https://onlinelibrary.wiley.com/doi/10.1029/2020GL091187

[28] T. Elfouhaily, B. Chapron, K. Katsaros, and D. Vandemark, "A unified directional spectrum for long and short wind-driven waves," *Journal of Geophysical Research: Oceans*, vol. 102, no. C7, pp. 15781-15796, Jul. 1997. [Online]. Available: http://doi.wiley.com/10.1029/97JC00467

[29] H. Li, B. Chapron, A. Mouche, and J. E. Stopa, "A New Ocean SAR Cross-Spectral Parameter: Definition and Directional Property Using the Global Sentinel-1 Measurements," *Journal of Geophysical Research: Oceans*, vol. 124, no. 3, pp. 1566-1577, Mar. 2019. [Online]. Available: https://onlinelibrary.wiley.com/doi/abs/10.1029/2018JC014638

[30] V. Kerbaol, B. Chapron, and P. W. Vachon, "Analysis of ERS-1/2 synthetic aperture radar wave mode imagettes," *Journal of Geophysical Research: Oceans*, vol. 103, no. C4, pp. 7833-7846, Apr. 1998. [Online]. Available: https://onlinelibrary.wiley.com/doi/10.1029/97JC01579

[31] D. Lyzenga, "Numerical Simulation of Synthetic Aperture Radar Image Spectra for Ocean Waves," *IEEE Transactions on Geoscience and Remote Sensing*, vol. GE-24, no. 6, pp. 863-872, Nov. 1986. [Online]. Available: http://ieeexplore.ieee.org/document/4072557/

[32] R. Husson and P. Vincent, "Sentinel-1 Ocean Swell Wave Spectra (OSW) Algorithm Definition," 2021.

[33] J. E. Stopa, F. Ardhuin, B. Chapron, and F. Collard, "Estimating wave orbital velocity through the azimuth cutoff from space-borne satellites," *Journal of Geophysical Research: Oceans*, vol. 120, no. 11, pp. 7616-7634, Nov. 2015. [Online]. Available: https://agupubs.onlinelibrary.wiley.com/doi/10.1002/2015JC011275

[34] V. Kudryavtsev, M. Yurovskaya, and B. Chapron, "2D Parametric Model for Surface Wave Development Under Varying Wind Field in Space and Time," *Journal of Geophysical Research: Oceans*, vol. 126, no. 4, Apr. 2021. [Online]. Available: https://onlinelibrary.wiley.com/doi/10.1029/2020JC016915

Marcel Kleinherenbrink obtained his master's degree at the faculty of Aerospace Engineering of Delft University of Technology, Delft, The Netherlands in 2013. He obtained his Ph.D. degree, on the observation of sea-level change and vertical land motion from space, at the faculty of Civil Engineering and Geoscience of Delft University of Technology in 2018. He specialises in the observation of the ocean surface using active radar instruments. In 2019 he became involved in the Earth Explorer 10 mission Harmony, for which he implements forward models, investigates inversion strategies and provides performance estimates. Since 2018, he has been involved in numerous (ESA) projects related to future missions and novel retrieval algorithms for ocean observations.

Paco Lopez-Dekker was born in Nijmegen, The Netherlands, in 1972. He received the Ingeniero degree in telecommunication engineering from the Universitat Politècnica de Catalunya (UPC), Barcelona, Spain, in 1997, the M.S. degree in electrical and computer engineering from the University of California, Irvine, CA, USA, in 1998, under the Balsells Fellowship, and the Ph.D. degree in clear-air imaging radar systems to study the atmospheric boundary layer from the University of Massachusetts, Amherst, MA, USA, in 2003. In 2003, he joined Starlab Barcelona, where he worked on the development of GNSS-R sensors and techniques. From 2004 to 2006, he was a Visiting Professor with the Department of Telecommunications and Systems Engineering, Universitat Autonoma de Barcelona. In March 2006, he was awarded a Ramon y Cajal Grant to conduct pioneering research on bistatic synthetic aperture radar (SAR) at Remote Sensing Laboratory, UPC. Between November 2009 and August 2016, he led the SAR Missions Group at the Microwaves and Radar Institute, German Aerospace Center, Wessling, Germany. Since September 2016, he has been an Associate Professor with the Geoscience and Remote Sensing Department, Faculty of Civil Engineering and Geosciences. He has been deeply involved in the development of several radar missions and mission proposals, and is the Lead Investigator of the Harmony ESA Earth Explorer 10 Mission. He has coauthored over 50 peer-reviewed journal articles and more than 125 conference contributions in a broad range of topics related to radar remote sensing.

Frédéric Nouguier received the Agrégation and M.S. degrees in applied physics from the Ecole Normale Supérieure de Cachan, Cachan, France, the M.S. degree in physical methods for remote sensing from the University of Paris-Diderot, Paris, France, and the Ph.D. degree in physics from Aix-Marseille University, Marseille, France, in 2009. He is currently a Research Scientist with the Laboratoire d'Océanographie Physique et Spatiale, Institut Français de Recherche pour l'Exploitation de la Mer, Plouzané, France. He has experience in applied mathematics, physical oceanography, and electromagnetic wave theory and its application to ocean remote sensing.

Bertrand Chapron was born in Paris, France, in 1962. He received the B.Eng. degree from the Institut National Polytechnique de Grenoble, Grenoble, France, in 1984, and the Doctorate National (Ph.D.) degree in fluid mechanics from the University of Aix-Marseille II, Marseille, France, in 1988. He was a Post-Doctoral Research Associate with the NASA/GSFC/Wallops Flight Facility, Wallops Island, VA, USA, for three years. He has experience in applied mathematics, physical oceanography, electromagnetic waves theory, and its application to ocean remote sensing. He is currently responsible for oceanography at the Space Laboratory, IFREMER, Plouzané, France.