A Probabilistic Approach for Predicting Vessel Motion

Dear Editor,

This letter addresses the challenge of forecasting the motion of real-world vessels over an extended period with a limited amount of available data. By employing stochastic differential equation (SDE) modeling, we integrate both deterministic and stochastic components of the available information. Subsequently, we establish a recursive prediction methodology based on Bayes’ rule to update the model state when new measurements are received. Furthermore, we develop a stochastic model tailored specifically to vessel dynamics and introduce an approximation method to tackle computational complexities. Finally, we present an application example and conduct a comparative experiment to validate the effectiveness and superiority of the proposed method.

Introduction: The emergence of the automatic identification system (AIS) has provided an accurate and reliable measurement source that complements traditional methods such as marine radar. This has facilitated various applications of the vast amount of measurement data available. The prediction of vessel motion, a cornerstone in numerous applications of AIS data, is recognized for the complexities inherent in vessels’ dynamics and the intricate nature of the marine environment. However, it remains relatively less explored compared to the existing volume of research in the fields of vehicle motion prediction and planning, as well as pedestrian trajectory prediction [1]-[4].

There are primarily three modeling methods to solve the prediction problem. The first method [5] is based on recurrent neural networks, such as long short-term memory (LSTM) networks and encoder-decoder models. This approach aims to learn a nonlinear recurrent model from a vast amount of historical measurement data and utilize the learned model to recursively predict vessels’ motion. The work presented in [6] utilizes a recurrent encoder-decoder neural network for both prediction and estimation of the prediction uncertainty. The second modeling method [7] is based on machine learning, with one of the most representative techniques being Gaussian processes, where the data is utilized to fit both the mean function and the covariance function. The study outlined in [8] tackles ship trajectory prediction uncertainty by characterizing lateral motion uncertainty using Gaussian process. The third method [9], such as the extended Kalman filter (EKF), derives a model with uncertainty based on physical laws, typically employing a kinematic model with Gaussian noise. Perera et al. [10] introduced an EKF for vessel state estimation, subsequently applying it to predict vessel trajectories.

While the first method heavily relies on abundant measurement data for effective modeling, its performance may suffer in scenarios with limited training data, leading to robustness and transferability issues. The second method may encounter difficulties in modeling high-dimensional data and is highly dependent on the choice of kernel function. The prediction outcomes of the third method, generated by a model comprising a kinematic model and white noise, frequently result in straight-line trajectories. This oversimplified trajectory fails to capture the complexity of real navigation trajectories.

In this letter, the Itô SDE is employed to leverage the limited available information and establish the model. The data includes sparse measurements of vessel motion, such as geodetic latitude, geodetic longitude, yaw angle, speed over ground (SOG), and course over ground (COG). Additionally, it provides details on navigation environments, including nautical charts, as well as vessel characteristics like length, width, and type. Subsequently, the Markov property of the Itô diffusion and Bayes’ rule are utilized to update the model state upon receiving new measurements. Following this, a specific demonstration model and approximation algorithm are presented. Finally, an application example is provided, and a comparative experiment is conducted to illustrate the effectiveness and superiority of the proposed method.

Notations: Throughout this letter, subscript N signifies the time $t_N$. Sometimes it is convenient to write $X(t)$ instead of $X_t$, and thus we will use these two notations interchangeably.

Proposed method: Consider a set of N temporal non-homogeneous measurements, denoted by $\mathbb{Y} = \{{\bf{y}}_1, \ldots, {\bf{y}}_N\}$, $N > 1$. These measurements represent realizations of the observable stochastic process $\{Y_1, \ldots, Y_N\}$, which is generated by a continuous-discrete state space model characterized by the following equations:

where $k = 1, \ldots, N$. Itô SDE (1) represents the dynamic model, while (2) represents the observation model. The continuous state variable $X_t \in \mathbb{R}^n$ is a latent random vector observed indirectly through its effect on the distribution of $Y_k$. The discrete observation variable $Y_k \in \mathbb{R}^p$ is assumed observable at time $t_k$. The function ${\bf{f}}$ denotes the drift coefficient determining the nominal dynamics of the model, while ${\bf{D}}(X_t) \in \mathbb{R}^{n \times m}$ is the dispersion matrix that determines how the m-dimensional standard Wiener process enters the model. ${\bf{h}}$ represents the observation model function, and $V_k \sim {\cal{N}}({\bf{0}}, {\bf{R}}_k)$ is a Gaussian observation noise.

The prediction problem aims to determine the probability density function (PDF) for the state variable $X_t$ at time $t_f > t_N$, which depends on the measurement set $\mathbb{Y}$ while simultaneously considering the navigation environment for this specific problem. This problem can be mathematically formulated as

where ${\bf{y}}_{1:N}$ represents the measurements from time $t_1$ to $t_N$. The knowledge of this PDF would allow prediction of the state for any performance criterion.

Navigating vessels typically adhere to predetermined waypoints established based on various criteria, such as mission objectives, collision avoidance, international regulations, and other factors. While reconstructing the exact waypoints followed by vessels in the real world based on a limited number of measurements is impractical, general experiences in the navigation environment allow for an estimation of potential waypoints. For example, a common guideline is for vessels to avoid grounding. We represent the waypoint that remains valid in the future time span $[t_N, t_p)$ by a random constant $S \in \mathbb{R}^{q}$, where $t_p$ is an unknown time depending on the upcoming measurements. The random constant S is defined such that ${d} S_t = 0$, $S(t_N) = S$ for $t \in [t_N, t_p)$. Therefore, (3) can be expanded as

The time-homogeneous Itô diffusion $\{X_t\}$ in (1) satisfies the Markov property [11] within the time horizon $[t_N, t_p)$. Thus, we have

where the augmented state variable $Z_t \stackrel{\triangle}{=} [X_t,\, S_t] \in \mathbb{R}^{n+q}$. The second equation in (4) leverages the Markov property of $\{X_t\}$: $p({\bf{x}}_f | {\bf{x}}_N, {\bf{y}}_{1:N}) = p({\bf{x}}_f | {\bf{x}}_N)$. Additionally, the assumption that $\{S_t\}$ remains a random constant for $t \in [t_N, t_p)$ is utilized.

The SDE for the augmented state variable $Z_t$ can be formulated as follows:

where ${\bf{g}}$ represents the new drift coefficient and ${\bf{L}}(Z_t) \in \mathbb{R}^{(n+q) \times m}$ denotes the augmented dispersion matrix. The Wiener process $\{W_{t}\}$ and the corresponding observation model remain consistent with those in (1) and (2), respectively. The distinction between ${\bf{f}}(X_t)$ and ${\bf{f}}(X_t, S_t)$ lies in the treatment of S: in the former, it is a deterministic constant, while in the latter, it is regarded as a random constant. This distinction also applies to ${\bf{D}}(X_t)$ and ${\bf{D}}(X_t, S_t)$.

Lemma 1 [12]: The PDF $p({\bf{z}})$ of the solution to (5) satisfies the Kolmogorov forward equation

with the initial condition $p({\bf{z}}(t_1))$ at time $t_1$.

When considering $p({\bf{z}}_N | {\bf{y}}_{1:N})$ as the initial PDF, $p({\bf{z}}_f | {\bf{y}}_{1:N})$ can be obtained by solving the partial differential eqaution (PDE) in (6).

Assumption 1: $p({\bf{x}}_N, {\bf{s}}_N | {\bf{y}}_{1:N}) = p({\bf{x}}_N| {\bf{y}}_{1:N}) p({\bf{s}}_N | {\bf{y}}_{1:N})$.

Remark 1: Assumption 1 implies that $S_N$ depends on the measurements up to time $t_N$ and is independent on the state variable $X_N$. The independence between $S_N$ and $X_N$ suggests that the waypoints, which reflect the intentions of the sailor in determining the vessel’s control actions, are generally unaffected by the vessel’s current state. Instead, they are influenced by more realistic features.

Accordingly, $p({\bf{z}}_N | {\bf{y}}_{1:N})$ can be expressed as

The second equation in (7) illustrates the general filtering process where Bayes’ rule is applied. The relationship between $p({\bf{x}}_N | {\bf{y}}_{1:(N-1)})$ and $p({\bf{x}}_f | {\bf{y}}_{1:N})$ is established by (4) and (7). This enables the prediction process to be carried out by recursively computing these two equations, where $p({\bf{x}}_1)$ and $p({\bf{s}}_1)$ need to be initialized, along with the specification of the PDF $p({\bf{s}}_k | {\bf{y}}_{1:k})$ each time a new measurement is obtained. Given that large vessels navigating in the ocean typically move slowly relative to the vast sea surface, $p({\bf{s}}_k | {\bf{y}}_{1:k})$ remains constant over a relatively long time span.

Experiments: The data acquired from AIS typically comprises five types of measurements: geodetic latitude (ϕ), geodetic longitude (λ), yaw angle (ψ), SOG (U), and COG (χ). The yaw angle, also referred to as heading, represents the Euler angle between the body-fixed frame and the local geodetic frame along the $z$-axis. The SOG and COG are computed from the vessel’s previous position, current position, and the time interval between them. The kinematic equation is directly expressed in geodetic coordinates [13]

where $u_n \stackrel{\triangle}{=} U \text{cos} (\psi + \beta)$ and $u_e \stackrel{\triangle}{=} U \text{sin} (\psi + \beta)$ represent the northward and eastward components of the vessel’s linear velocity, respectively. The angular velocity of the vessel is denoted as r. Here the sideslip angle β is defined as $\beta \stackrel{\triangle}{=} \chi - \psi$. The presence of the term $\tan(\phi)/[R_N(\phi)+h]$ distinguishes the geodesic on the reference ellipsoid from the misleading “straight line”. The units and range of measurement values are detailed in Table I (see Supplementary Material).

Assumption 2: ${d}{U} = d_1 {d}W^{(1)}_{t}$, where $\{W^{(1)}_{t}\}$ represents the one-dimensional standard Wiener process.

Assumption 3 [14]: ${d}{\beta} = 0$.

Remark 2: Assumption 2 implies that the SOG U changes slowly when the vessel is in a navigation state. Furthermore, Assumption 3 indicates that the sideslip angle β remains small and constant while the vessel follows the path to waypoints.

The dynamic model of r adopts the first-order Nomoto model [15]. The SDE for $[U, \beta, r]^{T}$ is given as follows:

where the one-dimensional standard Wiener process $\{W^{(2)}_{t}\}$ represents the noise stemming from parameter uncertainty and external disturbances. The rudder angle δ is the control input that needs to be designed for the model to follow the path to waypoints.

The random constant S is defined as the vector $S = [\phi_d, \lambda_d, \alpha_d]^{T}$, where the first two elements represent the desired latitude and longitude, respectively, while $\alpha_d$ denotes the desired azimuth angle. It follows a Gaussian distribution given by $S_N | Y_{1:N} \sim {\cal{N}}({\bf{m}}_s, {\bf{P}}_s)$. We consider the problem of following a straight-line path to reach the designated state S. The kinematic equations for line-of-sight (LOS) guidance systems are detailed in (10), where the equilibrium point of the cross-track error is proven to be uniformly semiglobally exponentially stable (USGES) [14]. The rudder angle command δ is generated using a heading autopilot of PD type

where the controller gains $K_p$ and $K_d$ are chosen to achieve the closed-loop dynamics with a damping ratio $\zeta \approx 0.53$ and a bandwidth $\omega_b \approx 0.05\;\text{rad/s}$. The control law in (10) is derived from the linearized version of kinematic model in (8). The variables x and y denote the coordinates of the vessel’s position in the local geodetic frame, where the frame’s origin is located at the waypoint S. The coordinates can be calculated as follows:

where the algorithm referred to as “Inverse” corresponds to the second algorithm detailed in [16]. The parameters occurred from (8) to (11) are listed in Table II (see Supplementary Material).

Given the complexity of solving (6), we prioritize the statistical moments of $X_t$ for prediction instead of directly handling its PDF. To address computational challenges, we employ the Gaussian assumed density approximation algorithm [17], detailed in Algorithm 1 (see Supplementary Material).

In the experimental simulation, the first measurement is used to determine the initial mean ${\bf{m}}^{-}_1$ with an assumption of angular velocity $r=0$. Given that the measurement value is accurate to the last decimal place, the initial covariance matrix ${\bf{P}}^{-}_1$ is constructed as follows: $\text{diag}[(\frac{10^{-4} \pi}{180})^2,\; (\frac{10^{-4} \pi}{180})^2, \;(\frac{10^{-1} \pi}{180})^2,\; (\frac{10^{-2}}{1.94384})^2, \; (\frac{10^{-1} \pi}{180})^2,\;(\frac{10^{-2} \pi}{180})^2]$. The diagonal elements in the measurement noise covariance matrix ${\bf{R}}_k$ are assigned as the first five diagonal elements of ${\bf{P}}^{-}_1$.

In Fig. 1, four selected measurements, each spanning an average time of over 10 minutes, are acquired from a cargo ship navigating on the sea surface near the Tsing Ma Bridge in Hong Kong, China. The distribution of random constants, representing the estimated waypoints based on the sea-continent constraint, is illustrated. The prediction time horizon is set to 15 minutes. Upon the arrival of measurements, the mean value of the Itô diffusion $\{X_t\}$, and the 1 standard deviational ellipse is displayed. Before obtaining future measurements, the confidence envelope curve for the ellipses is plotted. Fig. 1 illustrates that the measurements fall within the predicted ellipse during straight navigation; however, they lag behind the ellipse after the ship turns left. The discrepancy arises from the omission of the vessel’s speed reduction dynamics during the turn in the model represented by (9).

Figure  1.  Position prediction

DownLoad: Full-Size Img PowerPoint

To demonstrate the superiority of our method, we compare it with the EKF. In the comparative experiment, we maintain all parameters unchanged, with the sole adjustment of setting $\delta = 0$. We introduce a new metric called the geodesic error rate (GER), which is measured in meters per minute and defined as

where $\hat{p}_k$ and $p_k$ represent the geodetic coordinate tuples of the predicted position and measured position, respectively. The geodesic error (GE) $s_{12}$ is determined by the algorithm in (11). The GER assesses the prediction error in a temporal-spatial context, providing an accumulated rate of the error over time. Additionally, we introduce another metric called the ellipse area proportion index (EAPI), which is defined as the ratio of the area of the predicted ellipse to the product of the length and width of the ship. The GE, GER, and EAPI for each measurement are detailed in Table III (see Supplementary Material).

Table III shows that the proposed method outperforms the EKF across all metrics, especially noticeable after left turns. This emphasizes the substantial influence of representing waypoints as random constants. Moreover, the EAPI metric highlights that the predicted ellipses are considerably smaller in the proposed method compared to the EKF, indicating higher prediction precision.

Conclusion: We present a probabilistic approach for predicting vessel motion over an extended period, especially in scenarios with limited data, while taking into account the navigation environment. Our approach leverages stochastic modeling and Bayes’ rule to provide a systematic solution to the prediction problem. Through numerical simulations and comparative experiments using a customized stochastic vessel model, we validate the effectiveness and superiority of the proposed method. In future research, we aim to establish a more precise stochastic framework for modeling the constraints or incentives that influence sailors’ decision-making processes, moving beyond simple random constants for waypoints.

Acknowledgments: This work was supported by the National Natural Science Foundation of China (62073019), the Key R&D Program of Hebei Province (22340301D), China Postdoctoral Science Foundation (2021M703021), and Hebei Postdoctoral Science Foundation (B2021003031)

Supplementary Material: The supplementary material of this letter can be found in links https://drive.google.com/file/d/1Miumx4pkqaYKyl9v6Oj9HajWkOR2I2EB/view?usp=drive_link.