Identification of Ship Hydrodynamic Derivatives Based on LS-SVM with Wavelet Threshold Denoising

Abstract

Nowadays, system-based simulation is one of the main methods for ship manoeuvring prediction. Great efforts are usually devoted to the determination of hydrodynamic derivatives as required for the mathematical models used for such methods. System identification methods can be applied to determine hydrodynamic derivatives. The purpose of this work is to present a parameter identification study based on least-squares support-vector machines (LS-SVMs) to obtain hydrodynamic derivatives for an Abkowitz-type model. An approach for constructing training data is used to reduce parameter drift. In addition, wavelet threshold denoising is applied to filter out the noise from the sample data during data pre-processing. Most of the resulting derivatives are very close to the original ones—especially for linear derivatives. Although the errors of high-order derivatives seem large, the final predicted results of the turning circle and zigzag manoeuvres agree pretty well with the reference ones. This indicates that the used methods are effective in obtaining manoeuvring hydrodynamic derivatives.

1. Introduction

With the development of large-scale and professional ships over past decades, ships have become more difficult to control, and accidents may occur. Navigational safety issues related to manoeuvrability have attracted much attention worldwide recently. The International Maritime Organization (IMO) has formulated standards for ship manoeuvrability [1] to improve the quality of ship design and manufacturing. In addition, with the introduction of the concept of green ships, energy-saving and emission reduction have become vital goals. Accurate ship manoeuvring prediction can help control ships in an appropriate manner and reduce energy consumption.

There are usually two methods for evaluating ship manoeuvrability at the initial stage of design. One is to perform free-running tests, which can directly obtain the characteristic parameters including advance, transfer, tactical diameter, overshoot angles, etc. The other is to predict manoeuvring performances based on system simulations. For the latter, the Abkowitz-type model [2] is one of the most widely used ship manoeuvring mathematic models. However, many hydrodynamic derivatives need to be determined in advance. Therefore, it is essential to obtain hydrodynamic derivatives when an Abkowitz-type model is applied to predict ship manoeuvring motion.

Common methods for obtaining hydrodynamic derivatives may include databases, empirical formulae, captive model tests, computational fluid dynamics (CFD) calculations, and system identification (SI). The SI-based approach is an effective and convenient method that is less costly. Moreover, this method can be applied to full-scale trials to avoid scale effects. There have been many successful applications of SI methods for the determination of hydrodynamic derivatives over past decades. Traditional SI methods may include least-squares (LS) estimation [3], improved least-squares [4], extended Kalman filter (EKF) [5], and recursive prediction error (RPE) [6] methods. However, the application of these methods is limited due to some issues. For example, such methods are highly dependent on the selection of initial values and are sensitive to noise. With the development of artificial intelligence, some new SI methods are expected to overcome these defects, such as artificial neural networks (ANN) [7], support vector machines (SVMs), Bayesian optimisation (BO) [8], genetic algorithms and optimisation methods [9]. In this paper, SVM is used to identify hydrodynamic derivatives. SVMs have good generalisation performance and global optimal solutions. Furthermore, they can effectively avoid the problem of dimensionality, and there is no dependence on the initial parameter estimates. In recent years, Luo and Zou identified the hydrodynamic derivatives of the Abkowitz-type model by using a least-squares support-vector machine (LS-SVM) algorithm [10]. Then, particle swarm optimisation (PSO) [11] and the artificial bee colony algorithm (ABC) [12] were applied to find the optimal hyper-parameter in the LS-SVM algorithm. However, the LS-SVM algorithm still has a shortcoming, which is the lack of sparsity of the solution. To overcome this problem, $\epsilon -\mathrm{SVM}$ was applied to adjust the sparsity of the solution and to identify the hydrodynamic derivatives [13]. Besides this, some black-box modelling methods that do not need hydrodynamic derivatives were also applied to predict ship manoeuvring motion [14,15,16].

In the application of an identification method, the accuracy of the identification results is closely related to the training data. However, the training data obtained from experiments often contain noise. Meanwhile, the structure of the high-fidelity manoeuvring model is usually very complex and redundant. Such models are often prone to over-fitting [17]. In previous studies, most scholars directly used noisy data when applying identification methods. For example, Sutulo and Guedes Soares developed a classic genetic algorithm using the Hausdorff metric, and it was validated using simulated responses polluted with white noise [9]. Xue and Liu et al. used Gaussian process regression to identify a dynamic model with input noise [18]. Wang and Zou et al. applied an improved nu-SVM to predict ship manoeuvring motion with input noise [19]. However, few scholars consider applying denoising methods to the processing of training data. Moreover, parameter identification results are usually far from the true values [20]. This phenomenon may be caused by parameter drift. Parameter drift is a common problem when the SI method is applied to obtain model parameters. To diminish parameter drift, some scholars have tried several methods. Hwang applied “parallel processing” to process the training data but did not completely eliminate parameter drift [20]. Shenoi et al. simplified Son and Nomoto’s nonlinear model using sensitivity analysis [21]. Yoon and Rhee reduced parameter drift by modifying the input scenario [22]. Luo and Li reduced the multicollinearity in the input matrix by adding excitation to the training samples [23]. Luo reduced multicollinearity using a difference method and an additional signal method [24]. Most of these methods are mainly from the perspective of reconstructing training samples to reduce parameter drift. Their results indicate that parameter drift cannot be eliminated, but can be significantly reduced.

Recently, wavelet threshold denoising [25] has drawn attention in marine engineering. Wavelet threshold denoising is a theory based on wavelet analysis. As a new time-frequency analysis method, wavelet analysis is known as the “mathematical microscope”, because of its ability to perform multi-resolution analysis, which can focus on any details in the signal for performance of multi-resolution time-frequency analysis. Therefore, wavelet threshold denoising has been widely applied to data pre-processing. Zhang and Zou applied wavelet threshold denoising to filter out the noise in the training data and then used the denoised training data to identify the response model [26]. However, the response model contains only a few parameters, and the prediction accuracy of the response model is not very high.

In this study, an LS-SVM method with wavelet threshold denoising is introduced into the identification process of hydrodynamic derivatives for an Abkowitz-type model. The hydrodynamic derivatives, based on Reynolds-averaged Navier–Stokes (RANS) simulations, are used to obtain the original training data. The wavelet threshold denoising is applied to filter out the noise from the polluted training data during data pre-processing. The hydrodynamic derivatives in the Abkowitz-type model are identified by LS-SVM. A modified regression model is used to reduce parameter drift without adding additional excitation, which usually has to be considered in traditional approaches. Generalisation performance of the identified model is evaluated by predicting manoeuvring motions not included in the training data.

The structure of the paper is organized as follows: Section 2 describes the ship mathematical model. The algorithms of LS-SVM are presented in Section 3. In Section 4, the regression model and the wavelet threshold denoising algorithm are described. Moreover, a simulated experimental example is presented to demonstrate the validity of the identified model. Section 5 summarizes the study with conclusions.

2. Mathematical Model of Ship Manoeuvring Motion

To describe ship motion, two coordinate reference systems were used here, as shown in Figure 1. One was the earth-fixed coordinates $O_e-X_e{Y}_e{Z}_e$; the other was the body-fixed coordinates $O-XYZ$, where $u$, $v$, $r$ were the surge velocity, sway velocity and yaw rate, respectively, $\psi$ was the heading angle, and $\delta$ was the rudder angle.

According to Newton’s law, the equations for surge motion, sway motion and yaw motion can be written as follows:

$$\begin{array}{c}m(\dot{u}-vr-x_G{r}^2)=F_X\hfill \\ \\ m(\dot{v}+ur+x_G\dot{r})=F_Y\hfill \\ \\ I_Z\dot{r}+m{x}_G(\dot{v}+ur)=F_N\hfill \end{array}$$

(1)

where $m$ is ship mass; $x_G$ is the longitudinal coordinate of the ship gravity center in the body-fixed coordinate system; $I_Z$ is the moment of inertia about the Z axis; $F_X$, $F_Y$, $F_N$ are the longitudinal force, transverse force and yaw moment; according to the Abkowitz-type model [2], the hydrodynamic forces and moment on the right side of Equation (1) may be approximated by the following expressions:

$$\begin{array}{c}{F}_X=X_u\Delta u+X_{uu}{(\Delta u)}^2+X_{vv}{v}^2+X_{\dot{u}}\dot{u}+X_{rr}{r^{\prime }}^2+X_{vr}vr+X_{\delta \delta }{\delta }^2+X_{v\delta }v\delta +X_{r\delta }r\delta \hfill \\ \\ F_Y=Y_{v}v+Y_{vv}{v}^2+Y_{vvv}{v}^3+Y_{\dot{v}}\dot{v}+Y_{r}r+Y_{rr}{r}^2+Y_{rrr}{r}^3+Y_{\dot{r}}\dot{r}+Y_{\delta }\delta +Y_{\delta \delta }{\delta }^2+Y_{\delta \delta \delta }{\delta }^3\hfill \\ +Y_{vrr}v{r}^2+Y_{vvr}{v}^{2}r+Y_{u\delta }\Delta u\delta +Y_{v\delta \delta }v{\delta }^2+Y_{vv\delta }{v}^2\delta +Y_{r\delta \delta }r{\delta }^2+Y_{rr\delta }{r}^2\delta \hfill \\ \\ F_N=N_{v}v+N_{vv}{v}^2+N_{vvv}{v}^3+N_{\dot{v}}\dot{v}+N_{r}r+N_{rr}{r}^2+N_{rrr}{r}^3+N_{\dot{r}}\dot{r}+N_{\delta }\delta +N_{\delta \delta }{\delta }^2\hfill \\ +N_{\delta \delta \delta }{\delta }^3+N_{vrr}v{r}^2+N_{vvr}{v}^{2}r+N_{u\delta }\Delta u\delta +N_{v\delta \delta }v{\delta }^2+N_{vv\delta }{v}^2\delta +N_{r\delta \delta }r{\delta }^2+N_{rr\delta }{r}^2\delta \hfill \\ \end{array}$$

(2)

Equation (2) is substituted into Equation (1), and they are rewritten in a dimensionless form. The acceleration terms are then on the left side, and the other terms are on the right side. By doing this, the following equations are obtained:

$$\begin{array}{c}(m^{\prime }-{X_{\dot{u}}}^{\prime }){\dot{u}}^{\prime }={F_1}^{\prime }(u,v,r,\delta )\hfill \\ \\ (m^{\prime }-{Y_{\dot{v}}}^{\prime }){\dot{v}}^{\prime }+(m^{\prime }{x_G}^{\prime }-{Y_{\dot{r}}}^{\prime }){\dot{r}}^{\prime }={F_2}^{\prime }(u,v,r,\delta )\hfill \\ \\ (m^{\prime }{x_G}^{\prime }-{N_{\dot{v}}}^{\prime }){\dot{v}}^{\prime }+({I_Z}^{\prime }-{N_{\dot{r}}}^{\prime }){\dot{r}}^{\prime }={F_3}^{\prime }(u,v,r,\delta )\hfill \end{array}$$

(3)

$$\begin{array}{c}{F_1}^{\prime }={X_u}^{\prime }\Delta u^{\prime }+{X_{uu}}^{\prime }{(\Delta u^{\prime })}^2+{X_{vv}}^{\prime }{v^{\prime }}^2+\left(m^{\prime }{x_G}^{\prime }+{X_{rr}}^{\prime }\right){r^{\prime }}^2+\left(m^{\prime }+{X_{vr}}^{\prime }\right)v^{\prime }{r}^{\prime }+{X_{\delta \delta }}^{\prime }{{\delta }^{\prime }}^2+{X_{v\delta }}^{\prime }{v}^{\prime }{\delta }^{\prime }+{X_{r\delta }}^{\prime }{r}^{\prime }{\delta }^{\prime }\hfill \\ \\ {F_2}^{\prime }=-m^{\prime }{u}^{\prime }{r}^{\prime }+{Y_v}^{\prime }{v}^{\prime }+{Y_{vv}}^{\prime }{v^{\prime }}^2+{Y_{vvv}}^{\prime }{v^{\prime }}^3+{Y_r}^{\prime }{r}^{\prime }+{Y_{rr}}^{\prime }{r^{\prime }}^2+{Y_{rrr}}^{\prime }{r^{\prime }}^3+{Y_{\delta }}^{\prime }{\delta }^{\prime }+{Y_{\delta \delta }}^{\prime }{{\delta }^{\prime }}^2+{Y_{\delta \delta \delta }}^{\prime }{{\delta }^{\prime }}^3\hfill \\ +{Y_{vrr}}^{\prime }{v}^{\prime }{r^{\prime }}^2+{Y_{vvr}}^{\prime }{v^{\prime }}^2{r}^{\prime }+{Y_{u\delta }}^{\prime }(\Delta u{)}^{\prime }{\delta }^{\prime }+{Y_{v\delta \delta }}^{\prime }{v}^{\prime }{{\delta }^{\prime }}^2+{Y_{vv\delta }}^{\prime }{v^{\prime }}^2{\delta }^{\prime }+{Y_{r\delta \delta }}^{\prime }{r}^{\prime }{{\delta }^{\prime }}^2+{Y_{rr\delta }}^{\prime }{r^{\prime }}^2{\delta }^{\prime }\hfill \\ \\ {F_3}^{\prime }=-m^{\prime }{x_G}^{\prime }{u}^{\prime }{r}^{\prime }+{N_v}^{\prime }{v}^{\prime }+{N_{vv}}^{\prime }{v^{\prime }}^2+{N_{vvv}}^{\prime }{v^{\prime }}^3+{N_r}^{\prime }{r}^{\prime }+{N_{rr}}^{\prime }{r^{\prime }}^2+{N_{rrr}}^{\prime }{r^{\prime }}^3+{N_{\delta }}^{\prime }{\delta }^{\prime }+{N_{\delta \delta }}^{\prime }{{\delta }^{\prime }}^2\hfill \\ +{N_{\delta \delta \delta }}^{\prime }{{\delta }^{\prime }}^3+{N_{vrr}}^{\prime }{v}^{\prime }{r^{\prime }}^2+{N_{vvr}}^{\prime }{v^{\prime }}^2{r}^{\prime }+{N_{u\delta }}^{\prime }(\Delta u{)}^{\prime }{\delta }^{\prime }+{N_{v\delta \delta }}^{\prime }{v}^{\prime }{{\delta }^{\prime }}^2+{N_{vv\delta }}^{\prime }{v^{\prime }}^2{\delta }^{\prime }+{N_{r\delta \delta }}^{\prime }{r}^{\prime }{{\delta }^{\prime }}^2+{N_{rr\delta }}^{\prime }{r^{\prime }}^2{\delta }^{\prime }\hfill \end{array}$$

(4)

where ${F_1}^{\prime },{F_2}^{\prime },{F_3}^{\prime }$ are the polynomials related to the state information, as shown in Equation (4), the change of surge speed relative to the initial surge velocity $u_0$ is $\Delta u=u-u_0$.

3. Least-Squares Support-Vector Machine

Least-squares support-vector machine (LS-SVM) is an improved support-vector machine (SVM) based on statistical theory [27]. Traditional SVM needs to solve an inequality constraint, which increases the computing time. LS-SVM can transform the inequality constraints in traditional SVM to equality constraints, and uses the error squared loss function as the experience loss of the training set. This greatly simplifies the calculation process. Therefore, LS-SVM has been successfully applied to the field of system identification. Its specific algorithm is introduced in the following:

Suppose the training sample is given as $T={\left\{x_i,y_i\right\}}_{i=1}^k$, where $k$ is the number of training samples, $x_i\in R^k$ is the input variable, and $y_i\in R^k$ is the output variable. The nonlinear function $\Phi (x)={\left\{\phi (x_i)\right\}}_{i=1}^n$ is selected to map the input space to the feature space. The form of the nonlinear function modeling is:

$$y(x)=w^T\phi \left(x\right)+b$$

(5)

where $w$ is a weight matrix and $b$ is the bias term. According to the risk minimization principle, the constraint problem can be established as:

$$\{\begin{array}{c}\underset{w,b,e}{\min }J(w,e)=\frac{1}{2}{\Vert w\Vert }^2+\frac{\gamma }{2}{\displaystyle \sum _{i=1}^n{e}_i^2}\hfill \\ s.t.y_i=w^T\phi \left(x_i\right)+b+e_{i}i=1,2,\cdots ,n\hfill \end{array}$$

(6)

where ${\Vert w\Vert }^2$ is the regularized part; $\gamma$ is the regularization parameter; $e_i$ is the error between the model and the training sample; and $e$ is the error variable. By introducing the Lagrange multiplier $\alpha$, Equation (6) can be written as:

$$L\left(w,b,e_i,\alpha \right)=\frac{1}{2}{\Vert w\Vert }^2+\frac{\gamma }{2}{\displaystyle \sum _{i=1}^n{e}_i^2-{\displaystyle \sum _{i=1}^n\left[{\alpha }_i\left(w^T\phi \left(x_i\right)+b-y_i\right)\right]}}$$

(7)

According to the Karush–Kuhn–Tucker conditions (KKT) [28], the following formula can be obtained:

$$\{\begin{array}{c}\frac{\partial L}{\partial w}=0\to w={\displaystyle \sum _{i=1}^n{\alpha }_i\cdot \phi \left(x_i\right)}\hfill \\ \frac{\partial L}{\partial b}=0\to {\displaystyle \sum _{i=1}^n{\alpha }_i=0}\hfill \\ \frac{\partial L}{\partial e_i}=0\to \gamma e_i={\alpha }_{i}i=1,\cdots ,n\hfill \\ \frac{\partial L}{\partial {\alpha }_i}=0\to y_i=w^T\phi \left(x_i\right)+b+e_{i}i=1,\cdots ,n\hfill \end{array}$$

(8)

The ${\alpha }_i$ and $b$ in Equation (8) can be obtained using the following solution:

$$\left[\begin{array}{cc}0& s^T\\ s& K(x_i,x_j)+{\gamma }^{-1}I\end{array}\right]\left[\begin{array}{c}b\\ \alpha \end{array}\right]=\left[\begin{array}{c}0\\ y\end{array}\right]$$

(9)

where $\alpha ={\left[{\alpha }_1,{\alpha }_2,\cdots {\alpha }_n\right]}^T$, $y={\left[y_1,y_2,\cdots y_n\right]}^T$, $s={\left[1,1,\cdots ,1\right]}^T$, $I$ is the identity matrix, and $K(x_i,x_j)$ is the kernel function:

$$K\left(x_i,x_j\right)=\langle \phi \left(x_i\right),\phi \left(x_j\right)\rangle$$

(10)

where $\langle \cdot ,\cdot \rangle$ is the inner product of the eigenspace. The common kernel functions are the linear kernel function, polynomial kernel function, RBF kernel function, and sigmoid kernel function. Referring to the work of Luo and Zou [10], the linear kernel function is adopted in this study. The linear kernel is defined as follows:

$$K\left(x_i,x_j\right)={x_i}^T{x}_j$$

(11)

The final function model can be obtained as:

$$f(x)={\displaystyle \sum _{i=1}^n({\alpha }_i\cdot K(x,x_i))+b}$$

(12)

4. Case Study

4.1. Construction of Regression Model

As mentioned above, parameter drift is a significant problem when the regression method is applied to obtain hydrodynamic derivatives. Parameter drift is caused by multicollinearity in the regression model. Multicollinearity means that although the results of the function fitting are perfect, the results of parameter identification will still show significant deviations. Luo and Li reduced multicollinearity by adding additional excitation to the training samples [23]. In this study, a modified regression model was used to identify hydrodynamic derivatives without adding additional excitation. The detailed structure of the regression model is described below.

Equation (4) can be rewritten as:

$$\begin{array}{c}{F_1}^{\prime }-m^{\prime }{x_G}^{\prime }{r^{\prime }}^2-m^{\prime }{v}^{\prime }{r}^{\prime }=X\cdot A(i)\hfill \\ \\ {F_2}^{\prime }+m^{\prime }{u}^{\prime }{r}^{\prime }=Y\cdot B(i)\hfill \\ \\ {F_3}^{\prime }+m^{\prime }{x_G}^{\prime }{u}^{\prime }{r}^{\prime }=N\cdot C(i)\hfill \end{array}$$

(13)

where the hydrodynamic derivatives and speed state variables in Equation (13) are as follows:

$$\begin{array}{c}\begin{array}{c}X={\left[{X_u}^{\prime },{X_{uu}}^{\prime },{X_{vv}}^{\prime },{X_{rr}}^{\prime },{X_{vr}}^{\prime },{X_{\delta \delta }}^{\prime },{X_{v\delta }}^{\prime },{X_{r\delta }}^{\prime }\right]}_{1\times 8}\hfill \\ \\ Y={\left[{Y_v}^{\prime },{Y_{vv}}^{\prime },{Y_{vvv}}^{\prime },{Y_r}^{\prime },{Y_{rr}}^{\prime },{Y_{rrr}}^{\prime },{Y_{\delta }}^{\prime },{Y_{\delta \delta }}^{\prime },{Y_{\delta \delta \delta }}^{\prime },{Y_{vrr}}^{\prime },{Y_{vvr}}^{\prime },{Y_{u\delta }}^{\prime },{Y_{v\delta \delta }}^{\prime },{Y_{vv\delta }}^{\prime },{Y_{r\delta \delta }}^{\prime },{Y_{rr\delta }}^{\prime }\right]}_{1\times 16}\hfill \\ \\ N={\left[{N_v}^{\prime },{N_{vv}}^{\prime },{N_{vvv}}^{\prime },{N_r}^{\prime },{N_{rr}}^{\prime },{N_{rrr}}^{\prime },{N_{\delta }}^{\prime },{N_{\delta \delta }}^{\prime },{N_{\delta \delta \delta }}^{\prime },{N_{vrr}}^{\prime },{N_{vvr}}^{\prime },{N_{u\delta }}^{\prime },{N_{v\delta \delta }}^{\prime },{N_{vv\delta }}^{\prime },{N_{r\delta \delta }}^{\prime },{N_{rr\delta }}^{\prime }\right]}_{1\times 16}\hfill \\ \end{array}\\ \begin{array}{c}\\ A\left(i\right)={\left[\Delta u^{\prime }(i),\Delta {u^{\prime }}^2(i),{v^{\prime }}^2(i),{r^{\prime }}^2(i),v^{\prime }\left(i\right)r^{\prime }\left(i\right),{{\delta }^{\prime }}^2(i),v^{\prime }\left(i\right){\delta }^{\prime }\left(i\right),r^{\prime }\left(i\right){\delta }^{\prime }\left(i\right)\right]}_{1\times 8}^T\hfill \\ \\ B\left(i\right)={\left[\begin{array}{c}{v}^{\prime }\left(i\right),{v^{\prime }}^2\left(i\right),{v^{\prime }}^3\left(i\right),r^{\prime }\left(i\right),{r^{\prime }}^2\left(i\right),{r^{\prime }}^3\left(i\right),{\delta }^{\prime }\left(i\right),{{\delta }^{\prime }}^2\left(i\right),{{\delta }^{\prime }}^3\left(i\right),v^{\prime }\left(i\right){r^{\prime }}^2\left(i\right),\\ {v^{\prime }}^2\left(i\right)r^{\prime }\left(i\right),\Delta u^{\prime }\left(i\right){\delta }^{\prime }\left(i\right),v^{\prime }\left(i\right){{\delta }^{\prime }}^2\left(i\right),{v^{\prime }}^2\left(i\right){\delta }^{\prime }\left(i\right),r^{\prime }\left(i\right){{\delta }^{\prime }}^2\left(i\right),{r^{\prime }}^2\left(i\right){\delta }^{\prime }\left(i\right)\end{array}\right]}_{1\times 16}^T\hfill \\ \\ C\left(i\right)={\left[\begin{array}{c}{v}^{\prime }\left(i\right),{v^{\prime }}^2\left(i\right),{v^{\prime }}^3\left(i\right),r^{\prime }\left(i\right),{r^{\prime }}^2\left(i\right),{r^{\prime }}^3\left(i\right),{\delta }^{\prime }\left(i\right),{{\delta }^{\prime }}^2\left(i\right),{{\delta }^{\prime }}^3\left(i\right),v^{\prime }\left(i\right){r^{\prime }}^2\left(i\right),\\ {v^{\prime }}^2\left(i\right)r^{\prime }\left(i\right),\Delta u^{\prime }\left(i\right){\delta }^{\prime }\left(i\right),v^{\prime }\left(i\right){{\delta }^{\prime }}^2\left(i\right),{v^{\prime }}^2\left(i\right){\delta }^{\prime }\left(i\right),r^{\prime }\left(i\right){{\delta }^{\prime }}^2\left(i\right),{r^{\prime }}^2\left(i\right){\delta }^{\prime }\left(i\right)\end{array}\right]}_{1\times 16}^T\hfill \end{array}\end{array}$$

(14)

where $X,Y,N$ are the unknown terms containing the hydrodynamic derivatives. The structure of the regression model for identification can be written as follows:

Input variables: $\left[A\left(i\right),B\left(i\right),C\left(i\right)\right]$

Output response:

$$\left[\begin{array}{c}\begin{array}{c}\begin{array}{c}\left(m^{\prime }-{X_{\dot{u}}}^{\prime }\right)L\frac{\dot{u}}{u^2(i)}-m^{\prime }{x_G}^{\prime }{L}^2\frac{r^2(i)}{u^2(i)}-m^{\prime }L\frac{v(i)r(i)}{u^2(i)}\\ \end{array}\hfill \\ \left(m^{\prime }-{Y_{\dot{v}}}^{\prime }\right)L\frac{\dot{v}}{u^2(i)}+\left(m^{\prime }{x_G}^{\prime }-{Y_{\dot{r}}}^{\prime }\right)L^2\frac{\dot{r}}{u^2(i)}+m^{\prime }L\frac{u(i)r(i)}{u^2(i)}\hfill \\ \hfill \end{array}\hfill \\ \left(m^{\prime }{x_G}^{\prime }-{N_{\dot{v}}}^{\prime }\right)L\frac{\dot{v}}{u^2(i)}+\left({I_Z}^{\prime }-{N_{\dot{r}}}^{\prime }\right)L^2\frac{\dot{r}}{u^2(i)}+m^{\prime }{x_G}^{\prime }L\frac{u(i)r(i)}{u^2(i)}\hfill \end{array}\right]$$

(15)

$\dot{u}$, $\dot{v}$, $\dot{r}$ can be obtained by discretization using the explicit Euler scheme as follows:

$$\dot{u}=\frac{u(i+1)-u(i)}{\Delta t},\dot{v}=\frac{v(i+1)-v(i)}{\Delta t},\dot{r}=\frac{r(i+1)-r(i)}{\Delta t}$$

(16)

where $\Delta t$ is the time step and $i$ and $i+1$ are two consecutive time steps. Then, the input variables and output response were input into the LS-SVM to verify the validity of the modified regression model in the next section.

4.2. Identification

The second variant of Korea Very Large Crude Carrier (KVLCC2) was selected for the present case study. The main particulars of KVLCC2 are shown in

Table 1. Main particulars of KVLCC2.

Parameters	Values
Length between perpendiculars	320.00 m
Breadth	58.00 m
Design draft	20.80 m
Block coefficient	0.81
$x$ $\mathrm{coordinate} \mathrm{of} \mathrm{CG} (x_G$)	11.136 m
$\mathrm{Inertia} \mathrm{in} \mathrm{yaw} (I_Z$)	$2.00\times {10}^{12} \mathrm{kg}/{\mathrm{m}}^2$
Rudder turning rate	${2.34}^{°}/\mathrm{s}$
Approach speed	7.956 m/s

.

In order to obtain the training data, the hydrodynamic derivatives based on RANS simulations were used to simulate ship free-running motions. To this end, a series of captive model tests for KVLCC2, such as static drift motion, static circle motion, dynamic sway motion, dynamic yaw motion, etc., were simulated by the RANS method to obtain the hydrodynamic forces. Hydrodynamic derivatives were then determined by regression analysis of the RANS results. More details on RANS computations can be found in the work of Yao et al. [29]. A ${20}^{°}{/20}^{°}$ zigzag test was simulated by using the Abkowitz-type model, where the hydrodynamic derivatives were obtained by RANS computations. The sampling interval was 1 s, and the sampling time was 1000 s. The regularization parameter $\gamma ={10}^7$ was chosen. The identification results of the LS-SVM are shown in

Table 2. Comparison of the original and predicted values, without considering noise.

X-Coef.	Original	LS-SVM	Y-Coef.	Original	LS-SVM	N-Coef.	Original	LS-SVM
${X_u}^{\prime }$	−0.0022	−0.0022	${Y_v}^{\prime }$	−0.01902	−0.0190	${N_v}^{\prime }$	−0.007886	−0.00788
${X_{uu}}^{\prime }$	0.0015	0.0015	${Y_{vv}}^{\prime }$	0.000639	0.000637	${N_{vv}}^{\prime }$	−0.000308	−0.000308
${X_{vv}}^{\prime }$	0.00159	0.00159	${Y_{vvv}}^{\prime }$	−0.1287	−0.1232	${N_{vvv}}^{\prime }$	0.00175	0.0028
${X_{rr}}^{\prime }$	0.000338	0.000337	${Y_r}^{\prime }$	0.005719	0.0057	${N_r}^{\prime }$	−0.003701	−0.0037
${X_{vr}}^{\prime }$	0.01391	0.0139	${Y_{rr}}^{\prime }$	−0.000002	−0.000002	${N_{rr}}^{\prime }$	−0.000002	−0.000002
${X_{\delta \delta }}^{\prime }$	−0.00272	−0.0027	${Y_{rrr}}^{\prime }$	−0.000048	0.00061	${N_{rrr}}^{\prime }$	−0.000707	−0.00053
${X_{v\delta }}^{\prime }$	0.001609	0.0016	${Y_{vrr}}^{\prime }$	−0.02429	−0.0203	${N_{vrr}}^{\prime }$	0.003726	0.0047
${X_{r\delta }}^{\prime }$	−0.001034	−0.0010	${Y_{vvr}}^{\prime }$	0.0211	0.0291	${N_{vvr}}^{\prime }$	−0.019	−0.017
			${Y_{\delta }}^{\prime }$	0.00408	0.0041	${N_{\delta }}^{\prime }$	−0.001834	−0.0018
			${Y_{\delta \delta }}^{\prime }$	−0.000114	−0.000114	${N_{\delta \delta }}^{\prime }$	−0.000056	−0.000056
			${Y_{\delta \delta \delta }}^{\prime }$	−0.003059	−0.0031	${N_{\delta \delta \delta }}^{\prime }$	0.001426	0.0014
			${Y_{u\delta }}^{\prime }$	−0.00456	−0.0046	${N_{u\delta }}^{\prime }$	0.00232	0.0023
			${Y_{v\delta \delta }}^{\prime }$	0.00326	0.0029	${N_{v\delta \delta }}^{\prime }$	−0.001504	−0.0016
			${Y_{vv\delta }}^{\prime }$	0.003018	0.0032	${N_{vv\delta }}^{\prime }$	−0.001406	−0.0014
			${Y_{r\delta \delta }}^{\prime }$	−0.002597	−0.0028	${N_{r\delta \delta }}^{\prime }$	0.001191	0.0011
			${Y_{rr\delta }}^{\prime }$	0.000895	0.000833	${N_{rr\delta }}^{\prime }$	−0.000398	−0.000396

. The acceleration derivatives, including ${X_{\dot{u}}}^{\prime }$, ${Y_{\dot{v}}}^{\prime }$, ${Y_{\dot{r}}}^{\prime }$, ${N_{\dot{v}}}^{\prime }$ and ${N_{\dot{r}}}^{\prime }$, were calculated by RANS computations. The values of these derivatives were as follows: ${X_{\dot{u}}}^{\prime }=-0.001135$, ${Y_{\dot{v}}}^{\prime }=-0.014508$, ${Y_{\dot{r}}}^{\prime }=-0.001209$, ${N_{\dot{v}}}^{\prime }=-0.000588$, ${N_{\dot{r}}}^{\prime }=-0.000564$.

The results in Table 2 show that most of the hydrodynamic derivatives were accurately determined. These results indicate that the regression model is valid. However, due to an error of the sensor and various external uncertain disturbances, including wind, waves, and current, the training samples obtained from the actual experiment were noisy compared with the simulated ones. Therefore, for real situations, the results obtained will have a significant deviation when the LS-SVM directly identifies hydrodynamic derivatives. For this, wavelet threshold denoising was employed here to reduce interference in the data. Based on Sutulo and Soares’s method of noise generation [9], Gaussian white noise was added to the sample data to simulate the real situation. The detailed process was as follows:

$${\zeta }_i={\zeta }_{0i}+{\zeta }^{\max }{k}_0{k}_{\zeta }{\zeta }_i$$

(17)

where $\zeta =u,v,r,\delta$,${\zeta }_{0i}$ is the original “clean” reference data, ${\zeta }^{\max }$ is the maximum absolute value of the clean data, $k_0$ is the general reduction factor used to represent different noise levels—which is set to $5\%$ in this study—$k_{\zeta }$ is a specific reduction factor based on the values of different variables, and ${\zeta }_i$ is a random variable with a variance of 0.2 following Gaussian distribution. According to the records of previous experimental data, the noise of the ship surge velocity is significantly less than the rate of yaw and drift angle, and the noise of the rudder angle can be ignored. Therefore, the reduction factor $k_{\zeta }$ was set to 0.2 for the surge velocity and 1.0 for the remaining responses.

4.3. Wavelet Threshold Denoising

Commonly, a signal model with noise can be written as follows:

$$s\left(t\right)=f\left(t\right)+n\left(t\right)$$

(18)

where $s\left(t\right)$ is the signal with noise; $f\left(t\right)$ is the original effective signal; and $n\left(t\right)$ is the noise signal. The purpose of denoising is to reduce $n\left(t\right)$ in $s\left(t\right)$.

The basic idea of wavelet threshold denoising [25] is to transform the signal and gene-rate wavelet coefficients. These wavelet coefficients contain vital signal information. In general, an effective signal is a low-frequency signal, and a noise signal is a high-frequency signal. The effective signal has a certain continuity in the time domain. However, the noise signal is discontinuous in the time domain. Therefore, the wavelet coefficients of effective signals are usually large, whereas the wavelet coefficients of noise signals are generally small. An appropriate threshold is required here. When the wavelet coefficient of the signal is greater than the threshold value, it is retained. When the wavelet coefficient of the signal is less than the threshold, it is set to zero. Figure 2 is the flow chart for wavelet threshold denoising.

The implementation steps of wavelet threshold denoising are as follows:

Step 1: Wavelet decomposition. Wavelet basis function and decomposition levels need to be determined in this step. Wavelet basis function is chosen by observing the initial sign $s\left(t\right)$. Daubechies 4 wavelet [26] was selected in this study. The decomposition level is a vital parameter in the filtering process. On one hand, noise separation is more effective when the decomposition level is larger. On the other hand, excessive decomposition levels can result in the loss of effective signal. The decomposition level is usually determined by pre-tests. In this study, the decomposition level was five after testing.

Step 2: Threshold processing. This step includes the selection of wavelet threshold and threshold processing functions. The heuristic threshold [30] is an ordinary wavelet threshold. The specific algorithm of the heuristic threshold is as follows:

$$\begin{array}{c}f\left(k\right)={\left(sort\left|s\right|\right)}^2\left(k=1,2,\cdots N\right)\hfill \\ \\ r\left(k\right)=[N-2k+{\displaystyle \sum _{i=1}^{k}f\left(i\right)}+\left(N-k\right)f\left(N-k\right)]/N\hfill \\ \\ {\lambda }_1=\sqrt{\min \left(r\right)}{\lambda }_2=\sqrt{2\ln \left(N\right)}\hfill \end{array}$$

(19)

where $s$ is the initial signal; $N$ is the length of the signal; $sort$ means ascending order; $k$ is the moment; $r\left(k\right)$ is the value of risk at that moment; and $\lambda$ is the threshold value. The formula for the heuristic threshold can be written as:

$$\begin{array}{c}eta=\left[{\displaystyle \sum _{k=1}^N{\left|s_k\right|}^2-N}\right]/Ncrit=\sqrt{\frac{1}{N}{\left(\frac{\ln N}{\ln 2}\right)}^3}\hfill \\ \\ \lambda =\{\begin{array}{cc}{\lambda }_2\hfill & eta<crit\hfill \\ \min \left({\lambda }_1,{\lambda }_2\right)\hfill & eta\ge crit\hfill \end{array}\hfill \end{array}$$

(20)

The usual threshold processing functions are the hard threshold processing function and the soft threshold processing function [31]. In this study, the latter is adopted. The specific algorithm is as follows:

$$\stackrel{⌢}{w}=\{\begin{array}{cc}[\mathrm{sgn}\left(w\right)]\left(\left|w\right|-\lambda \right)\hfill & \left|w\right|\ge \lambda \hfill \\ 0\hfill & \left|w\right|<\lambda \hfill \end{array}$$

(21)

where $w$ is the original wavelet coefficient, $\stackrel{⌢}{w}$ is the denoised wavelet coefficient.

Step 3. Wavelet reconstruction. According to the denoised wavelet coefficients obtained above, the inverse transform of wavelet decomposition is applied to reconstruct the denoising signal.

According to Equation (17), polluted ${20}^{°}{/20}^{°}$ zigzag test data can be obtained. The above methods were applied to denoise the polluted data. The comparison of original data, polluted data and denoised data are shown in Figure 3. It can be seen from the enlarged surge velocity comparison graph that the wavelet threshold denoising had a good denoising performance. Following this, the denoised data was brought into the regression model for identification by the LS-SVM. The regularization parameter $\gamma ={10}^4$ was chosen here. The results of identification are shown in

Table 3. Comparison of the original and predicted values after noise reduction.

X-Coef.	Original	Denoising	Y-Coef.	Original	Denoising	N-Coef.	Original	Denoising
${X_u}^{\prime }$	−0.0022	−0.0023	${Y_v}^{\prime }$	−0.01902	−0.0199	${N_v}^{\prime }$	−0.007886	−0.0078
${X_{uu}}^{\prime }$	0.0015	0.0014	${Y_{vv}}^{\prime }$	0.000639	0.000455	${N_{vv}}^{\prime }$	−0.000308	−0.000416
${X_{vv}}^{\prime }$	0.00159	0.0014	${Y_{vvv}}^{\prime }$	−0.1287	−0.0526	${N_{vvv}}^{\prime }$	0.00175	0.0076
${X_{rr}}^{\prime }$	0.000338	0.0003	${Y_r}^{\prime }$	0.005719	0.0053	${N_r}^{\prime }$	−0.003701	−0.0038
${X_{vr}}^{\prime }$	0.01391	0.014	${Y_{rr}}^{\prime }$	−0.000002	0.00004	${N_{rr}}^{\prime }$	−0.000002	0.000022
${X_{\delta \delta }}^{\prime }$	−0.00272	−0.0021	${Y_{rrr}}^{\prime }$	−0.000048	−0.0027	${N_{rrr}}^{\prime }$	−0.000707	0.0012
${X_{v\delta }}^{\prime }$	0.001609	0.0037	${Y_{vrr}}^{\prime }$	−0.02429	−0.0349	${N_{vrr}}^{\prime }$	0.003726	0.0097
${X_{r\delta }}^{\prime }$	−0.001034	0.00009	${Y_{vvr}}^{\prime }$	0.0211	0.0553	${N_{vvr}}^{\prime }$	−0.019	−0.0095
			${Y_{\delta }}^{\prime }$	0.00408	0.0035	${N_{\delta }}^{\prime }$	−0.001834	−0.0017
			${Y_{\delta \delta }}^{\prime }$	−0.000114	−0.00024	${N_{\delta \delta }}^{\prime }$	−0.000056	−0.000047
			${Y_{\delta \delta \delta }}^{\prime }$	−0.003059	0.001	${N_{\delta \delta \delta }}^{\prime }$	0.001426	0.00074
			${Y_{u\delta }}^{\prime }$	−0.00456	−0.0047	${N_{u\delta }}^{\prime }$	0.00232	0.0022
			${Y_{v\delta \delta }}^{\prime }$	0.00326	0.0100	${N_{v\delta \delta }}^{\prime }$	−0.001504	−0.00028
			${Y_{vv\delta }}^{\prime }$	0.003018	0.0204	${N_{vv\delta }}^{\prime }$	−0.001406	−0.001
			${Y_{r\delta \delta }}^{\prime }$	−0.002597	−0.00046	${N_{r\delta \delta }}^{\prime }$	0.001191	0.0023
			${Y_{rr\delta }}^{\prime }$	0.000895	−0.00053	${N_{rr\delta }}^{\prime }$	−0.000398	−0.0002

.

It can be observed from Table 3 that most of the predicted hydrodynamic derivatives were relatively accurate. However, some high-order hydrodynamic derivatives had significant errors. One possible reason is the influence of the wavelet threshold denoising process. As mentioned in Section 4.3, the choice of wavelet decomposition levels and the wavelet threshold will affect the denoising process. Therefore, one possible explanation is that the noise in the data was not completely filtered out, or that the wavelet decomposition levels were too large to filter out part of the effective signal. Another possible reason is parameter drift. As mentioned in Section 4.1, the fundamental reason for parameter drift is multicollinearity in the regression model. Although the phenomenon of parameter drift is not obvious when using clean data for identification, the problem of multicollinearity still exists in the regression model. When polluted data is used for identification, the noise may aggravate the multicollinearity in the model and make the parameter drift more obvious. In addition, high-order hydrodynamic derivatives may be more sensitive. There may not be enough information to identify them in the training data. In general, although the predictions of some high-order hydrodynamic derivatives were not sufficiently accurate, it does not mean that the identified model cannot be used to predict the manoeuvring motion.

4.4. Model Validation

The generalisation performance of the identified model should be checked. This was essential for the identified model in practise, because it was possible for the identified model to perform well for the cases from which training data were obtained, but it performed poorly for other applications. Therefore, the identified model obtained above will be applied to predict other manoeuvring motions, including ${20}^{°}{/10}^{°}$ zigzag motion, ${15}^{°}{/15}^{°}$ zigzag motion, and ${35}^{°}$ turning circle motion. The RANS-based derivatives of KVLCC2 were used to simulate manoeuvring to get the original reference data. The comparison between the predicted results and the original data is shown in Figure 4, Figure 5 and Figure 6. Furthermore, the root mean square error (RMSE) was employed to evaluate the prediction accuracy of the identified model. The RMSE is defined by the formula:

$$RMSE=\sqrt{\frac{1}{n}{\displaystyle \sum _{i=1}^n{\left({\stackrel{\wedge }{y}}_i-y_i\right)}^2}}$$

(22)

where ${\stackrel{\wedge }{y}}_i$ is the predicted value, and $y_i$ is the original data. When the value of RMSE is smaller, it means that the predicted result is more accurate. The RMSEs of ${20}^{°}{/10}^{°}$ zigzag motion, ${15}^{°}{/15}^{°}$ zigzag motion and ${35}^{°}$ turning circle motion are shown in

Table 4. Estimation of prediction accuracy by RMSE.

	${20}^{°}{/10}^{°} \mathbf{Zigzag}$	${15}^{°}{/15}^{°} \mathbf{Zigzag}$	${35}^{°} \mathbf{Turning}$
Surge speed	0.0859	0.0744	0.2013
Sway speed	0.1163	0.1090	0.0527
Yaw rate	0.0462	0.0458	0.0160

.

As shown in Figure 4 and Figure 5, the prediction accuracy of ${15}^{°}{/15}^{°}$ zigzag motion and ${20}^{°}{/10}^{°}$ zigzag motion was generally satisfactory. In a previous study, nu-SVM has been used to identify data containing noise [21]. Comparing the RMSE of ${20}^{°}{/10}^{°}$ zigzag motion with the results of nu-SVM, it can be observed that the predicted results of this study were more accurate. This indicates that wavelet threshold denoising can effectively improve the identification ability of LS-SVM. However, there was still a cumulative deviation in the predicted data—especially in the ${35}^{°}$ turning circle motion, as shown in Figure 6. One possible reason is that the identification results of some high-order hydrodynamic derivatives have larger errors due to parameter drift. When the ship moves at a large rudder angle, the effects of high-order hydrodynamic derivatives also gradually increase in the model. So, there is an accumulation of errors. Future work should further study this regression model or adjust the wavelet threshold denoising algorithm to solve this problem. On the whole, the identified model could successfully predict the manoeuvering motion of the ship, and the prediction accuracy was satisfactory.

5. Conclusions

In this paper, LS-SVM with wavelet threshold denoising was applied to identify the hydrodynamic derivatives for an Abkowitz-type model. To reduce parameter drift, we tried to modify the structure of the regression model. The hydrodynamic derivatives based on RANS simulations were used to simulate the ${20}^{°}/{20}^{°}$ zigzag test to obtain the original training data. Then, the original training data was used to validate the modified regression model. The results show that most of the resulting hydrodynamic derivatives were very close to the original RANS-based ones. This indicates that the modified regression model can reduce parameter drift. However, for real situations, there will be noise in the training data. Directly identifying hydrodynamic derivatives will lead to poor results. Therefore, wavelet threshold denoising was applied to filter out the noise from the sample data during data pre-processing. The denoised data was identified by using the above modified regression model. The ${20}^{°}/{10}^{°}$ zigzag motion, ${15}^{°}/{15}^{°}$ zigzag motion, and ${35}^{°}$ turning circle tests were employed as a supplement to verify the generalisation performance of the model obtained by the identification. The identified results indicate that most of the hydrodynamic derivatives were quite accurate, and the model obtained by identification had a strong generalisation performance.

Although the predicted results so far look acceptable, some problems still need to be addressed in future works:

• More efforts are needed to reduce parameter drift, as the influence of parameter drift on identification accuracy remains considerable.
• The selection of some parameters for wavelet threshold denoising still depends on experience. How to better determine decomposition layer number and wavelet threshold during noise filtering is a key problem to be figured out.