Accurate Depth Inversion Method for Coastal Bathymetry: Introduction of Water Wave High-Order Dispersion Relation

: This paper proposes a wave model for the depth inversion of sea bathymetry based on a new high-order dispersion relation which is more suitable for intermediate water depth. The core of this model, a high-order dispersion relation is derived in this paper. First of all, new formulations of wave over generally varying seabed topography are derived using Fredholm’s alternative theorem (FAT). In the new formulations, the governing equation is coupled with wave number and varying seabed e ﬀ ects. A new high-order dispersion relation can be obtained from the coupling equation. It is worth mentioning that both the slope square and curvature terms ( ( ∇ h ) 2 , ∇ 2 h , ( ∇ k ) 2 , ∇ 2 k , ∇ h · ∇ k ) of water wavenumber and seabed bottom are explicitly expressed in high-order dispersion relation. Therefore, the proposed method of coastal bathymetry reversion using the higher-order dispersion relation model is more accurate, e ﬃ cient, and economic.


Introduction
In coastal engineering, real-time updated and accurate water depth data are critical to solve many marine problems [1][2][3][4][5][6][7]. Traditionally, marine bathymetry data are mainly obtained by using sonar to carry out field exploration. Usually, the sonar equipment and the positioning system are installed in the measuring ship. The mesh point is arranged in the sounding water area for measurement, and the data processing is performed later. This field measurement is not only costly to operate, but also takes long duration. Due to the limitation of the regional environment, some areas cannot be measured, resulting in the lack of water depth data in these areas [8]. In addition, the ship survey is greatly affected by the weather conditions. Since traditional measurement methods have so many limitations, it is urgent to explore alternative methods.
With the development of aerial photography technology, radar and remote sensing technology, satellite derived bathymetry (SDB) is being adopted as a cheaper and more spatially extensive method for bathymetric mapping than traditional acoustic surveys [9][10][11]. However, the acquisition of depth information is affected by chlorophytes, suspended sediment, dissolved organic matter in air colour and other information in water. Therefore, in water areas with large interference from human activities and muddy water, extremely high requirements are put forward for the spectral analysis of remote sensing [12]. In other words, when the visibility of water is limited, it will lead to low accuracy in field measurements. J. Mar. Sci. Eng. 2020, 8,153 3 of 16 fixed in space and z = 0 is located at the calm water level. The wave potential Φ(x, y, z, t) must satisfy the following equation and boundary conditions: Φ tt + gΦ z = 0, at z = 0, Φ z + ∇h · ∇Φ = 0, at z = −h.

Fredholm's Alternative Theorem (FAT)
Fredholm's alternative theorem (FAT) [33] is introduced in this section and briefly described as follows. Mathematically, a non-homogeneous problem can be solved whenever the corresponding homogeneous problem has unique trivial solution. However, in the alternate case, when uniqueness fails, it is solvable if and only if certain orthogonality conditions are satisfied. In other words, if the homogeneous equation has a non-trivial solution, then the corresponding non-homogeneous equation has a solution if and only if the orthogonality condition is satisfied. This theory is referred to as FAT (Fredholm's alternative theorem).
Introduce a differential operator L and its adjoint N, acting on function φ 1 and φ 2 , respectively. Rewrite equation of Friedman as defined as where Lφ 2 = ∂ 2 φ 2 ∂z 2 and Nφ 1 = ∂ 2 φ 1 ∂z 2 Friedman indicated that z 1 and z 2 are arbitrary constants [33]. And J(φ 1 , φ 2 ) is the conjunct of the functions φ 1 and φ 2 , defined as If L = N, the differential is said to be formally self-adjoint. If, in addition, boundary conditions for φ 1 and φ 2 are such that the conjunct J vanishes identically, the operator is said to be Hermitian.
Considering the Laplace differential operator for water wave problems, let thus, the Laplace Equation (1) can be re-written as: Since the operator N commutes with the operator L, we may treat N as a constant in the process of dealing with operator L. In other words, Equation (1) will be first considered as an ordinary differential equation for function of z, together with boundary conditions.
For the water wave problem, the operator L is only formally self-adjoint, i.e., N = L, but not Hermitian due to the non-homogeneous boundary conditions at free surface and bottom. z 1 and z 2 are defined as water surface and seabed, respectively. Applying the FAT to the homogeneous and non-homogeneous solutions of the wave problem, i.e., substituting φ 1 and φ 2 to be solved, the orthogonality condition for waves over a varying bottom can be obtained Assuming that both φ 1 and φ 2 are Lebesgue square integrable, this theorem is valid for arbitrary boundary conditions. For water wave problems, both non-homogeneous boundary conditions on the free surface and bottom will be applicable. However, if the functions are not square integrable, the Fredholm alternative theorem can be extended to such cases as Garabedian indicated [34].

Wavenumber and Velocity Potential Coupling Equation
Assume that the velocity potential has the following form of variable separation: where Q = k(z + h), q = kh, σ = tanh(q), f = cosh Q/ cosh q is the eigen-function of homogeneous solution, φ(x, y, t) is the 2-D wave velocity potential and Φ(x, y, z, t) is the 3-D wave potential.
Substituting φ 1 = f and φ 2 = Φ, using equation (8), we have: Substituting the linearized free surface boundary condition and the bottom boundary condition, and paying attention to f = 1 or 1/ cosh q, and f z = ktanh(q) and 0 when z = 0 and z = −h, respectively, the new coupling equation for surface waves propagating over generally varying seabed may be formulated by utilizing equations above based on FAT governing the velocity potential φ (x, y, t) and the wavenumber k (x, y): After a long derivation of integration along the vertical direction in the above equation which is shown in the Appendix A, the time-dependent equation governing the velocity potential φ(x, y, t) and the wave number k (x; y) is obtained: where The water depth variable h is the local water depth and is expressed as the function of x and y, i.e., h(x, y); ω is the wave frequency, g is the gravitational acceleration and k(x, y) is the wave number. For simplification of derivation, σ = tanh(q) is used, where q = kh is the product of wave number and water depth. (∇h) 2 represents the square term of the slope of the seabed, ∇ 2 h represents the seabed topography curvature term, (∇k) 2 represents the square term of the slope of the water wave number, and ∇ 2 k represents the water wave number curvature term.

New High-Order Dispersion Relation Formulation
A new high-order dispersion relation can be obtained from the derived seabed wave coupling equation, and the water depth can be directly calculated under medium and small wave conditions. Let the wave velocity potential function be: where A(x, y, t) is the amplitude of φ(x, y, t), while θ is the phase of φ(x, y, t), which are defined by k = ∇ θ, ω= − ∂ θ ∂t . Further, c.c is complex conjugate. The coupling equation can be decomposed into real and imaginary parts of corresponding dispersion relation and wave action conservation, respectively [35]. Thus, the real part is which leads to a new high-order dispersion relation: where In this paper, the change of wave number comes from the change of water depth under the ideal condition of single sinusoidal topography. The equation can be recovered from the coupling equation by using the relation of the wave-number slope and curvature with those of the bottom [29], which are obtained from: According to Equations (26)- (35), the change of wave number and the change of velocity potential function are proportional to the square of slope and curvature of seabed. By substituting β 1−8 and Equations (26)- (35) into Equation (17), the specific formula of the expression of the high-order dispersion relation of the seabed with varying depth is obtained as follows: where water depth variable h is the local water depth and is expressed as the function of x and y, i.e.,h(x, y), ω is the wave frequency, g is the gravitational acceleration, and k(x, y) is the wave number. For simplification of the derivation, σ = tanh(q) is used, where q = kh is the product of wave number and water depth. Finally, a new high-order dispersion relation is obtained, which contains the square of the slope terms and the curvature terms of the seabed.

Numerical Solution Procedure
For simplicity, the two-dimensional wave problem is studied in this section. For 2-D waves in steady state φ(x, y, t) = Re φ(x, y)e −iωt , and the reflection coefficient Cr is calculated by C r = Re φ(i) − 1 , where the symbol Re represents the real part of a complex value [36]. The coupling equation becomes where The boundary conditions for a patch of rippled beds are where φ I = e ikx is the incident wave of unit amplitude, x 1 and x 2 represent the up-wave and down-wave limits of the computational grid, and L is the length of computation domain. These boundary conditions have been given previously by Kirby [37].
Considering potential applications in the three-dimensional wave, a FEM model is developed for general topography. The FEM model has advantages in the three-dimensional problems with complex geometries, where it is desirable to use irregular meshes.
Multiplying the entire left-hand side of Equation (43) with a weight function w, and integrating over the domain (0, L) gives the weighted-residual statement: Mathematically, the above equation is a statement that the numerical error needs to be zero in the weighted-integral sense. The trading of differentiability from φ to w provides the weak form The trading of differentiability from φ to w can only be performed if it leads to boundary terms that are physically meaningful. The choice of the approximation for weight function gives the boundary term P φ x φ, which has physical meaning of energy flux through a section. It is easy to find that the primary variable and the secondary variable are φ and φ x respectively. Thus φ x L 0 is the natural boundary condition. Using the notation of Reddy [38], we have where and are bilinear and linear forms, respectively. Here, the bilinear form is asymmetric. For a typical element, φ is approximated by where N j are cubic shape functions and φ j are unknowns at the nodes. The water depth h(x), slope and curvature of both h and k at each element in the FEM scheme can readily be evaluated as In the FEM schemes, Bubnov-Galerkin method is adopted, thus, the solution shape functions are used as weighting functions.
As noted by Reddy, not all differential equations admit the functional formulation, and in order for the functional to exist, the bilinear form must be symmetric in its arguments [38]. Since the weak form statement is equivalent to the differential equation and the specified natural boundary condition of the problem, the weak form FEM is used in this study.

Numerical Results vs. Experiment Data
This new method of nearshore bathymetry is characterized by the nonlinear relation of wave number caused by the change in seabed slope and curvature, which is embodied in the square term of the slope and the curvature term of the seabed. In order to ensure the accuracy of the model, the numerical results are compared with the experimental values. In this paper, a type of seabed topography is selected for verification: single sinusoidal ripples.
The single sinusoidal seabed topography is expressed as follows: In the two-dimensional coordinate system, the high-order dispersion relation (36) becomes: where The new high-order dispersion relation can be applied to surface waves propagating over single sinusoidal ripples. The dispersion relation is the real part of the water wave governing equation.  [39], which is equivalent in that the high-order dispersion relation is compared with the linear dispersion relation. Considering the same cases of their papers, the single sinusoidal seabed parameters are listed in Table 1.  It can be seen that the present model and Miles and Chamberlain's model perform exceptionally well in predicting experimental data in Figure 1a-c [40]. However, there are still differences, especially at the subpeak position. In case (a,b), Miles and Chamberlain's model is the same as the MSE, which is closer to the experimental data. However, the two models cannot predict the trend from high to low to high of the experimental data around 2k/K = 2. In case (c), the present model shows better results when predicting the subpeak trend compared to the other two models despite a slight phase shift. In future research, we should explain the reasons for those situations, so that the present model will be closer to the experimental data.
The numerical results of the present model are presented in Figure 1a-c, which show that the simulated results of the single sinusoidal seabed agree with the experimental values very well in general. The fitting effect between the present model and laboratory is much better than that of the original mild-slope equation (MSE), which can be considered that the high-order dispersion relation is improved compared with the linear dispersion relation. It is more noteworthy that practical verification should be illustrated in future research, including comparisons with more complete models for water wave propagation over variable bathymetry, as the consistent coupled-model theory presented by Athanassoulis and Belibassakis [41].
Correlation coefficients are calculated between the numerical solutions of the high-order dispersion relation model and the experimental values in the three cases, which are shown in Table 2. The correlation , where S i and O i denote predicted and theoretical or observed data, respectively, S and O are mean values of S i and O i , and M is the number of evaluation point [42]. It can be observed that the correlation coefficients increase with the number of ripples. In other words, the high-order dispersion model is more adaptable to real undulating ocean bed forms.
In other models, it is about 0.75 of the maximum correlation coefficients between the solution generated by surface wave and seabed [42]. When considering the influence of bottom shape changes and the wave number changes caused by the topographic change, the high-order dispersion relation model shows a higher fitting degree and accuracy than the linear dispersion relation. Table 3 shows the wave number and wave wavelength (m) when the water depth (cm) is different under single sinusoidal seabed profile (the shape is consistent with case a-c) with a 1s wave period in the lab wave scale.  It can be seen that the present model and Miles and Chamberlain's model perform exceptionally well in predicting experimental data in Figure 1a-c [40]. However, there are still differences, especially at the subpeak position. In case (a,b), Miles and Chamberlain's model is the same as the MSE, which is closer to the experimental data. However, the two models cannot predict the trend from high to low to high of the experimental data around 2k/K = 2. In case (c), the present model shows better results when predicting the subpeak trend compared to the other two models despite a slight phase shift. In future research, we should explain the reasons for those situations, so that the present model will be closer to the experimental data.
The numerical results of the present model are presented in Figure 1a-c, which show that the simulated results of the single sinusoidal seabed agree with the experimental values very well in general. The fitting effect between the present model and laboratory is much better than that of the original mild-slope equation (MSE), which can be considered that the high-order dispersion relation is improved compared with the linear dispersion relation. It is more noteworthy that practical verification should be illustrated in future research, including comparisons with more complete models for water wave propagation over variable bathymetry, as the consistent coupled-model theory presented by Athanassoulis and Belibassakis [41].  Table   2. The correlation coefficient  [42]. It can be observed that the correlation coefficients increase with the number of ripples. In other words, the high-order dispersion model is more adaptable to real undulating ocean bed forms.  When the seabed profile is the same, the wavelength and wave number of the high-order dispersion relation model and linear dispersion model under different water depths are given in Table 3. When the water depth is shallow or deep, that is, the water depth is 10.0, 15.6, and 30.0 cm, the results calculated by the present dispersion relation are not much different from the wave number and wavelength. However, when the water depth is 20.0 cm, that is, at the intermediate water depth range, the accuracy difference of the high-order dispersion relation and linear dispersion relation is five times that of other water depths as mentioned above. This indicates that the high-order dispersion relation has better adaptability in the intermediate water depth. Therefore, the new high-order dispersion relation is used to calculate the intermediate water depth with higher accuracy. Considering Equation (55), we can see that the coefficient in front of the curvature term and slope square term is a function of q. Therefore, the high-order dispersion relation mentioned above has higher accuracy in the intermediate water depth, which can be observed from Figure 2. As shown in Figure 2, the solid line represents the coefficient in front of the curvature term of seabed topography, while the dotted line represents the coefficient in front of the square of slope term of seabed topography. The bottom form is the same shape as case a-c. It can be observed that the value of 1 P and 2 P fluctuate greatly when the value of q is in the interval (1,5). That is to say, the depth inversion has to consider the bottom topography in the intermediate because the nonlinear effect caused by terrain change on wave propagation is difficult to ignore. At the same time, it can be seen that the influence of the terrain curvature and slope square terms on the wave propagation cannot be ignore when the curvature or slope of seabed is very large. This can also explain why the error increases significantly when linear dispersion relation is used in the intermediate water depth.

Discussion
In this section, the future research topic will be envisaged, and a set of bathymetry inversion models with high precision is preconceived by using a high-order dispersion relation. The following is the design of the implementation steps of the water depth inversion: Firstly, a digital image or video image of the surface-temporal sequence of waves in the target sea area can be collected from UAVs (unmanned aerial vehicle) among other ways such as radar and remote sensing.
Then, the acquired image is processed with the maximum total coherence frequency band, and the corresponding frequency and wave number pairs are generated at each calculation point [43]. For the obtained digital image or video image, the Fourier transform is used to decompose time-varying It can be observed that the value of P 1 and P 2 fluctuate greatly when the value of q is in the interval (1,5). That is to say, the depth inversion has to consider the bottom topography in the intermediate because the nonlinear effect caused by terrain change on wave propagation is difficult to ignore. At the same time, it can be seen that the influence of the terrain curvature and slope square terms on the wave propagation cannot be ignore when the curvature or slope of seabed is very large. This can also explain why the error increases significantly when linear dispersion relation is used in the intermediate water depth.

Discussion
In this section, the future research topic will be envisaged, and a set of bathymetry inversion models with high precision is preconceived by using a high-order dispersion relation. The following is the design of the implementation steps of the water depth inversion: Firstly, a digital image or video image of the surface-temporal sequence of waves in the target sea area can be collected from UAVs (unmanned aerial vehicle) among other ways such as radar and remote sensing.
Then, the acquired image is processed with the maximum total coherence frequency band, and the corresponding frequency and wave number pairs are generated at each calculation point [43]. For the obtained digital image or video image, the Fourier transform is used to decompose time-varying pixels and the Fourier coefficients of the time-varying pixels are further normalized. Then, select the subset of normalized Fourier coefficients around the depth location to be determined, and calculate the intersection of all pixels in the subset density spectrum. Matching wave phases are determined for each selected frequency, and generate a set of frequencies and corresponding wave numbers at each calculation point.
Finally, the obtained wave numbers and frequencies are substituted into the high-order dispersion relation model for iterative calculation to obtain the target sea depth. The final obtained depth value of the target sea area is a continuous time process. The specific implementation process is shown in Figure 3.

Conclusions
The resulting reflection coefficient for a single sinusoidal bottom is compared with the results of numerical solutions of the mild-slope equation and modified mild-slope equation. The present model which includes five terms (       , h h k k h k ) can perform exceptionally well in predicting experimental data, especially predicting the trend from high to low to high of the experimental data around the subpeak (2k/K = 2). Although the capability of this model in the weak reflection area is not good enough, in future research, more attention will be paid to this aspect to improve the accuracy of the model. A new high-order dispersion relation emerges that offers significantly increased accuracy over a single sinusoidal bottom, which is derived from the real part of the coupling equation.
We have presented a model frame that space-time sequence digital images or videos combine with high-order dispersion relation for bathymetry mapping. The high-order dispersion relation model has higher adaptability in the intermediate water depth with complex seabed, such as bed with sand bar formation etc.
The reasons why selecting a single sinusoidal seabed form for model verification can be

Conclusions
The resulting reflection coefficient for a single sinusoidal bottom is compared with the results of numerical solutions of the mild-slope equation and modified mild-slope equation. The present model which includes five terms ((∇h) 2 , ∇ 2 h, (∇k) 2 , ∇ 2 k, ∇h · ∇k) can perform exceptionally well in predicting experimental data, especially predicting the trend from high to low to high of the experimental data around the subpeak (2k/K = 2). Although the capability of this model in the weak reflection area is not good enough, in future research, more attention will be paid to this aspect to improve the accuracy of the model. A new high-order dispersion relation emerges that offers significantly increased accuracy over a single sinusoidal bottom, which is derived from the real part of the coupling equation.
We have presented a model frame that space-time sequence digital images or videos combine with high-order dispersion relation for bathymetry mapping. The high-order dispersion relation model has higher adaptability in the intermediate water depth with complex seabed, such as bed with sand bar formation etc.
The reasons why selecting a single sinusoidal seabed form for model verification can be summarized as follows: (1) There are internationally recognized model experiments and valuable data from Davies, which are the basis of model validation. (2) It can be widely used in practical engineering. For example, the profiles of the most artificial sand dam and submerged dike are sinusoidal. (3) Mathematically, any seabed shape can be decomposed into sinusoidal form by the Fourier transform.
The acquisition of the space-time sequence digital images or video images of the surface wave in the target sea area is fast in time, low in cost, easy and secret to operate, eliminating the drawbacks and risks of conventional bathymetric inversion. Using the new high-order dispersion relation, the calculated water depth is a continuous time process with high precision. Combining the advantages mentioned above, the new method using high-order dispersion greatly makes the acquisition of water depth efficient and accurate.
Furthermore, this method not only achieves the real-time monitoring of the water depth of the target sea area, but also provides real-time water depth data for the construction of most marine engineering and waterway transportation projects. Therefore, the above method takes advantages of the remote sensing technology to obtain the synoptic bathymetry of coastal waters safely, economically and quickly. It is of great significance to bathymetric data related to coastal engineering.
Due to a lack of data right above the sinusoidal patches in the experiment of Davies, the present model is largely comparable to experiment data between wave gauge 1 and wave gauge 2 for wave propagation through the ripples [39]. While this is not the direct verification for the high-order dispersion relation, which considers the combined refraction diffraction and reflection that is induced by bathymetry, including slope square and curvature terms ((∇h) 2 , ∇ 2 h, (∇k) 2 , ∇ 2 k, ∇h · ∇k) of water wavenumber and seabed bottom, its application to simulate practical field experiments on coastal wave transformation that is induced by bathymetry will also improve the model. Future research of the present model will facilitate its extension to various directions such as, e.g., to free-surface elevation, to water wave potential, and to three-dimensional problems.