UBathy: A New Approach for Bathymetric Inversion from Video Imagery

: A new approach to infer the bathymetry from coastal video monitoring systems is presented. The methodology uses principal component analysis of the Hilbert transform of video images to obtain the components of the wave propagation ﬁeld and their corresponding frequency and wavenumber. Incident and reﬂected constituents and subharmonics components are also found. Local water depth is then successfully estimated through wave dispersion relationship. The method is ﬁrst applied to monochromatic and polychromatic synthetic wave trains propagated using linear wave theory over an alongshore uniform bathymetry in order to analyze the inﬂuence of different parameters on the results. To assess the ability of the approach to infer the bathymetry under more realistic conditions and to explore the inﬂuence of other parameters, nonlinear wave propagation is also performed using a fully nonlinear Boussinesq-type model over a complex bathymetry. In the synthetic cases, the relative root mean square error obtained in bathymetry recovery (for water depths 0.75 m (cid:54) h (cid:54) 8.0 m) ranges from ∼ 1% to ∼ 3% for inﬁnitesimal-amplitude wave cases (monochromatic or polychromatic) to ∼ 15% in the most complex case (nonlinear polychromatic waves). Finally, the new methodology is satisfactorily validated through a real ﬁeld site video.


Introduction
Decision making in coastal zone management requires a knowledge of the bathymetry [1,2]. Obtaining accurate bathymetries has interest on its own, since it allows, e.g., to know how the waves propagate to the shore and how the morphology has evolved in time (if several bathymetries are available) or to decide if it is necessary to dredge the mouth of a harbour. Measuring bathymetric time-series also allows to validate morphodynamic numerical models, especially if obtained with a relatively high frequency. The morphodynamic models, in turn, can be a helpful tool to predict future changes and to analyze the impact of potential human actions [3,4]. There is, consequently, a large interest in obtaining accurate bathymetries [5,6].
In situ bathymetric measurement techniques include a wide variety of approaches that range from swath-sounding sonar systems [7] to bottom-contacting vehicles such as the Coastal Research Amphibious Buggy (CRAB) [5]. While in situ techniques provide excellent bathymetries at high spatial resolution, they are both expensive and highly time consuming. Except during experimental field campaigns, they are usually obtained at most a few times per year. Alternatives to in situ methods include, among other, LiDAR techniques [8,9], satellite images [6,10], video images [11,12], and X-band radar images [13,14]. The reader is referred to the work in Reference [15] for a review. One technique will be preferable relative to the others depending on aspects such as the size of the area to be analyzed or the desired spatial and temporal resolution. For instance, satellite information is unbeatable for large domains but it has a low spatial resolution and a limited temporal resolution. On the other hand, video monitoring systems, the object of this work, provide excellent spatiotemporal resolution for domains up to a few kilometers.
Video monitoring systems [16][17][18], mainly developed after the relatively recent advent of digital cameras, have been shown to be a powerful and low-cost tool to monitor the coast. These systems are useful in a wide range of studies, such as, shoreline detection and coastal variability [2,19,20], intertidal bathymetry [21,22], or the study of the morphodynamics of beach systems [23][24][25][26]. Video systems have also been shown to be able to give estimates of the bathymetry through the wave propagation linear dispersion relationship, which relates the water depth to the wavenumber and the wave frequency. For this purpose, the wave frequency and (space-varying) wavenumber are to be obtained from a sequence of snap images. Reference [11] obtained the wavenumber from the wave speed (and the wave frequency), which, in turn, were obtained from timestack images of several transects normal and parallel to the coast. From then on, there has been a considerable effort to develop reliable techniques to obtain the bathymetry from video images, working both in 1D space, i.e., in transects, or by treating the whole images, i.e., 2D [12,27]. Some of the proposed approaches use video images combined with numerical models and/or radar images. Of mention, there have also been attempts to obtain the bathymetry without the dispersion equation by using the wave dissipation pattern observed in time-averaged (timex) images instead of a sequence of snaps [28].
The cBathy algorithm [12] is nowadays the most popular algorithm to obtain 2D bathymetries from video stations [29,30]. The code is made up of two main parts: first, a bathymetry estimate is obtained for each hourly video (typically of few minutes at ∼0.5 Hz); second, given a bunch of hourly (estimated) bathymetries, they are smoothed through a Kalman filter [31] to obtain the final hourly estimates. In regard to the first part (i.e., the estimate of the bathymetry from one video, which is the topic this work is focused on), cBathy first transforms each pixel intensity time history to the frequency domain. To obtain the water depth at any given point, it then considers a neighborhood of the point and obtains the wave spatial pattern through frequency-domain empirical orthogonal functions of the Cross Spectral Matrices (CSMs) averaged within frequency bands of interest. The latter is done to handle noisy data such as in the video images. The analysis of the set of CSMs allows to get the dominant frequencies (those of which the signal is more coherent) and their corresponding wavenumbers. The dispersion relationship then allows to obtain estimates of the local water depth for each dominant frequency and a weighted average of them. Some limitations and/or known problems of cBathy have been reported in the literature [30,32] (e.g., dealing with high wave heights, wet/dry tiles, or long waves).
In this work, we propose an alternative to cBathy to obtain nearshore bathymetries from videos, which is based on Principal Component Analysis (PCA) of the Hilbert transform of the video images. This methodology consists of retrieving wave patterns from a time-space complex-PCA analysis of video images, i.e., avoiding the frequency-domain analysis, to subsequently obtain nearshore bathymetries. The use of PCA in the wave propagation problem was actually introduced prior to the popularization of video monitoring stations [33]; in this work, PCA was applied to data recorded through tide gauges and the goal was detecting infragravity waves. The problems observed in that work were related to the fact that the distance between the gauges was relatively large compared to the wavelength. This is actually not expected to be a problem for video images detecting wind waves, for the pixel size projected in the space is much smaller than the wavelength in the area of interest. The presented algorithm, uBathy, takes its name from cBathy but with "u" referring to ULISES [34], the software it has been developed in.
The aim of the present contribution is to demonstrate the validity of this new methodology and to find the algorithms and parameter values that optimize the result. In order to fully control the conditions, most tests are performed with synthetic bathymetries and waves. To test the performance of uBathy in real conditions, the bathymetry inferred from a video captured in a field site is also included. The proposed methodology is first presented and illustrated with a simple 1D case (Section 2). The validity of the proposed approach is analyzed through synthetic wave fields, both linear and nonlinear, propagating over two different bathymetries (Section 3). Special emphasis is devoted to the recovery of the wave frequency and the wavenumber from the PCA analysis of the video. The influence of the most relevant parameters is further discussed (Section 4), such as the temporal and spatial discretization of the video or the parameters that influence the post-process of the PCA to recover the bathymetry. Additionally, a video from a field site is analyzed with the presented methodology to estimate a real bathymetry. The result is compared with the one obtained from cBathy. To conclude, the most important findings are listed (Section 5).

Methodology
The present algorithm relies on PCA, which is briefly described in Section 2.1 for completeness. From this algebraic decomposition, the frequency and wavenumber of the different wave components are determined (Sections 2.2 and 2.3), and henceforth, the bathymetry is obtained (Section 2.4). All the steps are illustrated with a simple 1D case (for ease of representation).

Principal Component Analysis
Consider any spatiotemporal real-valued function, f (x, t) , discretized in space and time into a matrix X so that X mn = f (x m , t n ) , with m = 1, . . . , M and n = 1, . . . , N. When working with real video images, x m will be the real-world coordinates corresponding to the image pixels, t n will be the time of each snapshot, and f will typically be the value of the pixel intensity for the snapshot. Each column of X corresponds to a stacked "snapshot", while each row is the time record at a given point.
Let Y be the demeaned matrix, built up by subtracting from each column of X the time average space-vector: Considering M time-depending variables, one for each row of the matrix Y, their covariance matrix is Y · Y T /N. Singular value decomposition, the key in PCA, allows to rewrite Y ∈ M M×N (i.e., an M × N matrix) as with U ∈ M M×Q , S ∈ M Q×Q , V ∈ M N×Q , and Q = min {M, N} . Further, U T · U = V T · V = I Q×Q , the identity, and S are a diagonal matrix of real positive values (by convention, in decreasing order in the diagonal). Following usual notation, the qth column of U, a spatial vector, is the "Empirical Orthogonal Function" EOF q , while the qth row of the product S · V T , a time vector, is the "Principal Component" PC q . From Equations (1) and (2), The EOFs (columns of U) are a set of orthonormal vectors in space, while the PCs (rows of S·V T ) are a set of orthogonal vectors in time. Each pair {EOF q , PC q } is a mode of the decomposition.
The PCs can be interpreted as a rewritten version of the information in Y using the EOFs as a basis in Equation (2). The covariance matrix of the data expressed in this new basis is i.e., a diagonal matrix. The contribution of each PC q and the corresponding EOF q to the total signal is quantified by the values of the diagonal matrix S · S T , which indeed represents the explained variance, σ 2 q = (S · S T ) qq . For complex-valued signals, the above results hold as long as transposed matrices are substituted by conjugate-transposed matrices. In that case, both EOF q and PC q become complex-valued vectors but σ 2 q remains real valued. Propagation of small-amplitude waves over a flat bed can be described as [35] where J is number of components of the wave field, a j is the amplitudes, k j = (k xj , k yj ) is the wave vectors (with corresponding wavenumbers k j = | k j |), ω j = 2π/T j is the angular frequencies (with T j being the wave periods), and ϕ j is the wave phase lags. The time-wise Hilbert transform is and the corresponding space and time discretized matrix reads For large time domains, it is f m ∼ 0 and Y mn ∼ X mn . The J time-wise complex vectors, i.e., exp (i ω j t n ) , tend to be orthonormal as the space-wise vectors are, i.e., exp (−i (k j · x m + ϕ j ) ) , and the Equation (5) is already the PCA decomposition of the signal into its monochromatic components (with ω j and k j ). In this case, the explained variance for the qth mode reads Hence, each component of the wave field can be linked to a mode of the PCA. Moreover, the angle of a given PC and EOF can be used to obtain the ω j and k j , respectively, of the corresponding wave component since Above, we have considered k j · x, while in wave propagation over uneven beds, k j is not constant and the spatial phase φ j (x) = ∫ k j · dx has to be used instead. As shown below, this will not actually be a limitation.
The usefulness of PCA to analyze wave propagation over uneven beds (i.e., with variable k) is illustrated with the linear 1D propagation of waves. The bathymetry considered is defined by a water depth, h, equal to To represent more realistic sea state situations at the coast, where waves with different frequencies and wavenumbers can coincide, the superposition of two wave trains has been considered. The free surface elevation η (waves travelling rightwards) is The amplitudes, a j (x) , are obtained from the amplitudes at x = 0, a 0 j , using fundamentals of linear wave propagation theory [35]. The wavenumbers k j are related to the angular frequencies and the water depth through the dispersion relationship Three synthetic linear wave propagation cases are considered (Table 1), including a case with only one wave train (monochromatic), a case with two wave trains of different frequencies (bichromatic), and a case with two wave trains of the same frequency but opposite propagation directions (reflective). This allows to describe different features of the PCA. The discretization considers x m = m∆x and t n = n∆t, with ∆x = 1 m, m = 1, . . . , M = 200 (so that x max = 200 m), ∆t = 0.25 s, and n = 1, . . . , N = 400 (so that t max = 100 s). Table 1. Wave conditions in the seaward boundary for the 1D examples: For each wave train (two at most), T j is the period, a 0 j is the wave amplitude at x = 0, and dir j is the direction of wave propagation (+, rightwards).
For the monochromatic case, the analysis of the signal yields one main mode that represents 99.5% of the variance ( Table 2). For the bichromatic case, two main modes are obtained (corresponding to the two components), with variances σ 2 1 = 87.6 % and σ 2 2 = 12.1 % (i.e., that account for the 99.7% of the total variance). Finally, when two waves with the same period and moving in opposite directions are superimposed (reflective case), only one mode that represents 99.7% of the variance is obtained out of the PCA. The angles of the PCs and EOFs of the three cases are represented in the top panels of Figures 1  and 2, respectively. In the monochromatic case, the slope of the angle α t of the PC, leaving aside the jumps from π to −π, is approximately constant. The slope of the angle α x of the EOF is gentler at x = 0 than at x = 200 m and acknowledges the dependence of the wavenumber on the water depth. In the bichromatic case, each mode corresponds to one of the superimposed waves: the slope of α t is gentler for the second mode than for the first one, consistently with the fact that ω 2 = 0.757 rad/s < 1.232 rad/s = ω 1 . In the reflective case, where there is only one mode, α t is similar to the monochromatic case but, now, α x has a wavy behavior, resulting from the superposition of two waves with the same frequency.

Wave Frequency and Wavenumber: Phase Fitting
The wave angular frequency (just frequency hereafter) at a certain time t 0 is to be determined by fitting the angle of the PC in a vicinity of radius R t around t 0 , making use of Equation (7a). However, to prevent the discontinuities in α t , the latter is first centered at t 0 , where complex conjugated values are denoted with an overline. Figure 3 illustrates howα t avoids the discontinues for t 0 = 10 s and considers ∆t = 1 s and R t = 3 s. In a second step, witht = t − t 0 ,α is fitted through the expressionα so that the frequency at t 0 is estimated as ω = p 1 . This fitting method will be referred to as phase fitting. An analogous phase fitting procedure can be applied to determine the wavenumber, k, from α x in Equation (7b).  . The observed overshoots in ω (t) at the domain boundaries are related to the discrete Hilbert transform. For each PCA mode corresponding to a travelling wave, the frequency must be constant. This constant value is estimated by averaging ω (t) but by skipping the values at a time distance of ≈ T near the boundaries. This averaged frequency will hereafter be referred to as ω, and it is plotted in red in Figure 1. For the monochromatic and reflective cases, ω∼ω 1 = 1.232 rad/s ( Table 2, expressed as period), and for the bichromatic case, the first mode gives ω∼ω 1 and the second mode gives ω∼ω 2 = 0.757 rad/s. In all cases, the relative errors in ω 1 and ω 2 are less than 0.05 %. The standard deviation, σ ω , of the values of ω (t) used to obtain ω can be regarded as a measure of the quality of the recovered frequency. Above, σ ω /ω < 1 % for the first modes and σ ω /ω∼2 % for the second mode in the bichromatic case.
The second row of Figure 2 shows the values of k obtained in the 1D examples through phase fitting of α x with R x = 2 m (= 2∆x). The behavior of the wavenumber is smooth (with small oscillations) except for the reflective case (right), which is unrealistic but consistent with the corresponding α x . In order to cope with reflected waves, which can occur in some real conditions, an alternative method to obtain k must be developed.

Wavenumber: Function Fitting
If non-negligible reflection occurs, the EOF contains information on both the incident and reflected waves, i.e., with the condition | k a | = | k b | = k. To obtain k and the rest of the parameters involved and following the procedure described before for the frequency (Section 2.2) but now including a normalization, we first define , with x 0 being the point where the wavenumber is being estimated. The expression (10) is then wherex = x − x 0 . The function of Equation (11) is therefore fitted to the normalized EOF by optimization in the neighborhood of x = x 0 . The seven optimization parameters areÂ a ,Â b , ϕ a ,φ b , | k a | = | k b | = k, and the two angles (directions) of the wavenumbers. The wavenumbers obtained with this method match the analytical values in all three cases of the 1D example (not shown). This method to get k, that will be referred to as function fitting, is computationally much more expensive than phase fitting (around two orders of magnitude).

Depth Inversion
Once ω (constant) and k (space-varying) have been estimated, the local water depth, h, can be inferred from the dispersion equation [12]. With γ = ω 2 /gk, it reads Figure 2 (third and fourth rows) shows the results of h in the three 1D cases using phase fitting and function fitting approaches for k (ω is always computed with phase fitting). The exact bathymetry of Equation (8) is also included in the figure (blue lines). Most interestingly, the function fitting method allows to properly recover the bathymetry in the reflective case, where the phase fitting method fails. From Figure 2, for the nonreflective cases and for the first mode, the function fitting method avoids the small oscillations observed in the phase fitting method. The function fitting can also produce spurious peaks, which are due to the difficulties in the optimization to obtain the parameters in Equation (10). This errors can be avoided by further increasing the computational time.
To reduce the errors associated to the small oscillations observed for the phase fitting method with a low computational cost and to handle real videos that might include time intervals with large noise, an extension of the above scheme is proposed. First, instead of performing one unique PCA, all the sub-videos obtained from a moving time window of width w t are analyzed and only the dominant PCA mode of each sub-video is considered for the analysis. Second, for each time window, t i , a pair {ω i , σ ωi } is obtained from the first PC and sub-videos for which σ ωi /ω i > 15% are discarded by assuming that the recovered ω i is not good enough. Third, for each time window, k is obtained from the first EOF using the phase fitting method. Since the recovery of k might be unfeasible at some points, those points where the correlation coefficient of the fitted α x are below 0.70 are filtered out as well. Finally, in a fourth step and following Reference [12], the bathymetry at each point is the result of the best fit of the dispersion relationship (12) using all the pairs {ω i , k i } obtained in a neighborhood R x of the point (R x = 0 meaning that only the values at the point are considered). This extension will be referred to as the windowing method. The last row in Figure 2 shows, for the same examples, the results of the windowing method with w t = 40 s and R x = 0. While the reflective case cannot be recovered, as expected, it is seen how the recovered bathymetry fits the exact one in the monochromatic and bichromatic cases.
The oscillations and related errors, observed using the phase fitting method can also be reduced by increasing the total time of the video (for nonreflective cases only). Following the above example, Figure 4 shows the evolution of the root mean square error of the obtained bathymetry, E RMS h , as a function of time t max using phase fitting (without windowing) for the monochromatic case. As a general trend, the error (from the oscillatory pattern in h) reduces as t max increases. The oscillatory behavior in Figure 4 is related to the time domain adjusting (or not) to a multiple of the wave period. For t max = 300 s, the error is below 5.5 cm. In Figure 2, corresponding to t max = 100 s, the errors E RMS h for the monochromatic case are 10.5 cm (phase fitting), 4.3 cm (function fitting), and 2.8 cm (phase fitting + windowing).

Results
The proposed methodology introduces parameters such as the resolution of the spatiotemporal discretization (∆t and ∆x), the radius (R t and R x ) of the neighborhoods to recover the wave frequency and the wavenumber from the PCs and EOFs, the video duration (t max , assuming that the video starts at t = 0), or the parameters defining the windowing (w t and R x ). The influence of these parameters is studied using synthetic linear and nonlinear 2D wave fields. Linear wave propagation equations assume that the wave height is infinitesimally small and that, for alongshore uniform bathymetries, simple analytical solutions can be computed. The synthetic linear wave fields will be used to analyze the influences of ∆t, ∆x, R t , and R x . These linear solutions, however, dismiss wave reflections at the shore, and are unable to represent other important features of real wave fields such as wave decomposition (energy transfer from one frequency to other) or wave-wave interactions. For this reason, realistic nonlinear numerical models are also used to examine other phenomena: propagation over complex bathymetries, generation of subharmonics, and wave reflection. The synthetic nonlinear 2D wave fields will also allow to investigate the influence of t max and w t (time windowing).

Analysis for Synthetic Linear Waves: Phase Fitting
Three monochromatic wave trains (W1, W2, and W3; Table 3) and their superposition (WS = W1 + W2 + W3) are propagated over an alongshore uniform bathymetry ( Figure 5A). This bathymetry corresponds to a beach profile with a shore-parallel bar located 80 m from the shore. Wave conditions are meant to provide different wave periods in a realistic range and different propagation directions to check whether the obtained bathymetry is mainly independent of the wave characteristics. Notice that, since linear wave theory assumes infinitesimal amplitudes, in this section, the amplitudes of Table 3 are only to show the influence of the relative strength of the different wave trains in the WS case. Table 3. Wave conditions in the seaward boundary for the analysis of synthetic 2D cases: For each wave train, T is the period, A is the wave amplitude in deep waters, θ is the angle with respect to the shore normal in deep waters, and ϕ is a phase lag.  Knowing the wave conditions in the seaward boundary (Table 3), the waves are propagated towards the coast using the linear wave theory [35]. The initial snapshots for W1, W2, W3, and WS are shown in Figure 6. In all cases, the spatial domain is 200 m × 300 m (in the alongshore and cross-shore directions) and the time domain is of 90 s. The video snapshots, X mn , are obtained by assigning to each pixel an intensity that is a linear function of the free surface elevation. In this section, only the phase fitting method is applied. Within this section, windowing is not considered so as to focus on the influences of ∆t, ∆x, R t , and R x . The influence of the spatiotemporal discretization (∆t and ∆x = ∆y) of the signal is analyzed considering all combinations of ∆x = ∆y in {1 m, 2 m, 4 m, 10 m} and ∆t in {0.25 s, 0.5 s, 1 s, 2 s, 4 s}. A major result of the PCA for monochromatic waves is to obtain one main mode that explains, in all cases, more than 98% of the total variance of the signal and that matches the corresponding wave train. The results found for the three monochromatic wave trains are qualitatively similar, and the focus hereafter is on W1.

Wave Train T [s]
For illustrative purposes, Figure 7 shows an example of the angles α t and α x , capturing the refraction of the propagation as the wave train travels to the shore (compare α x with W1 in Figure 6). For ∆t = 4 s, the PCA decomposition yields useless results for any ∆x, and hereafter, results with ∆t = 4 s have been disregarded.  Table 4 for ∆x = 2 m. It turns out that results are independent of ∆x. The case R t = 4 s, giving large errors, is not considered in the following. The wavenumber k from the EOFs has been computed for the 2D spatial domain using values of R x in {2 m, 4 m, 8 m, 12 m, 16 m} and the different combinations of ∆x and ∆t. Figure 8 (top panels) shows an example of the recovered k as well as the local relative error, k , obtained from the first (and only) EOF of the wave train W1. As depicted in the figure, k is below ∼ 1%. The global relative RMSE, RMS k , obtained for the full set of exploration results is shown in Table 5 for the domain restricted to depths h 0.75 m (see also Figure 5 where h = 0.75 m is highlighted). As occurred with ∆x in the frequency recovery, the influence of ∆t is minor when computing k (only ∆t = 0.5 s is shown in Table 5). Table 4. Relative errors for ω, ω (in %), as a function of ∆t and R t for ∆x = 2 m, corresponding to linear propagation of W1.    Once ω and k are obtained from the only PCA mode of each monochromatic wave train, the corresponding bathymetry can be derived. Figure 8 (bottom panels) shows the bathymetry and the error obtained from the same case above. The inversion produces small errors, the largest errors being located in the shallower area. As in the 1D example, very small oscillations appear manifestly in k and especially in h. A summary of the results from the PCA for each monochromatic wave train and the corresponding global relative errors in ω, k, and h is presented in Table 6 (upper  half) for W1, W2, and W3. The errors in k and h are significantly larger for W3, which corresponds to the wave field with the smallest period. The PCA of the linear polychromatic wave field (WS in Figure 6) has been performed only for the default values ∆t = 0.5 s, ∆x = 2.0 m, R t = 1.0 s, and R x = 8 m, so as to explore the ability of the method to identify the different wave components and to infer the bathymetry from each one of them. In this case, three main modes are obtained (see Table 6 (lower half)) that accumulate a 99.2 % variance. The periods of each of these modes coincide with those of the corresponding constitutive linear waves. Further, the variance of each of them corresponds to the one predicted by Equation (6) from each of the amplitudes (which are 71.1 %, 25.7 %, and 2.8 %). The angle α x of these three modes (Figure 9, upper panels) reproduce well each respective linear wave train. A summary of the errors in k and h is included in Table 6 (lower half). The errors in k and h increase with the mode number, and again, the largest errors in h occur for the third mode, which corresponds to the smallest period (W3). However, the bathymetry can be successfully retrieved from the first mode, with a relative error of only 3 %.

Analysis for Synthetic Nonlinear Waves: Function Fitting and Windowing
Nonlinear wave propagation over a nonuniform bathymetry has also been analyzed. The goal is to test the capability of the new bathymetry inversion method under more challenging conditions and to gain a better understanding of function fitting and windowing methods. The proposed bathymetry is based on the one for linear waves but with the addition of three sand banks: two in the region of the bar, simulating a crescentic bar, and another one near the shore, simulating a transverse bar ( Figure 5B). Nonlinear waves are modelled with the fully nonlinear Boussinesq-type model FUNWAVE [36] (see Appendix A for details). Unlike in the linear propagation model used above, wave reflection is allowed to happen in FUNWAVE. Only two wave fields are considered in this section: W1 (monochromatic) and WS (polychromatic). The influence of the wave height is analyzed through three multiplying factors for the wave amplitudes in Table 3 (which now represents the real wave amplitudes). The three multiplying factors are F = 1.0 as reference case, F = 0.25 for smaller nonlinearities (i.e., wave amplitudes A 1 = 2.5 cm, A 2 = 1.5 cm, and A 3 = 0.5 cm), and F = 2.5 (A 1 = 25.0 cm, A 2 = 15.0 cm, and A 3 = 5.0 cm). Figure 10 shows snapshots of the wave fields. In the forthcoming analysis, ∆t = 0.532 s, ∆x = 2.0 m, R t = 1.1 s, and R x = 8 m. The computations for the first wave field case, W1, are presented to show how the nonlinear nature of the waves is revealed in the PCA and the corresponding EOFs ( Table 7). The first mode for F = 2.5 corresponds to the wave field W1 (period T ≈ 7.95 s); the second one has half the period and reflects the nonlinear transfer of the dominant wave field W1 into its first harmonic (of twice the frequency of the main one). Further on this, the third EOF has a period which is 1/3 of that of the dominant wave field. As wave amplitudes decrease, contributions to higher harmonics are reduced and, for F = 0.25, this nonlinear transfer is so small that only one EOF is obtained. For the polychromatic wave field, WS, the complete analysis is performed and the results obtained using phase and function fitting are summarized in Table 8. For F = 0.25, the three modes retrieve the three wave fields W1, W2, and W3 (compare the periods in Table 8 with the ones in Table 3). For F = 1.0 and F = 2.5, the modes retrieve only W1 and W2. The first mode provides better results in all cases, and function fitting method improves the results for F = 0.25 and F = 1.0. The upper half of Figure 11 shows, for F = 1.0, how the first two modes correspond to the first two components of the wave field (W1 and W2), while the third one seems to mix characteristics from a harmonic of W1 with W3 (given that T 1 /2 ≈ T 3 ). Figure 11 also includes the bathymetry obtained using phase fitting. Again, as in the linear case, the bathymetries retrieved by modes 2 and 3 have larger errors than that corresponding to mode 1, and therefore, only the first mode will be considered for windowing.  The results obtained by windowing the polychromatic case WS are shown in Table 9. The three multiplying factors F, two video lengths (t max = 90 s, as above, and t max = 150 s), and five different values of w t are considered. The number of sub-videos, that is a function of t max and w t for given ∆t, is included in Table 9 in parentheses. The results have been obtained for R x = R x = 8 m, thus introducing some spatial filtering. The results for R x = 0 have larger errors (∼ 50% higher, not shown). Applying windowing method improves the results of phase fitting method (for certain optimal values of w t ), and the obtained relative RMSE in h can remain below 15 % for the three values of F. Table 9. Results for the first mode for nonlinear wave propagation of the polychromatic WS case with different F factors, total video length t max , and window width w t : Relative RMSE, RMS h , is given for h 0.75 m. Here, "p" and "f " stand for results using phase and function fitting for t max , respectively. The number of sub-videos are included in parentheses. As an example, Figure 12 shows, for the three inversion methods, the relative errors, h , obtained for F = 1.0, t max = 150 s (and using w t = 60 s for windowing). In this case, both function fitting and windowing improve the result of phase fitting in all domains. However, as in the 1D example, function fitting method shows, at some points, peaks due to an optimization failing.

Error Sources
As already seen in the 1D linear example (Figure 2 monochromatic and bichromatic cases), while the errors in k are similar both in deep and shallow waters (they are almost negligible in the first mode), the errors in h are much larger in the deeper region. This is a consequence of how the errors, both in ω and in k, propagate to the inverted water depth. Recalling the dispersion relationship, we can write for the relative errors h , ω , and k , with Here, γ = ω 2 /gk. Figure 13 shows δ k and δ ω as functions of γ. This figure is similar to that by Reference [11] for δ k but, here, includes δ ω . For γ 0.8 (which is equivalent to kh 1.10 or ω 2 h/g 0.88), the propagation errors are large and increase rapidly. The physical reason is that waves do not feel the bottom if water depth is much larger than their wavelength. Figure 13. Propagation of the errors in k and ω to water depth h when using the dispersion relation.
Note that the number γ increases both if the water depth increases or if the wave period diminishes. A critical analysis of the results in the view of the values of γ is crucial. In the 1D linear example shown in Section 2.4, γ goes from 0.53 in the shallow area to 0.93 in the deeper area, with the corresponding observed amplification of the errors in h (Figure 2). In the 2D examples, the bathymetry retrieved from W3 gives the largest errors (e.g., Table 6) because the period is the smallest (T = 5.0 s) and γ 0.8 for h 5.5 m, i.e., inside the studied domain (for W1, T = 7.95 s, so that γ 0.8 for h 13.5 m, outside the studied domain).

Sensitivity to ∆t, ∆x, R t , and R x
The influence of the spatiotemporal discretization of the signal has been analyzed through the linear propagation of W1. One major result of the PCA for monochromatic waves is to obtain one main mode that explains more than 98%. The error of the wave frequency, ω , depends mainly on ∆t and R t (not on ∆x). The radius R t (that has to be ∆t) is required to be smaller than T/2 (Table 4) to avoid the jumps ofα t , i.e., it has to hold ∆t R t < T/2. Once this condition is satisfied, the errors in ω are small in all cases. Similarly, for the recovery of k through phase fitting, ∆t plays a minor role and the condition ∆x R x < L/2 must hold, where L is the wavelength (which will depend on h for a given ω). In this case, the error RMS k reduces as ∆x and R x reduce (Table 5), getting stable for ∆x 4 m and R x 8 m (the wavelength in that case ranged from ∼ 20 m to ∼ 65 m).

Sensitivity to the Inversion Method, T max , and W t
While phase fitting is the only alternative proposed for the recovery of ω from PCs, three different approaches are considered to recover k from the EOFs: phase fitting, function fitting, and windowing. Windowing for k relies on phase fitting; applying windowing on function fitting has been disregarded due to the high computational cost.
Function fitting turns out to be the best choice whenever wave reflection is not negligible or there are two wave fields with the same frequency (e.g., Figure 2, function fitting panels). Windowing is able, in general, to reduce the oscillatory patterns of phase fitting, but it does not work for cases with non-negligible reflection (e.g., Figure 2, windowing panels). Regarding the more realistic nonlinear 2D wave fields, windowing with a convenient w t can reduce the errors compared to both phase fitting (p) and even function fitting ( f ) for F in {1.0, 2.5} (Table 9). For F = 0.25, windowing can improve the phase fitting results (likely removing part of the oscillatory patterns) but not those of function fitting. Actually, wave reflection, though small, is not negligible for F = 0.25. In that case, there is little wave breaking, less energy dissipation, and therefore more wave reflection. Wave reflections are in general small in dissipative beaches but may play a role in reflective beaches (e.g., beaches with large slopes or short waves) or in areas with structures such as harbours. Whenever wave reflection is not negligible, function fitting should be considered. Figure 14 shows how windowing mitigates oscillatory patterns. The figure shows the relative errors of h recovered from three different sub-videos. The oscillatory pattern of the error propagates similar to the own wave field component. Therefore, provided that there are sufficient sub-videos, the errors are compensated to some extent in the averaging process. The more sub-videos there are available, the better the windowing will filter the oscillatory pattern as long as the quality of the ω and h recovered from each sub-video are sufficiently good. This fact can be seen in the Table 9 (especially for F = 1.0). For small values of w t , the error is greater than that obtained for the analysis of the complete video (t max duration); however, as w t increases, the error diminishes. For excessively large values of w t , the error rises again as the number of sub-videos is reduced. The dependency on w t becomes stronger for increasing F (more nonlinear waves).

Field Data
The present methodology has finally been applied to a video from a real field site in which waves are not only nonlinear but also affected by noise. The results have been compared to those obtained from cBathy. The cBathy code and the wave video snaps for test have been obtained from the GitHub distribution [37], managed by the Coastal Imaging Research Network. The cBathy code and images used (October 22, 2011 at 15h in Duck, NC) are part of the study presented by Reference [12]. The video consists of 2048 snaps at 2 Hz. The spatial domain of the video covers a region of 1000 m × 500 m with an irregular mesh of 8576 points (the spatial mesh not being regular does not make a difference for the present approach). The dominant wave period is around T∼15 s, the wavelengths are around 100 m, and wave breaking was present in some regions of the domain (lower left part). The ground truth was obtained through a Coastal Research Amphibious Buggy (CRAB) on 19 October 2011, and the results are known in a regular mesh (Figure 15), so that all the results from the video analysis were interpolated to the CRAB mesh for comparison purposes. The reader is referred to the work by Reference [12] for further details on the data set.
The video of around 16 minutes has been analyzed using the windowing method with w t = 80 s (i.e., around 5 periods). The results without windowing were unsatisfactory (not shown). For each sub-video, only the dominant mode is considered. Also, R t = 1.5 s (for frequency recovery), and given the long wavelengths, R x = R x = 30 m (for wavenumber recovery). Regarding cBathy, the code is applied without modifying any parameter and the Kalman filter is not used (since we consider only one video). The results obtained both from cBathy and the new methodology ("uBathy") are shown in Figure 15, together with the ground truth bathymetry provided by the "CRAB" (the plot is doubled for ease). In regions with observed wave breaking where the dispersion relation is not applicable, the errors increase with both methodologies. As shown in Table 10, uBathy improves the results obtained with cBathy. It not only recovers a higher amount of points (40 % more) but also provides smaller average error (bias) and RMSE. This proves that the new proposed methodology is also valid to handle the noisy wave conditions occurring in real beaches. The computational times to analyse the video with "cBathy" and "uBathy" were of the same order of magnitude.

Future Work
Several known issues require further analysis or remain still open. Such issues will be investigated in the future by applying systematically the present approach to real field site videos. First, the wavelength depends on h and, therefore, the values of R x (and R x ) could be a function of a previous estimate (if available) of the bathymetry. Second, when using real videos, there will be the possibility that some regions of the wave field are particularly noisy at given time intervals (e.g., due to passing of moving objects). This suggests the extension of the (time) windowing scheme proposed above to a space-time windowing scheme. Third, the window width w t can be chosen as a function of the wave period T (namely, w t ∝ T). The exact length of each sub-video is relevant for the quality of its result, according to Figure 4 (in that figure, t max plays the role of w t , since there is just one video). By properly choosing the windows width, the errors could probably be reduced. Fourth, when dealing with a series of hourly videos, following Reference [15], a Kalman filter should be used.
Finally, for real field site wave conditions, it is recommendable to retrieve the bathymetry for adequate conditions: monochromatic waves of small height (ideally), with an adequate wave period for the desired depths to be measured. In case of macro-tidal conditions, the method can be applied both in high tide (to obtain the bathymetry of the shallower area) and in low tide (to obtain the deeper area bathymetry). In dissipative beaches, where wave reflection is minor, the method using phase fitting with windowing should provide better results, whilst in reflective beaches, function fitting method should be applied.

Conclusions
A new methodology to retrieve the bathymetry out of wave propagation recorded by coastal video monitoring systems has been presented. It is based on Principal Component Analysis (PCA) of the Hilbert transform of video images. The method is first tested and validated with synthetic wave fields over known bathymetries. A first set of examples of wave fields are obtained with linear wave theory, which describes the propagation of waves of infinitesimal height over an alongshore uniform bathymetry. To generate more realistic conditions, a fully nonlinear Boussinesq-type model is also applied to propagate finite-amplitude waves over a more complex alongshore variable bathymetry. Finally, a field site video is also used to test the method under real wave conditions. A major result of the present contribution is that PCA successfully provides a decomposition of the videos into a set of modes associated to the different components of the wave propagation field, even when waves have large amplitudes (i.e., large nonlinearities). In the latter case, the PCA also allows to isolate the subharmonic components.
The frequency (ω) of the wave trains are obtained by phase fitting the angles of the PCs, which successfully works in all the studied conditions. Three different approaches have been developed to obtain the wavenumber (k). Performing a phase fitting of the angles of the Empirical Orthogonal Functions (EOFs) can only resolve well the linear wave cases, but it fails under more realistic conditions. Making a function fitting of the angles of the EOFs enables to accurately obtain the wavenumbers in most of the tested conditions and it can even identify incident and reflected constituents, but it has a high computational cost. Applying a time windowing to the phase fitting method greatly improves the results and provides accurate values of the bathymetry for all tested cases, being much more efficient computationally, but it fails when there are reflected waves. The latter is the recommended method to use in dissipative beaches, but in sites with significant wave reflection, the function fitting method is the only valid approach.
Once ω and k are obtained, the local water depth is successfully estimated by inverting the wave dispersion relation (for water depths 0.75 m h 8.0 m), after establishing the optimal values for all the parameters. When the methodology is applied to a field of monochromatic or polychromatic linear waves on the alongshore-uniform bathymetry, relative errors in h do not exceed 3.5 % (using the phase fitting method for k). For the more realistic case of polychromatic nonlinear waves over a complex bathymetry, the relative root mean square errors in h are around 15 % (using the windowing method for k). An application to a real video obtained in a field site confirms the capability of the present methodology to handle realistic wave conditions. Compared to a state-of-the-art bathymetry extraction code, the present new approach recovers a 40 % larger amount of points and the overall error is smaller.

Acknowledgments:
The authors would like to thank R. A. Holman for providing the cBathy code and for the stimulating discussions. The authors also acknowledge the useful suggestions from E. Abella.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

Appendix A. Nonlinear Wave Field Generation
To model nonlinear waves, FUNWAVE-TVD [36] has been used, which simulates wave propagation over a rectangular domain employing a set of fully nonlinear Boussinesq equations. The wave forcing input is modelled by a source term in the equations (the wavemaker), localized in some internal region of the domain. It needs to be placed over a constant depth and must have a thickness comparable to the wavelengths present in the domain (usually > 0.25 wavelengths). The wavemaker used is called WK_DATA2D and adds up monochromatic waves of frequency, orientation, and amplitude specified by the user. The source code has been modified to also accept the phases of the monochromatic waves (originally, it used random ones). The coastal and the offshore boundaries are modelled with sponges followed by reflective walls. Direct sponges have been used, which attenuate the value of variables over the sponged cells. The values used for the parameters verify that the sponge thickness is of the order of the wavelength and that most of wave energy is absorbed. In the laterals, periodic boundary conditions are implemented.
The bathymetry shown in Figure 5 (right panel) is implemented here in the so-called video zone, which has dimensions of 200 m × 300 m (alongshore and cross shore) as in the linear wave case. However, the FUNWAVE domain has been extended to fit the wavemaker and the sponges ( Figure A1). Instead of the original seaward depth limit of 7.3 m, the profile is extended up to 8 m depth and then becomes constant. For the frequencies used, this depth results in a maximum wavelength of about 100 m. Therefore, the thickness of the wavemaker is 25 m and that of the sponges is 100 m. The extension in the offshore boundary is of 250 m to fit the wavemaker, the sponge, and the separation spaces for safety. A coastal extension of 122 m that includes a sponge is also implemented. In order to limit reflection and to avoid breaking the bathymetry is clipped to 0.3 m. However, some wave reflection still occurs, although the reflected wave has a tiny amplitude. In the laterals, a space of 140 m at each side is used to minimize potential influences of the periodic boundary conditions into the domain ( Figure A1). The grid size is of 1 m both in the cross-shore and in the alongshore directions. The first 200 s are the warm-up time and the subsequent 150 s becomes the videos used for bathymetry inversion.