Sea state estimation based on vessel motion responses: Improved smoothness and robustness using B ´ ezier surface and L1 optimization

Floating structures oscillate in waves, where these wave-induced motions may be critical for various marine operations. An important consideration is thereby given to the sea states at the planning and operating stages for an offshore project. The most important information extracted from a sea state is the directional wave spectrum, indicating wave direction, significant wave height, and wave spectrum peak period. Among several available methods of measuring and estimating the directional wave spectrum, the wave buoy analogy technique based on vessel motion responses is an in situ and almost real-time solution without extra costs of devices. If the forms of the wave spectra are not predefined in the estimation, the method is called a nonparametric approach. Its most remarkable advantage is the flexible form, but the smoothness should be regulated. After the discrete Fourier transform has been applied to the measured vessel motions, smoothing is necessary. However, this process results in disturbed vessel cross-spectra and a lowpass characteristic of the windowing function. This paper presents a nonparametric approach for directional wave spectrum estimation based on vessel motion responses. It in- troduces novel smoothness constraints using B ´ ezier surface and includes a more robust estimate using L1 optimization. Both techniques are applied to the wave buoy analogy for the first time. Numerical simulations are conducted to verify the proposed algorithm.


Introduction
Environmental conditions play significant roles in the design and safety of marine operations. Among these conditions, waves are sometimes the most critical factors influencing the dynamics of floating structures [1][2][3]. Although long-term statistics and short-term weather forecasts can supply a large amount of knowledge regarding the operation site, in situ prevailing sea state information is often

Construction of the linear equations between the wave spectrum and response spectra
The translational and rotational motions (surge, sway, heave, roll, pitch, yaw) are indexed respectively by the index set J = {1, 2, 3, 4, 5, 6}. With a total number N d of the selected DOFs, the set of the selected DOFs in the calculation is I⊆ J. The DOFs that are not being actively controlled are the most common choices. Taking a vessel with a speed of advance as an example, the sway, heave, roll, and pitch motions are typically used, i.e., I = {2, 3, 4, 5}. In the proposed stationkeeping scenario, I = {3, 4, 5}. To estimate the wave direction, at least one DOF with asymmetric RAOs about the vessel centerline should be considered. Hence, either sway or roll motion, i.e., i = 2 or i = 4, must be included [6].
In this paper, the wave spectrum and motion response spectra are assumed to be linearly related by RAOs. The construction of the linear equations is based on several assumptions. First, the sea state is assumed to be stationary over the measurement period. In addition, the direction of incoming waves and vessel heading are assumed to be constant to simplify the problem, and the advancing speed is assumed to be zero so that the vessel motion spectra and RAOs are assumed to be consistent over the measurement period. All these assumptions ensure a stationary condition during the data collection period.
There are two variables in the directional wave spectrum E(ω,α), i.e., the frequency ω and incoming wave direction α with respect to North (see Fig. 1). The wave direction is defined as the direction to which the waves are traveling. To discretize the directional wave spectrum, the wave frequency and direction are divided into a network of wave components. The indices m and k are defined to indicate the specific wave component as M = {m|1, ⋯, N ω } and K = {k|1, ⋯, N β }, where N ω denotes the total number of wave frequency components and N β is the total number of wave directional components. When the vessel heading ψ and wave direction α remain constant in a short-term period, the wave heading β is defined as the angle between the longitudinal ship axis and the wave train approaching the vessel, i.e., β k : = α k − ψ. (1) For the sake of simplicity, β is considered hereafter, that is, the wave direction relative to the vessel heading is considered, e.g., β = π denotes head sea.
The directional wave spectrum E(ω, β) is assumed to be the superposition of products of long-crested wave spectra Q κ (ω) and spreading functions N κ (β), i.e., where Κ denotes the index of the wave components. Two components are enough to model most wave spectra, i.e., κ = 1 means the wind sea component and κ = 2 captures the swell component.
In the frequency domain, the cross-spectra for vessel motion response and the directional wave spectrum are related by RAOs. RAOs, which are calculated by software for seakeeping analysis based on strip theory or 3D panel model such as ShipX and WAMIT, are complex-valued linear transfer functions, i.e., Φ(ω, β) = R(Φ(ω, β)) + iI(Φ(ω, β)), where R(Φ(ω, β)) and I(Φ(ω, β)) are the real and imaginary parts of RAO Φ(ω, β), respectively. These transfer functions show both amplitudes and phase information. Note that the applied RAOs often deviate from the truth due to the model simplification and uncertainties from the loading conditions. However, we assume that the RAOs are perfectly known, as common in literature [19,22,23,28]. A cross-spectrum for vessel motion response at a specific frequency ω m is the discretized integral over the wave heading β, i.e., where S ij denotes the cross-spectrum for motion response of the ith and jth DOFs, Φ i denotes the transfer function of the ith DOF, Φ denotes the complex conjugate of Φ, The cutoff frequency ω Nω should be selected where the values of the spectrum are approximately zero, that is, S(ω) ≈ 0, when ω > ω Nω .
A linear equation can be constructed for both the real and the imaginary parts of S ij (ω m ). For i = j, S ii ∈ R is a real number. For i ∕ = j, S ij ∈ C is a complex number. Therefore, there are three parts in Eq. (3), i.e., S ii (ω), R(S ij (ω)), and I(S ij (ω)). For a specific frequency component ω m , there are a total of N 2 d linear equations. Writing in a compact form yields Combining all N ω frequencies in matrix form yields the final affine equation: where For a ship with an advancing speed, Eq. (3) should contain three parts since an encounter frequency may be mapped to three wave frequencies for following sea conditions. However, this only complicates the construction of matrix A in Eq. (7). Hereafter, we focus on solving Eq. (7).

Problem statement
The basis of the directional wave spectrum estimation problem is to reconstruct vector f from measured vector b and known matrix A by solving the linear equation (7). Since N 2 d N ω < N β N ω normally holds, the dimension of vector b is smaller than that of vector f, so the number of unknown variables is larger than the number of known measurements. This indicates an infinite number of solutions to Eq. (7). One should note that the least square (LS) methods, such as the pseudoinverse, are not feasible since the properties of a physical wave spectrum impose implicit constraints on the smoothness and non-negativeness.
To solve Eq. (7), a typical cost function for convex optimization is given by subject to f ≥ 0, (8b) where r 1 and r 2 are the powers. ||x|| p = (|x 1 | p + ⋯ + |x n | p ) 1/p denotes the p-norm of a vector x ∈ R n . In earlier studies, p = 2 is usually selected. A more general form using the p-norm is employed here for the first time. The diagonal elements in the diagonal positive definite matrix C are the weights of the constraints (penalty parameters). The directional wave spectrum is believed to be smooth and continuous. Matrix L contains the smoothness constraints. There are three smoothness conditions, i.e., smoothness of frequency,   and Af in the L p1 -norm sense, and the second term is a L p2 penalty representing the smoothness constraints. Constraints (8c) and (8 d) are optional; alternatively, they can be integrated into matrix L.
The proposed approach faces several challenges in early studies. First, the frequency smoothness and wave heading smoothness constraints are optimized separately. The smoothness condition at a discrete node depends on constant slopes to its two neighbor nodes in frequency-and heading-axes respectively, i.e., The smoothness of the 2D spectrum surface is not considered. Moreover, the effects of farther nodes in both axes are not considered. Calculation of cross-spectra can also be challenging due to the discrete sampling of sensor signals in time domain. The cross-spectra are calculated through discrete Fourier transforms, such as the fast Fourier transform (FFT), of discrete measurements. However, the FFT outputs are not smooth and thereby are not suitable inputs to the following algorithm directly since there exists noise in the resulting cross-spectra. Smoothing or multivariate autoregressive (MAR) modelling must be conducted. They are also subject to disturbances and a lowpass characteristic of the windowing function are left in the resulting cross-spectra.
Measurement noise and the involvement of observers are another problem. The vessel motions are measured by sensors, which also introduce measurement noise and bias. Consequently, sensor noise introduces additional power into all frequencies, especially the high-frequency part, of the long-crested wave spectrum. Because a vessel is not sensitive to high-frequency wave components, the power of high-frequency responses is relatively low. Hence, a small disturbance in the high-frequency part of the response spectrum may result in a relatively large estimation error. Besides, state estimators or observers are adopted to filter the sensor noise and estimate the unmeasured states as some states are expensive or impossible to be directly measured by a sensor, such as velocity. However, the observers cannot fully eliminate the measurement noise and bias caused by the sensors. An observer has the characteristics of a lowpass filter, which eliminates high-frequency information and causes the spectra of the vessel motions to be biased and shifted.
The forth challenge comes from the selection of weights in C, which require tuning and can greatly influence the estimate. A small C results in loose constraints and "noisy" estimates. However, an excessively large C has a negative effect on the fitted spectrum peak, which diminishes the amplitude of the estimated peaks. In practical applications, C is unlikely to change between different scenarios, but we assume that C stays within a range of reasonable values, although the range can be only roughly known. Therefore, a wave spectrum estimation approach providing stable and accurate estimates with different levels of disturbances in the cross-spectra is valuable.
All these disturbances degrade the precision and robustness of the estimate. Taking f as the signal vector and b as the measurement vector, we note that the nonparametric sea state estimation approach is very similar to the signal recovery problem in the realm of signal processing. Since the wave spectrum mainly stays in a directional range and the components in a wide range of wave directions are almost zero, the wave spectrum information b is assumed to be sparse.
The objective of this paper is to improve the smoothness and robustness in the reconstruction of the directional wave spectrum from measured vessel motions. The main contributions are the employment of a Bézier surface and L1 optimization for sea state estimation using vessel motions. We seek the proper ways to (1) improve the smoothness of the estimated spectrum and (2) limit the influences of disturbances in the vessel motion response spectra without focusing too much on fine-tuning the weight parameters in C.

Bézier curve and Bézier surface
A Bézier curve is a 2D curve with excellent linear interpolation properties. The curve between the starting point and endpoint is an affine combination of the positions of a series of control points, and the coefficients are simple polynomial functions. This technique has been widely used in various graphic design programs and computer-aided design programs.
An example of a Bézier curve is illustrated in Fig. 2. The position of a point on a Bézier curve is given by where n is the polynomial degree, γ ∈ {0, 1, ⋯, n} is the index of the control points, P γ denotes the position of the γ-th control point, and (caption on next page) Z. Ren et al. μ ∈ [0, 1] is a variable denoting the relative position between the starting point and endpoint. μ = 0 is the staring point, and μ = 1 is the endpoint [30]. The Bernstein polynomial for the point γ is given by where can be considered as the weights for the specific points on the curve. A control point gives higher weights to the closer points. Hence, the influence from distant control points is limited.
n if the control points are evenly distributed along the x-axis. To ensure local smoothness, the goal is to minimize the difference between P(μ γ ) and P γ . The constant-slope criterion in Ref. [6] is a special application of a Bézier curve with n = 2, i.e., P k− 1 − 2P k + P k+1 = 0. This condition is applied to either the wave frequency ω or the wave heading β individually.

Table 1
The environmental parameters in the sensitivity study.             The smoothness constraints therefore minimize the distances between the control points and each corresponding points on the fitted surface, i.e., When n 1 ≥ 2 and n 2 ≥ 2, the position of a point on the surface (but not on the edges) is influenced by the n 1 n 2 control points. There exist (n 1 n 2 − 4) constraints since the positions of the four corners are fixed.
For the sea state estimation algorithm presented in Section 2, the construction of the selected node is illustrated in Fig. 4. The variables ω and β are discretized, and the resulting discrete sequences are placed as the horizontal and vertical axes. The intersections of the discrete wave frequencies and directions form an N ω × N β mesh with white nodes. The mesh actually represents a side surface of a cylinder since β = 0 is the same as β = 2π. We consider the size of a Bézier surface to be n 1 × n 2 , as shown by the blue, yellow, and red nodes in Fig. 4. The top-left corner (m, k) can be started at k ∈ K and m ∈ M/{N ω − n 1 + 1,⋯, N ω }. No constraint is given to the blue nodes since the corner's positions are fixed. The yellow nodes are on the edges, and the red nodes are on the surface. Different weights c(γ 1 , γ 2 ) > 0 can be given to different types of points. The algorithm is summarized in Fig. 5. For each surface, the values of μ 1 and μ 2 are constant for a specific red node. Since the axes are equally divided, μ 1 = γ 1 n1 and μ 2 = γ 2 n2 . For one specific yellow and red node with the location (γ 1 ,γ 2 ), the corresponding line in matrix L can be calculated. The element with index (m +γ 1 − 1) . Since each point is considered several times on different Bézier surfaces, the weight c should be approximately selected.

L1 optimization
A sparse representation means that a small number of coefficients contribute a large proportion of the total energy. This is very relevant for a wave spectrum distribution. Taking the wave spectrum as an example and rewriting the corresponding vector f into Fig. 19. Wave spectrum estimates based on perfect and disturbed motion spectra (Sea state 6). matrix form, the element in the resulting matrix is set to zero when it is smaller than a small positive constant ε. The sparsity patterns are visualized in Fig. 6 with respect to different values of ε. The black squares in Fig. 6 correspond to the nonzero wave components, which constitute the powerful components in the wave spectrum; i.e., the power of each of these components is large enough to be noticed. The superposition of these wave components produces the wave field. Increasing the threshold ε increases the sparsity in f.
However, we note that increasing the sparsity only filters some less powerful wave components, and the dominating wave components are not seriously affected. Increasing the sparsity reduces the likelihood of the unnecessarily rising parts in the estimates, especially in the high-frequency part. Therefore, estimates can be achieved with greater robustness by increasing the sparsity.
Several optimization approaches are widely applied to enhance the sparsity of recovered signals [31]. The L0 measure, a traditional method, counts the number of the nonzero elements. However, the L0-norm is NP-hard with poor behavior in the presence of noise, and its derivatives are always zero. Hence, the L0-norm is not commonly used in optimization.
To improve the estimation robustness and reduce the influences of disturbance in the motion response spectra, the L1-norm is applied instead of the L2-norm, i.e., p 1 = r 1 = 1 and p 2 = r 2 = 1 in Eq. (8), namely sparse regression.
Sparse regression has several advantages. Foremost among them, L1 optimization is more robust to wildpoints; that is, L1-norm regularization is less sensitive to outliers in disturbed response spectra [24]. In addition, the influences of the imperfect and noisy motion response spectra are better limited. Taking the data-fitting term ||Af − b|| p1 into account, the L1-norm can avoid overfitting problem and increase sparsity.
Considering the smoothness constraints ‖CLf‖ p2 , the L1-norm reduces the effects of the constraints at the local peaks in the wave spectrum, resulting in smaller reduction in the local peaks of the estimated wave spectrum. For both the data-fitting term and the smoothness constraints in (8), the possible choices of the loss function are the Euclidean distance, the mean squared error, and the sum of absolute errors, i.e., (p 1 , r 1 ) ∈ {(2, 1), (2, 2), (1, 1)} and (p 2 , r 2 ) ∈ {(2, 1), (2, 2), (1, 1)}. The best combination of (p 1 , r 1 ) and (p 2 , r 2 ) will be discussed in the next section.
Although L1 optimization is not convex and may pose challenges to the solver, it can be changed into a classic programming problem. There are a number of methods to solve L1-norm regularized problems, such as changing the variables and using quadratic programming or interior point methods [32]. For instance, if (p 1 , r 1 , p 2 , r 2 ) = (2, 2, 1, 1), the problem is known as a least mixed norm (LMN) solution, or basis pursuit denoising [33,34]. The L2-norm can be expanded as ‖Af − b‖ 2 2 The cost function (8a) for the augmented unknown variables with additional inequality constraints becomes where 1 is a column vector of ones. Rewriting (14) into a compact form with , and I is the identity matrix. Since A is symmetric positive semidefinite, problem (15) is convex and can be solved by typical algorithms for convex optimization.

Overview
Numerical simulations are conducted to verify the algorithm. A double-peak spectrum E(ω, β) is adopted. The directional spectrum is a superposition of two spectra, i.e., κ = 2. There are three parameters for each Q κ , i.e., the wave height H s , peak wave frequency ω p = 2π Tp , wave heading β, and coefficient s which controls the sharpness of the corresponding spectrum [35]. The equation of a long-crested wave spectrum is given by where Γ( ·) is the gamma function. We choose λ 1 = λ 10 and λ 2 = λ 20 exp(λ 21 H s,2 ), where λ 10 = 3.0, λ 20 = 1.54, and λ 21 = − 0.062. The spreading function is given by [36].
where s κ is a coefficient controlling the spreading.
A 6360-tonne supply vessel with a length between perpendiculars of 82.8 m, a breadth of 19.2 m, and a draught of 6 m is adopted as an example [37]. The vessel is stabilized at the desired position and heading by a DP system. The ship motion cross-spectra are calculated based on motion RAOs calculated by ShipX [38], as shown in Fig. 7. We select I = {3, 4, 5}, N ω = 30, N β = 20, and n 1 = n 2 = 3.
The following simulated cases are conducted, (1) perfect response spectra, (2) disturbed response spectra, and (3) response spectra filtered by a lowpass filter. In the simulations, different levels of white noises and a lowpass filter are added directly to the cross-spectra calculated in the frequency domain, instead of generating the time series and calculating the cross-spectra from the time-domain measurements. This is because the selection of smoothing techniques couples with the proposed algorithm. It will influence the cross-spectra and consequently affect the final estimation. Hence, it would be difficult to validate the performance of the proposed algorithm quantitatively if we choose a specific smoothing technique.
The environmental parameters are tabulated in Table 1. The significant wave heights are randomly selected between 0 and 2 m. The wave spectral peak periods stay in the range of 6-16 s. The coefficients s 1 and s 2 are selected from 10 to 15 and from 25 to 65, respectively. Mean square errors (MSEs) of the estimation errors are used to analyze the results.
In the simulations, a number of cost functions are adopted by using different values of p 1 , r 1 , p 2 , and r 2 . The cost functions are indexed by (p 1 , r 1 , p 2 , r 2 ); see Table 2.

Simulation result 1: perfect response spectra
The perfectly measured response spectra and the estimation results, e.g., for Sea states 1-3, are presented in Figs. 8-11. The colorbar in each image ranges from 0 to the upper limits of the true wave spectrum. Hereafter, the colors of values larger than the colorbar upper limits remain the same as that of the upper limit (yellow in the figures). After tuning the coefficients C in different cost functions, the estimates from all cost functions agree with the in situ sea state well.
The smoothness constraints of the Bézier surface work well to regulate the wave spectrum. Because of the smoothness constraints, the peak of the spectrum estimate is slightly limited. When decreasing the value of C, the estimated peak approaches the real value. However, the increase is no longer remarkable when C becomes smaller than a certain level.
In [6], the constraint in Eq. (8d) are not considered, resulting in a faster computational speed and a more stable convex optimization problem. However, the lack of nonnegative constraints causes meaningless negative estimates at the low and high frequencies.
We note that the proposed method improves the estimation performance by only slightly slowing down the computational speed. The   computational time difference is not large due to the fast-developed hardware advancement unless the number of variables is extremely high. Using a computer with an Intel Core i7-4790 CPU @ 3.6 GHz, it took up to 12 s to execute the calculation. In the following simulations, the values of elements in matrix C remain the same as those in Simulation 1.

Simulation result 2: response spectra with disturbances
In this section, a series of simulations are conducted to evaluate the robustness of the cost functions using noisy motion cross- spectra. White noise is added to the motion response spectra to simulate the disturbance caused by the FFT. The amplitude of the white noise is a ratio of the corresponding spectral peaks; see Figs. 12, 14, and 16. The noise levels are adopted as 3%, 6%, and 10% of the corresponding spectral peaks. The black and blue nodes denote the perfectly and noisily sampled data, respectively. Since the noise is randomly generated for different noise levels, the estimates of a scenario with a lower noise level are, normally but not definitely, better than those with a higher noise level. The corresponding estimates are presented in Figs. 13, 15, and 17, respectively. Although all the cost functions can estimate the directional wave spectrum with perfect motion cross-spectra, their performance is very different with disturbed spectra. From the results, the estimates deviate by the disturbed cross-spectra. The cost functions giving the most accurate estimates are (p 1 ,r 1 ,p 2 ,r 2 ) ∈ {(1, 1, 1, 1), (2, 1, 1, 1), (2, 2, 2, 2)}. This conclusion holds for all the 20 sea states. More results can be found in Figs. 18-24. The deviations increase with the growing amplitude of disturbances.
The estimate given by the cost function with (p 1 , r 1 , p 2 , r 2 ) = (2, 2, 2, 2) is risky. When using the LS method, an accurate estimate appears only with a minor possibility (Figs. 23 and 24). Whereas most time, the estimates of (p 1 , r 1 , p 2 , r 2 ) = (2, 2, 2, 2) no longer agree with those of the other two cost functions (Figs. 18-22). Very small flaws in the motion spectrum leads to large errors in the estimated sea states. The estimates drift away from the real values, showing the LS method is not robust when addressing disturbances in motion response spectra. In contrast, the cost functions with (p 1 , r 1 , p 2 , r 2 ) ∈ {(1, 1, 1, 1), (2, 1, 1, 1)} still identify the main trend of the wave spectrum.
Although, the performance of cost function (p 1 , r 1 , p 2 , r 2 ) = (2, 2, 2, 2) using disturbed measured motion response spectra can be improved by increasing the value of C, a reduction in the estimated peaks was observed. Hence, increasing the value of C only limits the estimation deviation but cannot guarantee reasonable estimates. For some scenarios, such as the perfect and 3% noise level in Fig. 19, if the amplitude of the minor component is very small, L1 optimization still can estimate these components, while LS method can only estimate one major component.
The MSE of estimate error are listed in Table 4. The cost function (p 1 , r 1 , p 2 , r 2 ) = (1, 1, 1, 1) has a lower deviation than the other cost functions, as expected. To summarize, the L1-norm improves the solution robustness against disturbances in the measured motion response spectra.

Simulation result 3: response spectra filtered by a lowpass filter
Another simulation is conducted to evaluate the influence of the lowpass filter problem caused by observers. The response spectra are filtered by a first-order lowpass filter H(s) = 1 Ts+1 with T = 1.2 s. The filter response spectra and the corresponding estimates are, for example, presented in Figs. [25][26][27]. The lowpass filter only reduces the amplitudes of the response spectra, resulting in a scaled-down estimate of directional wave spectrum, but the effect is limited. To illustrate the performance of the algorithms in a more obvious way, the value of the cutoff frequency 2π/T is selected to be smaller than that in practical measurements purposely.
All three cost functions, to some degrees, manage to estimate of the directional wave spectra; however, the amplitudes of the peaks in the estimates are slightly lower than the real values. Greater reductions are observed near the high-frequency peaks. The amplitude reduction at the spectral peaks is smaller when applying the L1 optimization (p 1 ,r 1 ,p 2 ,r 2 ) = (1, 1, 1, 1). Compared with the LS method (p 1 ,r 1 ,p 2 ,r 2 ) = (2, 2, 2, 2), the estimate of cost function (p 1 , r 1 , p 2 , r 2 ) = (1, 1, 1, 1) shows a better concentration near the peak of each wave component. In summary, the proposed algorithms work stably with the lowpass filter characteristics caused by observers.

Conclusion
In this paper, the Bézier surfaces and L1 optimization are applied to estimate the in situ directional wave spectra using measured vessel motion responses. A stationkeeping scenario is taken as an example. Simulation results show that the Bézier surface guarantees the bidirectional smoothness. The L1-norms in Eq. (8) improves the robustness of the optimization problem, as compared with LSbased L2 optimization. When the L1-norm is applied to both the data-fitting term and the smoothness constraints, accurate estimates are received with good resilience with respect to disturbances in the response spectra and subsequently the lowpass-filtered responses.
The performance of the wave buoy analogy is limited by the accuracy of the selected RAOs, as well as the fact that the vessel acts as a lowpass filter. It is almost impossible to perfectly calculate the RAOs, due to the uncertainties of loading conditions, system linearization, and simplification. Future studies therefore will be focused on reducing the influences of errors in RAOs. In addition, further attention should be paid to the fact that vessel responses, given the size of the vessel, may not be sensitive to the high-frequency tails of wave spectra. In addition, the performance of the proposed method to more general and complex applications, such as vessels with advancing speeds, will be studied in future.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.