Next Article in Journal
Global Land High-Resolution Cloud Climatology Based on an Improved MOD09 Cloud Mask
Next Article in Special Issue
Development of a Fully Convolutional Neural Network to Derive Surf-Zone Bathymetry from Close-Range Imagery of Waves in Duck, NC
Previous Article in Journal
Multi-Agent Deep Reinforcement Learning for Online 3D Human Poses Estimation
Previous Article in Special Issue
LIDAR Scanning as an Advanced Technology in Physical Hydraulic Modelling: The Stilling Basin Example
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Updates to and Performance of the cBathy Algorithm for Estimating Nearshore Bathymetry from Remote Sensing Imagery

1
College of Earth, Ocean and Atmospheric Sciences, Oregon State University, 104 Ocean Admin Bldg, Corvallis, OR 97331-5503, USA
2
Earth Observation Lab, Centre National d’Etudes Spatiale (French Space Agency), 18 Avenue Edouard Belin, CEDEX 9, 31401 Toulouse, France
*
Author to whom correspondence should be addressed.
Remote Sens. 2021, 13(19), 3996; https://0-doi-org.brum.beds.ac.uk/10.3390/rs13193996
Submission received: 27 July 2021 / Revised: 29 September 2021 / Accepted: 1 October 2021 / Published: 6 October 2021
(This article belongs to the Special Issue Advances in Remote Sensing in Coastal and Hydraulic Engineering)

Abstract

:
This manuscript describes and tests a set of improvements to the cBathy algorithm, published in 2013 by Holman et al. [hereafter HPH13], for the estimation of bathymetry based on optical observations of propagating nearshore waves. Three versions are considered, the original HPH13 algorithm (now labeled V1.0), an intermediate version that has seen moderate use but limited testing (V1.2), and a substantially updated version (V2.0). Important improvements from V1.0 include a new deep-water weighting scheme, removal of a spurious variable in the nonlinear fitting, an adaptive scheme for determining the optimum tile size based on the approximate wavelength, and a much-improved search seed algorithm. While V1.2 was tested and results listed, the primary interest is in comparing V1.0, the original code, with the new version V2.0. The three versions were tested against an updated dataset of 39 ground-truth surveys collected from 2015 to 2019 at the Field Research Facility in Duck, NC. In all, 624 cBathy collections were processed spanning a four-day period up to and including each survey date. Both the unfiltered phase 2 and the Kalman-filtered phase 3 bathymetry estimates were tested. For the Kalman-filtered estimates, only the estimate from mid-afternoon on the survey date was used for statistical measures. Of those 39 Kalman products, the bias, rms error, and 95% exceedance for V1.0 were 0.15, 0.47, and 0.96 m, respectively, while for V2.0, they were 0.08, 0.38, and 0.78 m. The mean observed coverage, the percentage of successful estimate locations in the map, were 99.1% for V1.0 and 99.9% for V2.0. Phase 2 (unfiltered) bathymetry estimates were also compared to ground truth for the 624 available data runs. The mean bias, rms error, and 95% exceedance statistics for V1.0 were 0.19, 0.64, and 1.27 m, respectively, and for V2.0 were 0.16, 0.56, and 1.19 m, an improvement in all cases. The coverage also increased from 78.8% for V1.0 to 84.7% for V2.0, about a 27% reduction in the number of failed estimates. The largest errors were associated with both large waves and poor imaging conditions such as fog, rain, or darkness that greatly reduced the percentage of successful coverage. As a practical mitigation of large errors, data runs for which the significant wave height was greater than 1.2 m or the coverage was less than 50% were omitted from the analysis, reducing the number of runs from 624 to 563. For this reduced dataset, the bias, rms error, and 95% exceedance errors for V1.0 were 0.15, 0.58, and 1.16 m and for V2.0 were 0.09, 0.41, and 0.85 m, respectively. Successful coverage for V1.0 was 82.8%, while for V2.0, it was 90.0%, a roughly 42% reduction in the number of failed estimates. Performance for V2.0 individual (non-filtered) estimates is slightly better than the Kalman results in the original HPH13 paper, and it is recommended that version 2.0 becomes the new standard algorithm.

1. Introduction

Coastal bathymetry and its change over time are of paramount importance to understanding and modeling coastal morphology and the health of natural coastal systems. In particular, bathymetry acts as the critical bottom boundary condition for numerical models of nearshore hydrodynamics and sediment transport and, hence, is vital to predictions of coastal evolution needed by coastal communities to understand their vulnerabilities and mitigation options and, thereby, to enable decision support and policy making [1,2,3]. In terms of drivers or boundary conditions, hydrodynamics such as wave conditions or tides are often measured frequently or can be approximated relatively accurately through wave or tidal models. Bathymetry, on the other hand, is often unknown or outdated, such that numerical models, that are all sensitive to bathymetry accuracy, will perform poorly and cannot be trusted [4,5]. This is a particular problem, as the focus shifts to coastal response of the world’s coastlines to climate change [6,7]. There are few solutions to this dilemma. Due to their expensive nature, boat surveys are rare and are limited to calm conditions, so bathymetry data are usually sparse and outdated and, therefore, not responsive to the continuing bathymetric changes that occur on the time scales of events, weeks, and seasons.
The most promising solutions to these limitations are through the increasingly popular use of remote-sensing techniques [6,8]. Different applications and problem scales are compatible with different methodologies. For example, airborne LiDAR technology can be quite accurate, O (cm), and can cover large zones, but remains relatively costly so is limited in temporal update and by turbid coastal waters [9,10]. Satellite-derived coastal bathymetry observations [11,12,13,14,15,16,17] or even sampling of the topo-bathy continuum [18] can be carried out globally but with limits on repeat times due to orbital constraints, issues of atmospheric dilution of quality, and overall lower accuracy, O (m), than in situ methods. However, these global methods are entering a new era of increasing accuracy and public availability [19,20]. By comparison, local bathymetry measurements using video cameras that are either land-based or mounted on inexpensive drones or using land-based radar observations can provide low cost, complementary, near-continuous, near real-time, bathymetric observations, O (dm), enabling the monitoring of high-frequency morphological changes such as storm impacts and providing information that is useful to coastal zone managers [4,17,21,22,23].
Video cameras at the coast were initially used to identify the morphology and cross-shore scales of submerged sand bars [24,25,26] and to track the shoreline in time to construct the intertidal beach [27]. Around the early 2000s, bathymetry estimation using land-based video cameras [28], airborne deployments [29,30], and X-Band radar [31] started to develop. Since the early 2000s, methods have improved significantly: estimations became more accurate, consistent, and more versatile, hence, accessible, enabling for example drone applications [17,32,33,34]. In general, there are three distinctive approaches, ones that use numerical models to match wave breaking patterns [35,36,37], others that use time-varying signatures and stay in the spatio-temporal domain [34,38] and still others that exploit the time signal but with analysis in the spectral domain [34,38]. The latter two approaches are linked to a great extent and use the linear dispersion relation for free surface waves to estimate a local water depth [39]. Here, we focus on the spectral approaches, and in particular, on the cBathy algorithm ([1]; this Matlab toolbox is available as open source at https://github.com/Coastal-Imaging-Research-Network/cBathy-Toolbox, accessed on 20 July 2021). While this manuscript is based on video imagery, it is recognized that this method and improvements are equally applicable to X-band radar e.g. [40,41]. cBathy has been incorporated into operational numerical modeling systems [42] and has been applied in mixed wave-current domains of a tidal estuary [43]. The growth of this community of practice has led to the formation of the Coastal Imaging Research Network (CIRN), a web-based group that shares skills and developments among worldwide researchers.
Since the original 2013 publication [1], cBathy has gone through a number of revisions but typically without extensive testing and documentation. The original form of the algorithm is now called version 1.0 and has been followed by versions 1.1, 1.2, and now a composite update, which we choose to call version 2.0. The goal of this paper is to document these revisions and test the different versions of the algorithm against a new testbed dataset consisting of 39 surveys spanning from 2015 to 2019 from the Field Research Facility at Duck, NC, USA. Below, we will start by summarizing the elements of the original 1.0 algorithm plus the various significant updates. We, then, introduce the test bed dataset, which can be used to test both cBathy and any other algorithm in nearshore remote sensing. Finally, we will describe and execute a testing protocol that was used to document performance changes with each algorithm update.

2. cBathy Versions—Version 1.0

Version 1.0 is fully described in HPH13, but is summarized herein. For more detail, the reader is referred to the original manuscript. The algorithm is based on the linear dispersion relationship
σ 2 = gk   tan h kh
that relates the wavenumber, k, and frequency, σ, of ocean waves to the water depth, h (k is 2π divided by the wavelength, L, while σ is 2π divided by the wave period, T).
The input data for cBathy are usually derived from cameras and consists of time series of image intensity, I (xp, yp, t), at a set of discrete pixel locations, xp, yp that span the domain of interest and adequately sample the typical ocean wave scales, while not oversampling them (across-shore and alongshore spacing is commonly 5 and 10 m, respectively, see blue dots in Figure 1). The analysis is carried out at a map array of model locations, xm, ym, (example red dot in Figure 1) and at each map location is based on the observed wave phases in a tile of observations within some user-specified cross-shore and longshore length scales, Lx and Ly, of the model location (green dots in Figure 1; see Figure 2 for an example phase map). At a set of dominant radial frequencies, σ, the two components of wavenumber, kx and ky, are derived from the phase ramp slopes in the x (cross-shore) and y (alongshore) directions and are combined to yield the magnitude of k.
Thus, the goal is to estimate the dominant wave frequencies and wavenumbers at each model location and merge that information into a spatially smooth estimated bathymetry using the dispersion relationship. The algorithm is divided into three phases: (1) estimation of frequency–wavenumber pairs, (2) estimation of bathymetry, h ^ , from suites of those estimates, and (3) Kalman smoothing of separate hourly estimates to create a running average bathymetry product, h ¯ , that is robust to noise and error. In the following, we sometimes use the terms f-combined and f-averaged to refer to phase 2 and 3 results, respectively. The updates since 2013 have all addressed improvements to phases 1 and 2, so the algorithm review below will omit phase 3 Kalman filtering. Note that the performance tests in HPH13 only considered phase 3 results, since phase 2 output was noisy.
The algorithm must adapt to incident waves with periods from almost 20 s to about 4 s (shorter waves than this are insensitive to depth in other than very shallow water so are not helpful) without knowing a priori what conditions will be. Because the analysis is performed in frequency space, the first step is to Fourier transform each pixel time series and normalize the Fourier magnitudes, allowing focus on Fourier phase only. Cross-spectra are computed between all pixels in a tile for a suite of candidate frequencies that are usually spaced to allow about 40 degrees of freedom (for a typical run length of 1024 s, frequencies are spaced by 0.02 Hz). The (usually four) dominant frequencies are those with the largest total coherence in the resulting cross-spectral matrix and a wavenumber is, then, estimated for each of these frequencies.
Wavenumbers are found from maps of Fourier phase for each tile (Figure 2). To guard against cases of waves from multiple directions within a single frequency band, as a first step, the cross-spectral matrix is decomposed into complex eigenvectors, v(xp, yp), and only the one with the largest eigenvalue is retained. This normalized eigenvector is modelled as a single dominant plane wave form
v = e x p i k c o s α x p + k s i n α y p + φ
where α is the wave angle, and ϕ is a scalar phase angle. In version 1.0, the three parameters k, α, and ϕ were found using a weighted nonlinear least-squares fit between the observed and modeled eigenvector phase maps. Details of the weighting are included in the original paper [1].
The goal of phase 2 of the algorithm is to use a suite of σ-k phase 1 estimates to estimate the depth at each individual model location. Each model point can yield up to four frequencies and phase 2 incorporates estimates from adjacent xm, ym locations from within the tile, weighted by the inverse distance from the current estimation point. In version 1.0, the weighting was taken to depend on the normalized eigenvalue and skill of the model fit, both from phase 1, as well as the inverse distance from the current estimate location in phase 2. The final phase 2 result is the single depth that is the nonlinear best fit to predicted depth using the input suite of frequency–wavenumber pairs and the dispersion relationship. Error estimates, herr, are produced for each depth (see [1], error estimation is unchanged in recent upgrades).

2.1. Version 1.1 Update

Versions 1.1 and 1.2 were consecutive steps in cBathy improvement but were never fully tested or used separately. Thus, we will describe each, then combine them into a single upgrade to cBathy version 1.2 for testing. We now describe the algorithm improvements.
As waves are measured in increasingly deep water, they become decreasingly sensitive to depth. Thus, small errors in the measured wavenumber become associated with large errors in estimated depth. To guard against this excessive sensitivity in version 1.0, values deeper than a user-specified maximum depth were neglected. However, this maximum depth should instead be wave-frequency dependent, and biases are introduced in deeper water results by simply truncating estimated values. Version 1.1 corrected this problem by accounting better for the dispersion relation sensitivity in phase 2 weighting. Equation (1) can be rewritten
h = 1 k t a n h 1 σ 2 g k = 1 k t a n h 1 k 0 k  
where k0 is the wavenumber in deep water. We can define k0/k = L/L0 = Γ, where Γ is a non-dimensional wavelength (or wavenumber), going from small in shallow water to a maximum value of 1.0 in deep water. We can define sensitivity to wavenumber error by the equation
Δ h h = μ Γ Δ k k  
where μ is the sensitivity, a function only of Γ, which can be found numerically and represents the fractional sensitivity of bathymetry estimates to fractional errors in estimated wavenumber. Figure 3 shows the sensitivity of wavenumber inversion as a function of non-dimensional wavelength. In shallow water, the value is 2.0 so that fractional bathymetry errors are twice the magnitude of fractional wavenumber errors. This is the best case. Closer to deep water (Γ approaching 1.0), the sensitivity rises rapidly, approaching 10 when wavelengths are 0.9 of their deep-water value.
This sensitivity provides a convenient weighting measure. If we denote the version 1.0 phase 2 weighting value as W for any σ-k pair, we can divide this weight by the sensitivity as a simple method of preferentially weighting σ-k pairs that are in shallower water and de-emphasizing those approaching the deep-water limit. Thus, the main change in the version 1.1 algorithm is this modified weighting in phase 2. A side benefit is that the user no longer needs to specify a maximum depth beyond which to no longer believe an estimate. Note that L0 = gT2/2π is frequency dependent so this weighting puts a strong preference on long-period waves.

2.2. Version 1.2 Update

Version 1.2 dealt mostly with a problem of poor bathymetry estimates along the seams between cameras, an issue that only applies to image collection systems whose sampling is distributed over multiple cameras. Camera image data that are a common input to cBathy analysis are often collected from multiple cameras and merged into a map of data coverage, for example, as shown in Figure 1. When these data were analyzed as a single array, estimated bathymetries along camera seams were often anomalous.
Two causes were identified. The most obvious issue is variations in camera geometry used to map from image to world coordinates [45]. The sampling pixel array is usually designed as a regular world grid, mapped to pixels for each camera using originally accurate camera geometries. If camera geometries shift slightly over time, the projected world spacing of those pixels on either side of a camera seam will begin to differ from the original spacing. Thus, waves will appear to move too quickly across a shortened sample gap or too slowly for a stretched gap, leading to errors in estimated depth.
The second plausible cause of cross-seam anomalies can come from a loss of synchronization of the cameras. Argus cameras are usually synchronized by either an electronic trigger or by computer bus synchronization. However, there can be rare frame slips that allow time shifts between cameras that would be interpreted as wave-speed anomalies for tiles that span camera seams.
Initial attempts were made to mitigate both of these problems. However, it soon became clear that the simple solution was to never mix pixels from multiple cameras in any tile [46,47]. Thus, the pixels from the camera with the most pixels in any tile are used, and the others are simply neglected. This results in fewer degrees of freedom and increased predicted error, but otherwise had no significant negative impact.
Version 1.2 also included a new method for finding the seed wavenumber and wave angle for the nonlinear search, but this change has been overtaken by a much better seed routine in version 2.0 (discussed below) and will not be discussed further here.

2.3. Version 2.0

The upgrade to version 2.0 involves three significant changes and a code reorganization, so is considered a major revision and is given a new leading version number. In order of importance, the changes are: (a) to change from tile sizes that are fixed by the user to those that are chosen automatically depending on expected wavelengths, (b) to change from solving for three variables, k, α, and ϕ0, at each map point to solving for only k and α, and (c) the introduction of a much better algorithm to find seed values for k and α before the nonlinear search for each tile. The cumulative consequence of these modifications is a major restructuring of the code. Each component will be described in turn.

2.3.1. Automatic Tile Sizes

In all earlier versions, cross-shore and alongshore tile sizes were set by the user using the parameters, Lx and Ly. Suggestions for best values were ad hoc with a belief that the search would work best if the tile was typically about one wavelength long, but an implicit faith that even mismatched tile sizes would solve well somewhere within the nonlinear fitting routine. Thus, the same tile size was used for, e.g., 4 s and 16 s waves despite at least a four-fold difference in their wavelengths. Although sub-optimal, this approach was still successful, since tiles that were poorly designed simply failed to converge and returned nan’s rather than a poor depth.
Issues with this approach became especially clear for short-period incident waves, for example, 4–5 s waves. These could have 4–5 wavelengths within the tile sampling region, so convergence of the nonlinear search usually behaved very poorly, since there was not a simple global cost function minimum unless you started with a very accurate seed. The problem was made worse by a feature of the code that decimates the original number of pixels in a tile down to a user-defined maximum, maxNPix, to help reduce processing time. Commonly, a standard tile was reduced from 250 collected pixels down to 80 that would be analyzed, enough to map the phase of a single typical wavelength but way too few to map out five short wavelengths in a tile. The resulting sparse sampling per wavelength often led to aliasing of the true signal.
Version 2.0 fixes these problems by defining the phase 1 tile sizes to be kL times the expected wavelength, where kL is an empirical scalar taken to be approximately 1.0. However, the expected wavelength depends of the frequency and depth, neither of which is known a priori. Thus, there is strong motivation to develop a seed-finding algorithm (see below) that can provide good initial estimates of k under all wave conditions. Thus, given an initial tile based on generic user input, version 2.0 feeds all of the available pixels to the routine to find the seed k and α, then crops the original tile size to kL times this wavelength, a size that varies considerably. The number of pixels in this truncated tile is, then, reduced, if necessary, to maxNPix to speed up processing. Because all tiles are roughly one wavelength no matter what the frequency, it is likely that fewer pixels are needed for search convergence, again speeding up processing.

2.3.2. Reduction in the Number of Search Variables

In all earlier versions, the model for wave phase, Equation (2), had three parameters, each of which was found using a standard Levenburg–Marquardt solver. However, only two of the parameters, k and α, are expected to be well behaved in a nonlinear gradient-descent search. The third variable, ϕ0, is a scalar offset between the measured and modeled phase maps and can jump around a great deal, in ways that are inconsistent with a search for a global cost function minimum. Thus, this variable is not well estimated and likely just confuses the search.
In a very early version of this algorithm known as Beach Wizard [35], the solution was sought not in x-y space phase maps as here, but in cross-spectral lag space. This had the advantage that a phase offset was not required (phase differences just increased with lag from zero), but it meant that if you had N pixels in your tile you were searching in an N2 space (each pixel compared to every other pixel). In addition, visualization of measured and modeled results (e.g., Figure 2) was not as clear in lag space. Thus, we wish to retain the simplicity of working in x-y space maps but using a method for estimating ϕ0 for each search iteration that will allow a sensible search for k and α. The solution is to force the measured and modeled phase to be the same at the tile center (the pixel closest to xm, ym). This is done by finding dϕ, the difference in phase at the middle pixel, and multiplying all modeled complex values of the eigenvector, v, by eidϕ. The nonlinear search is, then, reduced to two dimensions.

2.3.3. Improved Seed Algorithm

Both the success of adaptive tile sizes and of the nonlinear search depend on the accuracy of the initial seed search values for k and α. The search is complicated by the fact that Fourier phase has 2π jumps whose slight mispositioning adds to cost functions out of proportion to the error. It also rapidly became clear that the full resolution of the tile was needed during the seed search, i.e., decimation down to maxNPix could only occur after the search seeds were found and the full tile was sub-sampled to the adaptive tile size.
The solution for finding more robust seed values makes use of the Radon transform. First, the observed phases from the eigenvector were interpolated onto a regular grid, a somewhat noisy process due to the phase jumps. Next, the Radon transform was found to estimate the phase variance when projected along a suite of candidate wave directions (−45 to +45° in 2° increments; Figure 4). The angle that corresponded to the maximum variance was selected as the seed wave angle, α. Finally, the phase map was interpolated onto a new grid that is oriented in the direction of wave propagation, and the median of the phase gradient in that direction was used as the seed value of k. This search also returned the location of the pixel that lay closest to the center of the tile (xm, ym), to be used in finding dϕ.

2.4. Algorithm Organization Changes

The changes in version 2.0 have forced a major restructuring of the program logic flow. Figure 5 and Figure 6 show flow charts for versions 1.2 (also representing version 1.0) and 2.0, respectively. The reorganization only impacts phase 1 of the algorithm, shown in red boxes, with the partitioning into tiles happening early in version 1.2 (in subBathyProcess) but later in version 2.0 (in csmInvertKAlpha).
Instead of reducing the number of pixels prior to any analysis steps, it was clear that full pixel resolution was required for finding the k-α seeds and for establishing the required size of the tile. Thus, the full tile of size Lx, Ly is initially passed into the main analysis routine, csmInvertKAlpha, where the routine prepareTiles is called to (a) find the dominant frequencies, (b) for each frequency, find the dominant eigenvector and the k-α seeds, then c) reduce the tile to an adaptive size and to a maximum number of pixels. These outputs are, then, passed back to csmInvertKAlpha to carry out the nonlinear search for best-fit values and their errors and, then, to build the results structures and find the depths from the σ-k results.
Several new fields have also been added to the bathy output MATLAB structure. These provide a set of traditional Argus image proxies, computed from the input time stack data over the region of interest at the array of desired xm, ym points. These include time exposure, brightest and darkest images, all at the coarse resolution of the analysis array but adequate to determine, for example, the locations of wave breaking. Within the phase 1 fDependent sub-field, there are now maps that are included for diagnostic and performance improvement purposes. These include maps of the seed values of k and α used for the nonlinear search as well as maps of the number of pixels used in each tile and of the number of model calls during the nonlinear fitting routine, a key to algorithm speed. In addition, there is a map of which camera is used for each tile. Within the fCombined sub-field, there is now a map of the effective mean frequency, fBar, used in the phase 2 bathymetry calculation, found as the weighted mean of each frequency that contributed to the phase 2 final bathymetry estimate. Finally, the elapsed CPU time for each analysis is saved.

3. Bathymetry Test Bed Datasets

HPH13 tested the cBathy version 1.0 against a collection of 16 bathymetry surveys collected at the Field Research Facility (FRF) between 2009 and 2011 as well as one survey at Agate Beach, Oregon, on the US West Coast. In 2015, the Argus station at the FRF was upgraded with improved cameras and computer hardware, so it was decided to base the updated tests in this paper on more recent surveys.
The new dataset includes 39 FRF surveys collected from 20 May 2015 to 17 April 2019. The full dataset is described in the document “The2019cBathyDataTestBed”, which is located on the CIRN (Coastal Imaging Research Network) GitHub site (https://github.com/Coastal-Imaging-Research-Network) in the cBathy toolbox and at the date of publication is located in the version 2.0 branch in a folder called cBathyTestBed (this branch should become the master branch in the near future). The MATLAB data structure, bathyTestDataSet, includes the full set of actual survey data, the list of related cBathy stack files, which can be loaded directly from the Coastal Imaging Lab, and the environmental data on waves and tides. These are accompanied by a testbed description document and instructions on how to download cBathy time stacks. To be consistent with HPH13, the cBathy files for each survey include all collections for the day of the survey plus for each of the three preceding days, a total of 96 time stacks for each survey.
In HPH13, cBathy results were computed for each time stack and, then, smoothed into the phase 3 Kalman-filtered estimate—only the Kalman results were compared to ground truth, since the individual phase 2 results were viewed as too noisy. In this paper, we will examine Kalman results as did HPH13, but we will also compare individual phase 2 results from a large subset of runs within each four-day cBathy sampling window. Analysis runs were selected at 1200, 1400, 1600, and 1800, local standard time, for each of the four days, avoiding morning sun glare problems. This yielded a total of 624 phase 2 bathymetry map estimates in addition to the 39 Kalman smoothed results. Survey data, while only collected on the fourth day, are assumed to be representative for the full four-day sample periods.
Figure 7 shows a plot of the heights, periods, and directions of waves measured at both the 26 m depth waverider and the 8 m array seaward of the FRF property. Wave heights were typically fairly low, compatible with the safety needs of in situ surveying, although some of the surveys followed major storm events so include larger wave conditions in the collections from the three prior days. Tidal data are taken from the pier-end NOAA tide station and span a full range of almost 2 m.

4. cBathy Phase 2 Version Performance Statistics

Figure 8 compares cBathy results from the different versions for an example day, 19 November 2017, 1959 GMT, chosen as representative because the rms error (0.41 m) was nearly equal to the median of all 624 examples (0.46 m) and because there was some breaking over the offshore sand bar (wave breaking has been found to affect performance). The 26 m waverider wave height was 0.72 m and the period was 6.0 s, the wave direction at the 8 m array was 78° (shore normal is 72°, so this is from 6° south of normal), and the tide elevation was 0.05 m. Bathymetry estimates for each cBathy version (left three panels) are thresholded to only show those for which the predicted error is less than or equal to 0.5 m. For comparison, the right-hand panel shows the ground-truth survey, taken two days later. Ground-truth values are only used if the loess interpolation error in the gridded data is less than 0.15, omitting some larger interpolation errors where there were occasional gaps in survey coverage. Similarly, ground-truth values are only used if the true water depth (survey depth plus tide) was greater than zero, i.e., we only consider locations that are wet. Consistent with HPH13, performance statistics excluded the region of the FRF pier, i.e., the region between y = 400 and 600 m (not inclusive), which interferes with both the camera views and the waves.
Similar to HPH13, we characterize the performance of the three algorithms in terms of the bias, the rms error (rmse), the 95th percentile absolute error exceedance value (Δh95), and the percentage of successful coverage (estimated error <0.5 m) for estimation locations with survey depths greater than 0 m, including tide (called coverage). The results from this example run show a steady improvement in all statistics for versions 1.0 to 1.2 to 2.0 of: bias = [−0.36, −0.05, −0.04] m, rmse = [0.90, 0.59, 0.41] m, Δh95 = [2.58, 1.30, 0.86] m, and coverage = [55%, 77%, 91%].
Previous papers have observed that cBathy estimates are worst in the surf zone and particularly at the onset of breaking [21]. The onset of breaking can be detected based on coarse-resolution time exposure image (one of the new fields in version 2.0), by defining the breaking onset as any region whose intensity gradient, averaged over 20 m in the cross-shore, is less than −2.0 intensity units per cross-shore meter. Neither V1.0 nor V1.2 had any successful estimates in the breaker region (dark red regions landward of roughly x = 200 m). For V2.0, bias and rms errors were 1.20 m and 1.24 m, respectively, reinforcing the poor cBathy performance statistics in this region (for this example run).
HPH13 also found that performance varied with true depth. cBathy estimates were partitioned by depth as 0–1 m, 1–2 m, …, 7–8 m and performance estimates found for each depth bin (Figure 9). As expected, performance degrades in the shallowest two bins in the surf zone (note that for the shallowest bin, no estimation points succeeded for V1.0, and only 3 and 15 points were successfully estimated for V1.2 and V2.0, respectively). Both the bias and rmse are good for depths greater than 2 m (outside the surf zone) for versions 1.2 and 2.0 and are similar to those from Kalman-filtered estimates found by HPH13 for version 1.0. Version 2.0 performs best for all statistics, but particularly for the percentage of successful coverage.
Several new fields in version 2.0 can be explored. Most interesting is fBar, the weighted mean frequency of waves that contributes to the phase 2 f-Combined depth estimate. Figure 10 shows these results for the example day, showing that map areas to the north (large y-values) are more dominated by shorter northeast waves (f = 0.21 Hz, α ~ 30°, north of normal), while views looking offshore and south from the cameras in the center region saw longer period waves from the southeast (f = 0.17 Hz, α ~ −15 to −20°, south of normal). In general, fBar favors lower frequencies in a broad spectrum. This result is consistent with predictions by Walker [48] and observations by Holman et al. [32] that the strongest optical signals occur when viewing directly into oncoming waves.
Finally, we can compare the seed values of k and α versus final nonlinear fit values as a test of the effectiveness of the new seed algorithm (Figure 11). Eighty-four percent of wave angle fit values are within 10° of the seed, while k seeds are less accurate with only 60% of best-fit values being with 50% of the k seeds (which seems to be good enough for successful nonlinear searches, given good seed angles).

4.1. Bulk Analysis of cBathy Version Statistics

In all, 624 cBathy analyses were considered including 39 surveys and 16 cBathy collections for each survey (a subset of the full set of collections). Because each of the three cBathy versions returned results over different regions (i.e., had different coverages), it is harder to directly compare statistics such as bias and rms error. Thus, we will first discuss the performances of the algorithms in terms of coverage, then we will discuss other error statistics considering both the full comparison as well as comparisons based on coverage regions that were in common between algorithms.
Figure 12 shows histograms of the bulk coverage statistics for each cBathy version. It is apparent that V2.0 does much better than earlier versions; for example, the median coverage values are 86, 87, and 95%, for V1.0, V1.2, and V2.0, respectively. All algorithms provide decent coverage with only about 41 of the 624 runs having less the 50% success. There were several reasons for low successful coverage. Figure 13 shows the dependence of coverage on the significant wave height at the 8-meter array. Performance clearly declines with wave height with a rough rule of thumb that for Hs > 1.2 m, expected performance will become poor (consistent with Brodie et al. [21]). Snapshots from the 41 cases of coverage worse than 50% were examined manually and showed that poor performance was also caused by foggy or rainy days (although many rainy days were also successful) as well as low evening light.
We next wish to compare the basic performance statistics between the three versions. Figure 14 shows histograms of each statistic for each version. Surprisingly, the first three statistics, bias, rmse, and Δh95, look almost indistinguishable between versions, despite what we had hoped were significant improvements in the algorithm as shown in the example run above. This result is a consequence of the fact that for each version, results with estimated errors, herr > 0.5 m were rejected from the analysis. Thus, the statistics are computed for different regions of coverage. To make a more direct comparison (apples versus apples), the map of successful coverage for V2.0 was saved for each collection and used as the basis for computing performance statistics for V1.0 and V1.2. This usually larger region, thus, includes poorer performing regions from the earlier algorithm versions. This is confirmed by Figure 15, which shows the same histograms of statistics but is now based on a common area of sampling defined by the V2.0 error estimates. For all three statistics, bias, rmse, and Δh95, the performance of V2.0 is superior. The fact that statistics are roughly the same when herr is determined by each individual version is a testament to the robustness of that error estimate. In the following, all statistics will be based on the common area of coverage determined for V2.0.
The bias statistics are seen themselves to be slightly positively biased (estimates too deep), with mean biases for the three versions of 0.19, 0.14, and 0.16 m, respectively. There appear to be two causes for larger biases: performance issues with large waves and breaking and cases of low coverage, so poor imaging conditions or noisy statistics. Removing cases for which Hs is greater than 1.2 m reduces the mean bias to 0.15, 0.11, and 0.12 m, respectively. Removal of estimates for which coverage was less than 50% reduced the original mean bias to 0.15, 0.1, and 0.1, respectively, while removal of both issues (consider only cases with Hs ≤ 1.2 m and coverage ≥ 50%, reducing the dataset from 624 to 563 runs) led to mean biases for the three algorithm versions of 0.15, 0.09, and 0.09 m, respectively.
Ninety-five percent of the bias statistics had magnitudes less than 0.61, 0.58, and 0.59 m for the three cBathy versions when the full dataset was considered. When large waves (>1.2 m) and low coverage (<50%) data were removed, the 95% exceedance of bias was reduced to 0.46, 0.43, and 0.34 m, respectively. As a further step, the most extreme cases of error were examined visually. Most cases were associated with fog or rain (low coverage), storms (large waves), or an apparent absence of longer period (>4 s) waves.
Figure 15 shows that the rms and 95% exceedance errors also improved with later versions of the cBathy algorithm. Table 1 summarizes the mean values for the four basic performance statistics for each of the three cBathy versions. These statistics are shown for the full dataset (top of Table) as well as for the reduced dataset (bottom of Table). It is clear that version 2.0 is superior in all categories of performance measurement.
Figure 16 shows the four bulk statistics against significant wave height, plotted by color for each cBathy version, for all 624 data runs. Points are plotted in order of V1.0 (red), V1.2 (green), then V2.0 (blue), so later values hide earlier ones. That said, it is clear that V2.0 values are better than those of V1.0 and V1.2. This is particularly apparent for rmse statistics where there are fewer cases of large error for V2.0 (fewer blue dots at large values). For comparison, Figure 17 shows the same statistics for the reduced dataset of 563 points for which wave height is less than or equal to 1.2 m, and successful coverage is greater than 50%. The y-scales are the same as for Figure 16 to ease comparison, although the range of wave heights (x-axis) is reduced.
As was shown by the example run in the previous section, statistics vary with true depth. The full dataset was partitioned by depth (0–1 m, 1–2 m, …, 7–8 m) and the mean of each of the core statistics computed for each depth bin. The results are shown in Figure 18. It is apparent that the improvement in V2.0 performance is concentrated in the shallow (<2 m) and deep (>6 m) regions with much smaller differences in intermediate depths. It is unclear why the performance differences are greater in deeper water. However, in shallow water, the worsening performance is likely linked to wave breaking in the surf zone. To test this hypothesis, we partitioned the dataset into breaking and non-breaking regions for each data run, based on the mean intensity brightness from the time exposure image (calculated from the time stack in version 2.0 but applicable to all versions). Time exposure intensities greater than 150 were considered to be surf zone, while values less than 150 were considered non-breaking. Figure 19 shows the bias and rms statistics for the breaking (dashed line) and non-breaking (solid line) partitions, reinforcing the increased errors when waves break, which always occurs preferentially in shallow depths. Note that while errors are largest in shallow depths for all versions of the algorithm (Figure 18), they are significantly better for version 2.0.
Finally, the quality of the new k-α seed algorithm was tested for version 2.0 (seeds are not saved for earlier versions). Over the 624 data runs, the standard deviation of the residual between the seed and the best fit value of wave angle, α, was typically between 5 and 10 degrees with a mean value of 6.2 degrees. The errors of the k seeds are relatively larger with a mean of 0.05 m−1, but with improved wave angle seeds, this was sufficiently close to yield greater success in the nonlinear search.

4.2. Kalman-Filtered Results

In the first paper describing cBathy (HPH13), only Kalman-filtered results were presented because individual results from version 1.0 were felt to be too noisy or to have too many analysis gaps. Kalman filtering allowed these gappy areas to be merged with data from other runs based on the estimated error, herr, using a Kalman filter. Thus, regions with bad results due to, for example, wave breaking, can be filled with results from runs in which there was no breaking at that location. Kalman results accumulate over a series of prior data collections, with each new run improving the estimate until the results stabilize after roughly one or two days. In this case, we have included four cBathy bathymetry estimates per day for four days, the three days prior to each survey plus the actual survey day. The Kalman result from approximately 3:00 pm on the day of the survey is used to compare with the survey ground truth. Figure 20 shows the Kalman-filtered results for the example survey used in Figure 8, 21 November 2017. Results from the three versions are similar but with slight differences in coverage and depth estimates, for example, the region around x = 200, y = 830, which appears anomalous in version V1.0 but not in later versions. All statistics are better after Kalman filtering compared to the individual cBathy estimates. Bias for the three versions (in order V1.0, V1.2, V2.0) was 0.006, −0.07, and −0.03 m, while the rmse was 0.41, 0.37, and 0.35 m, and the 95% exceedance was 0.85, 0.83, and 0.77 m, respectively. Coverage of locations that were immersed at the final time were all 100%, a credit to averaging over the varying tide elevations.
Continuing the comparison to the non-filtered results in the previous section, Figure 21 shows the four basic statistics plotted by depth. The advantages of Kalman filtering become apparent. While there are still some anomalies in the shallowest depth bin, they are largely absent for depths greater than 1 m. Moreover, coverage is now 100% for all versions.
These results are consistent with the bulk statistics for the 39 surveys, shown in Figure 22 and plotted versus wave height for the survey date data run. Bias, rmse, and Δh95 are all very good with V2.0 results outperforming earlier versions. Note that the coverage (of submerged depths) is now always greater than 95% due to the Kalman filtering, with V2.0 performing the best (the lowest V2.0 coverage was 99.3%). Again, performance is best for lower waves, although the Kalman filtering averages over multiple prior runs. Mean statistics for the three algorithms averaged over the 39 surveys are shown in Table 2. Again, version 2.0 outperforms earlier versions.
Finally, the effectiveness of predicted errors from V2.0 was compared to errors that were measured in comparison to ground truth. For each selected Kalman dataset (3:00 pm on the survey date for each of the 39 surveys), the ratio (labelled Rerr) of the measured absolute error to both the Kalman and non-Kalman (phase 2) predicted error was taken, averaged over the domain for which Kalman and non-Kalman predicted errors were less than or equal to 1.0 and 0.5 m, respectively (this removes the influence of extreme predicted errors). In HPH13, this value was only tested for Kalman results and was found to be approximately 7. For this more extensive dataset and the version 2.0 algorithm, the values are smaller, as shown in Figure 23. The mean of the Kalman-filtered results is 4.47, smaller than the value of 7.0 noted by HPH13. The mean of the non-Kalman-filtered results was 2.0. During this comparison, it was also noted that the process of Kalman filtering reduced the actual error of the Kalman estimates by an average factor of 4.6 compared to the non-Kalman single run estimates. The same analysis, carried out for version 1.0 results, showed that the mean of the predicted to measured errors for Kalman- and non-Kalman-filtered results were 5.22 and 2.29, with the Kalman ratio being a bit smaller than the factor of 7.0 found by HPH13 for their dataset. The reduction in the algorithm error for Kalman- and non-Kalman-filtered predictions was 5.3, slightly larger than with V2.0, likely due to the improvements in the non-Kalman estimates.

5. Discussion

Improvements to the cBathy algorithm contained in versions 1.2 and, then, 2.0 have led to significant improvements in cBathy performance. One of the main improvements is the larger region of successful coverage, defined as the fraction of the domain that yielded bathymetry estimated with a predicted error ≤0.5 m. The typical improvement in coverage of around 5–10% may seem small, but when viewed in terms of the fraction of locations for which estimates fail, the improvement is considerable, with a mean reduction in failed estimates of 28% for the phase 2 full dataset (42% for the reduced dataset) and 89% for the Kalman result (although failed estimates are rare with Kalman filtering).
Interestingly, the other performance statistics were not greatly changed by the improvement in the algorithm. However, this observation is only true when comparing dissimilar populations, i.e., the reduced dataset that is acceptable for V1.0 and V1.2 compared to the successful region for V2.0. When the comparison is made over similar regions (that which was acceptable to V2.0), the bias, rmse, and Δh95 all showed improvements with V2.0. In some ways, it is reassuring that the quality measurements (predicted herr) are performing sensibly.
With all of the algorithms, it was found that performance degraded when the waves were large or when the data were of poor quality for one or more reasons such as fog, rain, or sun glitter. The latter issues can be hard to objectively identify but can be represented by a simple test that if coverage is less than 50%, there are likely to be problems with the data collection, and it should be ignored completely. To remove issues of large wave breaking in the surf zone, we follow the recommendations of Brodie et al. [21] that data runs for which the wave height was greater than 1.2 m should be omitted from subsequent analysis.
This yields a potential concept of operations that could apply to a beach that is similar to Duck, NC. Analysis can be automated for all data collections, but for those with either a wave height greater than 1.2 m or an estimated data coverage less than 50%, the results should be excluded from subsequent Kalman filtering or from further analysis consideration. For other beaches that have no wave height measurements, it may be possible to introduce a quality filter that uses the coarse time exposure image (estimated in V2.0) to detect wave breaking and exclude those regions. For the tests above, a threshold of image intensity >150 was used to detect breaking, but this was ad hoc, and image intensity is affected by automated gain control so is not a fundamental variable. Additionally, there are beaches such as on the northwest coast of the US where waves are always large and surf zones are usually wide. Waiting for small waves at such sites may not be practical and further investigation is needed on a useful concept of operations. For beaches with wide surf zones, it may be necessary to use either a more complete numerical model approach to address finite-amplitude increases in the dispersion relation, or to use an approximate compensation that might be based on detected breaking in the time exposure image. Either is beyond the scope of this paper.
The data used in this study are excellent for cBathy testing in that the cameras are fixed to a tower (image geometry is fixed), the data runs are all 1024 s, a standard length for incident wave analysis, and many consecutive runs are available for Kalman filtering. However, there is increasing interest in less structured sampling methods, particularly using small commercial drones (e.g., [32,33]). For these applications, Kalman filtering is likely not available, so the improvements in phase 1 and 2 components of the algorithm that are embodied in V2.0 are important.
Finally, the CPU time was found to be influenced by the algorithm changes. All analyses were done on a modest Linux desktop machine with an Intel i7 CPU using MATLAB R2016b running four parallel threads. The three versions took an average of 92, 91, and 118 s to compute for each data run. The improved seed algorithm for k and α is likely responsible for the increased CPU time in V2.0.
It is recommended that V2.0 should supersede previous versions for all practical uses.

6. Conclusions

In 2013, Holman et al. ([1]; HPH13) published the new cBathy algorithm for the estimation of bathymetry from optical observations of wave celerity. The current manuscript describes a set of algorithm improvements that have been made since that time, labeling the original version as V1.0 and comparing it to an intermediate version, V1.2, that has seen considerable use but limited testing, and a current version, V2.0, that is significantly changed compared to the original. Tests of the performance improvements associated with each version were made using an updated dataset of 39 bathymetric surveys collected from 2015 to 2019 at the Field Research Facility in Duck, NC, each spanning a region of 500 and 1000 m in the cross-shore and alongshore, respectively. cBathy returns products of two types, a non-filtered, phase 2 bathymetry and a Kalman-filtered, phase 3 product that averages intelligently over a suite of phase 2 results that have variable quality. Since phase 2 results were viewed by HPH13 as too noisy or having too many gaps, they only presented Kalman results. For the 39 bathymetries of the new dataset, the performance of Kalman stage results improved between V1.0 and V2.0 with bias, rms error, and 95% exceedance error improving from 0.15, 0.47, and 0.96 m, respectively, for V1.0, to 0.08, 0.38, and 078 m for V2.0. The percentage of successful returns (predicted error ≤0.5 m) was 99.1 for V1.0 and 99.9% for V2.0.
With the increasing use of cBathy for isolated data collections, for example from drones or cameras of opportunity for which Kalman filtering is not possible, there has been increasing interest in the quality of individual phase 2 results. Phase 2 performance statistics were not published by HPH13 but have been computed for this new dataset including 624 phase 2 bathymetry estimate maps. For this full dataset, the bias, rms error, and 95% exceedance errors for V1.0 were 0.19, 0.64, and 1.27 m, respectively, and for V2.0, they were 0.16, 0.56, and 1.19 m. The average successful coverage was 78.8% for V1.0 and 84.7% for V2.0 (about a 28% reduction in the number of failed estimates). Thus, V2.0 performs almost as well as the Kalman predictions of V1.0
Larger errors were associated with both large waves (a wide surf zone) as well as a variety of poor imaging conditions such as fog, rain, or darkness. As a practical mitigation, it is recommended that at this site, conditions for which the significant wave height is greater than 1.2 m or for which successful coverage is less than 50% (a proxy for image quality problems) should be excluded from consideration. With this reduced dataset of 563 runs, the bias, rms error, and 95% exceedance errors for V1.0 were 0.15, 0.58, and 1.16 m and for V2.0 were 0.09, 0.41, and 0.85 m, respectively. Successful coverage for V1.0 was 82.8%, while for V2.0, it was 90.0%, an approximately 42% reduction in failed estimate locations. It is noted that this concept of operations for removing poor estimates is specific to Duck, NC. Finally, predicted errors for Kalman and non-Kalman results are found to be too small by a factor of 5.22 and 2.29, respectively, for V1.0 and 4.47 and 2.0 for V2.0.
It is recommended that V2.0 becomes the new standard for cBathy bathymetry estimation.

Author Contributions

The development and improvement of cBathy has occurred over decades. R.H. has done much of the work while E.W.J.B. was particularly involved with development and testing of the seam problem solution and the updated seed algorithm. Writing was led by R.H. with E.W.J.B. contributing particularly to the introduction material and in helpful proofreading. Both authors have read and agreed to the published version of the manuscript.

Funding

RAH would like to thank the Office of Naval Research, Coastal Geoscience Program, the Naval Research Lab, as well as the US Geological Survey for support over those many years. Over the years, EWJB was funded through the Quality Research/Research Excellence Framework CarePlan of the University of Plymouth, IRD and CNES post-doc funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The full dataset is described in the document “The2019cBathyDataTestBed”, which is located on the CIRN (Coastal Imaging Research Network) GitHub site (https://github.com/Coastal-Imaging-Research-Network) in the cBathy toolbox and at the date of publication is located in the version 2.0 branch in a folder called cBathyTestBed (this branch should become the master branch in the near future).

Acknowledgments

The development of cBathy from its roots has spanned decades and multiple research grants. RAH would like to thank the Office of Naval Research, Coastal Geoscience Program, the Naval Research Lab, as well as the US Geological Survey for support over those many years. Over the years, EWJB was funded through the Quality Research/Research Excellence Framework CarePlan of the University of Plymouth, IRD and CNES post-doc funding. cBathy is also a result of a community of researchers around the world. As always, we are grateful to the US Army Corps of Engineers for providing survey and wave data and for their continuing collaborations over many decades. Thanks to all.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Holman, R.A.; Plant, N.G.; Holland, K.T. cBathy: A robust algorithm for estimating nearshore bathymetry. J. Geophys. Res. 2013, 118, 1–15. [Google Scholar] [CrossRef]
  2. Davidson, M.A.; Van Koningsveld, M.; de Kruif, A.; Rawson, J.; Holman, R.; Lamberti, A.; Medina, R.; Kroon, A.; Aarninkhof, S.G.J. The CoastView project: Developing video-derived Coastal State Indicators in support of coastal zone management. Coast. Eng. 2007, 54, 463–475. [Google Scholar] [CrossRef]
  3. Kroon, A.; Davidson, M.A.; Aarninkhof, S.G.J.; Archetti, R.; Armaroli, C.; Gonzales, M.; Medri, S.; Osorio, A.; Aagaard, T.; Holman, R.; et al. Application of remote sensing video systems to coastline management problems. Coast. Eng. 2007, 54, 493–505. [Google Scholar] [CrossRef]
  4. Radermacher, M.A.; De Schipper, M.A.; Reniers, A. Sensitivity of rip current forecasts to errors in remotely-sensed bathymetry. Coast. Eng. 2018, 135, 66–76. [Google Scholar] [CrossRef] [Green Version]
  5. Sembiring, L.; van Dongeren, A.; WInter, G.; Ormondt, V.; Briere, C.; Roelvink, D. Nearshore bathymetry from video the application to rip current predictions for the Dutch Coast. J. Coast. Res. 2014, 70, 354–359. [Google Scholar] [CrossRef]
  6. Benveniste, J.; Cazenave, A.; Vignudelli, S.; Fenoglio-Marc, L.; Shah, R.; Almar, R.; Andersen, O.; Birol, F.; Bonnefond, P.; Bouffard, J.; et al. Requirements of a coastal hazard observing system. Front. Mar. Sci. 2019, 6, 348. [Google Scholar] [CrossRef] [Green Version]
  7. Melet, A.; Teatini, P.; Le Cozannet, G.; Jamet, C. Earth observations for monitoring marine coastal hazards and their drivers. Surv. Geophys. 2020, 41, 1489–1534. [Google Scholar] [CrossRef]
  8. Caballero, I.; Stumpf, R.P. On the use of Sentinel-2 satellites and lidar surveys for the change detection of shallow bathymetry: The case study of North Carolina inlets. Coast. Eng. 2021, 169, 103936. [Google Scholar] [CrossRef]
  9. Johnson, S.Y.; Cochrane, G.R.; Golden, N.A.; Dartnell, P.; Hartwell, S.; Chochran, S. The California Seafloor and Coastal Mapping Program—Providing science and geospatial data for California’s State Waters. Ocean Coast. Manag. 2017, 140, 88–104. [Google Scholar] [CrossRef] [Green Version]
  10. Lyzenga, D.R. Passive remote-sensing techniques for mapping water depth and bottom features. Appl. Opt. 1978, 17, 379–383. [Google Scholar] [CrossRef] [PubMed]
  11. Stumpf, R.P.; Holderied, K.; Sinclair, M. Determination of water depth with high-resolution satellite imagery over variable bottom types. Limnol. Oceanogr. 2003, 48, 547–556. [Google Scholar] [CrossRef]
  12. Poupardin, A.; Heinrich, P.; Hebert, H.; Schindele, F.; Jamelot, A.; Reymond, D.; Sugioka, H. On the role of frequency dispersion on the transw-Pacific tsunamis: Study of the 2010 and 2015 Chilean tsunamis. In Proceedings of the Oceans 2016 MTS/IEEE Monterey, Monterey, CA, USA, 19–23 September 2016. [Google Scholar]
  13. Poursanidis, D.; Traganos, D.; Chrysoulakis, N.; Reinartz, P. Cubesats allow high spatiotemporal estimates of satellite-derived bathymetry. Remote Sens. 2019, 11, 1299. [Google Scholar] [CrossRef] [Green Version]
  14. Caballero, I.; Stumpf, R.P.; Meredith, A. Preliminary Assessment of Turbidity and Chlorophyll Impact on Bathymetry Derived from Sentinel-2A and Sentinel-3A Satellites in South Florida. Remote Sens. 2019, 11, 645. [Google Scholar] [CrossRef] [Green Version]
  15. Pike, S.; Traganos, D.; Poursanidis, D.; Williams, J.; Medcalf, K.; Reinartz, P.; Chrysoulakis, N. Leveraging Commercial High-Resolution Multispectral Satellite and Multibeam Sonar Data to Estimate Bathymetry: The Case Study of the Caribbean Sea. Remote Sens. 2019, 11, 1830. [Google Scholar] [CrossRef] [Green Version]
  16. Bergsma, E.W.; Almar, R.; Maisongrande, P. Radon-Augmented Sentinel-2 Satellite Imagery to Derive Wave-Patterns and Regional Bathymetry. Remote Sens. 2019, 11, 1918. [Google Scholar] [CrossRef] [Green Version]
  17. Bergsma, E.W.; Almar, R.; Rolland, A.; Binet, R.; Brodie, K.L.; Bak, A.S. Coastal morpholgy from space: A showcase of monitoring the topography-bathymetry continuum. Remote Sens. Environ. 2021, 261, 112469. [Google Scholar] [CrossRef]
  18. Turner, I.L.; Harley, M.D.; Almar, R.; Bergsma, E.W. Satellite optical imagery in Coastal Engineering. Coast. Eng. 2021, 167, 103919. [Google Scholar] [CrossRef]
  19. Bergsma, E.W.; Almar, R. Coastal coverage of ESA’ Sentinel 2 mission. Adv. Space Res. 2020, 65, 2636–2644. [Google Scholar] [CrossRef]
  20. Brodie, K.L.; Palmsten, M.L.; Hesser, T.J.; Dickhudt, P.J.; Raubenheimer, B.; Ladner, H.; Elgar, S. Evaluation of video-based linear depth inversion performance and applications using altimeters and hydrographic surveys in a wide range of environmental conditions. Coast. Eng. 2018, 136, 147–160. [Google Scholar] [CrossRef] [Green Version]
  21. Gawehn, M.; van Dongeren, A.; de Vries, S.; Swinkels, C.; Hoekstra, R.; Aarninkhof, S.G.J.; Friedman, J.I. The application of a radar-based depth inversion method to monitor near-shore nourishments on an open sandy coast and an ebb-tidal delta. Coast. Eng. 2020, 159, 103716. [Google Scholar] [CrossRef]
  22. Bouvier, C.; Balouin, Y.; Castelle, B. Video monitoring of sandbar-shoreline response to an offshore submerged structure at a microtidal beach. Geomorphology 2017, 295, 297–305. [Google Scholar] [CrossRef]
  23. Lippmann, T.C.; Holman, R.A. Quantification of sand bar morphology: A video technique based on wave dissipation. J. Geophys. Res. 1989, 94, 995–1011. [Google Scholar] [CrossRef]
  24. van Enckevort, I.M.J.; Ruessink, B.G. Video observations of nearshore bar behaviour. Part 1: Alongshore uniform variability. Cont. Shelf Res. 2003, 23, 501–512. [Google Scholar] [CrossRef]
  25. Van Enckevort, I.M.J.; Ruessink, B.G. Video observations of nearshore bar behaviour. Part 2: Alongshore uniform variability. Cont. Shelf Res. 2003, 23, 513–532. [Google Scholar] [CrossRef]
  26. Plant, N.G.; Holman, R.A. Intertidal beach profile estimation using video images. Mar. Geol. 1997, 140, 1–24. [Google Scholar] [CrossRef]
  27. Stockdon, H.F.; Holman, R.A. Estimation of wave phase speed and nearshore bathymetry from video imagery. J. Geophys. Res. 2000, 105, 22015–22033. [Google Scholar] [CrossRef]
  28. Piotrowski, C.C.; Dugan, J.P. Accuracy of bathymetry and current retrievals from airborne optical time-series imaging of shoaling waves. IEEE Trans. Geosci. Remote Sens. 2002, 40, 2602–2612. [Google Scholar] [CrossRef]
  29. Dugan, J.P.; Piotrowski, C.C.; Williams, J.Z. Water depth and surface current retrievals from airborne optical measurements of surface gravity wave dispersion. J. Geophys. Res. 2001, 106, 16903–16915. [Google Scholar] [CrossRef] [Green Version]
  30. Trizna, D.B. Errors in bathymetric retrievals using linear dispersion in 3-D FFT analysis of marine radar ocean wave imagery. IEEE Trans. Geosci. Remote Sens. 2001, 39, 2465–2469. [Google Scholar] [CrossRef]
  31. Holman, R.A.; Brodie, K.L.; Spore, N.J. Surf zone characterization using a small quadcopter: Technical issues and procedures. IEEE Trans. Geosci. Remote Sens. 2017, 55, 2017–2027. [Google Scholar] [CrossRef]
  32. Brodie, K.L.; Bruder, B.L.; Slocum, R.K.; Spore, N.J. Simultaneous Mapping of Coastal Topography and Bathymetry From a Lightweight Multicamera UAS. IEEE Trans. Geosci. Remote Sens. 2019, 57, 6844–6864. [Google Scholar] [CrossRef]
  33. Hashimoto, K.; Shimozono, T.; Matsuba, Y.; Okabe, T. Unmanned Aerial Vehicle depth inversion to monitor river-mouth bar dynamics. Remote Sens. 2021, 13, 412. [Google Scholar] [CrossRef]
  34. Van Dongeren, A.; Plant, N.; Cohen, A.; Roelvink, D.; Haller, M.; Catalan, P. Beach Wizard: Nearshore bathymetry estimation through assimilation of model computations and remote observations. Coast. Eng. 2008, 55, 1016–1027. [Google Scholar] [CrossRef]
  35. Aarninkhof, S.G.J.; Turner, I.L.; Dronkers, T.D.T.; Caljouw, M.; Nipius, L. A video technique for mapping intertidal beach bathymetry. Coast. Eng. 2003, 49, 275–289. [Google Scholar] [CrossRef]
  36. Soto, F.; Catalan, P. Bathymetry inversion in the surf zone via assimilation of remotely-sensed wave breaking energy dissipation. Coast. Eng. Proc. 2020, 26. [Google Scholar] [CrossRef]
  37. Almar, R.; Bonneton, P.; Senechal, N.; Roelvink, D. Wave celerity from video imaging. In Proceedings of the 31st Conference on Coastal Engineering, Hamburg, Germany, 31 August–5 September 2008; pp. 661–673. [Google Scholar]
  38. Plant, N.G.; Holland, K.T.; Haller, M. Ocean wavenumber estimation from wave-resolving time series imagery. IEEE Trans. Geosci. Remote Sens. 2008, 46, 2644–2658. [Google Scholar] [CrossRef]
  39. Simarro, G.; Calvete, D.; Luque, P.; Orfila, A.; Ribas, F. UBathy: A new approach for bathymetry inversion from video imagery. Remote Sens. 2019, 11, 2722. [Google Scholar] [CrossRef] [Green Version]
  40. Bergsma, E.W.; Almar, R. Video-based depth inversion techniques, a method comparison with synthetic cases. Coast. Eng. 2018, 138, 199–209. [Google Scholar] [CrossRef]
  41. Rutten, J.; de Jong, D.; Ruessink, B.G. Accuracy of nearshore bathymetry inverted from X-band radar and optical video data. IEEE Trans. Geosci. Remote Sens. 2017, 55, 1106–1116. [Google Scholar] [CrossRef]
  42. Honegger, D.; Haller, M.; Holman, R. High-resolution bathymetry estimates via X-band marine radar: 1. beaches. Coast. Eng. 2019, 149, 39–48. [Google Scholar] [CrossRef]
  43. Bak, A.S.; Brodie, K.L.; Hesser, T.J.; Smith, J.M. Applying dynamically updated nearshore bathymetry estimates to operational nearshore wave modeling. Coast. Eng. 2019, 145, 53–64. [Google Scholar] [CrossRef]
  44. Holman, R.A.; Stanley, J. cBathy Bathymetry Estimation in the Mixed Wave-Current Domain of a Tidal Estuary. J. Coast. Res. 2013, 65, 1391–1396. [Google Scholar] [CrossRef]
  45. Holland, K.T.; Holman, R.A.; Lippmann, T.C.; Stanley, J.; Plant, N. Practical use of video imagery in nearshore oceanographic field studies. IEEE J. Ocean Eng. 1997, 22, 81–92. [Google Scholar] [CrossRef]
  46. Holman, R.A.; Stanley, J. The history and technical capabilities of Argus. Coast. Eng. 2007, 54, 477–491. [Google Scholar] [CrossRef]
  47. Bergsma, E.W.; Conley, D.C.; Davidson, M.; O’Hare, T.J. Video-based nearshore bathymetry estimation in macro-tidal environments. Mar. Geol. 2016, 374, 31–41. [Google Scholar] [CrossRef] [Green Version]
  48. Walker, R.E. Marine Light Field Statistics; John Wiley and Sons, Inc.: New York, NY, USA, 1994; p. 675. [Google Scholar]
Figure 1. cBathy pixel array used for the analyses in this paper, overlain on a merged snapshot from 17 May 2015. For clarity, only ¼ of the pixels are shown (decimated by two in both x and y directions). The wavenumber for each analysis point, for example, xm = 250; ym = 750 shown above by the red asterisk, is found using phase map data from a surrounding region, shown by green dots above. Imagery is derived from six oblique-viewing cameras and merged into rectified images such as this [44,45].
Figure 1. cBathy pixel array used for the analyses in this paper, overlain on a merged snapshot from 17 May 2015. For clarity, only ¼ of the pixels are shown (decimated by two in both x and y directions). The wavenumber for each analysis point, for example, xm = 250; ym = 750 shown above by the red asterisk, is found using phase map data from a surrounding region, shown by green dots above. Imagery is derived from six oblique-viewing cameras and merged into rectified images such as this [44,45].
Remotesensing 13 03996 g001
Figure 2. Phase map for the example tile shown in Figure 1 for the dominant frequency, f = 0.0956 Hz. Observed phase is on the left, best-fit modeled (Equation (2)) phase on the right. The x and y components of wavenumber are derived from the components of slope of the phase ramps (colors from seaward to landward going from blue through green, yellow then red with 360° phase jumps as blue–red transitions).
Figure 2. Phase map for the example tile shown in Figure 1 for the dominant frequency, f = 0.0956 Hz. Observed phase is on the left, best-fit modeled (Equation (2)) phase on the right. The x and y components of wavenumber are derived from the components of slope of the phase ramps (colors from seaward to landward going from blue through green, yellow then red with 360° phase jumps as blue–red transitions).
Remotesensing 13 03996 g002
Figure 3. Sensitivity of the dispersion relationship showing the increasing sensitivity of the inversion process as deep water is approached.
Figure 3. Sensitivity of the dispersion relationship showing the increasing sensitivity of the inversion process as deep water is approached.
Remotesensing 13 03996 g003
Figure 4. Radon transform results for the example tile and model point shown in Figure 1 and Figure 2. Left panel shows the actual Radon transform versus candidate wave angle (x-axis) and projected distance (y-axis). The right panel shows the variance of the Radon transform as a function of wave angle. The angle with maximum variance is chosen as the wave angle seed, in this case 10.8°.
Figure 4. Radon transform results for the example tile and model point shown in Figure 1 and Figure 2. Left panel shows the actual Radon transform versus candidate wave angle (x-axis) and projected distance (y-axis). The right panel shows the variance of the Radon transform as a function of wave angle. The angle with maximum variance is chosen as the wave angle seed, in this case 10.8°.
Remotesensing 13 03996 g004
Figure 5. Organization of cBathy version 1.2 (and 1.0) showing the three phases of the algorithm by color.
Figure 5. Organization of cBathy version 1.2 (and 1.0) showing the three phases of the algorithm by color.
Remotesensing 13 03996 g005
Figure 6. New organization of cBathy version 2.0 with the partitioning into tiles happening later. “PrepImages” is added to create timex and other image types from the time stack and is not integral to the cBathy calculations.
Figure 6. New organization of cBathy version 2.0 with the partitioning into tiles happening later. “PrepImages” is added to create timex and other image types from the time stack and is not integral to the cBathy calculations.
Remotesensing 13 03996 g006
Figure 7. Time series of wave height, Hs, wave period, Tm, peak wave direction, and tide elevation for each of the 624 cBathy estimates collected over almost four years. Blue dots for the upper three panels correspond to estimates from the 26 m depth waverider buoy. These are overplotted with red dots from the 8 m array. Shore normal wave direction is 72°, indicated by the blue horizontal line in the wave direction panel.
Figure 7. Time series of wave height, Hs, wave period, Tm, peak wave direction, and tide elevation for each of the 624 cBathy estimates collected over almost four years. Blue dots for the upper three panels correspond to estimates from the 26 m depth waverider buoy. These are overplotted with red dots from the 8 m array. Shore normal wave direction is 72°, indicated by the blue horizontal line in the wave direction panel.
Remotesensing 13 03996 g007
Figure 8. Comparison of bathymetry from the three versions of cBathy, V1.0, V1.2, and V2.0 (left three panels) for 19 November 2017, 1959 GMT, with the ground-truth survey from two days following (right panel). All cBathy results have been filtered to remove any values with an estimated error greater than 0.5 m (dark red color).
Figure 8. Comparison of bathymetry from the three versions of cBathy, V1.0, V1.2, and V2.0 (left three panels) for 19 November 2017, 1959 GMT, with the ground-truth survey from two days following (right panel). All cBathy results have been filtered to remove any values with an estimated error greater than 0.5 m (dark red color).
Remotesensing 13 03996 g008
Figure 9. cBathy performance statistics for 19 November 2017, partitioned by depth. Values are plotted at the mean depth for each bin. No survey depths greater than 6 m were observed for this run.
Figure 9. cBathy performance statistics for 19 November 2017, partitioned by depth. Values are plotted at the mean depth for each bin. No survey depths greater than 6 m were observed for this run.
Remotesensing 13 03996 g009
Figure 10. Map of fBar, the weighted mean frequency that contributed to phase 2 depth estimates in version 2.0.
Figure 10. Map of fBar, the weighted mean frequency that contributed to phase 2 depth estimates in version 2.0.
Remotesensing 13 03996 g010
Figure 11. Comparison of seed wave angle and wavenumber to the final values from the nonlinear fitting routine. Example shows only values for the dominant frequency.
Figure 11. Comparison of seed wave angle and wavenumber to the final values from the nonlinear fitting routine. Example shows only values for the dominant frequency.
Remotesensing 13 03996 g011
Figure 12. Histograms of the bulk coverage statistics of each of the three cBathy versions.
Figure 12. Histograms of the bulk coverage statistics of each of the three cBathy versions.
Remotesensing 13 03996 g012
Figure 13. Percentage of successful coverage for sub-aqueous pixels versus significant wave height at the 8 m array for each of the three versions. Versions were plotted in order from versions 1.0 through 2.0 so the results from version 2.0 (blue dots) overplot earlier results.
Figure 13. Percentage of successful coverage for sub-aqueous pixels versus significant wave height at the 8 m array for each of the three versions. Versions were plotted in order from versions 1.0 through 2.0 so the results from version 2.0 (blue dots) overplot earlier results.
Remotesensing 13 03996 g013
Figure 14. Histograms of the four basic performance statistics for each of the three algorithm versions (colors).
Figure 14. Histograms of the four basic performance statistics for each of the three algorithm versions (colors).
Remotesensing 13 03996 g014
Figure 15. Histograms of the same statistics as Figure 14, but now based on common regions of sampling determined by the regions of successful estimates from the V2.0 version of cBathy.
Figure 15. Histograms of the same statistics as Figure 14, but now based on common regions of sampling determined by the regions of successful estimates from the V2.0 version of cBathy.
Remotesensing 13 03996 g015
Figure 16. Scatter plots of the four performance statistics for each of the three algorithm versions (see legend at top right) for the full dataset of 624 runs. Points are plotted in order, so blue points (V2.0) overplot earlier versions.
Figure 16. Scatter plots of the four performance statistics for each of the three algorithm versions (see legend at top right) for the full dataset of 624 runs. Points are plotted in order, so blue points (V2.0) overplot earlier versions.
Remotesensing 13 03996 g016
Figure 17. Comparison scatter plots of the four performance statistics for each of the three algorithm versions (see legend at top right) for the reduced dataset of 563 runs. Points are plotted in order, so blue points (V2.0) overplot earlier versions.
Figure 17. Comparison scatter plots of the four performance statistics for each of the three algorithm versions (see legend at top right) for the reduced dataset of 563 runs. Points are plotted in order, so blue points (V2.0) overplot earlier versions.
Remotesensing 13 03996 g017
Figure 18. Error bulk statistics binned by depth. Bin depths are plotted at the mean bin depth. Algorithm versions are indicated in the legend. Coverage is omitted, since all versions have been forced to use the V2.0 coverage.
Figure 18. Error bulk statistics binned by depth. Bin depths are plotted at the mean bin depth. Algorithm versions are indicated in the legend. Coverage is omitted, since all versions have been forced to use the V2.0 coverage.
Remotesensing 13 03996 g018
Figure 19. Histograms of the bias (left) and rmse (right) statistics for each of the three version of the cBathy algorithm (see legend). The data have been partitioned into non-breaking regions (solid lines) and surf zone regions (dashed lines; legend has “SZ” appended).
Figure 19. Histograms of the bias (left) and rmse (right) statistics for each of the three version of the cBathy algorithm (see legend). The data have been partitioned into non-breaking regions (solid lines) and surf zone regions (dashed lines; legend has “SZ” appended).
Remotesensing 13 03996 g019
Figure 20. Comparison of Kalman-filtered bathymetry from the three version of cBathy, V1.0, V1.2, and V2.0 (left three panels) for 21 November 2017, with the ground-truth survey from that day (right panel). This is the same survey used in Figure 8. All cBathy results have been filtered to remove any values with an estimated error greater than 0.5 m.
Figure 20. Comparison of Kalman-filtered bathymetry from the three version of cBathy, V1.0, V1.2, and V2.0 (left three panels) for 21 November 2017, with the ground-truth survey from that day (right panel). This is the same survey used in Figure 8. All cBathy results have been filtered to remove any values with an estimated error greater than 0.5 m.
Remotesensing 13 03996 g020
Figure 21. The four basic performance statistics plotted versus surveyed depth for the example survey date of 19 November 2017. Note that each plot includes three lines, one for each algorithm, with the third (blue) often overplotting the others.
Figure 21. The four basic performance statistics plotted versus surveyed depth for the example survey date of 19 November 2017. Note that each plot includes three lines, one for each algorithm, with the third (blue) often overplotting the others.
Remotesensing 13 03996 g021
Figure 22. The four performance statistics describing the results from the 39 Kalman results. Version colors are shown in the legend.
Figure 22. The four performance statistics describing the results from the 39 Kalman results. Version colors are shown in the legend.
Remotesensing 13 03996 g022
Figure 23. Histograms of the mean ratio of observed to estimated version 2.0 errors, averaged over each of the 39 surveys. Red and blue lines indicate Kalman and non-Kalman results, respectively.
Figure 23. Histograms of the mean ratio of observed to estimated version 2.0 errors, averaged over each of the 39 surveys. Red and blue lines indicate Kalman and non-Kalman results, respectively.
Remotesensing 13 03996 g023
Table 1. Means of the four basic performance statistics (rows) comparing across the three versions of the cBathy algorithm (columns). Statistics are computed for the full dataset of 624 runs (upper section) as well as for the reduced dataset of 563 runs for which successful coverage was greater than 50% (based on V2.0 coverage) and significant wave height was less than 1.2 m (lower section).
Table 1. Means of the four basic performance statistics (rows) comparing across the three versions of the cBathy algorithm (columns). Statistics are computed for the full dataset of 624 runs (upper section) as well as for the reduced dataset of 563 runs for which successful coverage was greater than 50% (based on V2.0 coverage) and significant wave height was less than 1.2 m (lower section).
Full Dataset
StatisticV1.0V1.2V2.0
Bias (m)0.190.140.16
rmse (m)0.640.820.56
Δh95 (m)1.271.251.19
Coverage (%)78.778.084.7
Reduced Dataset
Bias (m)0.150.090.09
rmse (m)0.580.660.41
Δh95 (m)1.161.010.85
Coverage (%)82.882.390.0
Table 2. Means for Kalman-filtered bathymetries of the four basic performance statistics (rows) comparing across the three versions of the cBathy algorithm (columns). Statistics are computed for the full dataset of 39 surveys, comparing the mid-afternoon Kalman bathymetry for each survey date to ground truth.
Table 2. Means for Kalman-filtered bathymetries of the four basic performance statistics (rows) comparing across the three versions of the cBathy algorithm (columns). Statistics are computed for the full dataset of 39 surveys, comparing the mid-afternoon Kalman bathymetry for each survey date to ground truth.
Full Dataset
StatisticV1.0V1.2V2.0
Bias (m)0.150.080.08
rmse (m)0.470.660.38
Δh95 (m)0.960.860.78
Coverage (%)99.199.299.9
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Holman, R.; Bergsma, E.W.J. Updates to and Performance of the cBathy Algorithm for Estimating Nearshore Bathymetry from Remote Sensing Imagery. Remote Sens. 2021, 13, 3996. https://0-doi-org.brum.beds.ac.uk/10.3390/rs13193996

AMA Style

Holman R, Bergsma EWJ. Updates to and Performance of the cBathy Algorithm for Estimating Nearshore Bathymetry from Remote Sensing Imagery. Remote Sensing. 2021; 13(19):3996. https://0-doi-org.brum.beds.ac.uk/10.3390/rs13193996

Chicago/Turabian Style

Holman, Rob, and Erwin W. J. Bergsma. 2021. "Updates to and Performance of the cBathy Algorithm for Estimating Nearshore Bathymetry from Remote Sensing Imagery" Remote Sensing 13, no. 19: 3996. https://0-doi-org.brum.beds.ac.uk/10.3390/rs13193996

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop