Task-Driven Learning Downsampling Network Based Phase-Resolved Wave Fields Reconstruction with Remote Optical Observations

Abstract

We develop a phase-resolved wave field reconstruction method by the learning-based downsampling network for processing large amounts of inhomogeneous data from non-contact wave optical observations. The Waves Acquisition Stereo System (WASS) extracts dense point clouds from ocean wave snapshots. We couple learning-based downsampling networks with the phase-resolved wave reconstruction algorithm, and the training task is to improve the wave reconstruction completeness ratio $C_R$. The algorithm first achieves initial convergence and task-optimized performance on numerical ocean waves built by the linear wave theory model. Results show that the trained sampling network can lead to a more uniform spatial distribution of sampling points and improve $C_R$ at the observed edge regions far from the optical camera. Finally, we apply our algorithm to a natural ocean wave dataset. The average completeness ratio is improved over 30% at low sampling ratios ($S_R\in [2^{-9},2^{-7}]$) compared to the traditional FPS method and Random sampling method. Moreover, the relative residual between the final reconstructed wave and the natural wave is less than 15%, which provides an efficient tool for wave reconstruction in ocean engineering.

1. Introduction

The rapid and real-time reconstruction of ocean wave elevation has critical applications in marine engineering, such as providing decision support for unmanned ships, optimizing marine equipment, and conducting safety analyses of offshore operations [1,2]. Most research on wave surface reconstruction involves fitting the wave surface by interpolating observed discrete wave points, which is limited to specific application scenarios [3,4,5]. In contrast, phase-resolved modeling offers an instantaneous representation of ocean wave motion and statistical properties related to surface wave dynamics [6]. This method is extensively employed for linear deterministic wave prediction and the wave load calculation [7,8,9,10,11]. Constructing a phase-resolved wave field and inverting the wave dynamic models necessitates high spatial resolution of wave observation points [12,13]. Additionally, the wave shadow effect and the uneven distribution of measurement points hinder the development of inversion algorithms and the accuracy of reconstruction [14,15,16]. To address these challenges, we propose a phase-resolved wave field reconstruction method via the learning-based downsampling network. This framework processes substantial amounts of inhomogeneous data from non-contact wave observations, and provides an optimal sampling strategy for downstream wave inversion tasks.

To reconstruct the phase-resolved wave field, observations must encompass a sufficiently large area and possess a high spatial resolution [17]. Traditional methods, such as wave staffs/buoys and X-band radars, face issues regarding the number of measured points and the spatial resolution, respectively [18,19,20,21]. LiDAR is a powerful remote sensing tool, and has been widely used in the marine, offshore and coastal fields [22,23]. Although LiDAR cameras offer higher spatial resolution and can directly measure wave elevation via the time-of-flight method, they are constrained by wavefront conditions, requiring an aerated and turbulent water surface [24,25]. Stereo imaging technology, employing high-resolution cameras, provides high-precision observations of the sea surface, and is relatively inexpensive and widely applicable [26]. Consequently, we utilize this method for data acquisition in this study. However, when stereo imaging technology operates at a grazing angle, it suffers from a spatially inhomogeneous distribution of measurement points. Since the reconstruction relies on the observation points, an uneven distribution of these points can lead to the loss of wave information, thereby negatively affecting the reconstruction accuracy at the edges of the observed region. Additionally, the extensive volume of point cloud data poses significant challenges for the rapid reconstruction of the phase-resolved wave field directly from the acquired points [27,28].

Point cloud sampling methods can reduce the data size while altering the distribution of observation points to enhance the performance and computational efficiency of reconstruction tasks. Traditional sampling methods include Random Sampling (RS), Grid Sampling, Kd-tree Sampling, and Farthest Point Sampling (FPS) [29,30]. The FPS method is widely used because it effectively preserves the overall features of the point cloud, though it is limited in processing large-scale point cloud data [31]. Recently, many deep learning-based point cloud downsampling methods have been proposed. Geometric deep learning-based point cloud compression employs techniques, such as local feature extraction and local feature aggregation, to achieve point cloud compression [32]. K-means clustering was first used to select representative points and to remove redundant points [33]. Furthermore, a combination of clustering and coarse-to-fine methods has been proposed to achieve fast point cloud simplification [34]. PointNetLK-CPD utilizes the PointNetLK and coherent point drift algorithms to perform unsupervised point cloud compression [35]. Neural Point Cloud Compression extracts point cloud feature information by learning optical flow representations, enabling high-ratio point cloud simplification while preserving geometric features [35]. These deep learning-based point cloud downsampling methods predominantly focus on the geometric structure information of the point cloud, achieving outstanding performance in various point cloud processing tasks. However, they cannot dynamically adjust the sampling strategy according to the requirements of different tasks, thereby affecting the performance of downstream tasks.

Task-driven point cloud downsampling methods have emerged in recent years, dynamically learning sampling strategies based on task requirements to enhance the performance of downstream tasks [36,37,38]. S-Net was the first to propose the idea of learning-based sampling, making the sampled points more adaptable to various tasks and application needs [39]. APSNet is an attention-based point cloud sampling network that focuses on all points by exploiting the history of previously sampled points to find the optimal sampling strategy to optimize the performance of various downstream tasks [40]. LighTN introduced a novel single-head self-correlation module and a sampling loss function, achieving state-of-the-art performance in preserving geometric features across various downstream tasks [41]. Skeleton-Guided Learning proposes a sampling method to preserve the object geometry and topology information during sampling [42,43]. REPS introduced the Global-Local Fusion Attention (GLFA) module, effectively integrating local and global features of the point cloud, thereby ensuring the quality of reconstruction tasks and downsampling [44]. The results of these studies indicate that task-driven point cloud downsampling methods can significantly enhance the performance of downstream tasks.

For wave reconstruction in ocean engineering, embedding physics-based numerical methods within data-driven approaches can enhance the performance of wave model acquisition and wave reconstruction tasks. When predicting waves at specific locations, artificial neural networks (ANNs) serve as effective supplementary tools to physical wave models, boosting their efficiency and reliability. By training ANNs to learn the discrepancies between wave model outputs and measured values, the accuracy of wave models can be improved [45]. Additionally, machine learning can reduce subjectivity and improve modeling accuracy by understanding the sources of uncertainty and quantifying parameter selection in coastal ocean modeling [46]. Bias correction techniques can effectively align wave model outputs with observational data, addressing model inaccuracies [47]. Bagged regression trees (BRT) have been used to predict and describe systematic errors in significant wave height model forecasts, helping to identify areas in the model output space that require refinement [48]. The combination of physical models like Simulating Waves Nearshore (SWAN) with machine learning algorithms (such as BRT and ANN) can be employed to estimate wave parameters in shallow estuaries and identify sources of error in physical models [49]. The physics-informed neural network model NWnets integrates prior knowledge of wave mechanics into soft computing deep learning algorithms, enabling efficient reconstruction of nearshore wave fields with limited measurements [50]. However, these studies neglect the effect of downsampling on wave reconstruction tasks.

Overall, existing research primarily focuses on leveraging networks and models to reduce the inherent errors within the models themselves. However, to the best of the authors’ knowledge, errors stemming from data input have not received adequate attention. Building on the concept of integrating data-driven approaches with physical models from prior studies, this paper proposes a learning-based downsampling strategy to identify the optimal input point cloud for downstream reconstruction tasks. Addressing errors from the perspective of data input minimizes the phase-resolved wave field reconstruction errors and ensures a uniform distribution of the downsampled point cloud. This approach opens up new possibilities for developing novel wave inversion algorithms. This manuscript is organized as follows: In Section 2, the proposed algorithm is described in detail. Section 3 reviews the background of phase-resolved wave field reconstruction, specifically focusing on linear wave theory and the numerical methods for solving wave model parameters. Additionally, we detail the construction of the dataset. Section 4 presents and discusses the test results of the proposed methodology. Conclusions are provided in the final section.

2. Proposed Methodology

To enhance the accuracy of real-time wave field reconstruction based on remote optical measurements, in this study, we developed a learning-based downsampling network. We defined a total loss function that considers task loss and sampling loss as the optimization objective to obtain the optimal input point cloud for downstream reconstruction tasks. The downsampling network comprises a point cloud simplification network and a soft projection operation, enabling sampling at any predetermined ratio. The resulting point cloud satisfies task optimality, similarity to the original point cloud, and distribution uniformity at the specified sampling ratio. The overall optimization framework for the task is illustrated in Figure 1. The dense wave point cloud $P_{n\times 3}$ is obtained by initializing the input wave stereo image pairs through the WASS pipeline [51], which contains 3D position information of n wave observation points, and is used as input data for the network. The data volume of $P_{n\times 3}$ is compressed to m by the simplified network to form a point cloud $Q_{m\times 3}$ with a smaller data volume, and the sampled output point cloud $R_{m\times 3}$ is obtained by soft-projecting $Q_{m\times 3}$ onto $P_{m\times 3}$. This step is used to further ensure the proximity between the simplified point cloud Q and the original point cloud P. This output serves as an input to the wave reconstruction algorithm for the downstream task. The sampling network is trained to obtain the optimal sampling strategy by minimizing the task loss based on the phase-resolved wave reconstruction. After converging, reconstructed sea-waves are exported. Next, we introduce the construction of learning-based sampling network and the algorithm of phase-resolved wave reconstruction.

Our objective is to extract an optimal point cloud subset $R^{*}\in {\mathbb{R}}^{m\times 3}$ from the original point cloud $P\in {\mathbb{R}}^{n\times 3}$ to enhance the performance of the downstream wave reconstruction task T. We define the optimization problem as follows:

$$R^{*}=\underset{R}{argmin}\mathcal{F}\left(T\left(R\right)\right),\mathrm{with}R\subseteq P,m\le n,$$

(1)

where $R\in {\mathbb{R}}^{m\times 3}$ represents the variation in the point cloud during training. To address this optimization problem, we employ a simplified network model based on multilayer perceptrons (MLP). The network comprises two components: the point cloud simplification network and the soft projection mechanism. This approach optimizes the point cloud for the given task while ensuring the sampled output point cloud closely resembles the source point cloud, thereby minimizing the loss of wave features. We aim to generate an optimal R for a given downstream task T from the original point cloud P, ensuring that the points in R are as close and similar as possible to the points in P.

The performance of the simplified network mainly depends on the definition of the loss function expressed as follows:

$${\mathcal{L}}_{total}=\alpha {\mathcal{L}}_{sampling}(R,P)+\beta {\mathcal{L}}_{task}\left(R\right),$$

(2)

where ${\mathcal{L}}_{sampling}$ and ${\mathcal{L}}_{task}$ denote the sampling loss and task loss, respectively, which together constitute the total loss ${\mathcal{L}}_{total}$. Here, $\alpha$ and $\beta$ are weighting factors, with $\alpha$ being set to 0.12 and $\beta$ to 0.88. The loss function of the downsampling network comprises the point cloud simplification network loss ${\mathcal{L}}_{simplify}(Q,P)$ and the projection loss ${\mathcal{L}}_{project}$:

$${\mathcal{L}}_{sampling}(R,P)={\mathcal{L}}_{simplify}(Q,P)+{\mathcal{L}}_{project},$$

(3)

where ${\mathcal{L}}_{simplify}(Q,P)$ represents the feature loss difference between the sampled points and the original point cloud, defined as:

$$\begin{array}{c}\hfill {\mathcal{L}}_{simplify}(Q,P)={\mathcal{L}}_a(Q,P)+\zeta {\mathcal{L}}_m(Q,P)+(\gamma +\delta |Q\left|\right){\mathcal{L}}_a(P,Q),\end{array}$$

(4)

where ${\mathcal{L}}_a(X,Y)$ and ${\mathcal{L}}_m(X,Y)$ are the average nearest neighbor loss and maximal nearest neighbor loss, respectively, to ensure the similarity between Q and P. The average nearest neighbor loss measures the average distance from each point in X to its nearest neighbor in Y, ensuring that points in Q are close to P, which is defined as:

$${\mathcal{L}}_a(X,Y)=\frac{1}{\left|X\right|}\sum _{\mathbf{x}\in X}\underset{\mathbf{y}\in Y}{min}\left|\right|\mathbf{x}-{\mathbf{y}\left|\right|}_2^2,$$

(5)

The maximal nearest neighbor loss measures the maximum distance from each point in X to its nearest neighbor in Y, thereby preventing any point in Q from deviating significantly from its corresponding points in P. It is defined as:

$${\mathcal{L}}_m(X,Y)=\underset{\mathbf{x}\in X}{max}\underset{\mathbf{y}\in Y}{min}\left|\right|\mathbf{x}-{\mathbf{y}\left|\right|}_2^2,$$

(6)

To further ensure the similarity between the sampled point cloud and the original point cloud, a soft projection operation is employed. In this operation, each point $q\in Q$ is “softly” projected onto the original point cloud P to generate a new projected point. The objective of the projection is to find a neighborhood of q in P. The projection point $\mathbf{r}\in R$ of q is computed as a weighted average of the k nearest neighbors in P:

$$\mathbf{r}=\sum _{i\in {\mathcal{N}}_P\left(q\right)}{\omega }_i{\mathbf{p}}_i,$$

(7)

where ${\mathcal{N}}_P\left(q\right)$ denotes the neighborhood of q in p, comprising the k nearest neighbors, and ${\omega }_i$ represents the corresponding weights. The specific calculation of these weights is as follows:

$${\omega }_i=\frac{e^{-d_i^2/{\beta }^2}}{{\sum }_{j\in {\mathcal{N}}_P\left(\mathbf{q}\right)}{e}^{-d_j^2/{\beta }^2}},$$

(8)

where $\beta$ is a learnable shape parameter that controls the shape of this probability distribution, $d_i$ is computed as the Euclidean distance $d_i=\left|\right|\mathbf{q}-{\mathbf{p}}_i{\left|\right|}_2$ between the points q and $p_i$. The weights ${\omega }_i$ can be interpreted as a probability distribution function defined over the neighboring points $p_i$ in P, with the projection point r as the expected value. As $\beta \to 0$, the probability distribution function ${\omega }_i$ will converge to the Kronecker delta function, which is defined as:

$$\omega \left(i\right)=\left\{\begin{array}{cc}1\hfill & i=i_{nearest}\hfill \\ 0\hfill & i\ne i_{nearest}\hfill \end{array}\right..$$

(9)

Therefore, when $\beta \to 0$, the entire probability mass becomes wholly concentrated at the nearest neighbor point $p{i}_{nearest}$, and the probability density at the other points will be 0. This makes the projected point r equivalent to the value of that nearest neighbor point. The projection loss used to optimize the shape parameter $\beta$ is defined as:

$${\mathcal{L}}_{project}={\beta }^2.$$

(10)

For each point r in R, we select the point with maximum projection weight ${\omega }_i$ from the original point cloud P. Similar to the work of Dovrat et al. [39], the unique set of sampled points is taken if multiple r correspond to the same $p_i$. This set is then complemented to m points using the Farthest Point Sampling (FPS) algorithm and, finally, the task performance is evaluated.

Next, we introduce the task loss ${\mathcal{L}}_{task}$, which is defined to minimize the reconstruction loss as:

$${\mathcal{L}}_{task}=1-C_R,$$

(11)

where $C_R$ is the completeness ratio, a metric to evaluate the quality of the 3D reconstruction. It can be computed as follows:

$$C_R=\frac{N_c}{N},\mathrm{with}{N}_c=\sum (\Delta z\le \delta ),$$

(12)

where N is the total number of points in the transformed point cloud. The value of $N_c$ depends on $\Delta z$. For each corresponding point, the absolute value of the residual is calculated as:

$$\Delta z=|z_t-z_r|,$$

(13)

where $z_t$ represents the elevation of the actual sea wave, and $z_r$ is the height of the reconstructed wave by the phase-resolved wave reconstruction algorithm discussed in the next section. The number of coverage points that satisfy the threshold $\delta$ to calculate the coverage rate. The threshold is set to $\Delta \varphi =0.001$.

3. Experimental Design and Construction of Data Sets

3.1. Phase-Resolved Wave Reconstruction Algorithms

With the methodology clearly outlined, we now focus on the experimental design and the processes involved in constructing the datasets for our study. To invert and reconstruct the phase-resolved wave field from the sampling point $R^{m\times 3}$, we employ the linear wave theory (LWT) model. The overall flow of the reconstruction algorithm is shown in Figure 2, where $A_0$ is added as the random initial condition. Then, we calculate the coefficient matrix for the LWT model through iteration, with a convergence criterion $\delta =1\times {10}^{-6}$.

We construct the LWT model by the Eulerian method for modeling rotation-free, noncohesive, incompressible fluids. A Cartesian coordinate system $(x,y,z)$ is considered, with the x-axis and y-axis located at the horizontal water surface and the z-axis along the vertical direction. The linear superposition of $n=1,...,N$ individual wave harmonics of amplitude ${\mathit{A}}_n$, angular frequency ${\omega }_n$, and direction of propagation ${\mathbf{\Theta }}_n$ yields the linear ocean surface representation in the horizontal plane of coordinates $\mathit{r}=(x,y)$:

$${\eta }^{lim}\left(\mathit{r},t\right)=\sum _{n=1}^N{\mathit{A}}_{\mathit{n}}cos\left({\mathit{k}}_n·\mathit{r}-{\omega }_{n}t-{\phi }_n\right),$$

(14)

where t is time, ${\phi }_n=2\pi {\mathcal{R}}_n$ is mutually independent (i.e., random) phases, with ${\mathcal{R}}_n\in [0,1]$ being a set of uniformly distributed random numbers. ${\mathit{k}}_{\mathit{n}}=k_n{\widehat{\mathit{k}}}_{\mathit{n}}=\left(k_{n}cos{\theta }_n,k_{n}sin{\theta }_n\right)$ and $k_n=2\pi /{\lambda }_n=\left|{\mathit{k}}_{\mathit{n}}\right|$ are wave number vectors and wavenumbers. With the deep water dispersion relationship, wavenumbers are found as $k_n={\omega }_n^2/g$, where g is the acceleration of gravity.

To simplify the following mathematical and algorithm developments related to free surface reconstruction, it is more convenient (and numerically accurate) to use the equivalent linear representation:

$${\eta }^{\mathrm{lin}}\left(\mathit{r},t\right)=\sum _{n=1}^N\left(a_{n}cos{\psi }_n+b_{n}sin{\psi }_n\right),$$

(15)

where $\{a_n,b_n;n=1,\dots ,N\}$ are $2N$ wave harmonic parameters describing the ocean surface, with $(a_n,b_n)=\left(A_{n}cos{\phi }_n,A_{n}sin{\phi }_n\right)$, and ${\psi }_n={\mathit{k}}_n·\mathit{r}-{\omega }_{n}t$ are spatio-temporal phases.

We invert the wave parameters based on a linear wave model to perform the sea surface reconstruction. The standard approach to wave reconstruction based on observed data is the variational approach [14], which optimizes the values of $2N$ unknown parameters $(a_n,b_n)$ by minimizing the root-mean-square (RMS) difference between the wave observations and their model representation. Thus, the cost function is defined as

$$F\left(\mathit{p}\right)=\frac{1}{2}\sum _{l=1}^L{(\eta \left(\mathit{p}\right)-{\eta }_l)}^2,$$

(16)

where $\mathit{p}=\left\{a_n,b_n\right\}$$(n=1,\dots ,N)$ is the control vector of $2N$ unknown model parameters, $\eta \left(\mathit{p}\right)$ is the unknown reconstructed sea surface elevation computed by the linear wave model (an given in Equation (14)), and ${\eta }_l$$(l=1,\dots ,L)$ is the point cloud elevation obtained from downsampled observations.

The model parameters are obtained by minimizing the cost function and solving the system of equations:

$$\frac{\partial F}{\partial a_m}=0,\frac{\partial F}{\partial b_m}=0,$$

(17)

where $m\in [1,...,N]$.

Developing these equations yields a linear system of $2N$ equations for $2N$ unknown parameters:

$$\begin{array}{cc}\hfill \sum _{l=1}^L\sum _{n=1}^N& \{a_{n}cos{\psi }_{ml}cos{\psi }_{nl}+b_{n}cos{\psi }_{ml}sin{\psi }_{nl}\}=\sum _{l=1}^L{\eta }_{l}cos{\psi }_{ml},\hfill \\ \hfill \sum _{l=1}^L\sum _{n=1}^N& \{a_{n}sin{\psi }_{ml}cos{\psi }_{nl}+b_{n}sin{\psi }_{ml}sin{\psi }_{nl}\}=\sum _{l=1}^L{\eta }_{l}sin{\psi }_{ml},\hfill \end{array}$$

(18)

where wave harmonic phases are defined as ${\psi }_{ml}={\mathit{k}}_m·{\mathit{r}}_l-{\omega }_m{t}_l$.

Rewrite Equation (18) in matrix form as:

$$A_{mn}{p}_n=B_m,$$

(19)

where $p_n$ is a vector made of the $2N$ unknown parameters $p_n=a_n,p_{N+n}=b_n$, and A is the $2N\times 2N$ matrix given by:

$$\begin{array}{c}A=[A_{mn},A_{m,N+n};A_{N+m,n},A_{N+m,N+n}],A_{mn}={\displaystyle \sum _{\ell =1}^L}cos{\psi }_{n\ell }cos{\psi }_{m\ell },\\ A_{m,N+n}={\displaystyle \sum _{\ell =1}^L}sin{\psi }_{n\ell }cos{\psi }_{m\ell },\\ A_{N+m,n}={\displaystyle \sum _{\ell =1}^L}cos{\psi }_{n\ell }sin{\psi }_{m\ell },A_{N+m,N+n}={\displaystyle \sum _{\ell =1}^L}sin{\psi }_{n\ell }sin{\psi }_{m\ell },\end{array}$$

(20)

and B is the $1\times 2N$ vector given by:

$$\begin{array}{c}B=[B_m,B_{N+m}],B_m={\displaystyle \sum _{\ell =1}^L}{\eta }_{\ell }cos{\psi }_{m\ell },B_{N+m}={\displaystyle \sum _{\ell =1}^L}{\eta }_{\ell }sin{\psi }_{m\ell }.\end{array}$$

(21)

When inverting waves in the observation area, especially short-peak waves, the condition number of the A matrix is often huge, making the inversion problem often pathological. To address this, we apply Tikhonov regularization, which yields consistent solutions by minimizing the approximated error function:

$$min\left({\left|\left|{\mathsf{A}}_{mn}{p}_n-{\mathsf{B}}_n\right|\right|}^2-{\xi }^2{\left|\left|p_n\right|\right|}^2\right),$$

(22)

where $\xi$ is the regularization parameter, which is selected by generalized cross-validation (GCV), which minimizes the error in the post-test unit weights, i.e., maximizes the confidence in the observations from which the GCV function is constructed. When the GCV function achieves the minimum value, the corresponding $\xi$ is the regularization parameter. After solving for the wave parameters, we can reconstruct the phase-resolved wave field $z_r$, which is used to calculate the task loss in Equation (13).

3.2. Constructing Datasets

In this work, the origin point cloud is measured by stereo cameras installed on the ship. To illustrate this process and verify the effectiveness of our framework, we first construct a data set from the numerical sea wave based on the LWT model. Figure 3 shows the sampling process by stereo cameras, which face downward at an angle of 45° along the negative direction of the z-axis. The horizontal aperture angle of the cameras is 70° and the vertical aperture angle is 55°. The stereo camera captures 32,468 pixel pairs in the viewing range. Three-dimensional coordinates of the point cloud are intersection points between the rays and waves’ surface.

We construct random waves in the horizontal coordinate range $(x,y)\in [-128,128]$ meters and the PM wave spectrum. The specific form of the spectrum is as follows:

$$S\left(\omega \right)=\frac{0.78}{{\omega }^5}\exp [-\frac{5}{4}{\left(\frac{0.795}{\omega }\right)}^4]\frac{2}{\pi }{cos}^2\theta ,$$

(23)

where $\omega$ is the frequency of the component wave and $\theta$ is the direction of the component wave. The coordinates are the wave’s three-dimensional coordinates $(x,y,z)$ in meters at the observation point, where z represents the wave elevation. We generate 2000 random wave samples in Matlab, using 1400 for training, 350 for validation, and 350 for testing. A single point cloud consists of 32,468 points.

Next, we obtain dense point cloud data from the stereo observation images of waves using the WASS method. This wave stereo 3D reconstruction pipeline is shown in Figure 4. Sextant’s publicly available YS 01 wave stereo image data were used as the input for the open-source project [52]. The dataset consists of three-dimensional wavelengths observed during the passage of atmospheric fronts, which is applicable to this study since the phenomenon results in a broad directional propagation of wave energy.

The reconstruction of the wave point cloud using the WASS method in this work is divided into three main steps: (i) Feature extraction and matching: The input stereo image is subjected to feature extraction and matching. Following this, the rotation and translation matrices are computed. (ii) Rectification: the extrinsic and intrinsic calibration data are used to rectify the two stereo images. (iii) Disparity map and point cloud generation: After rectification, a dense stereo algorithm obtains the disparity map between the two frames. Preliminary morphological filters are then applied directly to the disparity map. The initial unfiltered 3D point cloud is generated by triangulating the corresponding pixels. Points with outlier depths and those not belonging to the most significant connected component are filtered out.

An initial robust RANSAC estimation of the mean plane is performed [53]. The algorithm fits the model by continuously randomly sampling a subset of the data and then determining the inlier points based on a threshold of the residuals. This process is repeated to obtain the optimal planar model with the maximum number of inlier points:

$$ax+by+cz+d=0,$$

(24)

where $(a,b,c)$ is the normal vector $n_1$ to the plane and d is the intercept to the origin.

The goal is to rotate the initial plane to be parallel to the z-axis (i.e., horizontal) with an average vector of $n_2=[0,0,1]$. We must compute a rotation matrix $\Gamma$ such that $R\times n_1=n_2$. Compute the axis of rotation as the cross product of the two vectors $u=n_1{n}_2$. Compute the angle of rotation $\theta$ as the angle between the two vectors. Use the Rodrigues’ rotation formula to calculate the rotation matrix $\Gamma$:

$$\Gamma =I+sin\left(\theta \right)·{\left[u\right]}_{\times }+(1-cos\left(\theta \right))·{\left[u\right]}_{\times }^2,$$

(25)

where ${\left[u\right]}_{\times }$ is the antisymmetric matrix of vector u. The translation vector is calculated by $T=n_1-n_2$. The final coordinate transformation of the original point cloud P is based on the rotation matrix $\Gamma$ and the translation vector T. The transformed coordinate $P^{\prime }$ can be calculated by the following equation:

$$P^{\prime }=\Gamma ·P+T,$$

(26)

The WASS method acquires a 3D point cloud of ocean waves at the ${10}^5$ level of data, which is too large for a neural network. After experiments, compressing the number of point clouds to one-tenth can guarantee a more than $95\%$ reconstruction integrity ratio. To balance performance and computational efficiency, we used the farthest point sampling (FPS) method to sample the data of a single point cloud up to 32,468 points. After performing the above pre-processing operations on 2000 WASS-extracted 3D point cloud data of actual ocean waves, we took 1400 point cloud samples as a training set, 350 as a validation set, and 350 as a test set.

Our proposed downsampling network is first trained and evaluated using a simulated ocean wave dataset. We use a server with an NVIDIA GeForce RTX 2080 TI GPU for training, with a batch size set to 8, epochs set to 400, bottleneck size to 128, and output points set to 512. We use an Adam optimizer with an initial learning rate of 0.001. To ensure the uniformity of sampling, we compress each point cloud into a unit voxel for sampling and recover the original and sampled point cloud sizes to calculate the task loss when reconstructing the task.

4. Experimental Results and Discussion

With the experimental design and dataset construction methods in place, we first verify the effectiveness of our algorithm by using the simulated wave and then evaluate the performance on actual sea wave cases in this section. The objective is to obtain an optimal sampling strategy with a low sampling ratio and to analyze the underlying mechanisms.

4.1. Training and Evaluation Based on Simulated Wave Datasets

Figure 5a displays the evolution of total loss $L_{total}$ during the training. $L_{total}$ converges with an approximately exponential decay trend as the number of training epochs increases. After 200 training epochs, the losses approximately converge. Total loss decreased from an initial value of 0.37 to 0.05 at $Epoch=400$, which improves the performance of the downstream reconstruction task by 32%. To illustrate the sampling results and reconstruction effects of the network, we select two snapshots during training, that is CaseA ($Epoch=30$) and CaseB ($Epoch=200$) in Figure 5a. We present the one-dimensional wave cross-section form $(x,y)=(0,0)$ to $(x,y)=(0,$$-120)$ (red line shown in Figure 3).

Figure 5b shows observation points are not uniform in spatial distribution, which is denser for closer to the location of the camera and relatively sparse for those far away. Before optimizing the sampling strategy, the distribution of the sampling points was similar to that of the original observation point. Since the wave reconstruction effect of the downstream task is determined by the distribution of sampling points, it can be intuitively deduced from the figure that this sampling result would lead to an increase in reconstruction error with increasing distance from the camera. Figure 5c shows that the optimized sampling points are more uniformly distributed, and no longer densely clustered near the camera but evenly spread across the entire sampling area. The uniformly distributed sampling points can better reproduce the overall wave shape, avoiding the previous situation where the reconstruction quality was higher near the camera and lower in farther areas. The reconstruction effect curve shows that the trained network can predict the wavefront more accurately in the observation area, and the reconstruction error remains stable across the entire sampling area without exhibiting a rapid increase with increasing distance. Figure 6 shows the inversion results after training at different wavelengths and directions. We can reach the same conclusion as shown in Figure 5c, which validates that our method is robust for wave inversion on different wavelengths and angle spreads.

Then, we evaluate the impact of the proposed sampling strategy optimization method on the reconstructed wave contour map. To evolute the performance of our method, we define the spatial residuals between the reconstructed results and origin wave data as:

$$Res=\frac{|z_r-z_t|}{z_r},\mathrm{with}(z_r,z_t)\in (x,y),$$

(27)

where $z_r$ and $z_t$ are the height of reconstructed wave and origin wave fields, respectively. Figure 7a shows the original natural wave contour map as a benchmark comparison. Figure 7b presents the reconstruction results before training the sampling network. The reconstruction results are relatively accurate when close to the observation area (i.e., the location of the camera). However, the reconstruction accuracy decreases sharply with the increasing distance. There are apparent distortions and aberrations in the waveforms far from the observation area. This phenomenon is consistent with our previous analysis of the problem caused by the uneven distribution of sampling points.

In contrast, Figure 7c shows the reconstruction results after training with the optimized sampling strategy. After learning, the network output sampling points became more uniform, significantly improving the quality of the reconstructed contour maps over the entire observation area (both near and far from the observation area). The distribution of the peaks and troughs is more precise, with no severe distortion or flattening. To quantitatively evaluate the errors, the relative errors between the reconstruction results in Figure 7c and the original wave contours in Figure 7a are visualized in Figure 7d. From the error distribution, the reconstruction errors are controlled at a low level in $90\%$ of the observed area, verifying that our method effectively improves the problem of degradation in reconstruction accuracy caused by the uneven distribution of sampling points.

4.2. Training and Evaluation Based on Real Wave Data Set

The above experimental results validate the effectiveness of our proposed deep learning-based adaptive sampling policy network. To address the more complex features and details of the natural wave dataset, we improve and optimize the PointNet-based network architecture. We increase the number of channels/neurons in the convolutional and fully connected layers. For example, the output channels of the first convolutional layer increase from 64 to 128, and the output dimensions of the last fully connected layer increase from 256 to 512. Moreover, a batch normalization layer is added after each convolutional and fully connected layer to accelerate the training convergence and improve the model’s generalization ability. Although we retain the basic encoder–decoder framework of the PointNet, the depth of the network is appropriately adjusted by increasing the number of fully connected layers from 2 to 3. This adjustment enables the network to perform deeper feature extraction and fusion.

Our propose is obtaining a high completeness ratio $C_R$ between the reconstructed wave and real wave by a sampling strategy with as few as sampling points. Then, we define sampling ratio as:

$$S_R={\log }_2\left(\frac{N_s}{N_t}\right),$$

(28)

where $N_t=2^{15}$ is the number of total sampling point from WASS, and $N_s$ is the number we used. We introduce two other methods as a comparison, the FPS method and random sampling algorithms, to evaluate the performance of the present method. As shown in Figure 8, $C_R$ increases with the increases in $S_R$ for all three methods, and the sampling network keeps the leading position followed by the FPS method. For low sampling ratios $S_R\in [$$-9,$$-7]$, FPS and random sampling algorithms have poorer reconstruction results ($C_R<0.5$). However, our sampling network shows excellent ability in these ranges; for example, the completeness ratio can reaches 0.78 more than $37\%$ compared the FPS method ($C_R=0.43$) at $S_R=$ −7. When $S_R=$ −6, our method can reach $C_R=0.91$, which leads the reconstruction task performance by more than $30\%$ compared other two methods. In the higher sampling interval $S_R\in [$$-4,$$-2]$, the performance gap between the three methods narrows. However, our method maintains a slight lead throughout the interval, consistently achieving a 1–2% higher $C_R$ than the traditional methods. This highlights the superior performance of the algorithm in various sampling density scenarios for the downstream reconstruction task.

We chose three different sampling ratios $S_R=$$-8,$−6, and −5 and observed the convergence process of their training. As shown in Figure 9, the total loss in all three cases shows a decaying trend with the number of training epochs. For the case of $S_R=$ −8, the reconstruction task performance improves by 24% at the end of training. However, the number of points is too small to capture the small-scale wave features, leading to only the main wave information being extracted. When increasing the sampling ratios to $S_R=$ −6 and −5, we observe that these two cases share similar convergence curves. The main difference is that the error for $S_R=$ −5 converges slightly faster in the initial training phase. To balance the sampling density and the computational efficiency, we select the case of $S_R=$ −6 for subsequent analyses.

Figure 9 illustrates the spatial distribution of sampling points from the global views (the first row) and the horizontal projections (the second row). To better present the distribution of the point cloud during the training process, we chose one point for every five neighboring points to show the distribution. The distribution of sampling points during the pre-training phase (Phase-I in Figure 9) is notably uneven, as shown in Figure 10a1,a2. Most sampling points are concentrated in the central region of the observation frame and in proximity to the observation viewpoint. Conversely, sampling points are relatively sparse in the peripheral regions, away from the viewpoint along the frame’s edges. This inhomogeneous distribution arises from inherent variations in the density of the original point cloud data across different regions. Figure 10b1,b2. depicts the distribution of sampling points during Phase-II. Compared to Phase-I, the sampled points exhibit a more uniform distribution throughout the observation area, with a reduction in the density of the central region and an increase in the number of points in the edge regions, albeit with some residual sparse areas. Figure 10c1,c2. presents the final distribution of sampling points in the post-training stage (Phase-III). Following comprehensive network optimization, the sampling points demonstrate a uniform and random distribution in three-dimensional space, no longer exhibiting excessive concentration in specific regions. Both central and peripheral areas exhibit balanced densities without any conspicuous blank or overdense regions. The adaptive sampling algorithm’s gradual optimization of the sampling strategy during the training process transforms the initially uneven distribution, influenced by the observation viewpoint, into a final optimized sampling distribution that uniformly and efficiently covers the entire area.

Figure 11a illustrates the spatial residuals $Res$ in Phase-I, where dark regions are apparent in the residual plot, indicating significant discrepancies between the reconstruction results and the actual data. The highest residual values and poorest reconstruction quality are observed in the peripheral regions furthest from the observing system. In these regions, the density of the original point cloud data are low, resulting in sparse sampling points and posing substantial challenges for the subsequent reconstruction process. In contrast, the central region near the observation viewpoint exhibits superior reconstruction quality due to the high density and more uniform distribution of the sampling points. Figure 11b corresponds to the residual map during the mid-training stage. Compared to the initial state, the residual distribution exhibits a significant improvement. Nonetheless, high residual values persist in the edge regions distant from the observation viewpoint, indicating a need for further enhancement of the reconstruction quality. Figure 11c presents the final residual distribution at the post-training stage. Notably, both the central and peripheral areas show an almost complete absence of dark regions in the residual map, with the entire area displaying a light gray hue and residual values being exceptionally close to zero.

This observation indicates that, after sufficient training, our algorithm can accurately reconstruct the wave structure, closely approximating the natural wave structure and effectively mitigating the reconstruction errors induced by the initial uneven distribution of the sampling points. The residuals have been minimized. By comparing the residual plots across different training stages, we observe a substantial improvement in the algorithm’s reconstruction accuracy after undergoing learning and optimization, particularly in the edge areas where reconstruction quality was initially compromised due to the sparse point cloud. Post-training, these areas exhibit significant enhancement, validating the high reconstruction quality and reliability of our algorithm.

5. Conclusions

In this work, we develop a phase-resolved wave field reconstruction method by the learning-based downsampling network to handle ocean wave data obtained by optical systems. To do this, we first extract dense point clouds from ocean wave snapshots using WASS. Then, we couple learning-based downsampling networks with the phase-resolved wave reconstruction algorithm, training them to improve the wave reconstruction completeness ratio $C_R$. Next, we test the performance of our algorithm on numerical ocean waves modeled by the liner wave theory model. Results indicate that the trained sampling network leads a more uniform spatial distribution of sampling points and improves $C_R$ at the observed edge regions far from the camera. Finally, we apply our algorithm to a natural ocean wave dataset. It can improve the averaged completeness ratio over 30% at low sampling ratios ($S_R\in [2^{-9},2^{-7}]$) compared to the traditional FPS method and Random sampling method. Additionally, the relative residual between the final reconstructed wave and the natural wave is less than 15%.

Our algorithm enables fast and real-time phase-resolved wave field reconstruction in partial ocean engineering applications. In the future, this algorithm could play a significant role in short-term deterministic sea wave prediction (DSWP) and real-time analysis of wave loads. Moreover, it paves the way for developing new inversion algorithms, such as discrete Fourier transform, for wave reconstruction. This could further be applied to autonomous vessel trajectory tracking and path optimization, optimization of marine equipment, and safety analysis of maritime operations, providing possibilities for the intelligent development of the ocean engineering field.