Impact Statement
Ship wakes are known to influence coastal sediment transport and affect the habitats of marine life. Monitoring ships is crucial for enforcing drug and fishery regulations as well as for controlling port traffic. In this study, we apply a machine learning-based technique to analyse surface disturbances that can serve as indicators of real ships. Our findings demonstrate the effectiveness of convolutional neural networks in characterising both the velocity and length scale of ships from their wakes. Consequently, this technique has the potential to greatly enhance ship monitoring in operational settings.
1. Introduction
A distinct wave pattern is generated when a disturbance moves at a constant speed at the water’s surface. This phenomenon, also known as the Kelvin wake (Thomson, Reference Thomson1887), has been extensively studied because the wake structure is closely related to the kinematics and geometry of the sources of the disturbance. Note that Kelvin wakes should not be mistaken for Kelvin waves, which are produced by the geostrophic balance between the Coriolis force and gravity in the direction transverse to wave propagation (Kundu et al., Reference Kundu, Cohen and Dowling2016). Ship-induced Kelvin wakes represent an important anthropogenic impact in the coastal environment. For instance, a recent study measured the shipping events in the Venice Lagoon and identified ship wakes as one of the causes of shoreline erosion (Scarpa et al., Reference Scarpa, Zaggia, Manfè, Lorenzetti, Parnell, Soomere, Rapaglia and Molinaroli2019). Hence, understanding the hydrodynamics of these ship-related motions is critical for reducing their adverse environmental impact. On the other hand, ship wakes have enabled the remote sensing-based monitoring of maritime traffic, in which the ship’s location and velocity data prove valuable for ship collision prevention.
Because of the complexity in surface wave–ship structure interaction, the study of Kelvin wakes is often facilitated by applying a moving pressure to the free surface of an ideal fluid. For selected disturbance velocity and geometry, one can then derive an analytical solution for the wake angle,
$ \varphi$
max
, defined as the angle between trendlines consisting of local maxima of the surface elevation. In linear theory (Thomson, Reference Thomson1887),
$ \varphi$
max
was found to be a constant value of arcsin (1/3) = 19.5∘, independent of the speed of the disturbance. However, this was found to contradict observations. Rabaud and Moisy (Reference Rabaud and Moisy2013) examined satellite images of ship wakes and obtained a scaling law for large Froude numbers
$ \varphi_{max} $
∼ Fr
−1 by arguing that the waves generated by the surface disturbance may not exceed the ship length L, which imposes a lower limit on the characteristic wavenumber of the wake (here,
$Fr=U/\sqrt {gL}$
, where U is the ship velocity and g is the acceleration due to gravity). This hypothesis was later proved unnecessary as Darmon et al. (Reference Darmon, Benzaquen and Raphaël2014) derived the asymptotic value of the angle corresponding to the maximum amplitude of the waves and reproduced the same scaling law. In a following study, Benzaquen et al. (Reference Benzaquen, Darmon and Raphaël2014) extended their result to an elliptical pressure disturbance and discovered the relation
$\varphi _{max}\sim \sqrt {W}/\textit {Fr}$
, where W is the ellipse aspect ratio. Miao and Liu (Reference Miao and Liu2015) studied the wave patterns formed when the pressure disturbance is arbitrary. By decomposing the pressure into Fourier components and integrating the surface response corresponding to each individual component, they found that
$ \varphi$
max
may be a complex function of the Froude number. For example, they showed that, in the presence of a rectangular uniform pressure, the scaling law becomes
$ \varphi_{max} $
∼ Fr
−2. It is worth noting that
$ \varphi$
max
is a nonlinear function of the disturbance parameter (e.g. the speed U), even for an infinitesimally small pressure disturbance, and the resulting linear wave response (see the scaling laws mentioned above). Besides the aforementioned transition with Froude number, the wake angle is also known to be a function of the disturbance geometry. In general, a theoretical expression derived from asymptotic analysis is applicable only to highly idealised free surface disturbances.
Ship detection and characterisation from remote sensing images is an active area of research, where the objectives may include identification of a ship and the prediction of the ship’s velocity. A widely used technique for wake identification is the Radon transform (e.g. Copeland et al., Reference Copeland, Ravichandran and Trivedi1995), which integrates a raw two-dimensional (2-D) signal along all possible lines. Hence, the Radon transform performs best when used for identifying linear patterns in an image, such as the asymptotic wakes mentioned above. However, the wake angle is not a well-defined variable, as the definition of a local maximum of the surface elevation can be ambiguous as a result of the contamination by background surface waves and the observational errors in remote sensing. Recent years have seen the application of emerging machine learning techniques to the inverse problems in ship wake characterisation. For example, dynamic mode decomposition, an unsupervised learning technique, has been applied to the separation of the ship wake pattern from the background wind waves (Zhang and Jiang, Reference Zhang and Jiang2020). Supervised learning methods based on the convolutional neural network (CNN) have also been used for detecting ship wakes from remote sensing images. Kang and Kim (Reference Kang and Kim2019) applied CNN to satellite images, where the output data of the neural network are the latitude and longitude coordinates of a subregion marked as the ship, wake or sea. The ship velocity was then estimated from the azimuth offset by performing edge detection and Radon transform of the subregion. In another study (Xue et al., Reference Xue, Jin, Qiu and Yang2022), the output of the CNN not only includes the coordinates of a subregion, but also the ship’s location and the angle of the wakes. A majority of these studies focus on the feasibility and accuracy in terms of operational applications while neglecting the hydrodynamics constraints.
In this study, we investigate the inverse problem of inferring the properties of a moving disturbance from Kelvin wakes using CNNs. Unlike existing studies that typically focus solely on predicting velocity, we explore the potential of CNNs to simultaneously predict both velocity and the disturbance length scale. Additionally, we aim to interpret the CNN’s prediction capabilities and assess the impact of wave nonlinearity on accuracy. The remainder of this paper is organised as follows. In § 2, we introduce the techniques needed for building the prediction framework, which include a phase-resolved wave model for obtaining the training and testing data, the mathematical function for a free disturbance and the neural networks used. We outline the simulation set-up and training data generation in § 3. The training results are provided in § 4. Conclusions are given in § 5.
2. Methodology
2.1 Wave model
In this study, we assume that the water is incompressible, inviscid and irrotational, while the ship is modelled as a surface pressure disturbance. This strategy allows for an acceptable computational cost and has been widely used in theoretical and numerical studies of Kelvin wakes (see e.g. Benzaquen et al., Reference Benzaquen, Darmon and Raphaël2014; Miao and Liu, Reference Miao and Liu2015; Pethiyagoda et al., Reference Pethiyagoda, Moroney, Lustri and McCue2021). In the present study, the deep-water assumption is adopted, as the focus is on the generation of Kelvin wakes. In future studies, their impact on sediment transport can be evaluated as these waves propagate into shallower waters. Under the potential flow assumption, it can be shown that the surface gravity waves in deep water are uniquely determined by the surface elevation η(x, y, t) and the surface velocity potential
$ \Phi $
S
(x, y, t) ≡
$ \Phi $
(x, y, η(x, y, t),t) (Zakharov, Reference Zakharov1968).
For a steady pressure distribution moving at a constant speed U, it is convenient to rewrite the Zakharov equations in the moving frame of reference (Dommermuth and Yue, Reference Dommermuth and Yue1988)
\begin{align} \eta _t - U\eta _x+\nabla \eta \cdot \nabla \Phi ^S -(1+\nabla \eta \cdot \nabla \eta )W_S=0, \\[-24pt] \nonumber \end{align}
\begin{align} \Phi _t^S -U\Phi ^S_x+g\eta +\frac {1}{2}\nabla \Phi ^S\cdot \nabla \Phi ^S -\frac {1}{2}(1+\nabla \eta \cdot \nabla \eta )W_S^2+\frac {p_a(\boldsymbol {x})}{\rho _w}=0, \\[6pt] \nonumber \end{align}
where the gradient operator in the horizontal plane is defined as ∇ ≡ (∂/∂x, ∂/∂y),
x
= (x, y) is the horizontal coordinate vector,
$ W_{S} = \Phi_{z} $
is the surface vertical velocity, p
a
(
x
, t) is the external pressure disturbance imposed at the surface and ρ
w
is the density of water.
The above equations can be solved with the high-order spectral (HOS) method (Alam et al., Reference Alam, Liu and Yue2009; Alam, Reference Alam2012; Dommermuth and Yue, Reference Dommermuth and Yue1987) in a periodic domain by assuming that the velocity potential can be expanded in a perturbation series
$\Phi (\boldsymbol {x},\,z,\, t)=\sum _1^M\Phi ^{(m)}(\boldsymbol {x},\, z,\, t)$
with respect to a characteristic wave steepness, where M is the maximum perturbation order. At each order, it is assumed that
$ \Phi^{(m)} $
, as
$ \Phi $
, satisfies the Laplace equation, ∇2
$ \Phi^{(m)} = 0 $
, and the deep-water boundary condition, ∇
$ \Phi^{(m)} $
(
x
, z→ −∞, t) → 0. Further assuming periodic boundary conditions in the horizontal directions,
$ \Phi^{(m)} $
can be written as the eigenfunction of the Laplace equation
\begin{align} \Phi ^{(m)}(\boldsymbol {x}, z, t)=\sum _{j=1}^{\infty }\sum _{l=1}^{\infty } \Phi _{jl}^{(m)}(t) \exp \left (|\boldsymbol {k_{jl}}|z+i\boldsymbol {k_{jl}}\cdot \boldsymbol {x}\right ), \end{align}
where
$ \Phi^{(m)}_{jl} (t) $
are the Fourier coefficients to be calculated from
$ \Phi^{(S)} $
, and
k
jl
= (2jπ/L
x
, 2lπ/L
y
), with L
x
and L
y
being the computational domain size in the x and y directions, respectively. To determine
$ \Phi^{(m)}_{jl}(t) $
, the Dirichlet problem is first solved
\begin{align} \left \{ \begin{aligned} & \psi ^{(1)}=\Phi ^S(x,y,t),\quad m=1\\ & \psi ^{(m)}=-\sum _{l=1}^{m-1}\frac {\eta ^l}{l!}\frac {\partial ^l\psi ^{(m-l)}}{\partial z^l}.\quad m=2,3,\cdots ,M ,\end{aligned}\right . \end{align}
where ψ
(m) ≡
$ \Phi^{(m)} $
(
x
, 0, t), m = 1, 2, ⋯, M. Subsequently, the coefficients
$ \Phi^{(m)}_{jl}(t) $
can be computed in an efficient manner by taking the 2-D Fourier transform of ψ
(m).
The vertical velocity at the free surface is calculated by
\begin{align} W_S=\sum _{m=1}^{M}\sum _{l=0}^{M-m}\frac {\eta ^l}{l!}\frac {\partial ^{l+1}\psi ^{(m)}}{\partial z^{l+1}}. \end{align}
The system can then be integrated in time from η and
$\Phi^{S}$
prescribed at t = 0. While the HOS method was originally developed for solving the Zakharov equations, it can also be directly applied to the system in the moving frame of reference (Dommermuth and Yue, Reference Dommermuth and Yue1988; Miao and Liu, Reference Miao and Liu2015). To avoid numerical instabilities caused by excessive local steepness, an adaptive low-pass filter may be applied to η and
$ \Phi^{S}$
at each time step when M > 1 (Xiao et al., Reference Xiao, Liu, Wu and Yue2013).
2.2 Pressure disturbance
In our simulations, the pressure disturbance is placed near the end of the domain to allow for the wake pattern formation (see figure 1a). To reduce wave reflection, we impose a relaxation zone at the boundaries of the computational domain (Mei et al., Reference Mei, Stiassnie and Yue2005). In the present study, we consider three types of pressure distributions that have been widely used in previous studies for generating Kelvin wakes. For each type of pressure distribution, an example of their spatial distribution is shown in figure 1(b).
Figure 1. Schematic of (a) the computational domain, (b) the pressure distributions representing different ship types and (c) the CNN-based prediction framework. In (a), the region for extracting the training data is marked by the shaded square in red. In (b), the coordinates are shifted such that the geometric centres of the pressure distributions are at the origin (i.e. X = x − x c and Y = y − y c ), and the normalised pressure is defined as p* = (p−p min )/(p max −p min ), where p min and p max denote the minimum and maximum pressures.
For an elliptical pressure distribution, we follow the definition of Dommermuth and Yue (Reference Dommermuth and Yue1988)
\begin{align} P_e(x, y)=P_0 e^{-\Pi (r/R)} ,\end{align}
where r = (X
2+Y
2/W
2)1/2, X = x − x
c
, Y = y − y
c
, P
0 is the peak magnitude of the pressure and
$ \Pi $
is a smooth function. Here, x
c
and y
c
are the coordinates of the geometric centre, R is the radius of the ellipse and W is the ellipse aspect ratio. While the exact form of P
e
varies slightly in different studies, the impact on the far-field Kelvin-wake pattern is expected to be small.
The pressure distribution corresponding to a monohull is defined as the sum of two circular distributions
\begin{align} P_m(x, y)=P_0 \left [e^{-\gamma \left [(X-L_h / 2)^2+Y^2\right ]}- e^{-\gamma \left [(X+L_h / 2)^2+Y^2\right ]}\right ], \end{align}
where L h is the hull length and γ is the bow and stern geometry parameter.
The catamaran is defined as the sum of two monohull distributions at a distance of B
\begin{align} P_c(x,y)=P_m(x,y-B/2)+P_m(x,y+B/2). \end{align}
2.3 Convolutional neural network
The basis of a CNN is the convolutional layer, which can be written as
\begin{align} \boldsymbol {B}_{i,j}=\sum _{l=-\Delta }^{\Delta }\sum _{m=-\Delta }^{\Delta }\boldsymbol {A}_{l,m}\boldsymbol {w}_{i+l,j+m}+\boldsymbol {u}, \end{align}
where
A
and
B
denote the input and output data of the layer, respectively,
w
is a weighting matrix,
u
is an optional bias matrix and
$ \Delta $
is the filter size. Both
w
and
u
are determined in the training process. Hence,
B
can be seen as the discretised convolution (more accurately, cross-correlation) between a 2-D array A and
w
when
u
= 0. Hence, compared with a fully connected layer used in a multi-layer perceptron, a convolutional layer can identify the spatial structures of the data and is therefore suitable for image-like input.
LeNet-5 is a classical CNN architecture originally designed for handwritten digit recognition (Lecun et al., Reference Lecun, Bottou, Bengio and Haffner1998). It comprises three convolutional layers followed by pooling layers and two fully connected layers, which proved to be sufficiently accurate in classifying handwritten digits. While LeNet-5 demonstrated the power of CNNs in extracting spatial structures from 2-D data, it has a modest complexity compared with contemporary models. A more widely used modern CNN is the residual network (ResNet) (He et al., Reference He, Zhang, Ren and Sun2016), which utilises the residual connections and addresses the vanishing gradient problem during training. The residual blocks enable the network to focus on the difference between the input and output of each layer, which significantly improves the network’s ability to learn deeper representations and enhances its overall performance. In the present study, LeNet-5 and ResNet-18, a residual neural network with 18 layers, are used. The prediction framework is shown in figure 1(c) (see the Supplementary Material for the architectures of the CNNs).
2.4 Wave model validation
To demonstrate the accuracy and numerical capability of the HOS method, we perform simulations of Kelvin wakes produced by an elliptical pressure forcing applied to the free surface. We first vary the speed U (equivalently Fr) of the disturbance, while keeping the magnitude of the ellipse aspect ratio W and the semi-minor axis R constant (see dataset V
Fr in table 1). For the dataset V
W, we change the values of W instead. Figure 2(a) shows an example of the Kelvin wakes generated. Here, the normalised surface elevation is defined as
$ \eta^{\ast} = (\eta- \eta_{min})/(\eta_{max} - \eta_{min}) $
, where
$ \eta_{min}) $
and
$ \eta_{max}) $
denote the minimum and maximum elevations, respectively. From the instantaneous surface wave field, we identify the maximum elevation along each y direction and perform a linear regression to obtain the wake lines (see also the white dashed lines in figure 2a) and the wake angle
$ \varphi_{max} $
. The Froude number is defined as
$\textit {Fr}=U/\sqrt {gR}$
. The computed wake angles are then plotted as a function of Fr (figure 2b) and W (figure 2c), respectively. For both datasets, our simulation results agree well with the theoretical scaling law of
$\varphi _{max}\sim \sqrt {W}/\textit {Fr}$
, derived by Benzaquen et al. (Reference Benzaquen, Darmon and Raphaël2014).
Table 1. Parameters used in the HOS-based simulations for asymptotic analysis and model validation, where N x and N y are the grid numbers in the x and y directions, respectively, L is a characteristic length scale and U is the speed of pressure disturbance
Figure 2. Simulation results for validating the response of the wave system: (a) an example of the (normalised) surface elevation; (b) the maximum surface elevation angle as a function of Froude number; (c) the maximum surface elevation angle as a function of the ellipse aspect ratio. In (a), the locations of the local maximum surface elevation are denoted by the white dashed lines. The theoretical scaling law
$\varphi _{max}\sim \sqrt {W}/\textit {Fr}$
(Benzaquen et al., Reference Benzaquen, Darmon and Raphaël2014) is denoted by the black dashed lines in (b) and (c).
We then examine the nonlinear error related to the perturbation expansion in the HOS method. Let
$ \epsilon = \text{max}(\eta^{2}_{x} + \eta^{2}_{y})^{1/2} $
denote the characteristic steepness of the surface wave field, the error
$ e = \| \eta_{M} - \eta_{true} \| $
then should scale as
$ e ~ \epsilon_{(M+1)} $
, where
$ \eta_{M} $
is the surface elevation obtained from HOS simulations performed at order M and η
true
is the corresponding ground truth. This strategy is equivalent to the validation of the HOS method by examining the error in
$ \Phi_{z} $
, which was adopted in the simulations of Stokes wave (Dommermuth and Yue, Reference Dommermuth and Yue1987) and nonlinear surface capillary waves (Pan, Reference Pan2020). Ideally,
$ \eta_{true} $
needs to be calculated from analytical solutions, such as the Crapper wave used by Pan (Reference Pan2020). However, no such solutions exist for the nonlinear waves generated by the moving pressure of an arbitrary shape and we treat
$ \eta_{M}(M = 5) $
as
$ \eta_{true} $
instead.
Simulations are performed for M = 2, 3, 4, 5 (see datasets V
M in table 1). For each M, we repeat the simulations at different characteristic wave steepnesses by adjusting the strength of the peak pressure P
0. Figure 3 shows the variations of the normalised error
$ e / e_{0} $
with the steepness at different maximum perturbation orders. As expected, the error growth with increasing M agrees well with the theoretical prediction.
Table 2. Parameters used in the HOS-based simulations for training data generation. Note that each case includes multiple simulation runs at varying parameters, and the subscripts ‘e’, ‘m’ and ‘c’ in the case names denote the ellipse, monohull and catamaran types of disturbance, respectively
Figure 3. Simulation error as a function of the characteristic wave steepness. The errors are normalised by e 0, their corresponding values at the smallest steepness. The black-dashed lines correspond to the theoretical scaling of the error with the steepness
$ e \sim \epsilon^{M+1}$
, where M is the maximum perturbation order used in the simulations.
3. Training data generation
We perform 787 simulations to generate the training data across a wide range of parameters (see datasets T e , T m and T c in table 2). Since the length scales of the far-field Kelvin-wake patterns are often much larger than the characteristic size of the pressure distribution (see the validation cases in figure 2 for example), we choose larger grid numbers and smaller domain lengths compared with the validation cases to further reduce the impact of the (possible) Gibbs phenomenon. In all simulations, the surface elevation data are collected at a time instant when the Kelvin wakes are sufficiently close to the end of the domain. To reduce the computational cost, we set the pressure to an infinitesimally small amplitude and perform the simulations using the linear model by setting M = 1. We note that a novel machine learning-based technique has been developed to extend the HOS method, enabling predictions for ocean surface gravity waves to be made more than two orders of magnitude faster than with the original method (Mohaghegh et al., Reference Mohaghegh, Murthy and Alam2021). This strategy offers an efficient way to generate data for future studies.
To obtain the input to the neural network, the far-field wake pattern η in the region of (L x /4 ≤ x ≤ 3L x /4, L y /4 ≤ y ≤ 3L y /4) is first extracted from the raw simulation data such that the neural network does not memorise the surface pattern near the location of the pressure forcing (see figure 1a). A downsampling operation is then applied to the resulting 513 × 513 matrix to create the training data. Two spatial resolutions are considered: (L x /64, L y /64) and (L x /128, L y /128). The corresponding matrix sizes of the training data are 33 × 33 and 65 × 65, respectively. To increase the size of training samples, we then apply an affine transformation to the surface elevation data, including rotating between 0∘ and 360∘ at an interval of 10∘ and scale up to 160 % of the original image. The details on data augmentation are provided in the Supplementary Material. To improve the training stability, the augmented surface elevation is normalised, η* = (η−η min )/(η max −η min ). In the final step, different levels of random error up to 30 % of the maximum wave elevation are added to η* to address the overfitting issue commonly faced in deep learning, resulting in 4.53 × 105 images in total as the input data (training data available in UC San Diego’s data archive https://doi.org/10.6075/J0GH9JCX).
Figure 4. An example of normalised surface elevation (a) computed using the HOS method (i.e. raw data); and (b) obtained from the raw data using augmentation technique.
To demonstrate the effect of affine transformation, we show an example of the raw surface elevation data obtained from the HOS simulation (figure 4a) and the transformed data (figure 4b), obtained by performing a 50∘ rotation, scaling at a ratio of 1.6 uniformly in both directions, and adding a random noise of up to 30 % of (η max −η min ). The output of the neural network includes the speed U and the length scale of the moving pressure L, which corresponds to R for the elliptical pressure distributions, and L h for the monohull and catamaran distributions. A similar normalisation is performed to obtain the actual output of the neural network v = [U * L *] T where U * = (U−U min )/(U max −U min ) and L * = (L−L min )/(L max −L min ).
In the training process, we define the loss function as the mean squared error between the neural network prediction and the ground truth J = ( v pred − v true )2, where v pred and v true correspond to the CNN prediction and the ground truth of the vector v = [U * L *] T , respectively. The Adam optimiser is adopted and the initial learning rate is 0.001 (Kingma and Ba, Reference Kingma, Ba, Bengio and LeCun2015). All input images are first randomly shuffled and approximately 4.1 × 105 samples are then selected for training, while the other 4.5 × 104 samples are reserved for validation. Both models are used for training, and for the 65 × 65 input resolution, ResNet-18 has approximately 11 million parameters while LeNet-5 has 0.33 million.
4. Training results
4.1 Loss function and accuracy
For both models and input resolutions, the loss function converges within approximately 400 epochs (figure 5). At both resolutions, the LeNet-5 model shows a faster convergence rate than the ResNet-18 model because of its smaller parameter size. Meanwhile, the loss function decreases by approximately three orders of magnitude for the ResNet-18 model, which is significantly better than LeNet-5. We remark that using images at a different resolution may affect the performance of the ResNet-18 model.
Figure 5. Convergence of the loss function J normalised by its initial value J 0.
Figure 6. The CNN prediction and ground truth of the normalised (a) velocity U* and (b) length scale L* of the disturbance. Plotted are results for ResNet-18 trained using images with a resolution of 65 × 65. The black dashed lines correspond to a perfect prediction. Here, the size of the markers increases with the number of data points they denote.
To examine the performance of each model, we reconstruct the output vector
v
when the input data are selected from the validation data. The CNN predictions are plotted against the ground truth in figure 6 for resolution 65 × 65. For the velocity, the neural network shows an excellent performance with few outliers, whereas for the length scale, more outliers exist across all ranges of the length scale. We further define the accuracy function A(x) as the percentage of validation data that satisfies
$ \|x_{pred} - x_{true} \|_{2} / \|x_{true}\|_{2} \lt 5\% $
, and evaluate its value for U and L, respectively. As shown in table 3, increasing the input resolution from 33 × 33 to 65 × 65 appears to have limited impact on the prediction accuracy, while using the ResNet-18 model significantly increases the prediction performance. When ResNet-18 is used, the accuracy for predicting U is close to perfect. Despite the outliers observed in figure 6(b), the corresponding accuracy A(L) is 92.4 %, suggesting that the ResNet-18 model can provide a reasonable prediction of the length scale. We have also examined the performance of ResNet-18 on local data, which generally shows a robust performance (see details in the Supplementary Material).
Table 3. Summary of training configuration, loss and accuracy. Here, the normalised loss at convergence is calculated by taking the average of their values in the last 100 epochs
Figure 7. Radon transform analysis for selected examples of pressure disturbances of type: (a) ellipse, (b) monohull and (c) catamaran. The left, middle and right columns correspond to the input data (i.e. the normalised surface elevation η*), the Radon space of the data,
$\mathbf {R}[\eta ](\theta ,s)$
and the reconstructed Kelvin wakes,
$\hat{\eta}$
, respectively. The reconstruction,
$\hat{\eta}$
, is computed using
$\mathbf {R}[\eta ](\theta _1,s)$
and
$\mathbf {R}[\eta ](\theta _2,s)$
, where θ
1 and θ
2 correspond to the two local peaks (denoted by the white dots) in the Radon space.
For comparison, an analysis based on the Radon transform (Radon, Reference Radon1986) is performed for three selected examples of pressure disturbances (figure 7). These three examples are chosen such that the rotation angles used in the affine transformation for producing the training data are identical. Note that, for monohull and catamaran types of pressure disturbance, the wakes produced often show an interference pattern (for more details, see the visualisations of the raw surface elevation at https://doi.org/10.6075/J0RN386Z). The mathematical expression for the Radon transform is given by
\begin{align} \mathbf {R}[\eta ](\theta ,s)=\int _{-\infty }^{\infty } \eta (z \sin \theta +s \cos \theta ,-z \cos \theta +s \sin \theta ) \mathrm {d} z, \end{align}
where θ is the angle of the line pattern, and s is the arc length. In practice, the integration limits should be replaced by the size of an image. Similar to the study by Zilman et al. (Reference Zilman, Zapolski and Marom2004), the local peaks at angles θ
1 and θ
2 corresponding to the diverging wake patterns are identified (see the white dots plotted on
$\mathbf {R}[\eta ]$
in figure 7). Additionally, the surface elevation reconstructed from
$\mathbf {R}[\eta ](\theta _1,s)$
and
$\mathbf {R}[\eta ](\theta _2,s)$
reproduces the line patterns of the asymptotic wakes.
However, a direct comparison between the Radon transform and the CNN prediction is challenging for several reasons. Technically, these local peaks must be identified manually, which prevents automatic operations on large datasets. Furthermore, the local peaks in the Radon space shown in figure 7 are less distinct than the clusters identified from the Radon transform of satellite data used by Zilman et al. (Reference Zilman, Zapolski and Marom2004) (an example where the local peaks cannot be identified from the training data is shown in the Supplementary Material). This is likely because the relatively high spatial resolution and small domain size of the training data (compared with satellite data) reduce the signal to noise ratio of the asymptotic line patterns.
4.2 Result interpretation
One caveat of pure data-driven machine learning models is that the training process is a ‘black box’, which poses challenges to the interpretation of the result. To address this issue, we resort to the Shapley value, a concept originally developed for quantifying the contribution of each player in game theory (Shapley, Reference Shapley, Kuhn and Tucker1953) for interpreting the prediction capability of neural network. The classical Shapley value is defined as
\begin{align} \phi _i=\sum _{S \subseteq F \backslash \{i\}} \frac {|S|!(|F|-|S|-1)!}{|F|!}\left [f_{S \cup \{i\}}\left (x_{S \cup \{i\}}\right )-f_S\left (x_S\right )\right ], \end{align}
where f(x) is the model to be evaluated, with x being the input vector, F the (full) set of all features and
$ \phi_{i}$
measures the contribution of the ith component of x. Here, the summation loops over all subsets S of F when the ith component is excluded in the prediction.
Figure 8. The SHAP values computed for selected examples of pressure disturbance of type: (a) catamaran, (b) ellipse, (c) monohull. Here, the left, middle and right columns correspond to the input data (i.e. the normalised surface elevation), SHAP value corresponding to the prediction of U and the prediction of L, respectively.
Directly computing (11) is computationally too demanding, because the number of subsets |S| grows fast with the dimension of the input data. Therefore, we compute the SHAP (SHapley Additive exPlanation) values that are the Shapley values of a conditional expectation function of the original model (Lundberg and Lee, Reference Lundberg and Lee2017). The results for three example inputs are shown in figure 8. The spatial distributions of large SHAP value show a strong correlation with both the diverging and transverse waves in the Kelvin wakes, suggesting that the CNN prediction is based on the wave kinematics identified. As a comparison, the asymptotic wake angle (see figure 2) only identifies the local maximum of surface elevation that is heuristically determined.
4.3 Impact of nonlinearity
To evaluate the impact of nonlinearity on the prediction performance of CNN, we perform nonlinear simulations of Kelvin wakes generated by an elliptical pressure. The parameters of the surface disturbance are
$U/\sqrt {gL}=0.6$
, R/L = 5 and W = 0.5, belonging to the dataset T
e
. Six simulations are performed with varying pressure amplitude such that the characteristic wave steepness
$ \epsilon $
varies between 2.5 × 10−3 (weakly nonlinear) and 0.29 (strongly nonlinear). To demonstrate the nonlinear effect on the wave profiles, we plot the surface elevation along the transverse (y) direction at three locations behind the geometric centre of the surface disturbance (figure 9a). As shown, the broadening in the wake angle with increased nonlinearity agrees well with the findings of Dommermuth and Yue (Reference Dommermuth and Yue1988).
Figure 9. (a) Normalised surface elevation along various transverse locations in the weakly (
$ \epsilon = 0.0025 $
) and strongly (
$ \epsilon = 0.29 $
) nonlinear cases and (b) prediction accuracy as a function of the characteristic wave steepness.
Using the same data augmentation for generating the training data, the raw data in these nonlinear simulations are transformed into a validation dataset, which is then fed into the pretrained CNN model for predicting the velocity and the length scale. In figure 9(d), the prediction accuracy A(U) and A(L) are plotted as function of the characteristic wave steepness
$ \epsilon $
, an indicator of nonlinearity. As
$ \epsilon $
increases, A(U) is largely unchanged. Therefore, to train a neural network model for predicting the velocity of the disturbance, it is sufficient to generate Kelvin-wake patterns from a linear wave model. In contrast, A(L) sees a significant decrease with increased nonlinearity. The sensitivity of A(L) to nonlinearity suggests that the CNN model better captures the wave dynamics compared with the velocity.
5. Conclusions
We have presented a CNN-based framework for predicting ship characteristics. By leveraging a phase-resolved wave model, we have simulated Kelvin wakes produced by surface disturbances of varying types. It is found that the CNN model trained from the wave data shows an excellent performance in predicting the disturbance speed and a satisfactory performance when predicting the disturbance length scale. We have revealed two important features of the prediction performance of CNN in inferring the disturbance properties that were rarely explored in pure data-driven based studies. First, it is found that the SHAP values are strongly correlated with the transverse and diverging wave patterns, which not only substantiates the capability of CNN as a reliable prediction tool but also paves the way for building a further connection between CNN and Kelvin-wake dynamics in future studies. Additionally, we have, for the first time, demonstrated that the nonlinearity impact on the prediction accuracy of the pretrained CNN model varies significantly depending on the prediction quantity.
The present study has demonstrated several advantages of supervised machine learning based on CNNs over unsupervised algorithms. For unsupervised algorithms, such as the Radon transform, an empirical relation is required to quantify the disturbance velocity and length scale from the features extracted from the raw data (see discussions by Zilman et al., Reference Zilman, Zapolski and Marom2004). In contrast, this step is unnecessary for CNNs, as the labelled training data are generated from physics-based simulations. In other words, the CNN model itself serves as a reduced-order model of the underlying wave physics. Additionally, given the universal approximation capability of deep neural networks (Hornik et al., Reference Hornik, Stinchcombe and White1989), the CNN-based prediction framework can be readily adapted to predict other quantities, such as the type of disturbance. From an efficiency perspective, while producing the training data and training a CNN model are computationally expensive, making predictions using a trained ResNet-18 model can be done in real time on modern CPUs and GPUs.
We remark that the present study also sees several promising extensions in future studies. For instance, the performance of the CNN can be further examined against a physics-based data assimilation framework that leverages the HOS method to incorporate the surface wave dynamics (Wang and Pan, Reference Wang and Pan2021; Wang et al., Reference Wang, Zhang, Ma, Zhang, Li and Pan2022; Wu et al., Reference Wu, Hao and Shen2022; Wu et al., Reference Wu, Hao, Li and Shen2023). While we have only considered the Kelvin-wake part of ship wakes in this study, the same approach can be applied to non-Kelvin-wake parts, such as breaking waves, bubbly flows and viscous wakes (Reed and Milgram, Reference Reed and Milgram2002; Spedding, Reference Spedding2014) provided that the training data are generated by turbulent simulations based on the Navier–Stokes equations. Finally, the process for extending this framework to realistic remote sensing data is outlined. The training data (i.e. the synthetic normalised radar cross-section) can be generated by applying a transfer function to the surface elevation data (see review by Rizaev et al., Reference Rizaev, Karakuş, Hogan and Achim2022). A CNN model can then be trained using the synthetic satellite data as input to predict the ship characteristics.
Supplementary material
Supplementary information are available at https://doi.org/10.1017/flo.2025.4.
Acknowledgements
X. H. gratefully acknowledges the anonymous referees for their constructive and valuable comments. The author also thanks L. Lenain and N. Statom for providing the ASV C-Worker 5 ship geometry used in the graphical abstract, and S. Pramod Anumula for producing the graphical abstract.
Data availability
Training data and visualisations of the raw surface elevation are available in UC San Diego’s data archive https://doi.org/10.6075/J0GH9JCX.
Funding
This work was supported by Naval Innovation, Science, and Engineering Center under Office of Naval Research (ONR) grant N000142312831.
Competing interests
The author declares no conflict of interest.
Ethical standards
The research meets all ethical guidelines, including adherence to the legal requirements of the study country.
Open access
