A Self-Adaptive Method for Mapping Coastal Bathymetry On-The-Fly from Wave Field Video

Abstract

Mapping coastal bathymetry from remote sensing becomes increasingly more attractive for the coastal community. It is facilitated by a rising availability of drone and satellite data, advances in data science, and an open-source mindset. Coastal bathymetry, but also wave directions, celerity and near-surface currents can simultaneously be derived from aerial video of a wave field. However, the required video processing is usually extensive, requires skilled supervision, and is tailored to a fieldsite. This study proposes a video-processing algorithm that resolves these issues. It automatically adapts to the video data and continuously returns mapping updates and thereby aims to make wave-based remote sensing more inclusive to the coastal community. The code architecture for the first time includes the dynamic mode decomposition (DMD) to reduce the data complexity of wavefield video. The DMD is paired with loss-functions to handle spectral noise and a novel spectral storage system and Kalman filter to achieve fast converging measurements. The algorithm is showcased for fieldsites in the USA, the UK, the Netherlands, and Australia. The performance with respect to mapping bathymetry was validated using ground truth data. It was demonstrated that merely 32 s of video footage is needed for a first mapping update with average depth errors of 0.9–2.6 m. These further reduced to 0.5–1.4 m as the videos continued and more mapping updates were returned. Simultaneously, coherent maps for wave direction and celerity were achieved as well as maps of local near-surface currents. The algorithm is capable of mapping the coastal parameters on-the-fly and thereby offers analysis of video feeds, such as from drones or operational camera installations. Hence, the innovative application of analysis techniques like the DMD enables both accurate and unprecedentedly fast coastal reconnaissance. The source code and data of this article are openly available.

1. Introduction

Observations of near-shore hydrodynamics and bathymetry are used for various purposes: to study and manage the coast [1,2,3], to update early warning systems [4], to monitor swimmer safety [5,6,7], for dredging-and-dumping surveillance [8], and military landing operations [9]. An observation of hydrodynamics or bathymetry with areal coverage, a map, is thereby often beneficial if not a prerequisite to recognize relevant spatial details. A straightforward approach to map coastal hydrodynamics and bathymetry is via video-based remote sensing of a wave-field. In comparison to in-situ measurements, video-based remote sensing is less accurate; however, data acquisition is also less labor-intensive, and measurements have high spatial coverage by default.

Different instruments and video-processing methods are used to map hydrodynamics and bathymetry. In terms of instruments, videos may be recorded with stationary cameras [10], aircrafts [11], UAVs/drones [12,13,14,15,16], (navigational) X-Band radars [17,18,19,20,21], or satellites [22,23]. In terms of video-processing methods, various types exist, e.g., [17,20,24,25,26,27,28,29,30]. Most widely used among the coastal remote-sensing community are video-processing methods that analyse wave-dispersion properties. They prove applicable to different instruments [31,32] and allow to estimate several coastal parameters. Consecutive frames of a wave-field recording are scanned [20,33] to extract dominant wave frequencies, $\omega$ (rad/s), and associated wave lengths and directions via wave-number vectors, $\mathit{k}$ with [$k_x$, $k_y$] (rad/m) [18,21,34,35,36], but also wave celerity vectors, $\mathit{c}$ with [$c_x$, $c_y$] (m/s) [37,38,39], near-surface current vectors, $\mathit{U}$ with [u, v] (m/s) [11,20,40,41], and the apparent depth, d (m) [9,11,17,42,43,44]. While $\omega$, $\mathit{k}$, and $\mathit{c}$ are retrieved directly from wave patterns [45,46] or their corresponding frequency–wave-number spectra, e.g., [37], (Figure 1 left), $\mathit{U}$ and d are retrieved indirectly, albeit simultaneously [2,30,47,48,49], by matching frequency–wave-number spectra with a theoretical wave model (Figure 1 right). This study focused on mapping $\mathit{c}$, d, and $\mathit{U}$. Vector fields of $\mathit{k}$ are directly coupled to vector fields of $\mathit{c}$ via corresponding $\omega$.

The wave model is typically given by the Doppler-shifted linear dispersion relationship, Equation (1):

$${\omega }_{model}(d,\mathit{U})=\sqrt{g\left|\mathit{k}\right|tanh\left(\right|\mathit{k}\left|d\right)}+\mathit{U}·\mathit{k},$$

(1)

where ${\omega }_{model}$ (rad/s) is the wave model frequency, which is adjusted with unknown parameters d and $\mathit{U}$ until it optimally matches (and therewith explains) the observed $\omega$. The gravitational acceleration is given by g (m/s²). Note that $\mathit{U}$ is derived as the depth-averaged velocity of a depth-uniform current profile; however, in natural settings with depth-varying current profiles, $\mathit{U}$ represents the weighted average of velocities in the upper layer and is therefore typically referred to as the near-surface current vector [20,41]. Advanced techniques also allow to resolve the underlying current profile [50,51]. Yet, the advantage of condensing the profile to a bulk vector $\mathit{U}$ is a comparably simple expression for the Doppler shift ($+\mathit{U}·\mathit{k}$ in Equation (1)), by which $\mathit{U}$ can be estimated together with d [47,49]. Moreover, near-surface currents can then directly be visualized through maps of vector fields, e.g., [40,52].

Two aspects need to be considered for estimating coastal parameters from video: first, the optimal extraction of wave length and wave frequency characteristics from optical spectra, and second, the optimal estimation of $\mathit{c}$, d, $\mathit{U}$ from the found wave-length and -frequency data. Both aspects are successively treated in Section 1.1 and Section 1.2.

1.1. Optimizing Wave-Number–Frequency Extraction from Optical Spectra

One of the core challenges of wave dispersion-based video processing is to robustly identify the wave-number–frequency signature of gravity waves (${\mathit{k}}^{\uparrow },{\omega }^{\uparrow }$, where ↑ denotes the “gravity wave signature”). Different strategies [30,37,43,53] follow different pathways to do this (Figure 2). All strategies use greyscale video as a basis and inspect local subdomains to capture spatial variation in ${\mathit{k}}^{\uparrow }$. Searching for an optimal pathway to retrieve ${\mathit{k}}^{\uparrow },{\omega }^{\uparrow }$, the alternatives are briefly presented (Figure 2). Appendix A contains a more detailed review of how successive transformations in different strategies lead to the retrieval of ${\mathit{k}}^{\uparrow },{\omega }^{\uparrow }$ (Figure 2, arrows) and summarizes some benefits and drawbacks in

Table A1. Pros and cons of different strategies used to retrieve local gravity wave ${\mathit{k}}^{\uparrow },{\omega }^{\uparrow }$ signatures from video of a moving wave-field.

Strategy	Pros	Cons
e.g., [37]	Simple approach Spectral data: 2 × Nyquist limit	Comp. load: full-dimensional video
[30]	LOCPI: sugg. higher localisation Spectral data: 2 × Nyquist limit	comp. load: full-dimensional video
[43]	LOCPI: sugg. higher localisation Freedom local pixel sampling	Comp. load: full-dimensional video approximate ${\omega }^{\uparrow }$
[53]	G$\widetilde{\mathrm{OC}}$PI: sugg. higher localisation Freedom global pixel sampling Globally coherent wave patterns Predict subdomain sizes from ${\omega }^{\uparrow }$ Comp. load: reduce dimensionality	Sensitive to mixing Fourier modes
proposed	GOCPI: sugg. higher localisation Freedom global pixel sampling Globally coherent wave patterns Predict subdomain sizes from ${\omega }^{\uparrow }$ Comp. load: reduce dimensionality

.

Local cut-outs from the video form cubes in space time ($x,y,t$; in short $\mathit{x},t$) and can directly be used as analysis subdomains. Decomposing a cube reveals the wave constituents of different frequency, either directly in wave-number–frequency space (in $\mathit{k},\omega$) [37] (Figure 2, left pathway) or via their complex-valued phase images (in $\mathit{x},\omega$) [30] (Figure 2, centre pathway). Separation from the spectral noise floor finally yields ${\mathit{k}}^{\uparrow },{\omega }^{\uparrow }$ (Figure 2 bottom). Associated to one frequency component, a phase image is referred to as a one-component phase image [30] and shows the distinct wave pattern at a certain wave frequency. The suggested benefit of using one-component phase images is better localization of ${\mathit{k}}^{\uparrow }$. Instead of deriving local one-component phase images (LOCPI) from video cut-outs (Figure 2, centre pathway), global one-component phase images (GOCPI) of the full video domain can be derived [53] (Figure 2, right pathway). These GOCPI are then further analysed in local subdomains. Independent of using LOCPI [30], GOCPI [53], or both [43], the final phase structure should be spatially coherent to get ${\mathit{k}}^{\uparrow }$.

The benefits of GOCPI are that they can be generated taking global spatial coherence into account [53], while additionally the dimensionality of the video data can be reduced [53]. Reducing dimensionality is an important asset as it offers a reduction in the required computational working memory and thereby an increase in computational speed. For this purpose, a singular value decomposition (svd)(Appendix A, Equation (A1)) can be employed, which is a dimensionality reduction technique that is controlled by the variance in the video. The svd splits the video into modes, consisting of spatial and temporal structures. The first $r\sim O\left(1\right)-O\left(10\right)$ modes typically describe most of the video’s variance (see ${\sigma }_j$ in Appendix A, Equation (A1)) and thereby capture its essence. The remaining modes confine noise and can be discarded, which reduces the data load. Simarro et al. [53] directly use the spatial structures from svd as GOCPI. However, a spatial structure from svd may contain mixed wave patterns of different frequencies (see [54], Figure 6) such that its corresponding temporal structure is not a pure oscillation. Hence, using spatial structures from svd as GOCPI is suboptimal.

Aiming to map coastal parameters on-the-fly, ${\mathit{k}}^{\uparrow },{\omega }^{\uparrow }$ must be retrieved at minimal computational cost. A strategy similar to Simarro et al. [53] is therefore desirable but with pre-knowledge of oscillating spatial structures when generating GOCPI. This is what the dynamic mode decomposition (DMD) aims to do (Appendix A, Table A1).

The DMD is a mathematical procedure designed to reduce the dimensionality of video data by identifying dominant oscillating spatial structures. Different variants of DMD have been developed, e.g., [54,55,56,57,58,59,60], employed in various fields of science amongst which electrical engineering [61,62], system and control applications [57], neuroscience [63], but also physical oceanography [58] and coastal engineering [64].

Here, the spatial structures from DMD represent phase images of waves with different periods. The DMDs’ foundation lies in the assumption that snapshots in a video (sequence of frames) are related to each other through a linear dynamic system of oscillatory components. This means that the benefits of Simarro et al. [53], offering GOCPI and the possibility for dimensionality reduction can be complemented with a guarantee that GOCPI oscillate in time. What is more, the oscillation frequencies are inherently found.

To summarize, the DMD can be used to reduce video data to a set of GOCPI (Figure 2, right pathway) but with guaranteed oscillatory time behaviour. These GOCPI are the basis to find local ${\mathit{k}}^{\uparrow },{\omega }^{\uparrow }$ signatures (Figure 2 left, purple dots): while (global) ${\omega }^{\uparrow }$ are known from the DMD, ${\mathit{k}}^{\uparrow }$ still need to be locally deduced. This can be done in various ways, such as via characteristic spatial phase differences in sub-domains [30,43,53] or alternatively via (computationally cheap) spatial 2D-FFTs or particle image velocimetry (PIV) [23,65]. Derived local ${\mathit{k}}^{\uparrow },{\omega }^{\uparrow }$ form the basis to retrieve maps of wave celerity per GOCPI frequency, via $\mathit{c}\left({\omega }^{\uparrow }\right)={\omega }^{\uparrow }/{\mathit{k}}^{\uparrow }$[66], but also maps of depths and near surface currents d, $\mathit{U}$ through Equation (1).

1.2. Optimizing Depth and Near-Surface Current Estimates

Pursuing the best approximations to d and $\mathit{U}$, optimization problems using Equation (1) (with/without Doppler shift $+\mathit{U}·\mathit{k}$) have been stated differently, with some as (non-linear) least-squares (LS) minimization problems [20,30,43] and others as maximization problems of a normalized scalar product (NSP) [19,49]. These schemes have been formulated to handle an abundance of spectral information from Fourier decompositions. However, applying dimensionality reduction techniques to video data concentrates spectral information to its essence offering far fewer data points for an optimization of d, $\mathit{U}$. As such, inaccuracies and outliers in the spectral data gain negative influence on the solution. Standard LS minimization is not fit to handle this issue and with little data to fit, NSP maximization using a Heaviside step function is a crude approach. Outlier contamination is a common issue in applications such as the training of neural networks [67], because reducing the relatively large residual of an outlier is more attractive in minimizing the cost function than reducing the small residual of an inlier. In standard LS optimization, residuals are squared, whereby the impact of outliers on the solution is disproportionately large. However, LS problem statements have the benefit that many different methods have been developed to solve them [68], among which the Levenberg–Marquardt method is often applied to invert d, $\mathit{U}$ e.g., [2,43]. It is therefore attractive to adapt LS formulations such that they can handle outliers. This is achieved using various kinds of loss functions [69], among which the Cauchy (also called Lorentzian) loss function is an effective type.

Higher accuracy in d, $\mathit{U}$ can also be achieved by successively producing estimates and converging the results with a Kalman filter [29,43,70]. The quality of an estimate is thereby typically judged based on the sensitivity of the fit to the data. If the data are plagued by inaccuracies and outliers, this sensitivity increases, making the Kalman filter a suitable tool to mitigate their influence.

To summarize, while loss-functions can increase the quality of an individual d, $\mathit{U}$ estimate, the application of a Kalman filter increases the quality over several d, $\mathit{U}$ estimates. Both techniques can be used simultaneously.

1.3. Outlook of This Study

This study is about a self-adaptive, robust method to map $\mathit{c}$, d and $\mathit{U}$ on-the-fly from video of a wave-field. The DMD plays a key role in making the video-processing algorithm self-adaptive to the data and computationally fast: it reduces (video) data complexity, finds the dominant wave-components, and allows self-adaptive sampling schemes, which cause the standardly used computational cubes to instead become pyramids for optimal localization. For the optimization of d, $\mathit{U}$, the algorithm implements a loss-function to handle spectral outliers, which seemingly counteracts the error of overestimating $\mathit{U}$ and also reduces depth errors. The two errors are commonly interlinked [2,47,48,70]. The algorithm temporarily stores spectral data and employs Kalman filtering for quick convergence of measurements, and it comprises a set of rules and filters for autonomous decision making such that the algorithm does not need to be tuned to the field site. In summary, the algorithm’s main elements include:

• Data reduction and retrieval of wave components through DMD.
• Wave-component-dependent subdomains using pyramid cells.
• Counteracting spectral outliers and errors in d,$\mathit{U}$ with loss functions.
• Fast convergence of d,$\mathit{U}$ and recognition of current refraction through temporary spectral data storage.
• Additional fast convergence of d,$\mathit{U}$ using Kalman filtering.

Section 2 describes the workflow of the proposed video-processing algorithm, which includes explanations of the core principles of the DMD and the workflow of the algorithm. Section 3 describes the field sites and data. In Section 4, the algorithm is validated for measuring near-shore bathymetry (i.e., d) based on four field sites in the USA, the UK, the Netherlands, and Australia. A qualitative discussion on the algorithm’s ability to measure hydrodynamics (i.e., $\mathit{c}$, $\mathit{U}$) follows in Section 5, together with a discussion on the algorithm’s potential for mapping on-the-fly. The findings are concluded in Section 6.

2. Method

Ahead of presenting the mapping algorithm in Section 2.2, first the implemented dynamic mode decomposition (DMD) is explained and demonstrated in Section 2.1.

2.1. Dynamic Mode Decomposition

The DMD forces oscillatory time dynamics through a set of discrete linear differential equations whose solution consists of complex eigenvalues and eigenvectors. The eigenvalues represent the temporal oscillations, which may include a real part denoting growth or decay. The corresponding eigenvectors are the dynamic modes and represent spatial structures, which after entry-wise normalization to unity represent global one-component phase images (GOCPI).

Suppose a linear model $\mathit{A}$ can advance a (squeezed) frame ${\mathit{x}}_n$ at time $t_n$ to the next frame at time $t_{n+1}$: ${\mathit{x}}_{n+1}={\mathit{Ax}}_n$. Extending the model to advance a (time) series of N frames simultaneously, $\mathit{Y}=\mathit{AX}$, where matrices $\mathit{X}=[{\mathit{x}}_1,{\mathit{x}}_2,\dots ,{\mathit{x}}_{N-1}]$ and $\mathit{Y}=[{\mathit{x}}_2,{\mathit{x}}_3,\dots ,{\mathit{x}}_N]$ pair each frame ${\mathit{x}}_{n+1}$ with the previous frame ${\mathit{x}}_n$. It essentially requires $\mathit{A}$ to propagate a frame through time since $[{\mathit{x}}_1,{\mathit{x}}_2,{\mathit{x}}_3\dots ]=[{\mathit{x}}_1,{\mathit{Ax}}_1,{\mathit{A}}^2{\mathit{x}}_1\dots ]$, which in mathematical terms is referred to as a Krylov sequence. It indicates that $\mathit{Y}=\mathit{AX}$ is a discrete formulation that is closely tied to a system of continuous differential equations $d\mathit{x}/dt=\mathit{Ax}\left(t\right)$ and therewith an eigenvalue problem $\lambda \phi e^{\lambda t}=\mathit{A}\phi e^{\lambda t}$. Conceptually, $\mathit{A}$ hence describes a dynamical process, whose eigenvectors $\phi$ present the dynamic modes, with the associated frequencies captured in the complex eigenvalues $\lambda$. Note that, in contrast to the modes retrieved from an svd (Section 1, [53]), dynamic modes do not have to be orthogonal to each other, which implies that they do not have to be fully independent of each other.

The goal of the DMD is to find the eigenvalues and eigenvectors of $\mathit{A}$ at minimal computational cost. Finding the eigenvalues and eigenvectors straightforwardly by first calculating $\mathit{A}\approx {\mathit{YX}}^{†}$ ($†\equiv$ pseudo-inverse) poses a problem for computer memory. Say video footage has a frame size of $1000\times 1000$ px, then X and Y have a row size $m={10}^6$, resulting in an A matrix of size $m\times m={10}^6\times {10}^6$. Even in case a computer can handle such data loads, the calculation is slow, and the matrix size suggests a large amount of redundancy. The conceptual idea is to convert the eigenvalue problem from m into a lower dimension r. Typically, r is in the range 1–100 such that $r\ll m$, expressing a severe dimensionality reduction. By deflating the eigenvalue problem $m\to r$, eigenvectors and eigenvalues can be quickly calculated. Subsequently, the eigenvectors are inflated again $r\to m$ to yield the dynamic modes. Although the details of the conversions differ between DMD algorithms (e.g., standard DMD vs. optimized DMD [60,71], respectively), they share the common principle of using r modes from the svd of $\mathit{X}$ (Appendix A, Equation (A1)) for dimensionality reduction. There exists no general strategy to make a choice for r. It may be based on existing knowledge about the observed system [60], but it can also be based on the amount of singular values needed to capture a certain percentage of the variance in the video, for example, 99% [54], or on an algorithmic truncation of noise [72].

Note that $\mathit{X}$ (and $\mathit{Y}$) may first be converted to a time-analytic signal [53]. It is not necessary, but it has the benefit that the DMD, similar to an FFT, does not produce conjugate dynamic mode pairs. Then, r frequency components associate to r modes instead of $2r$ modes. Preventing the generation of conjugate modes also prevents the DMD to produce them imperfectly (i.e., with slightly different frequency to their counterpart). Analytic extension comes at a computational cost; however, since r can be halved to achieve the same number of frequency components, successive matrix multiplications within the DMD become computationally cheaper.

An important factor when applying a DMD is whether the recorded data contain standing-wave behaviour. Without adjustments, a DMD cannot capture standing waves [55]. For a fixed framerate, Tu et al. [55] offer a straightforward solution through augmentation of the video matrix $\mathit{X}$ with a time-shifted version of itself, by which the DMD acquires the skill to detect standing waves.

From the existing DMD algorithms, the optimized DMD based on variable projections [60] is an elaborate variant. Formulated as a least-squares optimization problem, it is more accurate than other DMD algorithms and theoretically allows video frames to be spaced non-equidistantly in time. Instead of splitting the video into two video matrices $\mathit{X}$ and $\mathit{Y}$, it uses a single, transposed video matrix ${\mathit{X}}^T={[{\mathit{x}}_1,...{\mathit{x}}_N]}^T$. A detailed explanation can be found in Appendix B. The methods skill for complex harmonic analysis was recently demonstrated in modelling rotating detonation waves [73]. Note that other DMD algorithms could be potent alternatives, for example, by allowing elaborate forms of regularization to handle noisy or lower-quality images [74].

Wave components in a wave field can be accurately extracted using the DMD [54]. To illustrate this, a wave field consisting of six known wave components is considered and recorded over a period of 32 s at 2 fps. The corresponding video matrix hence comprises 64 frames (${\mathit{X}}^T$ has 64 rows, see Appendix B, Equation (A2)). Subsequently, the signal is made time-analytic such that r modes associate to r frequency components. The DMD with $r=6$ modes identifies the underlying components precisely (Figure 3a, orange stars) (Figure 3b). It demonstrates that the DMD is not only powerful in analysing a short frame sequence but therein superior compared to a standard time-domain fast Fourier transform (FFT) (Figure 3a, green squares): the FFT returns (much) more spectral data and with significant redundancy, and yet, the pre-set frequency resolution restricts its ability to capture the six intrinsic wave components to a mere rough spectral representation. Half of the components (Figure 3a, first, third, and fifth component) are poorly captured in frequency and amplitude.

For a real wave field, the amount of wave components is not known a-priori and a choice needs to be made for the amount of modes r returned by the DMD. The question arises what the DMD returns if the choice for r is smaller or larger than the actual number of components in the wave field: if r is smaller, the DMD simply returns fewer components, but those components are still correctly represented (Figure 3c, $r=3$). If r is larger, the DMD still identifies the intrinsic components; however, it also adds spurious modes and these modes may be energetic, which indicates that r being too large is a scenario that should be avoided (Figure 3c, $r=12$). For observations of a real wave-field this scenario is unlikely, as waves have a stochastic nature and typically dense spectra [66]. In fact, for $r\sim O\left(1\right)-O\left(10\right)$, the number of modes is certainly less than the actual number of wave components (cf. Figure 3c, $r=3$), and the DMD’s representation of the wave field simplifies to the governing wave components.

A real wave field may also experience standing-wave dynamics as waves reflect back from the shoreline. Although, at flat dissipative beaches, sea-swell waves are typically progressive [75], at steeper beaches and longer sea-swell periods, incident wave reflection can be significant with reflection coefficients up to $K=0.4-0.45$ [76,77,78,79]. In the presence of hard structures like sea-walls even up to $K\approx 0.9$ [80]. Here, $K=[$0, 1] denotes [no, full] reflection [81]. High wave reflection signals (partially) standing-wave characteristics, prompting adaption of the DMD as of [55]. Doing so enables the DMD to cope with any degree of incident wave reflection (Figure 3c, $K=0$, $K=0.5$, and $K=1$) and thereby accommodates application to a broad range of wave field scenarios.

2.2. Mapping Algorithm

The proposed workflow for a self-adaptive and on-the-fly mapping algorithm of coastal hydrodynamics and bathymetry follows a series of steps (Figure 4, labels ①...⑬). The workflow requires video in greyscale and orthorectified format as basis input. If the video is in colour, standard BGR-to-grey conversions can be used to prepare the video. The orthorectification differs per field site and is briefly treated in Section 3. The workflow is set-up in such a way that a video feed could be processed. For that purpose, the algorithm marches forward in time by consecutively analysing short sequences of N video frames (e.g., $N=64$) and then updating the maps of $\mathit{c}$, d and $\mathit{U}$ after each sequence, which is referred to as mapping updates. Consecutive frame sequences may overlap, as for on-the-fly application, which is explained in Section 5 (see also [82], Figure 3). Default settings for the algorithm as used in this study are listed in

Table A2. Default parameter values of the mapping algorithm.

Parameter	Value
N	64 frames
Overlap	16 s (results, Section 4) Variable (on-the-fly, Section 5)
Analytic extension of$\mathit{X}$	True
Mode frequencies (= r if $\mathit{X}$ is analytic)	16
$[T_{min}$,$T_{max}]$	[3,15]s
Subdomain size per${\omega }^{\uparrow }$	$2\times L_{off}\left({\omega }^{\uparrow }\right)$ (maximum) $1\times L_{off}\left({\omega }^{\uparrow }\right)$ (minimum)
$[{\Gamma }_{min}$,${\Gamma }_{max}]$	[0.3,1.0]
$Rad$	75 m
Stationary time (temp. spectral storage)	60 s
$\alpha$	0.012
${\left|\mathit{U}\right|}_{max}$	0.75 m/s
$[d_{min}$,$d_{max}]$	[0.1,50] m
$Q_c$	0.0005 m${}^2/$s${}^3$
$Q_U$	0.0005 m${}^2/$s${}^3$

(Appendix C).

The workflow commences with the global analysis of a video frame sequence. The first step is to retrieve r global one-component phase images (GOCPI) through DMD (Figure 4, ①). GOCPI linked to frequencies outside the ocean wave band (Appendix C, Table A2, $[T_{min}$, $T_{max}]$) are discarded (Figure 4, ②). (Note that in case the video is not extended to a time-analytic signal, conjugate modes are also discarded; see Section 2.1). By construction of the DMD, the remaining number p of GOCPI with frequencies ${\omega }_1^{\uparrow }\dots {\omega }_p^{\uparrow }$ describe dominant wave field components. These are treated as equally important (regardless of their spectral weight b; see Appendix B).

Knowing ${\omega }^{\uparrow }$, the maximum wavelengths in each GOCPI are predictable. The size of the subdomains that are used to determine local ${\mathit{k}}^{\uparrow }$, can thereby be automatically tailored to the individual GOCPI (Figure 4, ③). This is done ahead of analysing any location. A conservative rule is to set the subdomain size to one or two offshore wave lengths, $L_{off}\left({\omega }^{\uparrow }\right)=2g\pi /{\left({\omega }^{\uparrow }\right)}^2$, where g denotes the gravitational acceleration. Here, $2\times L_{off}\left({\omega }^{\uparrow }\right)$ is used, unless this is too large, for example, near a boundary, where the size is reduced up to a minimum of $1\times L_{off}\left({\omega }^{\uparrow }\right)$. High-frequency GOCPI are analysed with smaller subdomains than low-frequency GOCPI, such that stacking the subdomains in layers on top of each other creates a pyramid-shaped cell at a certain location (see Figure 4, ③, yellow pyramid). Note that this contrasts with the usual cube shape, whose size is typically predefined manually, e.g., [41,83].

Since the pixel resolution is constant between cell layers, large layers used for lower frequencies likely encompass plentiful pixels, capturing the underlying waves in unnecessarily high spatial resolution. It is therefore computationally attractive to subsample larger cell layers. The rules for subsampling require a minimum resolution of 8 samples per $L_{off}$. Instead of skipping pixels for subsampling, the current set-up lowers the resolution through fast bilinear interpolation. It preserves more image information and makes the algorithm robust to videos with different pixel resolutions. For proper analysis, a cell layer is set to hold at least $24\times 24$ samples. Note that demanding eight samples per offshore wave length may exclude higher-frequency GOCPI from being analysed, if the pixel resolution of the video is not high enough to capture the correspondingly short wavelengths.

After subdomain sizes and resolutions have been determined, the local analysis of the GOCPI commences. The local analysis occurs around points in a grid. The processing of different grid point locations can be done in parallel, which increases computational speed. Since the algorithm aims to map results on-the-fly, the grid resolution is determined based on computational speed (see also Section 5). Note that computation speeds are different for different processing machines. A solution to find an optimal number of gridpoints for different processing machines could be to start with a low number of gridpoints and then increase the number of gridpoints for consecutive frame sequences until an optimum is reached.

To analyse a certain grid point, first a pyramid-shaped cell is built around it (Figure 4, ④) using GOCPI subdomains as cell layers with formats as determined under step ③. Each cell layer is first autocorrelated in space to accentuate the waveform and tapered with a two-dimensional Hanning-window, which focusses wave information to the centre point and prepares for analysis with two-dimensional fast Fourier transforms (2D-FFTs). For robustness, ${\mathit{k}}^{\uparrow }$ (Figure 4, ⑤) is estimated in two different ways: directly, by analysing spatial properties through 2D-FFT and indirectly, by analysing wave celerity through particle image velocimetry (PIV).

Performing a 2D-FFT on a cell layer yields its spectral content in space. Besides the footprint of a gravity wave component, ${\mathit{k}}^{\uparrow }$, this also includes spectral noise. Typically, an energy threshold aims to separate the two [41]. To avoid a search for an optimal energy threshold [2], here simply the spectral point with maximum energy is chosen as ${\mathit{k}}^{\uparrow }$. A different approach to acquire ${\mathit{k}}^{\uparrow }$ is via the wave celerity, ${\mathit{k}}^{\uparrow }={\omega }^{\uparrow }/{\mathit{c}}^{\uparrow }$. The wave celerity is encoded in the complex values of the cell layer. The difference between the real and imaginary parts denotes a temporal phase shift of ${90}^{\circ }$, hence resembling a quarter wave period. Performing PIV between the real and imaginary image of the cell layer yields the translation of (wave) patterns over a quarter wave period, which directly translates to ${\mathit{c}}^{\uparrow }$. Note that the temporal phase shift between real and imaginary parts is ${90}^{\circ }$ regardless of whether the video matrix $\mathit{X}$ is analytic or whether dealing with standing waves. While the 2D-FFT approach inherently presumes that the observed pattern in a cell layer describes a wave, the PIV approach does not, as it investigates movement of any pattern in the cell layer. Both approaches may have their benefits in case the pattern is noisy and are therefore used synchronically. If one or both estimates are unphysical, they are discarded in following the filter steps. Each estimate of ${\mathit{k}}^{\uparrow }$ gets a weight W assigned, whose value between $[0$, $1]$ gives an indication of the respectively [low, high] quality of the estimate (Figure 4, ⑤ red-coloured bar). A weight, W, is calculated from the correspondence between two images: the source image and the target image. For the 2D-FFT approach, this is the correspondence between the original pattern displayed by the cell-layer and the approximation of this pattern by its most energetic spectral wave component. For PIV, it is the correspondence between the translated real image of the cell-layer with the imaginary image of the cell-layer. The correspondence is calculated via $W=1-{\epsilon }_{im}$, where ${\epsilon }_{im}$ represents the normalized root-mean-square error between the source image relative to the target image [84].

Together, the estimates of ${\mathit{k}}^{\uparrow }$ over all ${\omega }^{\uparrow }$ layers in the pyramid cell form a sparse spectral point cloud (SSPC) of ${\mathit{k}}^{\uparrow },{\omega }^{\uparrow }$ pairs with assigned weights W. The SSPC expects to follow a cone shape as described by the dispersion relationship (Figure 1; Equation (1)). This pre-knowledge aids in the identification and removal of unphysical ${\mathit{k}}^{\uparrow },{\omega }^{\uparrow }$ points in the SSPC by means of a wavelength and direction filter (Figure 4, ⑥). Wave lengths are filtered using their ratio over the offshore wave length $\Gamma =|{\mathit{k}}_{off}|/|{\mathit{k}}^{\uparrow }|$, where $|{\mathit{k}}_{off}|=2\pi /L_{off}\left({\omega }^{\uparrow }\right)$. SSPC points with $\Gamma <0.3$ occupy the shallow water regime as they resemble wavelengths larger than twenty times the local depth [85] and are thereby too large to capture morphological detail. Such ${\mathit{k}}^{\uparrow }$ are deemed unsuited for a localized depth estimate, d. On the other hand, SSPC points with $\Gamma >1.0$ signalize an unphysical deep water regime if $\mathit{U}$ is small and are therefore also deemed unsuited. Note that in other algorithms this upper limit is often set lower, to $\Gamma >0.9$, which presents an approximate elbow value where the uncertainty in estimates of d becomes disproportionately large [39,53,86]; however, a higher limit $\Gamma >1.0$ preserves spectral points that are valuable for the estimation of $\mathit{U}$ and thereby also simultaneous estimates of d. Lastly, the wave direction filter excludes SSPC points, which do not align with the general direction by means of the svd filter of Gawehn et al. [2].

The idea behind the next two steps of the workflow (Figure 4, ⑦, ⑧) is to augment the SSPC left after ⑥ with additional spectral points to make it a dense spectral point cloud (DSPC), which captures directional spread in time and space and thereby allows for a solid inversion of d, $\mathit{U}$. As desired and designed, the SSPC as is holds the essence of the local hydrodynamics. The sparsity, however, brings along statistical uncertainty in time and space. An approach to acquire a DSPC is to combine several SSPCs into one. The additional SSPCs can be retrieved from preceding updates (i.e., previously analysed frame sequences) but also from surrounding grid locations within a collection radius $Rad$ (Figure 4, ⑧ attached green arrow; Appendix C, Table A2). Such an approach requires memorizing SSPCs in a designated short-term storage. Therein, the just-retrieved SSPC is stored (Figure 4, ⑦) and SSPCs from preceding updates and from surrounding grid locations are called up (Figure 4, ⑧). The size of the short-term storage for SSPCs is manageable since, by construction, each SSPC holds just a few essential spectral points. Note that accumulating spectral information over successive updates is only a valid approach if their computations occur (near) on-the-fly, because it assumes that the wave signal is stationary across several updates. Hence, SSPC data are discarded after a short stationary time period (e.g., 60 s; see Appendix C, Table A2) and replaced by new SSPCs.

An example illustrates the augmentation of an SSPC to a DSPC by means of stored spectral data: assuming a wave signal is stationary for 1 min, the algorithm stores SSPCs for 1 min. If the processing of a frame sequence takes 15 s, the storage contains stored SSPCs of all grid locations from the preceding 60 s/15 s = 4 updates. Now, as the algorithm starts to analyse a certain location, it retrieves one new SSPC for that location. By itself, this new SSPC might be sufficient to estimate d but not $\mathit{U}$. For that, a DSPC is required, which is achieved by augmenting the new SSPC with stored SSPCs from the same grid location and surrounding grid locations (Figure 4, ⑨), be it 10 within some radius $Rad$. The acquired DSPC then consists of $4_{stored}\times (1_{gridloc}+{10}_{surr.gridlocs})+1_{new}\times 1_{gridloc}=45$ SSPCs. In this example, the augmentation of the new SSPC of step ⑥ has hence produced a DSCP that is 45 times denser.

Using surrounding grid locations to augment spectral information has benefits and drawbacks. A drawback is that depth estimates d become less localized. A benefit is that wave refraction caused by near surface currents $\mathit{U}$ is captured in space, which is essential for estimating $\mathit{U}$ (and thereby also improves d). To minimize the loss in localization of d, but to still have improved $\mathit{U}$ estimates, the SSPCs that make up the DSPC are weighted differently. These SSPC weights are multiplied with the individual quality weights W of the spectral points and are subsequently normalized to the range $[0$, $1]$. In the current set-up, the new SSPC at a certain grid location weighs 50%. Stored SSPCs of the same grid location weigh together 25%, and stored SSPCs from surrounding grid locations weigh together 25%. In total, this means that 75% of the spectral information focuses on the currently analysed grid location, and 25% focuses on the surrounding area, thereby aiming to keep d estimates localized but still constructing a DSPC that includes enough current refraction for $\mathit{U}$ estimates.

The DSPC from step ⑨ may include a minority of spectral points from incorrect wave directions. This occurs because some of the SSPCs, which the DSPC consists of, contain too little spectral information to determine wave directions at step ⑥. Therefore, the direction filter repeats for the DSPC (Figure 4, ⑩). Now, the DSPC is ready for the algorithm to retrieve wave celerities per frequency, ${\mathit{c}}^{\uparrow }\left({\omega }^{\uparrow }\right)$ (Figure 4, ⑪) and depths and surface currents d, $\mathit{U}$ (Figure 4, ⑫). While ${\mathit{c}}^{\uparrow }$ are directly computed via ${\mathit{c}}^{\uparrow }={\omega }^{\uparrow }/{\mathit{k}}^{\uparrow }$, the inversion of d, $\mathit{U}$ is done by fitting a wave model—here, the Doppler-shifted dispersion relationship (Equation (1))—to the spectral points of the DSPC. The fit results from a nonlinear regression that aims to minimize the sum of residuals between the model frequencies ${\omega }_{j,model}$ (Equation (1)) and the observed frequencies ${\omega }_j^{\uparrow }$ per spectral point j using Equation (2). The implicit link of observed ${\mathit{k}}_j^{\uparrow }$ with ${\omega }_{j,model}$ and ${\omega }_j^{\uparrow }$ is assumed to be trivial.

$$\begin{array}{c}\mathrm{minimize}F(d,\mathit{U})=\frac{1}{2}\sum _j{W}_j\rho \left(f_j{(d,\mathit{U})}^2\right)\\ \mathrm{subject} \mathrm{to}\\ d_{min}<d<d_{max}\\ {\left|\mathit{U}\right|<\left|\mathit{U}\right|}_{max}\\ \mathrm{with}\\ \rho \left(f_j{(d,\mathit{U})}^2\right)={\alpha }^{2}ln\left(1+\frac{f_j{(d,\mathit{U})}^2}{{\alpha }^2}\right)\\ \mathrm{and}\\ f_j(d,\mathit{U})={\omega }_{j,model}(d,\mathit{U})-{\omega }_j^{\uparrow },\end{array}$$

(2)

where $F(d,\mathit{U})$ is the cost function, minimized by adjusting the regression parameters d and $\mathit{U}$. Typically, residuals per spectral point, j, are evaluated by the product of their weight $W_j$ with the difference $f_j(..)$ between ${\omega }_{j,model}$ and ${\omega }_j^{\uparrow }$. Here, $f_j(..)$ is first modulated by a Cauchy loss function $\rho (..)$ to penalize outliers. A predefined softmargin $\alpha$ tunes $\rho (..)$ to the optimization problem. Setting $\alpha$ such that the Cauchy loss function heavily penalizes residuals larger than $\sim 1/5$th of the average DMD frequency resolution produces accurate d, while counteracting overfitting of $\mathit{U}$. Lastly, bounds for depth $[d_{min}$, $d_{max}]$ and surface current magnitude ${\left|\mathit{U}\right|}_{max}$ set the range in which a solution for d,$\mathit{U}$ is sought. To improve estimates, first a regression is done without a Doppler shift to acquire a close estimate for d that can be used as an initializer for the follow-up regression including a Doppler shift. Adding the first regression has little impact on computation times, since its estimate of d closely approximates the local minimum of the second regression. The minimization of Equation (2) occurs using a sequential least-squares quadratic programming method (SLSQP) [87], which omits potentially expensive computations of Hessians and allows for straightforward implementation of loss-functions.

The final step of the mapping algorithm is a Kalman filter (Figure 4, ⑬) that judges the quality of d,$\mathit{U}$ results based on the sensitivity of the fit with regard to earlier updates. It causes d,$\mathit{U}$ estimates to quickly converge over successive updates. The implementation of the Kalman filter is identical to Gawehn et al. [2]. If the observation periods are short, say, less than 5 min, morphodynamics but also hydrodynamics can be assumed stationary, such that process variance is negligible. If observation periods are longer, say, 10–20 min, it may become important to capture changes in surface current direction (e.g., due to the formation of a rip current). To allow for such applications, the Kalman filter assumes small process variances $Q_c$ and $Q_U$ for phase celerities and near-surface currents, respectively (i.e., $Q_c=0.0005$ m${}^2$/s${}^3$, $Q_U=0.0005$ m${}^2$/s${}^3$; see Appendix C, Table A2).

3. Field Sites and Data

Videos of four different field sites around the world were used to test the algorithm’s performance (Figure 5): Duck (North Carolina, USA) [43], Porthtowan (UK) [88], Scheveningen (NL) [89], and Narrabeen (AU) [90]. Specifics on video collections are listed in

Table 1. Video collection.

Field Site	Instrument	Camera Height (m)	Camera Tilt (${}^{\circ }$)	Video Length (min)	Frame Rate (fps)	Pixel Size (m)	Reprojection Error (at Distance)
Duck	Argus	43	68–82	17	2	5	<1 m (<500 m)
Porthtowan	Argus	44	75–85	17	2	5	-
Scheveningen	UAV	110	61	9	2	2	∼ 1 m (<200 m)
							∼7 m (400–600 m)
Narrabeen	UAV	89	73	9	2	2	< 1 m (<250 m)
							∼7 m (1.5–2 km)

. Cameras were positioned at heights of 43–110 m to observe waves over a large distance without wave-shadowing effects. The recorded videos are available in orthorectified format. The geometries have been solved using ground control points (GCPs), by matching GCP image coordinates with world coordinates. The accuracy of the geometries differs per site and was quantified by the GCP reprojection error in world coordinates. These reprojection errors are generally in the (sub)meter range, but in case of Scheveningen and Narrabeen, they increased to on average $\sim 7$ m at distances of 0.5–2 km from the camera. At Porthtowan, the reprojection errors are unknown, yet slight errors in the geometry likely exist, particularly further afield [91]. The videos from Duck and Porthtowan (UK) are $\sim 17$ min long and were recorded with Argus stations. The videos of Scheveningen and Narrabeen are ∼ 9 min long and were recorded with UAVs. All videos have a framerate of $\Delta t=2$ fps. Pixel resolutions are $\Delta px=5$ m for Duck and Porthtowan, and $\Delta px=2$ m at Scheveningen and Narrabeen. The hydrodynamic conditions varied over the sites between $H_s=0.8-1.63$ m and $T_p=5.0-10.0$ s (

Table 2. Hydrodynamic conditions during field recordings.

Field Site	${\mathit{H}}_{\mathit{s}}$ (m)	${\mathit{T}}_{\mathit{p}}$ (s)	$\mathit{WL}$ (m)
Duck	0.79	5.0	0.08
Porthtowan	1.03	10.0	−0.96
Scheveningen	0.75	5.5	0.60
Narrabeen	1.63	8.5	0.67

).

The videos of all sites were analysed using sequences of 64 frames equalling 32 s. Marching forward in time, the next frame sequence overlapped 50% with the previous sequence (i.e., overlap of 32 frames at 2 fps), resulting in mapping updates every 16 s of video. This overlap was fixed for reproducibility, since computation times differ per processing machine. For on-the-fly applications, the overlap varies based on processing speed. This is simulated in Section 5 using a standard laptop. Each image sequence was decomposed into $r=16$ dynamic modes. The chosen r represents a balance between significantly reducing dimensionality, while retrieving enough pyramid cell layers to estimate d, $\mathit{U}$. Details on other settings are found in Appendix C, Table A2.

In Section 4, the results of the mapping algorithm are presented. The field-site of Duck, NC, USA was used as the lead case to elucidate the processing steps of the workflow for an arbitrary real case. The final results are presented for all field-sites: Duck (NC, USA), Narrabeen (Australia), Scheveningen (the Netherlands), and Porthtowan (UK). The quality of the final results for the depth estimates, d, was assessed by comparison with in-situ bathymetry data. Maps for wave direction and celerity and near-surface currents are discussed in Section 5.

4. Results

The field site of Duck gives an illustrative example (Figure 6) of the processing steps described in Section 3. After decomposing image sequences of Duck (Figure 6a) into 16 dynamic modes, the modes were normalized to global one-component phase images (GOCPI) and filtered for frequency and resolution (as of Figure 4). The frequency filter discarded eight GOCPI, and another four GOCPI failed the criterion of minimum 8 px$/L_{off}$, leaving four GOCPI for further analysis (Figure 6b). These remaining GOCPI reveal intricate wave patterns and finely capture wave refraction towards the coast. The associated frequencies ${\omega }^{\uparrow }$ (rad/s) $\approx \{0.59,0.78,0.97,1.17\}$ were quite constant across successive image sequences, showing that wave periods T(s) $=\{10.6,8.1,6.5,5.4\}$ govern the mixed wave field. On the basis of these wave components, the subdomain sizes for local analysis were determined. Stacking them in layers for some grid location revealed a pyramid-shaped cell, and tapering the layers put the focus toward its centre grid point (Figure 6c).

Performing PIV and 2D-FFTs to the autocorrelated layers revealed a sparse spectral point cloud (SSPC) with eight points, with ${\mathit{k}}^{\uparrow }$ values increasing for increasing ${\omega }^{\uparrow }$, sketching a typical wave-dispersion curve (Figure 6d). In the example, the PIV and 2D-FFT estimates for ${\mathit{k}}^{\uparrow }$ were almost indistinguishably close in three of the four spectral layers; however, for the uppermost layer, they lay further apart. In this case, the point that lays off-track corresponds to the 2D-FFT estimate. This is signalled by its quality weight W, which is lower than for the PIV counterpart (Figure 6d, blue vs. yellow point). Augmenting the SSPC with stored SSPCs from the past minute (i.e., from four previous updates) within a radius of $Rad=75$ m yields the dense spectral point cloud (DSPC) (Figure 6e). The choice of $Rad=75$ m was arbitrary and represents a balance between extracting information from close by, while capturing sufficient current refraction (especially of shorter period waves). Other choices for $Rad$ may be made, but the algorithm is not too sensitive to this parameter since the surrounding SSPCs within $Rad$ resemble just 25% of the total spectral weight. As desired, this lower weighting is apparent from the corresponding spectral points in the DSPC (Figure 6e, blue points). Note that in contrast to Duck, points in DSPCs of Porthtowan, Scheveningen, and Narrabeen were more dispersed over the frequency domain, because GOCPI frequencies vary more across successive updates (not shown). Finally, the DSPC was fitted with the Doppler-shifted linear dispersion relationship (Equation (1)) to produce characteristic cones corresponding to certain d,$\mathit{U}$ estimates (Figure 6f, magenta cone). Combining the local estimates from all grid locations yields global maps of d,$\mathit{U}$.

These resulting maps were not only quickly retrieved, but they also show that estimates of d are accurate (Figure 7). For all four field sites, depth maps compare with the ground truth. Recall that the algorithm was not tuned to the individual field sites. While the first depth update–after 32 s of video–is still rough, it already gives a clear overall picture of the bathymetry with shallower and deeper regions (Figure 7, 1st update). Nearshore sandbars at Duck (Figure 7a) are readily visible. With the fifth update–after 96 s of video–(Figure 7, 5th update) the depth maps approximated the ground truth (cf. Figure 7, 5th update and ground truth). Estimates quickly improved after the first update, indicating that the temporary spectral storage together with Kalman filtering effectively converged d estimates. Mapping updates continued to improve and become spatially more coherent towards the end of the videos (Figure 7, last update).

For all videos, differences with the ground truth are minimal over large parts of the observed area (Figure 7, difference, light areas). Regions with errors $|\Delta d|>0.5$ m are mostly found in shallow water or the deep water boundary of the observed domain. At Duck and Scheveningen, depth overestimation of $\Delta d\approx -0.5$ to $-1.0$ m occurs around the sandbars (Figure 7a,c, difference, blue areas), while in other shallow parts, the errors are smaller. At Porthtowan and Narrabeen, the overestimations are larger $\Delta d\approx -0.5$ to $-2.0$ m and generally occur in shallow water (Figure 7b,d, difference, blue areas). The reason is probably that wave heights at Porthtowan and Narrabeen are larger $H_s=1.0-1.6$ m (Table 2) and more nonlinear close to shore, compared to Duck and Scheveningen where $H_s<0.8$ m and wave breaking restricts to the sandbars, after which most wave nonlinearity is lost. Another, physical reason for near-shore depth estimates appearing larger is local wave set-up [26], which is not accounted for in the comparison with the ground truth. Regions and reasons for depth underestimation differ per site. At Duck, depths were only underestimated around the pier, which blocks the view to the underlying waves (Figure 7a, red patch). At Porthtowan and Scheveningen, depths were underestimated by $\Delta d\approx 0.5-1.5$ m near the offshore boundary (Figure 7b,c, red patches). At Porthtowan (Figure 7b), the underestimation was mainly caused by the fact that relevant lower-frequency cell layers in near-boundary pyramid cells were too large to be used. The underestimated region also lay >800 m from the camera, where inaccuracies in geometry influence depth estimates [91]. The size of errors in the video geometries (Table 1) suggests a limited effect on depth estimates in general, yet some depth error may be induced in regions further afield (e.g., for $T=8$ s and $d=10$ m, an error $\Delta L=7$ m in wave length causes a depth error of $\Delta d\approx 2$ m; see also Figure 1 of [39]). At Scheveningen (Figure 7b), the underestimated region begins further away from the offshore boundary, where boundary effects should be less pronounced. Here, the underestimation likely stems from the relatively small waves, $H_s=0.75$ m and $T_p=5.5$ s, who feel little of the underlying bottom. At Narrabeen, underestimation mainly occurs at the boundary farthest from the camera. This underestimation is likely caused by similar boundary effects as at Porthtowan, since here the underestimated region also lies more than 900 m from the camera. In conclusion, depth maps can show regions that are less accurate, yet all in all, the maps approximate local bathymetries: on average 80% of the mapped area had errors $\Delta d<1$ m.

Direct comparison of estimated depths against the ground truth confirms generally accurate depth maps (Figure 8). Median depth estimates per given depth are mostly close to ground truth (cf. Figure 8, double green curves and black 1:1 line). For the video of Duck (Figure 8a), the similarity is visible for the entire depth range from $d=1.5-5.0$ m, except for the scour hole where the pier obscures the underlying waves. At Porthtowan, median depth estimates deviate more from the ground truth but remain $\Delta d<2.0$ m for depths $d<10$ m (Figure 8b). Errors are largest in the breaking region where depths are $d=2-3$ m, which is a common observation in depth inversion studies e.g., [92]. This is similar for Narrabeen (Figure 8d). Here, however, differences between median depth estimates and the ground truth are minimal for the entire depth range beyond that, $d=3-17$ m. In contrast to Porthtowan ($H_s=1.0$ m), the wave heights at Narrabeen were larger ($H_s=1.6$ m); moreover, the few large boundary errors in the statistics were outnumbered by otherwise accurate estimates (Figure 8d, relatively few outliers). Median depth estimates at Scheveningen (Figure 8c) also confirm earlier spatial observations of slight depth overestimation near the shallow sandbar and underestimation for $d>7$ m.

Although for each field site the first mapping update after 32 s of video was considerably more scattered compared to the last update (cf. Figure 8, white and blue scattered dots), the median estimate is comparable (cf. Figure 8, double grey and double green curves). Hence, the first mapping updates already give a rudimentary impression of the bathymetry, albeit with more local uncertainty. It suggests that waiting for many mapping updates is superfluous and that it may be efficient to stop after a couple of updates, when a certain degree of accuracy in the depth map is reached.

Bulk errors decreased with the increasing number of mapping updates (Figure 9). Since depth errors are not always constant and linearly distributed over depth (Figure 8b) and may contain outliers (Figure 8d), the median bias was adopted to quantify structural over- or underestimation, and the interquartile range (IQR = 75th −25th percentile) was used to measure the scatter. In comparison to other common error measures, note that Duck, with quite linear and normally distributed depth errors (see Figure 8a), had a median bias that was almost identical to the commonly used mean bias, and its IQR was close to the root mean square error (not shown). Absolute median biases (Figure 9, dashed lines) start small $|\Delta d_{bias}|<0.5$ m with the first update and eventually become $|\Delta d_{bias}|<0.1$ m at Duck (with update 3), Scheveningen (with update 27), and Narrabeen (with update 10). The median bias at Porthtowan reached $\Delta d_{bias}\approx -0.25$ m. The IQRs (Figure 9, solid lines) decreased fast over the first couple of updates, signalling fast improvements in the accuracy of the depth map. After that, the convergence rates started to relax. The total IQR improvements were 0.9 m → 0.5 m (Duck), 1.7 m → 1.4 m (Porthtowan), 1.6 m → 0.6 m (Scheveningen), and 2.6 m → 1.1 m (Narrabeen). The elbow in the exponentially decreasing curves represents a compromise between the size of error and the number of updates. It is difficult to pinpoint the exact locations of the elbows, but for Duck and Porthtowan, they appeared to be somewhere between the 2nd–6th update, while for Scheveningen and Narrabeen they were more likely located between the 5th–10th update. The reason for the different elbow positions probably roots in the instrumentation. The videos of Duck and Porthtowan were recorded with stationary Argus stations, while the videos of Scheveningen and Narrabeen were recorded with UAVs. In contrast to the Argus stations, UAVs have some freedom to move and are not professionally tuned to the field site, which likely causes depth errors to be initially large; however, these errors decrease rapidly. Concluding the previous findings, a general rule of thumb might be to stop with the 5th update, after 1.5 min of video. Awaiting more updates may further improve results but also requires more computation time and a larger video length. A single estimate based on 32 s of video is enough to get a rough depth map. In the end, the use of the data, and desired accuracy, determine whether to stop the analysis sooner or later.

5. Discussion

Alongside the maps of depth, mapping updates for hydrodynamics, wave celerity, $\mathit{c}$, and local near-surface currents, $\mathit{U}$, were retrieved. Since no ground truth data are available for these parameters, they are discussed qualitatively in Section 5.1. The subsequent Section 5.2 discusses the algorithm’s ability to process video on-the-fly.

5.1. Maps of $\mathit{c}$ and $\mathit{U}$

While vector fields of $\mathit{c}$ remain quite stable after the first mapping update, vector fields of $\mathit{U}$ require more updates to converge. Here, approximately the 10th update (Figure 10) shows convergence of the local current patterns. Additionally, the current velocities and directions appear realistic. Note that local changes in $\mathit{c}$ and $\mathit{U}$ still occur after the 10th update, as allowed by the process variance in the Kalman filter (not shown).

Maps of $\mathit{c}$ reflect the general direction of wave propagation for all videos (Figure 10, left column). Local vectors are oriented normally to the recorded wave crests and trace characteristic wave refraction towards the coastline. Mean wave directions per wave period are quite similar for all field sites, except for Duck, which shows some directional spread between shorter and longer period waves (Figure 10a left, red vs. yellow arrows). Wave celerity decreases towards shallower water, being most apparent for Porthtowan and Narrabeen. This was expected as cross-shore differences in depth were the largest at these locations (Figure 10b,d, left). Moreover, the wave celerity of long-period waves started to decrease further offshore than for short-period waves (Figure 10d left, starting at offshore boundary, green arrows already get shorter while red arrows had constant length until close to the shoreline). This was also expected since longer waves feel relatively smaller depths.

Maps of $\mathit{U}$ show intricate patterns that are largely consistent in time and space (Figure 10 right). Coastal currents are typically tide or wave driven [85] and have magnitudes of decimetres per second (dm/s), e.g., [93]. They can reach meters per second (e.g., strong rip currents [94] such as the “Backpackers’ Express” at Bondi Beach [95]), but such conditions are not likely at the field sites analysed in this study. With mostly 0–5 dm/s, the maps of $\mathit{U}$ had the correct order of magnitude. At Duck, near-surface currents close to shore often pointed in offshore directions, while further from shore they also pointed in southerly directions (Figure 10a right). The near-shore offshore-directed flows probably represent the effect of undertow caused by normally incident waves (see Figure 10a left). It is interesting to note that surface current estimates using optical flow on wave-averaged images might predict opposite flow directions in shallow water, namely, shoreward directed currents, as recent data analyses for a similar wave situation at Duck suggest [96]. The hypothesis is that current estimates from optical flow are more indicative of Stokes drift, while wave-inverted current estimates capture the undertow. Especially around the thinner part of the sandbar in the south (see also Figure 7a, ground truth), streamlines from northerly and southerly directions often converge to form a common flow direction offshore, which could be indicative of a local rip current and which is not uncommon at this site [97]. Although it is visible throughout many updates (not shown), it is often hard to recognize among enlarged $\mathit{U}$ estimates on the sandbars. Note that $\mathit{U}$ estimates on the sandbars are influenced by distorted wave celerity estimates of breakers, which simultaneously leads to larger depth errors (see Figure 7a, difference). At Porthtowan, near-surface current estimates generally point in south and south-easterly directions, except for the deeper region where directions are less coherent. Additionally, at shallow parts, currents point away from the coast (Figure 10b right). The southerly directed flows may reflect some remaining tidal ebb flow consistent with the time of the video recording approaching low water. Similar to Duck, offshore-directed flows close to the shoreline are suggestive of a cross-shore-directed undertow under almost shore-normal wave incidence. Near-surface current estimates at Scheveningen were well in line with expectation. The video recording was taken just above a harbour (Figure 10c right). During flooding, the longshore tidal current accelerates around the harbour jetties forcing larger current velocities, and a characteristic region with eddies and opposite directed flows is formed on the lee side of the northern jetty. Both effects were visible, and current magnitudes of mostly 0.4–0.6 m/s off the coast also agree with typical current magnitudes off the nearby Rotterdam coast. Close to the coastline, $\mathit{U}$ estimates are again offshore-oriented, suggesting undertow. Similar to Scheveningen, coastal currents at Narrabeen run mostly alongshore. Here, however, they are likely not driven by tides but by waves instead. The direction of flow at this location is sensitive to the angle of wave incidence and alongshore differences in wave height and dissipation, e.g., [98]. For the recorded situation, the angle of wave incidence suggests northward flow; however, alongshore differences in wave height may still force a southward flow as mapped. The true current direction remains uncertain at this site. Summarizing, near-surface current estimates are coherent and can be explained, but the Scheveningen field-site is most relatable to expected tidal flow patterns (Figure 10c, right).

5.2. On-The-Fly Processing

Maps of d, $\mathit{c}$, and $\mathit{U}$ are realistic and ideally returned on-the-fly. To allow for an on-the-fly analysis during this study, specific choices were made to balance calculation time and spatiotemporal resolution of mapping updates. The algorithm marches forward in time by consecutively analysing small sequences of video frames, giving mapping updates after each sequence. The frame sequences partly overlap (see legend Figure 11a). For pre-recorded video, the overlap could be controlled. It was 50%, since frame sequences were 32 s long, and the algorithm was set to march forward in 16 s intervals. For on-the-fly analysis of a video feed, this approach does not work, as processing times start to play a role. A mapping update now needs to reflect the current situation and therefore is mostly based on the latest recorded video frames. Say the processing time of a frame sequence was 20 s, then during the processing, 20 s of new video frames have been recorded. The next frame sequence to analyse includes these new 20 s plus 12 s from the previous frame sequence representing the overlap. The amount of overlap hence depends on the processing time, which may vary per update. Statistically, the shorter a frame sequence, the more it represents a random sample of the wave field, where waves vary in height and direction, currents change, etc. [66]. Additionally, the video quality itself may vary due to changes in lighting, the accuracy of the orthorectification, etc. This impacts the processing time, which may hence be slower or faster causing smaller or larger overlap, respectively, between consecutive frame sequences.

For the results presented up until now, the processing time for an update—using a standard laptop with four CPUs and a working memory of 4 GB—was typically 30–60 s. If the fixed, 16 s overlap was changed to instead be variable-based on the processing time, this would be too slow. To enter the realm of on-the-fly computation, the processing time hence needs to be further reduced. The three main elements influencing processing time are: (i) the computational power, (ii) the grid resolution, and (iii) the amount of spectral data points stored and used to locally derive $\mathit{c}$, d, $\mathit{U}$. In this study, the effect of (ii) and (iii) was considered using the video of Duck as a testcase:

(ii) In a first step to increase computational speed, the grid resolution was reduced from 1260 to 720 grid points. As required to simulate on-the-fly analysis, the overlap was set to depend on the processing time. The resulting simulation returns mapping updates in variable time intervals with correspondingly variable overlap between analysed frame sequences (Figure 11a). Yet, on average, every 23 s, a mapping update is given. Although the grid resolution is lower, the bathymetry is still estimated in reasonable detail, as well as wave propagation and near-surface currents. At the 5th update, after 102 seconds, the results again converged to a large extent. The analysis could at this point be deemed on-the-fly; however, more frequent updates would further improve the user experience. A detailed view of the first $\approx 4$ updates reveals that they were returned at a faster pace than consecutive updates (Figure 11a, orange lines closer together, and more, darker green overlap between frame sequences). During the first updates, the amount of spectral data in storage still increased towards full capacity, which suggests that the rate of updates should increase for a smaller number of spectral data points.

(iii) Reducing the number of spectral data points that are stored and used reduces the number of spectral points in the dense spectral point clouds (DSPCs). This should be acceptable as long as the true local wave spectra are well represented by the DSPCs. By construction, a DSPC represents an aggregate of (mostly stored) sparse spectral points clouds (SSPCs) (Section 2). These SSPCs stem from surrounding grid points within a radius $Rad$ (Appendix C, Table A2, $Rad$ = 75 m) from the location under analysis and from previous updates within a short period where the wave signal is assumed to be stationary (see Section 2 and Figure 4). Now instead of using SSPCs from all grid points within $Rad$, a random subset can be selected. Here, a random subset of 12 surrounding grid points was selected from the full set of 48 grid points. The amount of SSPCs from this subset of grid points can be reduced even further by simply storing fewer SSPCs. The storage time, which is the time that the wave signal is assumed stationary (Section 2), was here reduced from 60 s to 30 s. The effect of both measures on processing speed was significant, while mapping results remained similar (not shown). The 5th update now occurred approx. 20 s earlier and the 10th update even 80 s earlier compared to the simulation with only a reduced grid resolution (cf. Figure 11a,b). Updates were, on average, returned every 13 s, which may be considered a quite continuous temporal output, as desired for on-the-fly application (Figure 11b). It is noteworthy to mention that convergence rates appear to depend rather on the number of updates than on the timing (not shown), which suggests that it may be favourable to return updates at a faster pace; however, a detailed analysis remains a subject for future study.

6. Conclusions

This study describes a fast and self-adaptive algorithm to map coastal parameters on-the-fly from aerial wave imagery. Updates of depth, d, wave propagation, $\mathit{c}$, and near-surface currents, $\mathit{U}$, are returned every few seconds, such that a video feed can theoretically be processed on-the-fly, and a user does not need to wait and engage in post-processing. The input requires orthorectified video with known pixel size (m) and frame rate (fps). Apart from that, the algorithm works unsupervised for the presented field sites. The basis for fast computational speed and increased automation lays in the use of the dynamic mode decomposition (DMD), which is a dimensionality reduction technique that disposes redundant video information. It reduces the video to a set of, here, 16, intrinsic wave patterns and autonomously finds the corresponding frequencies. Applying a DMD, the search for optimal spatial sampling schemes can also be automated, such that in the end, no manual choices underlie the local wave spectra, which are the basis to derive $\mathit{c}$ and invert d, $\mathit{U}$. Consecutive mapping updates are improved by taking specific measures, such as using an innovative system to temporarily store and call up spectral data but also by penalizing spectral outliers through a loss-function and Kalman filtering. The algorithm’s potential for all-round application was demonstrated using video data from stationary camera installations and UAVs, recorded at four different field sites, Duck (USA), Porthtowan (UK), Scheveningen (NL), and Narrabeen (AU), with hydrodynamic conditions ranging between $H_s=0.75-1.63$ m and $T_p=5-10$ s. Validating the depth maps, d, showed that they accurately reflected ground truth measurements. The maps quickly improved from the first update to typically the fifth update, at 1.5 min into the video, after which the rate of improvement relaxed. It suggests that for time-efficient coastal mapping of depths, 1.5 min of video suffice, such that the next location of interest can be recorded. The interquartile range (IQR) of depth errors decreased from the first update, with minimum and maximum values of 0.9 m (Duck) and 2.6 m (Narrabeen), respectively, to the last update, with values of 0.5 m (Duck) and 1.3 m (Porthtowan). Absolute depth biases were small throughout all updates, slightly improving from minimum and maximum values of respectively 0.1 m (Scheveningen) and 0.4 m (Porthtowan) at the first update, to 0.0 m (Duck) and 0.3 m (Porthtowan) at the last. Maps of wave celerity, $\mathit{c}$, and near-surface currents, $\mathit{U}$, could not be validated at this stage; however, qualitative assessment showed vector fields of $\mathit{c}$ to match observed wave propagation and vector fields of $\mathit{U}$ to match expected tide- and wave-induced currents. The algorithm finally demonstrated its potential for on-the-fly video feed analysis by taking computational processing times into account. Therewith the groundwork is laid for a fast and easy-to-use tool for coastal reconnaissance. Striving towards universal applicability, more field cases and community-driven development are desired, wherefore the code and data are openly available at https://doi.org/10.4121/c.5704333 (accessed on 22 November 2021).