*Article* **Identification and Prediction of Ship Maneuvering Motion Based on a Gaussian Process with Uncertainty Propagation**

**Yifan Xue <sup>1</sup> , Yanjun Liu 1,\*, Gang Xue 1,\* and Gang Chen <sup>2</sup>**


**Abstract:** Maritime transport plays a vital role in economic development. To establish a vessel scheduling model, accurate ship maneuvering models should be used to optimize the strategy and maximize the economic benefits. The use of nonparametric modeling techniques to identify ship maneuvering systems has attracted considerable attention. The Gaussian process has high precision and strong generalization ability in fitting nonlinear functions and requires less training data, which is suitable for ship dynamic model identification. Compared with other machine learning methods, the most obvious advantage of the Gaussian process is that it can provide the uncertainty of prediction. However, most studies on ship modeling and prediction do not consider the uncertainty propagation in Gaussian processes. In this paper, a moment-matching-based approach is applied to address the problem. The proposed identification scheme for ship maneuvering systems is verified by container ship simulation data and experimental data from the Workshop on Verification and Validation of Ship Maneuvering Simulation Methods (SIMMAN) database. The results indicate that the identified model is accurate and shows good generalization performance. The uncertainty of ship motion prediction is well considered based on the uncertainty propagation technology.

**Keywords:** system identification; ship maneuvering model; gaussian process; prediction uncertainty

## **1. Introduction**

Maritime transport plays a positive role in promoting the sustainable development of the country's economy [1], and it is also directly related to environmental pollution [2]. Accurate maritime traffic simulators (MTS) can provide an effective basis for the port and route planning and management [3], and can help liner shipping companies arrange vessel schedule efficiently [4,5]. Moreover, an accurate ship maneuvering model is of great practical value for ship trajectory prediction and controller design [6]. With the rapid development of maritime autonomous surface ships (MASSs) [7], autonomous navigation and collision avoidance systems require a more intelligent digital maneuvering model, which can predict the future dynamics of ships and estimate the uncertainty caused by the actions to be performed.

Modeling techniques for ship dynamic models involve parametric modeling and nonparametric modeling. Parametric modeling must define a complete mathematical structure in advance from a physical viewpoint and subsequently estimate the hydrodynamic derivatives through parameter identification techniques. Classic system identification methods are widely used for hydrodynamic parameter identification, such as least square estimation [8], the recursive prediction error (RPE) method [9]. However, the traditional methods are sensitive to noise, and the multicollinearity will significantly affect the identification accuracy [10]. Over the decades, a great number of new methods have been proposed to solve the above problems. Yoon and Rhee used ridge regression to suppress the parameter drift due to multicollinearity [11]. Revestido Herrero and Velasco Gonzalez proposed a

**Citation:** Xue, Y.; Liu, Y.; Xue, G.; Chen, G. Identification and Prediction of Ship Maneuvering Motion Based on a Gaussian Process with Uncertainty Propagation. *J. Mar. Sci. Eng.* **2021**, *9*, 804. https://doi.org/ 10.3390/jmse9080804

Academic Editors: Haitong Xu, Lúcia Moreira and Carlos Guedes Soares

Received: 4 July 2021 Accepted: 24 July 2021 Published: 27 July 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

two-step method based on extended Kalman filtering (EKF) to identify the parameters in the nonlinear model [12]. Sutulo and Guedes Soares adopted genetic algorithms (GA) with Hausdorff metric loss function to reduce the influence of white noise on parameter identification [13]. Least squares support vector machine (LS-SVM), with its good robustness and generalization ability, has been applied to various ship parametric model identification, and has been verified by simulation and experiment [14,15]. Recently, Xu et al. proposed an optimal truncated LS-SVM and validated this method by free-running tests [16,17] and planar motion mechanism tests [18]. The main advantage is that it can be successfully used for big data driven modeling or large-scale training set problems. However, parametric models have some inherent limitations. In the specified parametric framework, the unmodeled dynamics caused by external perturbations and noise [19] will greatly impact the parameter estimation. Moreover, the shapes of various unmanned surface vessels (USVs) are irregular, and the traditional parametric models obtained from classic ship types are not completely matched.

Unlike the parametric model, the nonparametric model does not require any predetermined equation framework constructed by prior knowledge [20]. Nonparametric modeling provides a wealth of techniques to extract information from measurement data, which can be translated into knowledge about hydrodynamic systems [21]. The typical representation of nonparametric modeling methods is neural networks. A recursive neural network (RNN) is first used to fit a maneuvering simulation model for surface ships [22]. Zhang and Zou presented the feed-forward neural network with Chebyshev orthogonal basis function for the black-box modeling of ship maneuvering motion [23]. Wang et al. proposed generalized ellipsoidal basis function fuzzy neural networks to identify the motion dynamics of a large tanker [24]. However, NNs require a considerable amount of training data, and the structure of NNs is difficult to determine. Long short-term memory (LSTM) NNs overcome these shortcomings with the transmission of long-term information and have been successfully used to identify USVs [25] and container ships [26]. The kernel-based method requires less training data and has a lower overfitting risk than the NN [27]. Locally weighted learning (LWL) with modified genetic optimization is presented to identify ship maneuvering systems with full-scale trials [28]. *ν*-SVM is proposed to establish the maneuvering motion model and validated by KVLCC2 ship experimental data [29]. In general, nonparametric modeling alleviates the drawbacks of parametric modeling, i.e., multicollinearity, parameter drifting and unmodeled dynamics.

Recently, the Gaussian process (GP) has drawn attention in nonparametric modeling in marine engineering. GP further strengthens the generalization ability of the kernel method with a priori introduction from a Bayesian perspective. GP is used to identify nonlinear wave forces [30], floating production storage and offloading (FPSO) vessel motion modeling [31], and ship trajectory prediction [32]. Ramire et al. first proposed using a multioutput GP to identify the dynamic model of a container ship [33]. Xue et al. presented a noisy input GP to improve the identification accuracy and verified it by using simulated ship motion data with artificial noise [34]. The experimental data of the KVLCC2 ship were used to construct the GP [35], but the accuracy of prediction in the experiment was not sufficiently high. In the prediction of ship motion based on GP, the prediction output of each time is used as the input to the next iteration, so uncertainty will accumulate. However, this ship dynamic modeling using GP does not consider the propagation of variance.

In this paper, to solve the problem of variance propagation in GP, an approximation method is applied. First, the input is assumed to follow a Gaussian distribution. The predictive distribution is approximated by a moment matching-based technique. To evaluate the effectiveness of the proposed scheme, the simulation case of a container ship and the experimental case of a KVLCC2 ship model from the Hamburg Ship Model Basin (HVSA) are taken as the study object. The identified models are assessed by the prediction error with other motion data not included in the training set.

The remainder of the paper is organized as follows. Section 2 describes the nonparametric ship dynamic model. The algorithms of GP with uncertain input are depicted in

Section 3. In Section 4, the identification scheme of the ship and experimental example are presented to demonstrate the applicability of the proposed method. Section 5 summarizes the study with conclusions. Section 3. In Section 4, the identification scheme of the ship and experimental example are presented to demonstrate the applicability of the proposed method. Section 5 summarizes the study with conclusions. metric ship dynamic model. The algorithms of GP with uncertain input are depicted in Section 3. In Section 4, the identification scheme of the ship and experimental example are presented to demonstrate the applicability of the proposed method. Section 5 summarizes the study with conclusions.

are taken as the study object. The identified models are assessed by the prediction error

The remainder of the paper is organized as follows. Section 2 describes the nonpara-

are taken as the study object. The identified models are assessed by the prediction error

The remainder of the paper is organized as follows. Section 2 describes the nonparametric ship dynamic model. The algorithms of GP with uncertain input are depicted in

#### **2. Ship Nonparametric Dynamic Model 2. Ship Nonparametric Dynamic Model**

heading angle.

*J. Mar. Sci. Eng.* **2021**, *9*, x FOR PEER REVIEW 3 of 15

with other motion data not included in the training set.

with other motion data not included in the training set.

*J. Mar. Sci. Eng.* **2021**, *9*, x FOR PEER REVIEW 3 of 15

For a surface ship, the dynamic model is usually described by a 3-DOF model, including the motion of surge, sway and yaw. Figure 1 shows the coordinate system of a surface ship maneuvering motion, including the Earth-fixed coordinates *O* − *X*0*Y*<sup>0</sup> and body-fixed coordinates *o* − *x*0*y*0. Here, *u*, *v*, and *r* are the state variables of surge velocity, sway velocity, and yaw rate, respectively, while *δ* is the rudder angle and *ψ* is the heading angle. For a surface ship, the dynamic model is usually described by a 3-DOF model, including the motion of surge, sway and yaw. Figure 1 shows the coordinate system of a surface ship maneuvering motion, including the Earth-fixed coordinates ܱ−ܻܺ and body-fixed coordinates −ݔݕ. Here, ݑ, ݒ, and ݎ are the state variables of surge velocity, sway velocity, and yaw rate, respectively, while ߜ is the rudder angle and ߰ is the heading angle. **2. Ship Nonparametric Dynamic Model**  For a surface ship, the dynamic model is usually described by a 3-DOF model, including the motion of surge, sway and yaw. Figure 1 shows the coordinate system of a surface ship maneuvering motion, including the Earth-fixed coordinates ܱ−ܻܺ and body-fixed coordinates −ݔݕ. Here, ݑ, ݒ, and ݎ are the state variables of surge velocity, sway velocity, and yaw rate, respectively, while ߜ is the rudder angle and ߰ is the

**Figure 1.** Reference frames for ships. **Figure 1.** Reference frames for ships. **Figure 1.** Reference frames for ships.

The ship maneuvering system is a nonlinear autoregressive model with an exogenous input (NARX) system [36], and the outputs at the next moment are based on the previous state variables. Figure 2 illustrates the modeling and prediction process of the ship dynamic model. In the first stage, ship motion data are collected by onboard sensors such as IMU and GPS. After data preprocessing, the machine learning technique is used to fit the surrogate time series model. Finally, other motions can be predicted through the learned model. The symbol "~" represents random variables. The ship maneuvering system is a nonlinear autoregressive model with an exogenous input (NARX) system [36], and the outputs at the next moment are based on the previous state variables. Figure 2 illustrates the modeling and prediction process of the ship dynamic model. In the first stage, ship motion data are collected by onboard sensors such as IMU and GPS. After data preprocessing, the machine learning technique is used to fit the surrogate time series model. Finally, other motions can be predicted through the learned model. The symbol "~" represents random variables. The ship maneuvering system is a nonlinear autoregressive model with an exogenous input (NARX) system [36], and the outputs at the next moment are based on the previous state variables. Figure 2 illustrates the modeling and prediction process of the ship dynamic model. In the first stage, ship motion data are collected by onboard sensors such as IMU and GPS. After data preprocessing, the machine learning technique is used to fit the surrogate time series model. Finally, other motions can be predicted through the learned model. The symbol "~" represents random variables.

**Figure 2.** Flowchart of the ship maneuvering identification and prediction using GP (speed output). **Figure 2.** Flowchart of the ship maneuvering identification and prediction using GP (speed output).

**Figure 2.** Flowchart of the ship maneuvering identification and prediction using GP (speed output).

According to the relevant studies of nonparametric ship dynamic modeling [25,34], the formulation of the ship discrete nonparametric model is as follows:

$$\begin{array}{l} u(k+1) = GP\_1(u(k), v(k), r(k), \delta(k)) \\ v(k+1) = GP\_2(u(k), v(k), r(k), \delta(k)) \\ r(k+1) = GP\_3(u(k), v(k), r(k), \delta(k)) \end{array} \tag{1}$$

The selected regressors of the GP are inspired by parametric models, including the Abkowitz [37] and Maneuvering Modeling Group (MMG) models [38]. The ship position variables can be obtained as follows:

$$\begin{aligned} \dot{x} &= u\cos(\psi) - v\sin(\psi) \\ \dot{y} &= u\sin(\psi) + v\cos(\psi) \end{aligned} \tag{2}$$

#### **3. Gaussian Process Regression Framework**

*3.1. Gaussian Process with Deterministic Input*

The following notation with a set of training data is defined:

$$\mathcal{D} = \left( [\mathfrak{x}\_t]\_{t=1\prime}^n, [y\_t]\_{t=1}^N \right) \tag{3}$$

where *x<sup>t</sup>* is an input vector, and output *y<sup>t</sup>* is given by

$$\begin{array}{l} y = f(\mathbf{x}) + \omega \\ w \sim \mathcal{N}(\mathbf{0}, \sigma\_{\omega}^{-2}) \end{array} \tag{4}$$

A standard GP is a collection of random variables and can be considered a collection of random variables with a joint Gaussian distribution for any finite subject [39]. GP is specified by the mean function *m*(*x*) and covariance function *k*(*x*, *x* 0 ) as

$$m(\mathbf{x}) = E[f(\mathbf{x})] \tag{5}$$

$$k(\mathbf{x}, \mathbf{x}') = \mathbb{E}\left[ (f(\mathbf{x}) - m(\mathbf{x})) \left( f(\mathbf{x}') - m(\mathbf{x}') \right) \right] \tag{6}$$

where *E* is the expectation operator. Then, the GP can be written as:

$$f(\mathbf{x}) \sim \text{GP}\left(\mathbf{m}(\mathbf{x}), k(\mathbf{x}, \mathbf{x}')\right) \tag{7}$$

The proposed model adopts the commonly used squared exponential covariance function:

$$k(\mathbf{x}\_i, \mathbf{x}\_j) = \sigma\_f^2 \exp(-\frac{1}{2}(\mathbf{x}\_i - \mathbf{x}\_j)^T \Lambda (\mathbf{x}\_i - \mathbf{x}\_j)) \tag{8}$$

where σ*<sup>f</sup>* and Λ are the amplitude and squared length-scale hypermeters, respectively.

Bayesian inference can be defined as the process of fitting a posterior probability model from a prior model with a set of training data D. The GP prior is given as:

$$p(f|X) = \mathcal{N}(m(X), k(X, X))\tag{9}$$

With these modeling assumptions in place, the likelihood function can be obtained,

$$p(y|f,X) = \prod\_{t=1}^{n} \mathcal{N}\left(y\_t; f\_{t\prime} \sigma\_\omega^2\right) \tag{10}$$

Then, combining the prior Equation (9) and the likelihood function Equation (10) with the Bayesian rule, the posterior probability distribution and predicted function values *f* ∗ can be calculated, at a set of deterministic inputs *X* ∗ .

$$
\begin{bmatrix} f^\* \\ y \end{bmatrix} \sim \mathcal{N} \left( \begin{bmatrix} m(X^\*) \\ m(X) \end{bmatrix}, \begin{bmatrix} K(X^\*, X^\*) & K(X^\*, X) \\ K(X, X^\*) & K + \sigma\_\omega^2 I \end{bmatrix} \right) \tag{11}
$$

which leads to the RGP regression predictive equations,

$$pf\*X\*\_{\prime}X\_{\prime}y = \mathcal{N}(\mathbb{Q}\_{\prime}f)\tag{12}$$

where the predictive mean and variance function are specified by:

$$m = m(X^\*) + K(X^\*, X)[K(X, X) + \sigma\_\omega^{-2}I]^{-1}(y - m(X))\tag{13}$$

$$\mathbf{s} = k(X^\*, X^\*) - \mathbf{K}(X^\*, \mathbf{X}) \left[ \mathbf{K}(\mathbf{X}, \mathbf{X}) + \sigma\_\omega^2 I \right]^{-1} \mathbf{K}(\mathbf{X}, \mathbf{X}^\*) \tag{14}$$

The hyperparameters in Equation (8) are usually obtained by maximizing the log of the marginal likelihood function. It is defined as [36]:

$$-\log p(y|X,\theta) = \frac{1}{2}(y - m(X))^T \left(\mathbf{K}(X,X) + \sigma\_\omega^2 I\right)^{-1} (y - m(X)) + \frac{1}{2}\log\left|\mathbf{K}(X,X) + \sigma\_\omega^2 I\right| + \frac{N}{2}\log 2\pi \tag{15}$$

This nonlinear and nonconvex optimization problem is usually solved by the gradient ascent-based methods, such as BFGS and the conjugate gradient (CG) algorithm [36].

#### *3.2. Prediction with Uncertain Inputs and Uncertainty Propagation*

In Equation (11), we assume that the input is deterministic, while the output is Gaussian distributed. This assumption is true in the one-step prediction. For multiple-step predictions, the traditional method is to recycle the one-step prediction. However, the uncertainties induced by each successive prediction cannot be ignored in the time-series tasks [40].

The uncertainty propagation problem can be addressed assuming that the input follows a Gaussian distribution [41]:

$$
\widetilde{X}^\* \sim \mathcal{N}\left(\widetilde{\mu}, \widetilde{\sum}\right) \tag{16}
$$

For convenience of expression and marking, the input variables are divided into speed variable *x* and control variable *u*, as *X* ∗ = [*x*, *u*]. The mean and variance are given as:

$$\widetilde{\Sigma} = \begin{bmatrix} \widetilde{\mu} = \begin{bmatrix} \mu, \operatorname{E}[u] \end{bmatrix}^T \\ \hline \operatorname{Cov}[u, \mathfrak{x}] & \operatorname{Var}[\mathfrak{x}, u] \\ \hline \operatorname{Cov}[u, \mathfrak{x}] & \operatorname{Var}[u] \end{bmatrix} \tag{17}$$

where *E*[*u*] and *Var*[*u*] are the mean and variance of the control variable; *Cov*[*x*, *u*] = *E*[*xu*] − *µE*[*u*].

The predictive distribution can be obtained by integrating over the input:

$$p(X^\* \middle| \widetilde{\mu\_\prime} \widetilde{\sum}) = \int p(f(\widetilde{X^\*}) \middle| \widetilde{X^\*}) p(\widetilde{X^\*} \middle| \widetilde{\mu\_\prime} \widetilde{\sum}) d\widetilde{X^\*} \tag{18}$$

However, this integration cannot be analytically computed, since the Gaussian distribution is mapped through a nonlinear function. Taylor approximation [40] or moment matching [42] is commonly used to approximate the integration. The computation of the Taylor approximation is expensive because it accounts for the gradient of the posterior mean and variance of the input. In this paper, moment matching is chosen. The moment matching method assumes that the unknown distribution only has two parameters: the

mean and the variance. The mean and variance at an uncertain input can be computed by the laws of iterated expectations and conditional variances [42]:

$$m\left(\widetilde{X^\*}^\*\right) = E\_{\widetilde{X^\*}}\left[E\left[X^\*\right]\right] \tag{19}$$

$$
\sigma^2 \left( \widetilde{X^\*} \right) = E\_{\widetilde{X^\*}} [Var[X^\*]] + Var\_{\widetilde{X^\*}} [E[X^\*]] \tag{20}
$$

Then, the predicted mean and variance at time t + 1 are given as:

$$
\mu\_{t+1} = \mu\_t + m\left(\widetilde{X}^\*\right) \tag{21}
$$

$$
\sum\_{t+1} = \sum\_{t} + \text{Cov}[\mathbf{x}\_t, y\_t] + \text{Cov}[y\_t, \mathbf{x}\_t] \tag{22}
$$

### **4. Case Study**

*4.1. Simulation Study on a Container Ship*

The first case uses the simulation maneuvers of a large container ship. The selected parametric numerical model is a nonlinear 4-degree-of-freedom (DOF) dynamic model [43]. The main particulars of the container ship are listed in Table 1. The model has been verified by experiments, which can well reflect the complex dynamic characteristics of the container ship. It is widely used in the testing of system identification algorithms [33,44].

**Table 1.** Particulars of the container ship.


The parametric maneuvering model is given as follows [43]:

$$\begin{cases} \left(m - X\_{\stackrel{\cdot}{u}}\right)\dot{u} - \left(m - Y\_{\stackrel{\cdot}{v}}\right)vr &= F\_{\mathcal{X}}\\ \left(m - Y\_{\stackrel{\cdot}{v}}\right)\dot{v} + \left(m - X\_{\stackrel{\cdot}{u}}\right)ur - Y\_{\stackrel{\cdot}{r}}\dot{r} = F\_{\mathcal{Y}}\\ \left(I\_{\mathcal{X}} - K\_{\stackrel{\cdot}{p}}\right)\dot{p} = F\_{\mathcal{K}} - WGM\varphi\\ \left(I\_{\mathcal{Z}} - N\_{\stackrel{\cdot}{r}}\right)\dot{r} + N\_{\stackrel{\cdot}{p}}\dot{v} = F\_{\mathcal{N}} \end{cases} \tag{23}$$

where *m* denotes the ship mass, *W* is the weight of the displaced water, *GM* is the metacenter height. *I<sup>x</sup>* and *I<sup>z</sup>* denote the moments of inertia of the ship about the x0, z<sup>0</sup> axes. *X* . *u* , *Y*. *v* , *Y*. *r* , *N*. *v* , and *N*. *r* are acceleration derivatives which can be determined using potential theory. *F<sup>X</sup> F<sup>Y</sup> FK*, and *F<sup>N</sup>* are the forces and moment disturbing quantity at x0-axis, y<sup>0</sup> -axis, and z0-axis, respectively. The nonlinear forces and moments are composed of Taylor expansion of hydrodynamic coefficient and speed.

With the 4-DOF model, 2 groups of maneuvers, including 10◦/10◦ and 20◦/20◦ zigzag tests, are undertaken under the following initial conditions: *u*<sup>0</sup> = 7 m/s, *v*<sup>0</sup> = 0 m/s, *r*<sup>0</sup> = 7 m/s, *δ*<sup>0</sup> = 0 rad and the propeller velocity is fixed at 70 rpm. Each maneuver lasted for 850 s, and the simulation interval was 0.1 s. The timeseries of the yaw velocity and rudder angle are shown in Figure 3.

**Figure 3.** 10°/10° and 20°20° zigzag tests for training. **Figure 3.** 10◦/10◦ and 20◦/20◦ zigzag tests for training.

With 2 s as the time interval of the training data, 850 points of the above training data are used to train the GP hyperparameters with maximum likelihood. All calculations are performed in MATLAB R2020a with 4.0 GHz CPU and 16 GB RAM. PILCO toolbox [45] with BFGS and CG algorithms is used to train the Gaussian process. The training process took 9 s in total and trained 3 GPs in Equation (1). The optimization parameter settings of each GP are listed in Table 2. With 2 s as the time interval of the training data, 850 points of the above training data are used to train the GP hyperparameters with maximum likelihood. All calculations are performed in MATLAB R2020a with 4.0 GHz CPU and 16 GB RAM. PILCO toolbox [45] with BFGS and CG algorithms is used to train the Gaussian process. The training process took 9 s in total and trained 3 GPs in Equation (1). The optimization parameter settings of each GP are listed in Table 2.

**Table 2.** Selection of the GP parameters for the container ship. **Table 2.** Selection of the GP parameters for the container ship.


Generalization verification is necessary for system identification. The ability of the identified model to predict other motions not included in the training data is called generalization. The generalization performance of the trained model is verified by predicting the motions, including the 15°/15° zigzag maneuver and port 30° turning circle test. The prediction results of the 15°/15° zigzag test are shown in Figure 4, where the predictions are compared with the raw data. The prediction results are consistent with the raw data. The prediction variance is also plotted in the figure and shows that the uncertainty is small. The prediction results of the 30° turning circle test are shown in Figure 5. The prediction results can well track the tendency of the raw data. The uncertainty in Figure 5 is much higher than that in Figure 4 because the training data, including 10°/10° and 20°/20° zigzag tests, completely reflect the dynamic characteristics of the ship when the rudder angle is less than 20°. However, the dynamic characteristics of a rudder angle of 30° are not included in the information range provided by the training data. Under this condition, the identified model can predict the motion of a large rudder angle with good generalization ability and maintains good accuracy while providing high uncertainty of prediction through the proposed method. Generalization verification is necessary for system identification. The ability of the identified model to predict other motions not included in the training data is called generalization. The generalization performance of the trained model is verified by predicting the motions, including the 15◦/15◦ zigzag maneuver and port 30◦ turning circle test. The prediction results of the 15◦/15◦ zigzag test are shown in Figure 4, where the predictions are compared with the raw data. The prediction results are consistent with the raw data. The prediction variance is also plotted in the figure and shows that the uncertainty is small. The prediction results of the 30◦ turning circle test are shown in Figure 5. The prediction results can well track the tendency of the raw data. The uncertainty in Figure 5 is much higher than that in Figure 4 because the training data, including 10◦/10◦ and 20◦/20◦ zigzag tests, completely reflect the dynamic characteristics of the ship when the rudder angle is less than 20◦ . However, the dynamic characteristics of a rudder angle of 30◦ are not included in the information range provided by the training data. Under this condition, the identified model can predict the motion of a large rudder angle with good generalization ability and maintains good accuracy while providing high uncertainty of prediction through the proposed method.

The root mean square error (RMSE) is applied to evaluate the model performance and is defined as: The root mean square error (*RMSE*) is applied to evaluate the model performance and is defined as:

$$RMSE = \sqrt{\frac{1}{n} \sum\_{i=1}^{n} (y\_i - y\_i)^2} \tag{24}$$
 
$$\dots \dots \dots \dots \dots \text{ is the numerable. The } BM \mathbb{C} \mathbb{F}\_2 \text{ of the } 1 \mathbb{F}^\diamond \text{ (15}^\diamond \text{) gives}$$

݊(ݕෝ−ݕ <sup>ప</sup> )ଶ ୀଵ (24) where *y*ˆ*<sup>i</sup>* is the prediction result, and *y<sup>i</sup>* is the raw value. The *RMSEs* of the 15◦/15◦ zigzag maneuver are listed in Table 3 with each DOF.

where ݕෝప is the prediction result, and ݕ is the raw value. The RMSEs of the 15°/15° zig-

zag maneuver are listed in Table 3 with each DOF.

*J. Mar. Sci. Eng.* **2021**, *9*, x FOR PEER REVIEW 8 of 15

*J. Mar. Sci. Eng.* **2021**, *9*, x FOR PEER REVIEW 8 of 15

**Figure 4.** Prediction of the 15°/15° zigzag test. **Figure 4.** Prediction of the 15◦/15◦ zigzag test. **Figure 4.** Prediction of the 15°/15° zigzag test.

**Figure 5.** Prediction of the 30° turning circle maneuver. **Figure 5.** Prediction of the 30° turning circle maneuver. **Figure 5.** Prediction of the 30◦ turning circle maneuver.

**Table 3.** Prediction accuracy assessed by the RMSE of the 15°/15° zigzag test.

**Table 3.** Prediction accuracy assessed by the RMSE of the 15°/15° zigzag test.

 **RMSE**  Surge speed 0.1130 Sway speed 0.0816 Yaw rate 0.0806

 **RMSE**  Surge speed 0.1130 Sway speed 0.0816 Yaw rate 0.0806


**Table 3.** Prediction accuracy assessed by the *RMSE* of the 15◦/15◦ zigzag text.

#### *4.2. Experimental Study of a Ship Scale Model 4.2. Experimental Study of a Ship Scale Model*

The experimental dataset of KCLCC2 from SIMMAN is used to further validate the proposed method. KVLCC2 is a scale model of large tankers. The main particulars of the scale ship model are detailed in Table 4. The model free-running tests are performed by the Hamburg Ship Model Basin (HVSA). The experimental dataset of KCLCC2 from SIMMAN is used to further validate the proposed method. KVLCC2 is a scale model of large tankers. The main particulars of the scale ship model are detailed in Table 4. The model free-running tests are performed by the Hamburg Ship Model Basin (HVSA).

**Table 4.** Particulars of KVLCC2. **Table 4.** Particulars of KVLCC2.


There is some interference and noise in the experimental dataset. Directly taking the speed as the input and output will reduce the identification accuracy. Using the acceleration obtained by the speed difference as the prediction output can effectively alleviate the influence of noise [29], and the new flowchart is shown in Figure 6. There is some interference and noise in the experimental dataset. Directly taking the speed as the input and output will reduce the identification accuracy. Using the acceleration obtained by the speed difference as the prediction output can effectively alleviate the influence of noise [29], and the new flowchart is shown in Figure 6.

**Figure 6.** Flowchart of the ship maneuvering identification and prediction using GP (acceleration output). **Figure 6.** Flowchart of the ship maneuvering identification and prediction using GP (acceleration output).

> To include more dynamic characteristic information of the ship, the experimental data of 10°/5°, 20°/5° and 30°/5° zigzag maneuvers are used to train the GP. Moreover, the empirical Bayes method [46] is applied here to reduce the noise of acceleration with 'wdenoise' MATLAB function. More details of the empirical Bayes denoising method for ship motion data can be found in our previous study [47]. Figure 7 shows the effect of noise reduction. There are 800 training points with a time interval of 0.6 s. It took 10.1 s to train 3 GPs for the ship maneuvering system, and the obtained hyperparameter settings are shown in Table 5. To include more dynamic characteristic information of the ship, the experimental data of 10◦/5◦ , 20◦/5◦ and 30◦/5◦ zigzag maneuvers are used to train the GP. Moreover, the empirical Bayes method [46] is applied here to reduce the noise of acceleration with 'wdenoise' MATLAB function. More details of the empirical Bayes denoising method for ship motion data can be found in our previous study [47]. Figure 7 shows the effect of noise reduction. There are 800 training points with a time interval of 0.6 s. It took 10.1 s to train 3 GPs for the ship maneuvering system, and the obtained hyperparameter settings are shown in Table 5.

*J. Mar. Sci. Eng.* **2021**, *9*, x FOR PEER REVIEW 10 of 15

**Figure 7.** Raw data and denoising data of the sway motion by the empirical Bayes method. **Figure 7.** Raw data and denoising data of the sway motion by the empirical Bayes method.

Then, the identified model is validated by comparing the experimental data with pre-

**Table 5.** Selection of the GP parameters for the container ship. **Table 5.** Selection of the GP parameters for the container ship.


dictions of 15°/5°, 35°/5° and 10°/10° zigzag maneuvers, as shown in Figures 8–10. The accuracy of the prediction speed assessed by RMSE is listed in Table 6. In the similar study [29], the same three sets of training data were used for training the nu-SVM. The prediction error of the proposed method can be compared with the nu-SVM in [29]. Figures 8 and 9 show that the experimental data and prediction follow similar trends, and the cumulative deviation is small. The comparison results between Table 6 and the error results in [29] indicate that the proposed method has stronger prediction ability than nu-SVM. However, in Figure 10, the deviation between prediction speed and experiment is obvious, especially in surge motion. The predicted acceleration is also plotted in Figure 11 to analyze the reason. There is a strong oscillation in the measurements of the surge speed. This oscillation in accretion causes a cumulative deviation in speed. As for the uncertainty, it can be observed that the variance of the predictions of the experiment is bigger than the simulation in the previous case. This is because there are more disturbances and noises in the experiment than the simulation. Then, the identified model is validated by comparing the experimental data with predictions of 15◦/5◦ , 35◦/5◦ and 10◦/10◦ zigzag maneuvers, as shown in Figures 8–10. The accuracy of the prediction speed assessed by *RMSE* is listed in Table 6. In the similar study [29], the same three sets of training data were used for training the nu-SVM. The prediction error of the proposed method can be compared with the nu-SVM in [29]. Figures 8 and 9 show that the experimental data and prediction follow similar trends, and the cumulative deviation is small. The comparison results between Table 6 and the error results in [29] indicate that the proposed method has stronger prediction ability than nu-SVM. However, in Figure 10, the deviation between prediction speed and experiment is obvious, especially in surge motion. The predicted acceleration is also plotted in Figure 11 to analyze the reason. There is a strong oscillation in the measurements of the surge speed. This oscillation in accretion causes a cumulative deviation in speed. As for the uncertainty, it can be observed that the variance of the predictions of the experiment is bigger than the simulation in the previous case. This is because there are more disturbances and noises in the experiment than the simulation.

**Table 6.** Prediction accuracy assessed by the *RMSE* of KVLCC2 with the proposed method.


*J. Mar. Sci. Eng.* **2021**, *9*, x FOR PEER REVIEW 11 of 15

**Figure 8.** Prediction of the 15°/5° zigzag test. **Figure 8.** Figure **8.** Prediction of the 15◦/5◦ zigzag test. **Figure 8.** Prediction of the 15°/5° zigzag test.

**Figure 9.** Prediction of the 35°/5° zigzag test. **Figure 9. Figure 9.**  Prediction of the 35 Prediction of the 35◦/5°/5◦° zigzag test. zigzag test.

*J. Mar. Sci. Eng.* **2021**, *9*, x FOR PEER REVIEW 12 of 15

**Figure 10.** Prediction of the 10°/10° zigzag test. **Figure 10.** Prediction of the 10◦/10◦ zigzag test. Yaw rate 0.0028 0.0099 0.0044

**Figure 11.** Prediction acceleration of the 10◦/10◦ zigzag test.

#### **5. Discussion and Conclusions**

In this work, a novel identification modeling and prediction scheme based on GP is proposed to identify the ship nonparametric maneuvering model. By introducing the moment matching approximation method, the multi-step prediction uncertainty of ship motion can be propagated. The performance of the proposed method has been tested with a large container ship and a scale ship model and shows good accuracy and generalization ability. Moreover, the uncertainty of propagation can help drivers or controllers make safe decisions. Through the simulation of the container ship, it is proven that the prediction uncertainty obtained by the proposed method is reliable enough. Where there is less dynamic information in the training data, the prediction uncertainty of turning circle motion is larger than that of zigzag maneuver. In addition, it has been demonstrated that the performance of the presented approach is superior to the nu-SVM method in the experimental case. There are also some limitations of this study: (1) The proposed method needs to spend more calculation time due to consider the uncertainty propagation compared with other methods. The sparse method can be used to improve computational efficiency. (2) Both the two verified cases in this paper are container ships. The applicability of the model to other ship types, especially new unmanned ships, needs further study.

Future work includes two main tasks: (1) Although the presented method has been verified by simulation and experimental data, full-scale trials with disturbances should be performed, including waves, currents, and wind. In this environment, the uncertainty prediction provided by this method will have great application value. (2) This method can be used in modern controllers such as model predictive control. The uncertainty of predictions can be introduced in the cost function to construct a cautious controller.

**Author Contributions:** Conceptualization, Y.X.; methodology, Y.X. and G.X.; software, Y.X. and G.C.; validation, G.C. and Y.L.; resources, Y.L.; data curation, Y.X. and G.C.; writing and editing, Y.X. and G.X. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by National Key R&D Program of China under grant number 2019YFB2005303, Shandong Provincial Key Research and Development Program Major Scientific and Technological Innovation Project under grant number 2019JZZY010802, National Key Research and Development Program of China under grant number 2017YFE0115000, National Natural Science Foundation of China under grant number 52001186.

**Institutional Review Board Statement:** This study does not involve humans or animals.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The authors would like to thank the Hamburg Ship Model Basin (HVSA) and SIMMAN for sharing the KVLCC2 experimental data.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

