#  Biases from spectral leakage in remote sensing of near-surface currents

Stefan Weichert, Benjamin K. Smeltzer, Simen Å. Ellingsen The work of S. Å.
Ellingsen was supported in part by the European Union (ERC, WaTurSheD, project
101045299) and the Research Council of Norway (project 325114). Views and
opinions expressed are however those of the authors only and do not
necessarily reflect those of the European Union or the European Research
Council. Neither the European Union nor the granting authority can be held
responsible for them.S. Weichert and S. Å. Ellingsen are with the Department
of Energy and Process Engineering, Norwegian University of Science and
Technology, 7491 Trondheim, Norway.B. K. Smeltzer is with SINTEF Ocean,
Marinteknisk senter, 7052 Trondheim, Norway

###### Abstract

Remotely measuring subsurface water currents from imagery of the wave field
has become a much-used technique. We study the biases and errors in such
measurements due to spectral leakage, and suggest mitigating procedures.
Deviations between peak values in the three-dimensional wave spectrum and the
known dispersion relation in quiescent water are extracted and interpreted as
current-induced Doppler shifts, from which the sub-surface current is
inferred. The use of discrete Fourier transforms, however, introduces spectral
leakage between nearby frequency bins. Analysing synthetically generated wave
data adhering to realistic input spectra we show that although no current is
in fact present, spurious currents can be “measured” which can amount to a
significant fraction of the phase speed at the spectral peak. We analyse the
effects of data tapering, method of Doppler shift extraction, limited
wavenumber and frequency resolution, peakedness and angular width of the input
spectrum, and average misalignment between waves and Doppler shift velocity
direction. The narrower the input wave spectrum in frequency and/or direction,
the greater the biases become. The use of a window function reduces the
severity in nearly all cases, yet mitigates the effects of limited resolution
more effectively in space than in time. When a current is present the absolute
biases remain essentially unchanged, when waves and currents are roughly
aligned tapering much alleviates the concomitant biases, whereas in the case
of a cross-current, biases remain significant even for tapered data.

###### Index Terms:

Remote sensing, Surface waves, Sea measurements, Dispersion, Sea surface,
Spectral Analysis, Signal processing algorithms

##  I Introduction

The prospect of measuring currents near the sea surface remotely from above is
a highly attractive one. Measuring depth-varying currents _in situ_ by
penetrating the surface requires the use of e.g. buoys, ships, gliders or
fixed instruments, all relatively expensive and able to measure a single point
or trajectory at a time, and often struggle to capture currents in the top few
metres. In comparison, remote measurement from above can be performed with
inexpensive equipment mounted on airborne platforms able to cover larger areas
in a short time (see, e.g., [1]).

By far the most common source of wave data for this purpose has been
measurements using HF or X-band radar [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,
14, 15, 16], primarily mounted on ships. Only observation of the wave phase
variation in time and space is required, however, meaning other methods are
equally applicable; the use of infrared [17] and polarimetric [18, 19] imaging
has been demonstrated, as has regular optical measurements (video) with
cameras mounted on quadcopter drones [20], and aircraft [21, 22, 23, 24].

A current varying with depth will affect the wave phase velocity differently
for different wavelengths, resulting in a measurable effective Doppler shift
which depends on wavenumber [2] and appears as a shift in the waves’ spectral
dispersion curve from that observed in quiescent water. Ever more advanced
methods have been developed in recent years for inferring the depth profile of
the sub-surface currents from such measured Doppler shifts [12, 25, 26], yet
the task of obtaining these shifts from an observation of the spatiotemporal
wave field, by reconstructing the current-modified dispersion surface in the
frequency-wavenumber spectrum, is itself a nontrivial task. Typically, the
spectrum is divided into wavenumber magnitude bins, whereby the Doppler shifts
are found for each bin separately. Perhaps the most frequently employed are
least-squares-based methods (e.g. [4, 8, 27]), whereas the alternative
Normalized Scalar Product (NSP) method is also in regular use (e.g. [10, 28,
13]). We compare these methods herein finding NSP to be unequivocally
favourable. A further method which we shall not consider here is the so-called
Polar Current Shell method (e.g. [13]), recently adapted for this purpose
[26], which has similar performance as NSP, but favourable in some
circumstances.

One should note that the same questions we seek to answer here, also apply to
bathymetry retrieval from the measured wave spectra (see e.g. [29] and
references therein). The extraction of the water depth from the measured
spectrum requires spectral intensity at low wavenumbers, which can be strongly
influenced by spectral leakage, as we will show. Since a finite water depth
enters the dispersion relation in a direction-independent and multiplicative
way in contrast to the additive and anisotropic term due to a current, it is
not obvious how the issues discussed herein affect the bathymetry retrieval
quantitatively, a question of potential importance which requires further
investigation.

###  I-A Outline

In this work we consider the effects of a strong spectral energy peak and the
corresponding effect it has on the extraction of Doppler shifts for
wavenumbers in the vicinity of the peak. The energy peak results in spectral
leakage to adjacent wavevectors in the wave spectrum and may cause errors in
the extracted Doppler shifts — as perhaps the clearest example we show that in
realistic seastates, spectral leakage can cause a significant spurious current
to be “measured” when none is in fact present.

The phenomenon of spectral leakage is briefly reviewed in section II-B, and we
discuss how it can be mitigated by data windowing. We then go on to describe a
numerical experiment where mock wave data is generated and analysed with the
methods in standard use for remote-sensing of currents from wave dispersion,
in section III. Results are reported studying how the spurious current
“measurements” depend on the method of Doppler shift extraction (section
III-C), the angular width and peakedness of a directional wave spectrum of the
Joint North Sea Wave Project type (JONSWAP) [30] (sections VI-B and VI-C), the
spatial and temporal resolution (sections VI-D and VI-D), as well as how the
situation changes when a uniform background current is present, in section
VI-E. We finally summarize and give a brief overview of ways whereby the
detrimental effect of spectral leakage can be reduced in practical
applications, in section VII.

##  II Background

When remote sensing of sub-surface currents from wave dispersion is performed,
the input spectrum is an observation of the motion of the water surface
resolved in time and space. The method is based on linear wave theory, so only
the phase of the waves is required, not the amplitude. The data is typically a
monotonic function of the sea surface elevation or its derivative, as a
function of position $(x,y)$ and time $t$. A three dimensional Fourier
transform is then applied in space and time to obtain a spectral signal as a
function of wave numbers $\bm{k}=(k_{x},k_{y})$ and frequency $\omega$. We
assume the spatial area and duration of the observation are $L\times L$ and
$T$, respectively, so that the resolution in wave-number and frequency are,
respectively, $\delta k=2\pi/L$ and $\delta\omega=2\pi/T$; we will refer to
these as a pixel or bin in wavenumber and frequency, respectively. (Note that
after we introduce nondimensional units in section III-A $\delta k$ and
$\delta\omega$ take the forms $1/L$ and $1/T$, respectively.) We will assume
infinitely deep water for simplicity herein.

###  II-A Theory

In a wave spectrum, the spectral signal is concentrated near the dispersion
relation $\omega=\omega_{\text{DR}}(\bm{k})$. The methods for sensing the sub-
surface current now extract a measured function $\omega(\bm{k})$. If a current
with moderately strong depth-dependence $\bm{U}(z)=(U(z),V(z))$ is present,
the dispersion relation for a wave with wave vector $\bm{k}$ is well
approximated as $\omega=\omega_{\text{DR}}(\bm{k};\tilde{\bm{c}})$ with the
dispersion function

$\omega_{\text{DR}}(\bm{k};\tilde{\bm{c}})=\omega_{0}(k)+\bm{k}\cdot\tilde{\bm{c}}$
(1)

where $k=|\bm{k}|$, $\tilde{\bm{c}}$ is the Doppler shift velocity (DSV) due
to the presence of a sub-surface current and the dispersion relation in deep,
quiescent water is

$\omega_{0}(k)=\sqrt{gk}.$ (2)

Inversion methods to infer depth-dependent velocity profiles from the measured
spectrum are based on the approximation [2, 31]

$\tilde{\bm{c}}=2k\int_{-\infty}^{0}\bm{U}(z)e^{2kz}dz.$ (3)

If a current is present which is uniform in depth, the resulting Doppler shift
will be independent of $\bm{k}$, while conversely, a Doppler shift which
varies with $\bm{k}$ implies the presence of a current which varies as a
function of $z$. We shall see that biases in the measured Doppler shift
$\tilde{\bm{c}}$ due to spectral leakage typically vary significantly with
$\bm{k}$, and hence the spurious currents which are “measured” will have a
nontrivial depth dependence. We do not pursue this question in detail.

The presence of the current thus introduces an observable Doppler shift
corresponding to the addition of a phase speed $\tilde{\bm{c}}(k)$ to the
phase velocity. Remote sensing of the depth-varying current $\bm{U}(z)$ is
then possible by measuring $\tilde{\bm{c}}$ and inverting equation (3) using
one of several methods available as reviewed in [1]. Methods for extracting
$\tilde{\bm{c}}$ from a measured spectrum are reviewed and compared in section
III-C.

Following Smeltzer et al. [26] we define the instructive quantities $\delta
c_{\delta\omega}$ and $\delta c_{\delta k}$ as

$\delta c_{\delta\omega}=\frac{\delta\omega}{k};~{}~{}~{}\delta c_{\delta
k}=\frac{\partial\omega}{\partial k}\frac{\delta k}{k},$ (4)

which estimate the change in the predicted phase velocity $c$ due to moving
the dispersion surface (1) by $\delta k$ along the wavenumber axis in the
spectrum, or by $\delta\omega$ along the frequency axis. These are thus
approximate measures of the uncertainty in velocity measurement introduced by
limited wavenumber and frequency resolution, respectively. Further discussion
may be found in section 4.2.1 of [26].

###  II-B Spectral leakage and windowing

Assume a continuous signal (in space or time) $\tilde{P}_{0}(t)$, the “true”
signal, is measured during a finite period of duration $T$. A sharp cut-off at
the beginning and end of the measurement is equivalent to multiplying
$\tilde{P}_{0}$ by a discontinuous top-hat function $w_{\text{box}}(t)$ which
is $1$ within a time interval of length $T$ and zero outside this “window”.
Multiplication with such a window is equivalent to a convolution of the
spectrum ${P}_{0}=\mathcal{F}\\{\tilde{P}_{0}\\}$ with a $\operatorname{sinc}$
function in frequency space (e.g. [32]),

$P_{\text{measured}}(f)=P_{0}(f)*T\operatorname{sinc}(fT),$ (5)

where $\operatorname{sinc}(a)=\sin(\pi a)/\pi a$, which is the Fourier
transform of the top hat function. The result is a blurring of the spectrum:
“true” spectral components will appear as spectral intensity not only in the
frequency bins closest to the actual frequency of said component, but also on
neighbouring ones that are more than one bin away, with an intensity
decreasing with distance. One can alleviate this by replacing the “box” window
with a window of choice, simply by multiplying the acquired data
$\tilde{P}(t)=\tilde{P}_{0}w_{\text{box}}(t)$ with said window function
$w(t)$. If the chosen $w(t)$ also vanishes outside the measurement interval,
the window replaces the top-hat. The choice of the optimal window function
much depends on the situation and data at hand.

Figure 1:  Fourier transform $\mathcal{F}\\{w\\}(f)$ of a selection of common
window functions $w(t)$. A frequency pixel is $\delta\omega=2\pi/T$,
normalized so that $\mathcal{F}\\{w\\}(0)=1$. This illustrates how long-range
spectral leakage can be mitigated by tapering the measured signal. The
drawback of increased short-range leakage (1-2 pixels) is often negligible.

The Fourier transforms $\mathcal{F}\\{w\\}(f)$ of some common window functions
are shown in figure 1 (See e.g. [33] for a more comprehensive comparison).
When taking a discrete Fourier transform (DFT) of $w(t)$ every real frequency
component $f$ will generate spectral intensity in frequency bins of width
$\delta f$ centered at $n\delta f$, $n\in Z$, with intensity
$\mathcal{F}\\{w\\}(n\delta f-f)$. For example, for the _Blackman window_
shown in figure 1, a frequency coinciding with a DFT frequency, i.e.
$f=m\delta f$, one obtains non-zero spectral intensity only for $n-m=0,1,2$,
leading to no long-range leakage whatsoever. Only the central lobe is sampled
in frequency space. For any real spectrum most frequencies will not coincide
with a DFT frequency, and so the side-lobes give rise to long-range leakage.
Now, the central lobe of the _No Window_ (or top hat) case is only one pixel
wide, but the leakage is very long-range, i.e. spectral intensity decays
slowly and stays well above $1\%$ for more than $20$ pixels. In contrast to
this, the _Blackman_ window decays rapidly to $10^{-3}$, at the cost of
increasing spectral leakage into the two nearest bins, 1 and 2. The best
choice of the windowing function is usually a compromise between suppressing
long-range leakage and blurring the spectrum (short-range leakage). In this
work, a Hann window, defined as

$w_{\text{H}}(t)=0.5-0.5\,\cos\left(2\pi\,t/T\right),$ (6)

is used, as it suppresses the long-range spectral leakage to less than 1% and
leaks significantly into frequencies less than two pixels away. In appendix -A
an example of extracted DSVs with different windowing functions is given to
illustrate how the choice affects the results.

The data in this work is surface elevation $\zeta(x,y,t)$ measured in three
dimensions $x$, $y$ and time $t$, which we pre-multiply by a 3D Hann window
constructed as $w(x,y,t)=w_{\text{H}}(x)w_{\text{H}}(y)w_{\text{H}}(t)$ prior
to subjecting it to a discrete 3-dimensional fast-Fourier transform (3DFFT)
(it is understood that $L$ replaces $T$ in equation (6) when the argument is
$x$ or $y$. The new signal $\zeta(x,y,t)w(x,y,t)$ goes smoothly to zero at the
edges of the domain of observation, in our case the square area with sides $L$
and time duration $T$.

##  III Methods

We proceed by producing synthetic surface elevation data $\zeta(x,y,t)$ by
superposing random linear plane waves of wave number $k$ and direction
$\theta$ from chosen spectra with varying properties, as detailed in section
III-B. Each wave is given a uniformly distributed random initial phase, the
frequency $\omega$ is found from equation (1), whereupon the waves are
propagated in time.

These wave “observations” are of course idealized, since different methods for
obtaining the actual surface elevation $\zeta(x,y,t)$ or true spectrum from
field data each come with their individual challenges and limitations, an
ongoing field of research in its own right. Taking such practical challenges
into account is beyond the scope of this work, and we use the ideal data to
isolate the effects of spectra leakage in the data analysis.

We mostly consider the case of quiescent water, i.e., there is no background
current and any Doppler shifts “measured” from the spectrum are spurious and
purely a consequence of spectral leakage. We also consider the case where a
constant background current $\bm{U}_{0}$ is present, in section VI-E. Biases
now manifest as deviations of the observed DSVs from the correct value.

###  III-A Nondimensional quantities

Dimensional quantities will be denoted by a superscript asterisk, all other
quantities are non-dimensionalized. The reference length-scale and time-scale
are defined based on a characteristic wavenumber $k_{0}^{*}$ and its
corresponding angular frequency in quiescent water
$\omega_{0}^{*}=\sqrt{gk_{0}^{*}}$. Thus, e.g. $T=T^{*}\omega_{0}^{*}/2\pi$,
$L=L^{*}k_{0}^{*}/2\pi$, $k=k^{*}/k^{*}_{0}$,
$\omega=\omega^{*}/\omega_{0}^{*}$, $U=U^{*}\sqrt{k_{0}^{*}/g}$. Unless
specified otherwise, $\omega_{0}^{*}(k_{0}^{*})$ is taken to be the location
of the peak of the energy spectrum, i.e., the spectral peak is at $\omega=1$
by definition, hence the peak in wavenumber space is close to to $k=1$.

###  III-B Wave spectrum

We generate wave fields from commonly used realistic model spectra with
varying directional broadness, assuming the form

$\hat{S}(\omega,\theta)=S(\omega)f(\theta)$ (7)

where $\theta$ is the angle between $\bm{k}$ and the $x$ axis,
$\cos\theta=\bm{k}\cdot\bm{e}_{x}/k$. We use the JONSWAP spectrum [30]

$S(\omega)=\tilde{N}\omega^{-5}\exp{\left[-\frac{5}{4}\omega^{-4}\right]\gamma^{r(\omega)}}$
(8)

with

$r(\omega)=\exp{\left[-\frac{1}{2}\left(\frac{\omega-1}{\sigma}\right)^{2}\right]}.$
(9)

For our purposes the value of $\tilde{N}$ is not of importance, we set
$\tilde{N}=1$. The parameter $\sigma$ is

$\sigma=\begin{cases}0.07,&{\text{if}}\ \omega\leq\omega_{p},\\\
{0.09,}&{\text{if}}\ {\omega>\omega_{p}}\end{cases}$ (10)

and the peakedness parameter $\gamma$ is varied (see section VI-C). The energy
spectrum for a selection of peakedness values $\gamma$ is depicted in figure
2.

Figure 2: Top: Energy spectrum $S(\omega)$ (equation 8) for a few values of
the peakedness parameter $\gamma$. The $\omega$ axis is scaled quadratically
to match the range of the bottom graph, as $k^{*}=\omega^{*2}/g$, i.e.,
$k=\omega^{2}$. Bottom: amplitudes $a(k,\theta=\theta_{0})$ (equation 12) of
superposed wave components. The spectra are normalized with respect to their
peak value.

The angular distribution is taken as a cosine-square with a full width
(distance of first roots) of $\Delta\theta$, i.e.

$f(\theta)=\left\\{\begin{matrix}\cos^{2}\left(\pi\frac{\theta-\theta_{0}}{\Delta\theta}\right)&,\,|\theta-\theta_{0}|\leq\Delta\theta/2\\\
0&\,\text{else}\end{matrix}\right.\,\quad\,$ (11)

We use this spectrum to prescribe the amplitudes $a(\bm{k})$ of the superposed
wave components in accordance with [34] as

$a(\bm{k})=\sqrt{2\,k^{-3/2}\,\hat{S}(\omega(k),\theta)\delta k}$ (12)

on an evenly spaced grid in $k_{x}$-$k_{y}$, with spacing $\delta k_{x}=\delta
k_{y}=\delta k$. Additionally, we set $a(\bm{k}$)=0 for $k/k_{p}>3.5$.

###  III-C Normalized scalar product

Two methods are in common use for extracting Doppler shifts from the measured
wave spectrum; least-squares (LS) methods (e.g. [4, 8, 27]) and the Normalized
Scalar Product (NSP) method (e.g. [10, 28, 13]). While there exists a wide
range of extensions and sophistications to LS methods, we focus our attention
on the NSP method. However, since two of the most famous works on spectrum
based current extraction ([4, 8]) use a simple form of the LS method we also
give a brief comparison of their performance in Appendix -C.

We employ the DSV extraction method as implemented by Smeltzer et al. [26];
see [10, 13, 28] for details on the NSP method more generally.

The starting point in either method is a measured, spatio-temporally resolved
free-surface $\eta(x,y,t)$ and its power spectrum obtained via a discrete
Fourier transform, $P(\bm{k},\omega)=|FFT[\eta(x,y,t)]|^{2}$. For each
wavenumber $k_{i}$ in a list, the spectral intensity on a cylindrical surface
with radius $k_{i}$ centered around the $\omega$-axis is defined

$F_{i}(\bm{k},\omega)=\sqrt{P(k_{i}\cos\theta,k_{i}\sin\theta,\omega)}.$ (13)

where $\theta$ is the azimuth angle in the $k_{x}$-$k_{y}$-plane. The
implementation of the algorithms used are formulated in Cartesian coordinates;
however, it is illustrative to use cylindrical coordinates for the following
conceptual considerations. More details of the implementation are given in
appendix -B.

To find the effective DSV $\tilde{\bm{c}}_{i}=\tilde{\bm{c}}(k_{i})$, first, a
characteristic function $G$ is defined that contains the components of
$\tilde{\bm{c}}_{i}$ as free parameters:

$G_{i}(\theta,\omega;\tilde{\bm{c}}_{i})=G^{+}_{i}(\theta,\omega;\tilde{\bm{c}}_{i})+G^{-}_{i}(\theta,\omega;\tilde{\bm{c}}_{i})$
(14)

where

$G^{\pm}_{i}(\theta,\omega;\tilde{\bm{c}}_{i})=\exp\left[-2\left(\frac{\omega\pm\omega_{\text{DR}}(\theta;k_{i},\tilde{\bm{c}}_{i})}{a}\right)^{2}\right].$
(15)

The normalized scalar product $N_{i}$ of the vectors $F_{i}$ and $G_{i}$ is
now maximized for each value of $i$ by varying the two components of
$\tilde{\bm{c}}_{i}$; it is calculated as

$N_{i}(\tilde{\bm{c}})=\frac{\left\langle
G_{i}F_{i}\right\rangle}{\left\langle G_{i}\right\rangle\left\langle
F_{i}\right\rangle}$ (16)

where $\left\langle...\right\rangle$ refers to an integral over all $\theta$
and $\omega$. In other words, $G_{i}$ can be thought of as a cosine in
$\theta$ (see figure 4) with offset $\omega_{0}(k_{i})$ amplitude
$k_{i}\tilde{c}$ and phase shift defined by the direction of $\tilde{\bm{c}}$;
its overlap with the measured intensity on the cylinder surface is maximized
to find the best DSV $\tilde{\bm{c}}(k_{i})$. This optimization step is
performed using the Nelder-Mead simplex method [35].

##  IV The normalized scalar product method and spectral leakage

For multidimensional data, spectral leakage is most prominent in the
directions parallel to the coordinate axes. Consider for example a 2D signal
$\eta(x,y)$ on a rectangular domain and its Fourier transform
$\tilde{\eta}(k_{x},k_{y})=\mathcal{F}\\{\eta\\}$. Since both $x$ and $y$ are
within a finite range, the effective window $w(x,y)$ is a product of top-hat
windows in the $x$ and $y$ directions, respectively, i.e.
$w=w_{1}(x)w_{2}(y)$. The Fourier transform of such a product is the product
of their respective Fourier transforms
$\tilde{w}(k_{x},k_{y})=\mathcal{F}\\{w_{1}(x)\\}\mathcal{F}\\{w_{2}(y)\\}=\tilde{w}_{1}(k_{x})\tilde{w}_{2}(k_{y})$.
Now, since the leakage for a top-hat window falls off as $1/k$, the product
$\tilde{w}(k_{x},k_{y})$ is smallest for a given $k$ when $|k_{x}|=|k_{y}|$
and largest when $k_{x}=0$ or $k_{y}=0$. An illustration of this can be seen
in figure 3. The extension to 3 or more dimensions is straightforward.

Figure 3: An illustration of spectral leakage in two dimensions. Left:
Logarithmic 2D FFT spectrum of a superposition of four sine-waves (i.e., the
spectrum is the sum of four delta functions in the $\bm{k}$ plane), with
frequencies denoted by triangles. The circles denote the DFT frequencies, i.e.
frequencies natural to the domain. From left to right, the leftmost spectral
peak coincides with a DFT frequency, while the rightmost peak lies in the
middle of DFT frequencies. Right: same as left, bit with a 2D Hann window
applied before doing the FFT. Note how the Hann window removes the background
and increases consistency between spectral peak, while broadening by roughly
one bin in all directions.

This (mostly) axes-parallel leakage is helpful in understanding how leakage
affects the DSV extraction using NSP.

A fundamental step in the NSP method is to pick out measured spectral
intensity on a cylinder surface with radius $k_{i}$. One can think of the
effect on the NSP as the spectral intensity projecting itself in the principal
directions onto the cylinder, with decreasing intensity the further the
surface is from the originating spectral intensity. Now, as the NSP method
essentially fits a function of the form
$\omega_{0}(k_{i})+A_{f}\cos(\theta-\theta_{f})$, with free parameters $A_{f}$
and $\theta_{f}$, to the spectral intensity on the cylinder surface defined by
$k_{i}$, as illustrated in figure 4, spurious intensity at
$\omega\neq\omega_{0}(k_{i})$ leads to non-zero values for
$A_{f}=\bm{k}\cdot\tilde{\bm{c}}$, implying a background current even in the
absence of one. In this work we consider spectra with a single peak at a given
wave vector $\bm{k}=(k_{x},k_{y})$.

Figure 4: Illustration of the dispersion relation
$\omega=\omega_{\text{DR}}(\bm{k};\tilde{\bm{c}})$ in connection within the
Nonlinear Scalar Product method. The top graphs show the dispersion relation
in absence of a current, $\tilde{\bm{c}}=0$ (a) and with a current,
$\tilde{\bm{c}}=\bm{U}_{0}$ (b). The intersections with a cylinder with radius
$k_{i}$ is highlighted. In graph (c) the spectral intensity is indicated in
greyscale, visible near the dispersion curve on the cylinder surface after
unrolling it, together with the lines of intersection from the top figures:
blue, dashed corresponds to quiescent water; red, solid is the best fit from
the NSP method resulting in a non-vanishing Doppler shift velocity.

We now consider the effect of spectral leakage for two distinct cases, one (1)
with the mean propagation direction $\theta_{0}$ being parallel to the
$k_{y}$-axis ($\theta_{0}=90^{\circ}$), and one (2) with $\theta_{0}$ between
$None$ and $None$.

###  IV-A Leakage for $\theta_{0}=$$$

Consider the extraction of the DSV for a wavenumber $k<1$, i.e. below the
$k$-value at the peak of the spectrum which is close to $k=1$ (see Section
III-A. The spectral leakage for the cylinder with radius $k$ (see Fig. 4) will
be dominated by short-range spectral leakage of wavenumbers close to $k$ and
long-range leakage from the spectral peak. the short-range spectral leakage
for the most part will merely blur the spectral intensity on the cylinder,
while the long-range leakage from the spectral peak leads to a new, possibly
dominating, spectral intensity spot on the cylinder, at the same angle, but a
higher frequency. In other words right above the real spectral intensity.
Depending on the ratio of intensities, distance, and width of the
characteristic function $G$, the NSP algorithm may find the correct frequency,
the frequency of the projected peak or a value in between these two to give
the best fit. Therefore the returned DSV will be parallel to the spectrum
propagation direction $\theta_{0}$. Similarly, if $k>1$ is considered, the
projection from the spectral peak will be below the real spectral intensity,
giving DSV antiparallel to the spectrum propagation direction $\theta_{0}$.

###  IV-B Leakage for $0<\theta_{0}<$$$

For $k<1$, it is now possible for spectral leakage from the peak to have no
axes-parallel projection onto the cylinder surface, thus not affecting the
extracted DSV for this wavenumber. Short-range spectral leakage and leakage
from lower wavenumbers now dominate, making the prediction of DSV direction
difficult. For $k>1$, one can still get an intuition for the effect of
leakage. As the peak of the spectrum lies in the first quadrant of the
$k_{x}-k_{y}$-plane, the shortest distance to a cylinder with radius $k$ is to
the quarter of the cylinder also lying in the first quadrant, and we can focus
our consideration on that. The leakage is now being projected onto the
cylinder at an angle relative to the mean propagation direction $\theta_{0}$.
Therefore, as the real spectral intensity now dominates for
$\theta=\theta_{0}$, the frequency for wavenumbers along $\theta_{0}$ is found
correctly. The perpendicular component of the DSV, however is determined
mostly by the position and relative magnitude of the spurious intensity
appearing at angles $\theta\neq\theta_{0}$. in quiescent water, for example,
one would obtain DSVs with directions $\theta\pm$$$.

Note that these considerations hold true for very narrow spectra, but are to
be understood as tendencies for spectra with considerable spectral width
$\Delta\theta$, because, for any angle $\theta$, real spectral intensity will
mostly dominate over spectral leakage, if present.

##  V Parameter Choice and Overview

For an overview and easier referencing table I contains the list of test cases
presented in this paper. Apart from the parameters stated therein, the
resolution of the input spectrum $\delta_{in}k$ and spatial resolution $\delta
x=\delta y$ as well as temporal resolution $\delta t$ had to be set. The
spatial resolution was chosen such that waves of the highest wavenumbers,
$k_{max}=4$, would be well resolved, $\delta x=1/28$. Similarly, time
resolution was set to $\delta t=1/14$. The resolution of the input spectrum
was set to $\delta_{in}k\approx 0.341/L$, with L being the domain size. This
was deemed small enough to mimic a continuous spectrum.

The domain size and video duration were chosen to be similar to typical
parameters found in airborne measurements. For instance, videos used in [36]
typically have a field of view of $L^{*}=$128512\text{\,}\mathrm{m}$$,
durations of $T^{*}=$2040\text{\,}\mathrm{s}$$ and wavenumbers with usable
spectral intensity from $0.22\text{\,}\mathrm{rad}\text{\,}{\mathrm{m}}^{-1}$
(see figure 1 in [36]). Taking a reference wavenumber of
$k^{*}_{0}=$0.4\text{\,}\mathrm{rad}\text{\,}{\mathrm{m}}^{-1}$$, this gives a
parameter range of $T\approx 5-10$, $L\approx 8-32$.

Equation 14 contains one free parameter, $a$, that determines the width of the
characteristic function. We set $a=4\,\delta\omega$, with $\,\delta\omega=1/T$
being the frequency resolution (see appendix -B for details).

The results shown in figures 5-9 display statistics of the extracted DSVs in
terms of the average $\left\langle|\tilde{\bm{c}}|\right\rangle$ and the
corresponding standard deviation $\sigma_{\tilde{c}}$ (times two for
illustrative purposes), calculated from 100 realizations. These quantities are
useful as they represent a mean bias and fluctuation that one has to expect
from a single measurement. Additionally, where possible, the implied velocity
resolutions $\delta c_{\delta\omega}$ and $\delta c_{\delta k}$ defined in
equation (4) are shown for reference.

| $\theta_{0}$ | $\Delta\theta$ | $\gamma$ | $L$ | $T$ | $U$ | figure | section
---|---|---|---|---|---|---|---|---
case-NSP-LS | $None$ | $None$ | 3.3 | 10 | 20 | 0 | 11 | -C
case-$\theta_{0}$-$\Delta\theta$ |  $None$, $None$ |  $None$…$None$ | 3.3 | 10 | 20 | 0 | 5 |  VI-A,VI-B
case-$\gamma$ | $None$ |  $None$…$None$ | 1…7 | 10 | 20 | 0 | 6 | VI-C
case-$L$-$T$ | $None$ | $None$ | 3.3 | 5, 10, 20 | 10, 20, 80 | 0 | 7 | VI-D
case-$U$ | $None$ | $None$ | 3.3 | 10 | 20 | 0…0.15 |  8, 9 | VI-E

TABLE I: Overview of parameter combinations. $\theta_{0}$: mean propagation
direction of input wavenumber spectrum, $\Delta\theta$: angular spread of
input wavenumber spectrum (see equation 11), $\gamma$: peakedness parameter
(see equation 8), $L$: domain spatial domain size, $T$: video duration, $U$:
background current. All parameters are non-dimensional, see section III-A

##  VI Results

In this section we consider the effects of wave-spectral properties,
resolution and data tapering on the spurious Doppler shifts “measured” when no
current is present, as well as the effect of a uniform current being present.
An overview of the parameter combinations used to obtain the following results
is given in section V.

###  VI-A Illustration of the effect of windowing (Hann window)

The mitigating effect of windowing/tapering is well illustrated when
considering the influence of the angular spread $\Delta\theta$ of the
wavenumber spectrum; simulations with variation both in $\Delta\theta$ and the
mean propagation direction $\theta_{0}$ were performed (see table I,
case-$\theta_{0}$-$\Delta\theta$, for all parameters). A selection of the
results are shown in figure 5 (upper row). A more detailed discussion of the
influence of $\Delta\theta$ and $\theta_{0}$ is given in section 5 below.

The propagation direction $\theta_{0}$ shows a strong influence on both the
mean bias and its variation. This behaviour is likely due to spectral leakage,
because spectral leakage appears as “streaks” in the spectrum along the
spectral axes $k_{x}$, $k_{y}$ and $\omega$ (See illustration in figure 3).
The second column in figure 5 shows the extracted DSVs of tapered surface
elevation data. The mean and random bias is strongly suppressed in all cases
and the dependency on $\theta_{0}$ is significantly reduced, compared to the
analysis without a Hann window. Note that $\Delta\theta=$$$ represents a rare,
very narrow spectrum, but the described effects are visible for all spectral
widths, decreasing with increasing $\Delta\theta$.

We observe that for smaller wavenumbers ($k<1$), using a Hann window does not
mitigate the effects of spectral leakage as strongly as for $k>1$. This is
because the Hann window suppresses long-range spectral leakage (more than one-
two pixels), while increasing the short-range spectra leakage (The central
lobe of the Hann window has a width of 2 bins instead of 1 for no-window). The
“steepness” $\partial\omega/\partial k$ grows with decreasing $k$, meaning
that spectral intensity that leaks in the $k$-direction appears as a strong
broadening in the $\omega$-direction at neighbouring $k$-values. In our
simulations the spectrum has a peak at $\omega=1$ (thus close to $k=1$) and
quickly decreases towards $k\approx 0.5$. Therefore, the combined effects of a
steep slope in the energy spectrum and the dispersion relation cause the
short-range spectral leakage to be most prominent for $k<1$.

The effect of a steep spectral slope — discussed further in section VI-C — can
be illustrated by considering that at, say, $k=0.6$ the spectrum shows
intensity at $\omega_{0}(0.6)$, but also intensity leaked from
$\omega_{0}(0.6+\delta k)$ and $\omega_{0}(0.6-\delta k)$, the latter of which
is much less than the former. The algorithm therefore finds a frequency
between $\omega_{0}(0.6)$ and $\omega_{0}(0.6+\delta k)$, which gives a non-
zero DSV (see also section IV).

Figure 5:  Doppler shift velocities (DSV) $\tilde{\bm{c}}$ in units of the
peak phase velocity for the parameter set of case-$\theta_{0}$-$\Delta\theta$
(see table I) extracted using the NSP method. The DSVs were extracted for 100
realizations with (right) and without (left) applying a Hann window
beforehand. The top and bottom graphs show the biases in DSVs in terms of the
average (top) and standard deviation (bottom) of the absolute value
$\tilde{c}=|\tilde{\bm{c}}|$. The solid and dashed line correspond to the
velocity resolutions implied by frequency and wavenumber resolution,
respectively.

###  VI-B Influence of the mean propagation direction and angular spread

To study the influence of the angular width on the mean and random bias in the
extracted DSVs, simulations were run with three propagation directions
$\theta_{0}\in\\{$$,$$,$$\\}$ and a range of spectral widths $\Delta\theta$,
of which two representative examples, ($$$,$$$), are shown in figure 5 (see
table I, case-$\theta_{0}$-$\Delta\theta$, for all parameters).

All cases show that an increase in spectral width $\Delta\theta$ significantly
decreases both mean and random biases across all wavenumbers. Results for
$\Delta\theta>$$$ are not shown here, as no notable further improvements were
observed.

The reduction in angular spread $\Delta\theta$ leads to an increase in
spurious DSVs $\tilde{\bm{c}}$, because the wave components of a narrow
spectrum are mostly influenced by the current component parallel to the mean
propagation direction $\theta_{0}$. The narrower a spectrum is, the smaller an
observable change in frequency due to the perpendicular current component
$\tilde{\bm{c}}_{\perp}$ becomes. Therefore, the influence of false spectral
intensity from any source (e.g. noise, aliasing, higher harmonics, spectral
leakage) outside the real spectrum can lead to strong perpendicular DSV
components when the spectrum is very narrow. For a more detailed explanation,
see Fig. 5 where we see the combined influence of a noisy spectrum and
spectral leakage.

Note, that all cases exhibit negligible spurious DSVs around $k\approx 1$.
This is simply a consequence of the spectral peak lying near this value.

The dependency of results on the mean propagation direction $\theta_{0}$ was
also investigated in a second manner: The videos were rotated numerically by
angles $\theta_{\text{rot}}$ up to $None$ before repeating the DSV extraction.
This is equivalent to rotating the camera in an experiment. Pairwise
comparison of results with same values of $\theta_{0}+\theta_{\text{rot}}$
show no significant difference.

Two conclusions can be drawn: First, the only effect $\theta_{\mathrm{rot}}$
can have on physical quantities is a trivial rotation of the coordinate
system, so any other $\theta_{\mathrm{rot}}$-dependencies can be traced back
to spectral leakage, which mainly occurs parallel to the video’s axes, thus
breaking rotational symmetry. Second, will see in section VI-E that when a
current is present, the measurements can indeed differ after camera rotation,
hence rotating the camera or data provides a useful check for spectral
leakage.

The use of a Hann window prior to DSV extraction greatly mitigates the effects
of spectral leakage, as discussed before, pushing even the narrow spectrum
cases down to sub-resolution (with respect to implied velocity resolutions
$\delta c_{\delta\omega}$ and $\delta c_{\delta k}$).

###  VI-C Influence of peakedness $\gamma$

Depending on how developed a sea is, the best fit for the frequency spectrum
uses a peak enhancement factor $\gamma$ between 1 and 7, where 3.3 is a
commonly used value for most applications [37]. Recently, Mazzaretto et al.
found that a global mean of $\gamma\approx 2.4$ is better suited [38].

Since the choice of $\gamma$ only affects the spectrum in a small range around
the peak, it also offers itself as a tool to examine the influence of steep
gradients in the spectrum. We therefore compare the DSVs for
$\gamma\in\\{1...7\\}$. (see table I, case-$\gamma$, for all parameters).

The results presented in figure 6 show that for spectra with small angular
spread, a more strongly peaked spectrum (higher value for $\gamma$) causes
increased biases (mean and random) across all wavenumbers except for a small
range around the spectral peak at $k\approx 1$, where the relative increase in
spectral intensity reduces the influence of spectral leakage from wavenumbers
$k\neq 1$. A wider angular spread of the spectrum, on the other hand,
mitigates this to a large degree, especially towards higher wavenumbers, where
the influence of spectral leakage from the peak is nearly eliminated and sub-
resolution biases are achieved. For angular spreads larger than
$\Delta\theta>$$$, the DSVs become independent of $\Delta\theta$, though this
also changes with spectral resolution.

Towards lower wavenumbers ($k<1$), without a Hann window, any extracted DSV is
unreliable, as spectral leakage dominates the result due to the exponential
decrease of ”real” spectral intensity towards $k=0.5$ and the steepness of the
dispersion relation (see previous discussion in section VI-A).

Once again the use of a Hann window improves the results so far that random
and mean biases are sub-resolution even for strongly peaked ($\gamma=7$),
narrow ($\Delta\theta=$$$) spectra, effectively eliminating the dependence on
$\Delta\theta$, and to a degree, on resolutions (see section VI-D).

Figure 6: Same as figure 5 but for a selection of parameter combinations from
the parameter set of case-$\gamma$ (see table I).

###  VI-D Influence of resolutions $\delta k$, $\delta\omega$

The range of spectral leakage is constant in terms of the number of pixels or
bins affected (see Fig. 1), and thus the leakage range in $k$ (and $\omega$)
is determined by the resolution $\delta k$ (and $\delta\omega$).

We therefore vary the spatial size $L$ and temporal duration $T$ of the videos
to change the frequency resolution $\delta\omega=1/T$ and wavenumber
resolution $\delta k=1/L$ to observe their influence on the DSV extraction
(see table I, case-$L$-$T$, for all parameters).

For a spectrum with an angular spread of $\Delta\theta=$$$ and a mean
propagation direction of $\theta_{0}=$$$ simulations were performed with a
range of video lengths $T=$580\text{\,}$$ and domain sizes $L=$520\text{\,}$$.
The resulting DSVs are shown in figure 7. (Note that the results for $T<10$
are not shown as they did not show relevant differences compared with results
for $T=10$)

Clearly, for untapered data, both improving $\delta k$ or $\delta\omega$
reduces the biases in DSVs, albeit improving the wavenumber resolution $\delta
k$ has a stronger influence.

Note how the maximum in biases in Fig. 7 moves to higher $k$ with decreasing
$T$. This is due to the characteristic function $G$ having a width $a$
proportional to frequency resolution, $a=4\delta\omega$. When this width is
reduced, the intensity that leaks from the spectral peak onto the cylinder at
$k>1$ can fall outside the reach of the characteristic function, thus reducing
or removing its influence.

When the data is tapered using a Hann window, this effect is mostly
eliminated, as the long-range spectral leakage is heavily suppressed.
Moreover, the influence of frequency resolution $\delta\omega$ is strongly
reduced (graphs of same color group together).

An exception to the rule of thumb that longer videos are always better can be
seen for the longest videos $T=80$ on the smallest domain $L=5$. Here, the
biases are actually greater than for the shorter cases with $T\leq 40$. This
is due to the characteristic function $G$ becoming narrow enough to not
encompass the width of spectral intensity in the $\omega$-direction. In one
dimension this would be proportional to $\delta\omega$, but in 2D or 3D can be
dominated by leakage in $k$, as this effective broadening in $\omega$ scales
with $\mathcal{O}(\partial\omega/\partial k\ \delta k)=\mathcal{O}(\delta
k/\sqrt{k})$.

Note, however, that the results presented in Fig. 7 all show biases well below
the implied velocity resolutions, provided the data is tapered using a Hann
window. If one has to choose between increasing the domain size or the video
duration, it is clear that an increase in domain size will give the most
benefit.

Figure 7: Same as figure 5 but with parameter set case-$L$-$T$ (see table I).
As in figure 5 The solid and dashed lines correspond to the implied velocity
resolutions from frequency and wavenumber resolution, respectively. From top
to bottom the resolution increases, i.e. $\delta w$ or $\delta k$ become
smaller.

###  VI-E Influence of a background current $U_{0}$

The cases shown above are without a background current. The presence of a
spatiotemporally uniform current $\bm{U}\neq 0$ breaks the rotational symmetry
of the dispersion relation, because of the additional, angular dependent term
$\bm{k}\cdot\bm{U}$ in equation (1). Like all velocities, $\bm{U}$ is given in
units of the phase velocity at the spectral peak (see Section III-A).

To see the effect of a current on the biases in the extracted DSVs, we assume
a (vertically and horizontally) constant current $\bm{U}$ with directions
$\theta_{U}=$$,$$,$$$ and current strengths in the range
$|\bm{U}|=$0.0250.15\text{\,}$$ (see table I, case-$U$, for all parameters).
The input wavenumber spectrum has an angular spread of $\Delta\theta=$$$, and
a peakedness of $\gamma=3.3$, representing a realistic scenario. The mean
propagation direction is held at $\theta_{0}=$$$, resulting in following,
crossing and opposing current, respectively.

The extracted DSVs are shown in figure 8. For $k>1$, we see a bias toward
wave-opposing DSVs, while for $k<1$, we see a bias toward wave-following DSVs.
This can be seen from over-/underestimations for $\bm{U}$
antiparallel/parallel to the wave propagation direction, respectively. The
absolute magnitude of errors is similar to the case of quiescent water (see
section VI-D). In the case of a pure cross-current, however, this also results
in a turning of the DSVs toward $\theta_{\tilde{\bm{c}}}=\pm$$$ for
wavenumbers around the spectral peak $k\lesssim 1.5$, as can be seen in the
third row of panels in figure 8. Particularly, where the spectral intensity
falls off to zero rapidly ($k\approx 0.5$), leakage from the peak of the
spectrum can dominate the signal, biasing the direction of DSVs towards the
wave propagation direction, i.e. $None$.

Tapering the data using a Hann window mitigates these biases in amplitude and
direction as in the case of quiescent water. For a cross-current, however, the
mean bias remains significant, i.e. on the order of 10% of the phase velocity
at the spectral peak and unusually also in the vicinity of the spectral peak.

The error in the direction $\theta_{\tilde{\bm{c}}}$ of DSVs for wavenumbers
below the peak is also reduced significantly for following and opposing
currents. In the case of a cross-current these errors remain significant,
($None$-$None$) if the current is relatively weak or the wavenumber goes
toward $k=0.5$.

As mentioned in section VI-B, turning the camera in an aerial measurement can
reveal the influence of spectral leakage (one could also rotate the resulting
images, bearing in mind that this is not in general a lossless operation and
might also reduce the field of view). Figure 9 shows the extracted DSVs for a
case with the same parameters as in figure 8, but with current and spectrum
rotated by $None$ (Similar results were obtained for
$1545\text{\,}\mathrm{\SIUnitSymbolDegree}$).

For the untapered data, we find even stronger mean biases after rotation,
especially for a cross-current. Here the DSVs around the spectral peak at
$k\approx 1$ drop below $50\text{\,}\mathrm{\char 37\relax}$ of the background
current, usually an unacceptable level of error. However, the use of a Hann
window again mitigates this, to the point where the results of different
rotation angles become virtually indistinguishable.

Figure 8: Same as figure 5 but with parameter set case-$U$ (see table I). The
dashed horizontal lines represent the three background velocity values.
Deviations from these are spurious. The third row shows the direction
$\theta_{\tilde{\bm{c}}}$ of the average Doppler shift velocity.

Figure 9: Same as figure 8 but with a $None$ rotated video. Note, how the case
of cross currents (blue right facing triangles) result in unusable DSVs
$\tilde{\bm{c}}$ especially around the spectral peak, when the data is
untapered. (top left)

##  VII Conclusions

We have investigated biases from spectral leakage in remote sensing of
currents from analysis of wave spectra. Apparent, spurious Doppler shifts in
the phase velocity are observed even in the absence of a current, and we
analyse how these depend on wave spectrum and the properties of the signal
processing procedure. Synthetically generated surface elevation data were used
to simulate a random sea state adhering to a JONSWAP wave spectrum [30] with a
$\cos^{2}$ directional distribution, resulting in cubes of data (videos).
These were subsequently analysed with different methods in common use to
extract the Doppler-shift velocities (DSV). Following reference [26], an
appropriate measure of random errors and biases are the “Doppler shift
resolutions” corresponding to the change in inferred velocity due to a shift
of one pixel in wavenumber or frequency.

A comparison of the normalized scalar product approach (e.g. [26]) and a
least-squares method for extracting DSV showed that the former is preferable
in all cases, and was therefore used for all subsequent analysis herein.

Assuming the simplest case of quiescent water (i.e. no current), a complex
interplay is found between wave-spectrum width and peakedness, and the wave-
vector and frequency resolution, together affecting the nature and extent of
spurious Doppler-shift “measurements”. Spectral leakage causes greater
problems when the wave spectrum is strongly peaked and highly directional so
that areas of the observed frequency-wave vector spectrum which are important
to Doppler-shift extraction have very low signal. Conversely, when the angular
spread is wide, $\Delta\theta>$$$ ($\Delta\theta$: full width of angular
distribution ), spurious DSVs are small, i.e., sub-resolution. For strongly
directional spectra $\Delta\theta<$$$, severe biases emerge, with amplitudes
on the order of the group velocity, depending on the JONSWAP peakedness
parameter $\gamma$. The biases are sensitive to resolution in frequency and
wave-number space, and especially in the absence of tapering (see below)
poorer resolution rapidly leads to unusable data for narrow and strongly
directional wave fields. Biases are most severe at low wavenumbers compared to
that at the spectral peak.

Tapering the video cubes with a 3D Hann window (e.g. [33]) lowers the biases
to the velocity resolution level implied by the wavevector and frequency
resolutions — $\delta c_{\delta k}$ and $\delta c_{\delta\omega}$ as defined
in equation (4) — or even below. Indeed, the effect of $\omega$-resolution is
mostly removed for tapered data (note that although biases are now sub-
resolution, the resolution itself will eventually be too poor for purpose).
This implies that in data acquisition one should prioritize large areas rather
than longer time series if spectral leakage is a problem.

The effect of spectral leakage is most pronounced in the $k_{x}$ and $k_{y}$
directions in the spectrum, and hence depends on the angle of propagation
relative to these. The biases increase towards an angle of $\theta_{0}=$$$ for
all wavenumbers outside a small range around the spectral peak $k=k_{p}$. This
dependence on $\theta_{0}$ is also mostly removed by tapering the data.
Clearly the observed current velocity cannot depend on which way the camera is
held, meaning that comparison with results when the video is rotated, say,
$None$ could give a simple indication of the severity of spectral leakage
problems.

For the case of a constant background current $\bm{U}$ with strengths up to
$0.15\,c_{0}(k_{p})$ a strong dependence on the angle between the current and
the waves is observed. While the magnitudes of the spurious DSVs are mostly
smaller than the implied velocity resolutions, we find a significant bias in
the direction of the DSVs around $k=k_{p}$ when the waves propagate
perpendicular to the current. This is also mitigated by tapering the data, but
not removed in the case of a cross current.

###  VII-A Recommendations for mitigation

Summarising the outcome of our analysis from a practical viewpoint we offer
the following considerations to mitigate the errors and biases related to
spectral leakage in remote sensing of currents from observed wave spectra. In
the extraction of Doppler shift velocities from the spectrum, the commonly
used least-squares method is not recommended except if calculation cost is a
severe restriction; a normalized scalar product procedure gives universally
better results (other methods are also in use, but were not tested). Tapering
the spatio-temporal data with a 3D Hann window greatly reduces the mean and
random biases and their dependence on spectral shape, spectral resolution, and
camera orientation. When errors due to spectral leakage are suspected,
rotating the camera (either the actual camera or the resulting images) by some
intermediate angle before analysis (we used $None$) and comparing results
could reveal whether long range spectral leakage is causing spurious results,
because leakage mainly occurs along the axes of the images. We find that the
wavenumber resolution plays a more important role in the DSV biases than
frequency resolution. The influence of the latter can be nearly eliminated by
tapering the data. Increasing the spatial domain size to improve wavenumber
resolution yields the largest improvement and should therefore be prioritized
over longer time series if spectral leakage is a concern. The Doppler shift
resolutions $\delta c_{\delta\omega}$ and $\delta c_{\delta k}$ defined in
equation (4) are useful as conservative measures of errors and biases due to
limited resolution and spectral leakage.

##  VIII Acknowledgments

The authors would like to thank Drs. Susanne Støle-Hentschel, Luc Lenain and
Nick Pizzo for discussions.

###  -A Comparison of window functions

There exists a wide variety of window functions employed across disciplines,
as discussed in e.g. [33]. The suppression of the long range spectral leakage
is the most important beneficial effect of using a window function, but
stronger suppression comes at the cost of increasing the short range leakage.
For example, the Blackman window has a central lobe that is about 20% wider
than that of the Hann window, as can be seen in figure 1. To examine its
detrimental effect, the extraction of DSVs has been repeated with a selection
of commonly used window functions. As the results are very similar we only
show the results for the Hann window

$w_{\text{}}(t)=0.5-0.5\,\cos\left(2\pi\,t/T\right),$ (17)

the approximate Blackman window,

$w_{\text{}}(t)=0.42-0.5\,\cos\left(2\pi\,t/T\right)+0.08\,\cos\left(4\pi\,t/T\right),$
(18)

and the Kaiser-Bessel window

$w_{\text{}}(t)=\frac{1}{T}I_{0}\left(b\sqrt{1-(2\pi\,t/T)^{2}}\right),$ (19)

with $I_{0}$ being the zeroth order modified Bessel function of the first kind
and $b$ a free parameter controlling the width of the central lobe. In figure
10 the extraced DSVs in terms of their main and random bias are shown. The
differences in the results are insignificant, except for the Blackman window
showing systematically higher biases due to its wider central lobe.

Figure 10: Doppler shift velocities (DSV) in units of the peak phase velocity
for a normal parameter combination ($\theta_{0}=$$$, $\Delta\theta=$$$,
$\gamma=3.3$, $L=10$, $T=20$, $U=0$) obtained using the NSP method after
tapering the data with different window functions. The top and bottom graph
show the DSVs $\tilde{\bm{c}}$ in terms of the mean (top) and standard
deviation (bottom), respectively.

###  -B Details of numerical implementation

Here follow further details on the implementation of the NSP method in section
III-C.

####  -B1 NSP, LS and discretized data

Because of discretization, obtaining the spectral intensity $P_{i}$ on a
cylinder surface needs to be replaced with the spectral intensity in a volume
around it, i.e. a cylinder shell, containing a wavenumber bin around $k_{i}$:

$F_{i}(\mathbf{k},\omega)=\begin{cases}\sqrt{P(\mathbf{k},\omega)},&\text{if
}|k-k_{i}|\leq\tilde{\Delta}k\\\ 0,&\text{otherwise}\end{cases}.$ (20)

This also renders the characteristic function $G_{i}$ (equation (14)) a
function of $k$. The integrals $\left\langle...\right\rangle$ in equation 16
then imply an additional integration over $k$ within the bin
$k_{i}\pm\tilde{\Delta}k$.

The concept of the NSP method stays the same, with two details added: First,
since the optimization parameter pair $\tilde{\bm{c}}$ is assumed constant
within a wavenumber bin, the extracted DSV is a weighted average within that
bin, effectively smoothing the function $\tilde{\bm{c}}(k)$. This also holds
for the LS method, but does not change equation (21), as it only increases the
number of triplets $(k_{x,j},k_{y,j},\omega_{j})$. Second, a new free
parameter is introduced with the bin width $2\tilde{\Delta}k$ that needs to be
chosen carefully. In this work, we use $\tilde{\Delta}k=2\delta k$, with
$\delta k=1/L$ the wavenumber resolution, which is a compromise between
increased smoothing (too large $\tilde{\Delta}k$) and strong noise, occurring
when too few pixels of the spectrum lie within a bin. (One could employ an
interpolation scheme to circumvent this; the algorithm used herein simply
masks the data, see equation 20)

Since the data are effectively averaged over $k$ with a running average of
width $2\tilde{\Delta}k$, the DSVs were extracted such that two consecutive
shells have an overlap of $3/4$, i.e. $k_{i+1}-k_{i}=\delta k$.

####  -B2 Width of the characteristic function

As mentioned in section V the characteristic function $G$ (equation 14) used
in fitting the dispersion relation to the measured spectrum contains one free
parameter, $a$,that determines its width. The choice of this parameter is
somewhat delicate, as too small a value causes single intensity pixels (often
outliers) or small high intensity regions to dominate the determination of the
best fit. This is especially problematic for small ($k<1$) wavenumbers, when
the spectral leakage from a noise-enhanced pixel near the spectral peak at
$k\approx 1$ causes a small high-intensity region on the cylinder shell,
whereas the “real” spectral intensity, not originating from spectral leakage,
is strongly broadened, thus having a larger total intensity but a smaller
maximum intensity. In this case, the algorithm effectively ignores the real
spectral intensity and leads to huge biases. On the other hand, too large a
value for $a$ also leads to an increased influence of spectral leakage, as a
too wide Gaussian is insensitive to shifts in its position. A shift with no,
or only small penalty to the overlap with the real spectral intensity, that
increases the overlap with intensity from leakage, is therefore more likely,
also leading to biases.

To decide on a good compromise we fit the characteristic function to the
Fourier transform of a Hann window, with $a$ as a free parameter. The result
is $a\approx 1.62\,\delta\omega$. To ensure the width is not too small this
value is roughly doubled to $a=4\,\delta\omega$. (Values for $a/\delta\omega$
between 2 and 6 were tried as well, but the results were most stable between 3
and 5.)

Note, that if spectral leakage occurred only in the $\omega$ direction,
$a=\delta\omega$ would be the ideal choice. However, as this is not the case,
broadening in the $\omega$ direction originates not only from leakage in the
$\omega$ direction, but also from leakage in the $k$-direction. The steeper
the dispersion relation $\omega_{\text{DR}}(k)$, the more this effective
$\omega$-leakage from $k$-leakage increases. A possible improvement of the
used NSP algorithm would therefore be to use an adaptive $a$, that increases
with $|\partial\omega/\partial k|$. In the interest of limiting the parameter
space, it was deemed necessary to stick to a single value for $a$.

###  -C Comparison NSP and a basic LS algorithm

Simple least squares based methods have by now been discarded in most
applications. However, their low computational cost and ease of implementation
are benefits to consider. We therefore repeated the DSV extraction using LS
and compare the results and performance. The least squares method used, like
the NSP method described in section III-C, first singles out the spectral
intensity on a cylinder surface with radius $k$. We discard data points with
an intensity $P<0.2$ (after normalization) and obtain a list of
$(k_{x,j},k_{y,j},\omega_{j})$. On these, the cost function

$C(\tilde{\bm{c}})=\sum_{j}(\omega_{DR}(\bm{k}_{j};\tilde{\bm{c}})-\omega_{j})^{2}$
(21)

is minimized to obtain the DSV $\tilde{\bm{c}}$. Repeating this for a list of
wavenumbers $k$ yields the desired value of $\tilde{\bm{c}}$ for each value of
$k$. As in the NSP method and, the minimization/optimization step is performed
using the Nelder-Mead simplex method [35].

We find that the LS method runs significantly faster (up to a factor of $3$);
clearly the difference in cost will depend on the hardware used as well as the
implementation. In terms of accuracy and precision, we find that the LS method
performed consistently worse, and never better, than the NSP method for
practical purposes. Figure 11 shows an example result for a test case in
quiescent water (see table I, case-NSP-LS, for all parameters). As can be
seen, the DSVs obtained via the LS method can show both a mean and random bias
that exceed the implied velocity resolutions $\delta c_{\delta k}$ and $\delta
c_{\delta\omega}$ for most wavenumbers, while the NSP method delivers sub-
resolution DSVs for all $k>1$. It is worth pointing out, that for wavenumbers
$k<1$ the NSP performs worse than the LS method. However, the random bias in
that range is so large as to make both methods unusable. Note that a more
advanced LS method may perform with similar accuracy and precision, as is
indicated in e.g. [10]. This would, presumably, lead to a computational cost
similar to that of the NSP method, eliminating the advantage. An iterative LS
method has been compared with NSP by Huang et al. (2012). For waves on
currents without vertical shear the methods performed similarly.

Figure 11: Doppler shift velocities (DSV) in units of peak phase velocity for
the parameter set case-NSP-LS (see table I) extracted using the NSP and a
least square method. The top and bottom graph show the DSVs $\tilde{\bm{c}}$
for quiescent water in terms of the mean (top) and standard deviation
(bottom). The solid and dashed line correspond to the velocity resolutions
implied by frequency and wavenumber resolution, respectively.

## Data availability

The data that support the findings of this study are available from the
authors upon reasonable request.

## References

  * [1] B. K. Smeltzer and S. Å. Ellingsen, _Current mapping from the wave spectrum_. Institution of Engineering and Technology (IET), 2021, ch. 15, pp. 357–378. [Online]. Available: https://arxiv.org/abs/2205.05373
  * [2] R. H. Stewart and J. W. Joy, “HF radio measurements of surface currents,” _Deep-Sea Res._ , vol. 21, pp. 1039–1049, 1974.
  * [3] E.-C. Ha, “Remote sensing of ocean surface current and current shear by HF backscatter radar,” Department of Electrical Engineering, Stanford University, Tech. Rep. no. D415-1, 1979.
  * [4] I. R. Young and W. Rosenthal, “A three-dimensional analysis of marine radar images for the determination of ocean wave directionality and surface currents,” _J. Geophys. Res._ , vol. 90, pp. 1049–1059, 1985.
  * [5] D. M. Fernandez, J. F. Vesecky, and C. Teague, “Measurements of upper ocean surface current shear with high-frequency radar,” _J. Geophys. Res._ , vol. 101, no. C12, pp. 28,615–28,625, 1996.
  * [6] K.-W. Gurgel and H.-H. Essen, “On the performance of a shipborne current mapping hf radar,” _IEEE journal of Oceanic Engineering_ , vol. 25, no. 1, pp. 183–191, 2000.
  * [7] C. C. Teague, J. F. Vesecky, and Z. R. Hallock, “A comparison of multifrequency HF radar and ADCP measurements of near-surface currents during COPE-3,” _IEEE J. Oceanic. Eng._ , vol. 26, pp. 399–405, 2001.
  * [8] R. Gangeskar, “Ocean current estimated from X-band radar sea surface images,” _IEEE Trans. Geosci. Remote Sens._ , vol. 40, pp. 783–792, 2002\.
  * [9] F. Ardhuin, L. Marié, N. Rascle, P. Forget, and A. Roland, “Observation and estimation of Lagrangian, Stokes, and Eulerian currents induced by wind and waves at the sea surface,” _Journal of Physical Oceanography_ , vol. 39, no. 11, pp. 2820–2838, 2009.
  * [10] W. Huang and E. Gill, “Surface current measurment under low sea state using dual polarized X-band nautical radar,” _IEEE J. Sel. Topics Appl. Earth Observ. in Remote Sens._ , vol. 5, pp. 1868–1873, 2012.
  * [11] K. Hessner, K. Reichert, J. C. N. Borge, C. L. Stevens, and M. J. Smith, “High-resolution X-Band radar measurements of currents, bathymetry and sea state in highly inhomogeneous coastal areas,” _Ocean Dyn._ , vol. 64, pp. 989–998, Jul. 2014.
  * [12] B. Lund, H. Graber, H. Tamura, C. C. III, and S. Varlamov, “A new technique for the retrieval of near-surface vertical current shear from marine X-band radar images,” _J. Geophys. Res.: Oceans_ , vol. 120, pp. 8466–8484, 2015\.
  * [13] W. Huang, R. Carrasco, C. Shen, E. W. Gill, and J. Horstmann, “Surface current measurements using x-band marine radar with vertical polarization,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 54, no. 5, pp. 2988–2997, 2016.
  * [14] B. Lund, B. K. Haus, J. Horstmann, H. C. Graber, R. Carrasco, N. J. M. Laxague, G. Novelli, C. M. Guigand, and T. M. Özgökmen, “Near-surface current mapping by shipboard marine X-band radar: a validation,” _J. Atmos. Oceanic Technol._ , vol. 35, pp. 1077–1090, 2018.
  * [15] R. Gangeskar, “Verifying high-accuracy ocean surface current measurements by X-band radar for fixed and moving installations,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 56, no. 8, pp. 4845–4855, 2018.
  * [16] S. Støle-Hentschel, B. K. Smeltzer, and Y. Toledo, “On the measurement of the dispersion relation by a radar and the implication on the current retrieval,” 2023, Submitted manuscript. [Online]. Available: https://arxiv.org/abs/2303.12181
  * [17] J. P. Dugan, H. H. Suzukawa, C. P. Forsyth, and M. S. Farber, “Ocean wave dispersion surface measured with airborne ir imaging system,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 34, no. 5, pp. 1282–1284, 1996.
  * [18] C. J. Zappa, M. L. Banner, H. Schultz, J. R. Gemmrich, R. P. Morison, D. A. LeBel, and T. Dickey, “An overview of sea state conditions and air‐sea fluxes during radyo,” _Journal of Geophysical Research_ , vol. 117, no. C7, 2012.
  * [19] N. J. M. Laxague, B. K. Haus, D. G. Ortiz-Suslow, C. J. Smith, G. Novelli, H. Dai, T. Özgökmen, and H. C. Graber, “Passive optical sensing of the near-surface wind-driven current profile,” _J. Atmos. Oceanic Technol._ , vol. 34, pp. 1097–1111, 2017.
  * [20] M. Streßer, R. Carrasco, and J. Horstmann, “Video-based estimation of surface currents using a low-cost quadcopter,” _IEEE Geoscience and Remote Sensing Letters_ , vol. 14, no. 11, pp. 2027–2031, 2017.
  * [21] J. P. Dugan, C. C. Piotrowski, and J. Z. Williams, “Water depth and surface current retrievals from airborne optical measurements of surface gravity wave dispersion,” _J. Geophys. Res._ , vol. 106, pp. 16 903–16 915, 2001.
  * [22] J. P. Dugan and C. C. Piotrowski, “Surface current measurements using airborne visible image time series,” _Remote Sens. Environ._ , vol. 84, pp. 309–319, 2003.
  * [23] M. Yurovskaya, N. Rascle, V. Kudryavtsev, B. Chapron, L. Marié, and J. Molemaker, “Wave spectrum retrieval from airborne sunglitter images,” _Remote sensing of Environment_ , vol. 217, pp. 61–71, 2018.
  * [24] L. Lenain, B. K. Smeltzer, N. Pizzo, M. Freilich, L. Colosi, S. A. Ellingsen, L. Grare, H. Peyriere, and N. M. Statom, “Airborne observations of surface winds, waves and currents from meso to submesoscales,” _Geophys. Res. Lett_ , vol. 50, p. e2022GL102468, 2023.
  * [25] J. Campana, E. J. Terrill, and T. de Paolo, “A new inversion method to obtain upper-ocean current-depth profiles using X-band observations of deep water waves,” _J. Atmos. Oceanic Technol._ , vol. 34, pp. 957–970, 2017.
  * [26] B. K. Smeltzer, E. Æsøy, A. Ådnøy, and S. Å. Ellingsen, “An improved method for determining near‐surface currents from wave dispersion measurements,” _J. Geophys. Res.: Oceans_ , vol. 124, pp. 8832–8851, 2019\.
  * [27] C. M. Senet, J. Seeman, and F. Ziemer, “The near-surface current velocity determined from image sequences of the sea surface,” _IEEE Trans. Geosci. Remote Sens._ , vol. 39, no. 3, pp. 492–505, 2001.
  * [28] F. Serafino, C. Lugni, and F. Soldovieri, “A novel strategy for the surface current determination from marine X-band radar data,” _IEEE Geosci. Remote Sens. Lett._ , vol. 7, pp. 231–235, 2010.
  * [29] B. Lund, B. Haus, H. Graber, J. Horstmann, R. Carrasco, G. Novelli, C. Guigand, S. Mehta, N. Laxague, and T. Özgökmen, “Marine x‐band radar currents and bathymetry: An argument for a wave number‐dependent retrieval method,” _Journal of Geophysical Research: Oceans_ , vol. 125, no. 2, p. e2019JC015618, 2020.
  * [30] K. F. Hasselmann, T. P. Barnett, E. Bouws, H. Carlson, D. E. Cartwright, K. Eake, J. Euring, A. Gicnapp, D. Hasselmann, and P. Kruseman, “Measurements of wind-wave growth and swell decay during the joint north sea wave project (JONSWAP),” _Ergänzungsheft zur Deutschen Hydrographischen Zeitschrift, Reihe A_ , vol. 12, pp. 1–95, 1973.
  * [31] S. A. Ellingsen and Y. Li, “Approximate dispersion relations for waves on arbitrary shear flows,” _J. Geophys. Res.: Oceans_ , vol. 122, pp. 9889–9905, 2017.
  * [32] D. A. Lyon, “The discrete Fourier transform, part 4: spectral leakage,” _Journal of object technology_ , vol. 8, no. 7, 2009.
  * [33] A. Nuttall, “Some windows with very good sidelobe behavior,” _IEEE Transactions on Acoustics, Speech, and Signal Processing_ , vol. 29, no. 1, pp. 84–91, 1981.
  * [34] H. Socquet-Juglard, K. Dysthe, K. Trulsen, H. E. Krogstad, and J. Liu, “Probability distributions of surface gravity waves during spectral changes,” _Journal of Fluid Mechanics_ , vol. 542, p. 195–216, 2005.
  * [35] J. C. Lagarias, J. A. Reeds, M. H. Wright, and P. E. Wright, “Convergence properties of the nelder-mead simplex method in low dimensions,” _SIAM J. Optim._ , vol. 9, pp. 112–147, 1998.
  * [36] L. Lenain, B. K. Smeltzer, N. Pizzo, M. Freilich, L. Colosi, S. A. Ellingsen, L. Grare, H. Peyriere, and N. Statom, “Airborne remote sensing of upper-ocean and surface properties, currents and their gradients from meso to submesoscales,” _Geophysical Research Letters_ , vol. 50, no. 8, p. e2022GL102468, 2023, e2022GL102468 2022GL102468. [Online]. Available: https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2022GL102468
  * [37] Y. Goda, _Random seas and design of maritime structures_ , 3rd ed. World Scientific, 2010.
  * [38] O. M. Mazzaretto, M. Menéndez, and H. Lobeto, “A global evaluation of the JONSWAP spectra suitability on coastal areas,” _Ocean Engn._ , vol. 266, p. 112756, 2022.