Long-Term Trajectory Prediction for Oil Tankers via Grid-Based Clustering

Abstract

Vessel trajectory prediction is an important step in route planning, which could help improve the efficiency of maritime transportation. In this article, a high-accuracy long-term trajectory prediction algorithm is proposed for oil tankers. The proposed algorithm extracts a set of waymark points that are representative of the key traveling patterns in an area of interest by applying DBSCAN clustering to historical AIS data. A novel path-finding algorithm is then developed to sequentially identify a subset of waymark points, from which the predicted trajectory to a fixed destination is produced. The proposed algorithm is tested using real data offered by the Danish Maritime Authority. Numerical results demonstrate that the proposed algorithm outperforms state-of-the-art vessel trajectory prediction algorithms and is able to make high-accuracy long-term trajectory predictions.

1. Introduction

About 90% of global trade is carried by maritime transportation [1]. With the continuing growth of international trade, modern maritime transportation calls for more intelligent methods for transportation management to achieve larger capacity, faster traveling speed and higher safety levels. To achieve these goals, accurate predictions of vessels’ future movement is important and can be used in many maritime applications, such as port management, anomaly detection and collision avoidance [2].

Despite the importance of vessel trajectory prediction, it remains a challenging task due to the diverse navigation environments and the stochastic nature of vessel movements. While most of the existing works on vessel trajectory prediction have focused on making short-term predictions [3,4,5,6,7], being particularly useful for collision warning and avoidance, this work investigates the problem of long-term vessel trajectory prediction. Long-term trajectory predictions can be used to guide the captains to operate the ship in a more fuel-efficient way and hence reduce the carbon dioxide emissions, with assistance also from additional information, including weather and current forecasts on the route. Long-term vessel trajectory predictions can also be used by the agencies of port management to obtain more accurate estimates of the remaining navigation distance and subsequently obtain more accurate predictions for vessel arrivals.

The prevalence of maritime transportation data makes it possible to take a data-driven approach to develop high-accuracy vessel trajectory algorithms [2,8]. For instance, the Automatic Identification System (AIS), a global autonomous tracking system that has been made compulsory for ships exceeding 300 tons, provides abundant and near real-time information about ships. Apart from static information such as MMSI, ship name and ship type, the AIS messages also contain the location information of the ship (longitude (Lon) and latitude (Lat)) and information about its traveling (e.g., speed over ground (SOG), course over ground (COG) and heading) [9]. This information, when gathered in large scale, can be used to characterize the voyage patterns in a corresponding area and can further be exploited to predict vessel trajectory.

This paper presents an AIS data based long-term vessel trajectory prediction algorithm, aiming to predict the trajectory of oil tankers from any location to a pre-defined destination, e.g., a port. The proposed algorithm takes a data-driven approach. First, the traveling patterns of tankers in an area of interest are extracted from the historical AIS data via key point clustering. The output of clustering is a set of waymark points, with each corresponding to a cluster. Each waymark point is characterized by the average Lon, Lat, COG and the number of key points in the cluster that it represents. In real-time trajectory prediction, the algorithm uses the Lon and Lat to filter the waymark points and uses the information of COG and the number of points to calculate a weighted distance between the filtered waymark points and the current reference point. A segment of the predicted trajectory is generated by connecting the reference point to the waymark point that has the smallest weighted distance. The complete long-term trajectory prediction is made by repeating this process until the reference point reaches the destination. For the proposed algorithm, the key point clustering only needs to be perform once. Once the set of waymark points is obtained, all subsequent oil tankers arriving at the same destination can use the obtained set of waymark points to achieve trajectory prediction. We verify the effectiveness of the proposed algorithm using real AIS data provided by the Danish Maritime Authority (DMA) [10].

The rest of this paper is organized as follows. Section 2 provides a brief review of related work in the area of trajectory prediction. Section 3 describes the proposed trajectory prediction algorithm. In Section 4, the proposed algorithm is applied to real AIS data and is compared to state-of-the-art algorithms. Finally, Section 5 draws conclusions and discusses possible future research directions. For ease of exposition,

Table 1. List of key notations.

Key Notation
${\lambda }_k$	Longitude
${\varphi }_k$	Latitude
${\theta }_i$	COG $\in [0^{°},{360}^{°})$
${\tilde{p}}_k=\left({\lambda }_k,{\varphi }_k,{\theta }_k\right)$	The k-th key point
$p_k=\left({\lambda }_k,{\varphi }_k\right)$	The location of the k-th key point
$q_m=\left({\overline{\lambda }}_m,{\overline{\varphi }}_m,{\overline{\theta }}_m,{\overline{\rho }}_m\right)$	The m-th waymark point, ${\overline{\lambda }}_m$, ${\overline{\varphi }}_m$, ${\overline{\theta }}_m$, ${\overline{\rho }}_m$ denotes the Lon, Lat, COG and the point density
$p_{DES}$	Destination location
$dire\left(p_m,p_n\right)$	The heading direction from position $p_m$ to position $p_n$
$dist\left(p_m,p_n\right)$	The spherical distance between position $p_m$ and $p_n$
$\overrightarrow{a}$	Direction $\overrightarrow{a}$, a unit vector of the origin of a two-dimensional coordinate system
$LDM\left\{\overrightarrow{a_1},\overrightarrow{a_2},...,\overrightarrow{a_m}\right\}$	The linear directional mean of direction $\overrightarrow{a_1},\overrightarrow{a_2},...,\overrightarrow{a_m}$, defined in Equation (2)
${\Delta }_{\overrightarrow{a},\overrightarrow{b}}$	The angular change from direction $\overrightarrow{a}$ to direction $\overrightarrow{b}$, defined in Equation (3)
${\mathcal{Q}}_n$	Candidate set identified in the n-th iteration, defined in Equation (4)
${\mathcal{Q}}_n^{\prime }$	Refined candidate set at the n-th iteration, defined in Equation (8), ${\mathcal{Q}}_n^{\prime }\subseteq {\mathcal{Q}}_n$
${\overrightarrow{d}}_{{cur}_n}$	The heading trend at the beginning of the n-th iteration, defined in Equation (5)
${\overrightarrow{d}}_{{\mathcal{Q}}_n^{\prime }}$	The mean COGs of the waymark points in ${\mathcal{Q}}_n^{\prime }$
${\overrightarrow{d}}_n$	The heading direction from position $p_{n-1}$ to position $p_n$
${\overrightarrow{d}}_{{DES}_n}$	The direction from position $p_n$ to the destination
${\delta }_m^{\prime }={\Delta }_{{\overrightarrow{d}}_{p_n,q_m},{\overrightarrow{d}}_{DE{S}_n}}$	Directional difference between ${\overrightarrow{d}}_{p_n,q_m}$ and ${\overrightarrow{d}}_{{DES}_n}$, defined in Equation (6)
${\delta }_m={\Delta }_{{\theta }_n,{\overline{\theta }}_m}$	Directional difference between the COG of $p_n$ and the COG of $q_m$, defined in Equation (7)

lists the key notation used in this paper.

2. Related Work

With the advances in high-precision positioning technologies such as the Global Positioning Systems and real-time radars, it is possible to acquire accurate information about positions for vehicles, aircraft, ships and pedestrians. The availability of high-accuracy positioning data, together with other information such as speed and acceleration, makes it possible to predict the trajectories of the targets of interest automatically, which could find applications in many areas, including terrestrial navigation [11,12], autonomous driving [13,14], and maritime traffic management [2,8,15,16]. This paper focuses on the trajectory prediction problem for vessels, for which the existing works can be classified into two categories, i.e., short-term trajectory prediction [3,4,5,6,7] and long-term trajectory prediction [17,18].

For short-term vessel trajectory prediction, the time horizon over which the predictions are made is usually from a few seconds to a few tens of minutes. Such vessel prediction algorithms are often developed for collision avoidance to ensure navigation safety. Hence, apart from the prediction accuracy, timeliness is also important [19]. In this category, classical algorithms follow from mathematical models of mobility and statistical techniques for accurate trajectory predictions. Examples of such methods include the Kalman filter-based algorithm [20], the Gaussian process-based algorithm [3], and the Markov process-based algorithm [21]. However, due to the complicated nature of vessel trajectories, these methods may not be able to capture the characteristics of the trajectories to predict and thus fail to provide accurate predictions.

The recent advances in machine learning techniques, along with the assistance of abundant AIS data, have stimulated a burst of research on machine-learning-based vessel trajectory algorithms [5,6,7]. For instance, Zhang et al. [5] used a hybrid method based on LSTM and KNN for short-term trajectory prediction. The algorithm switches between LSTM and KNN based on the trajectory densities of different areas. That is, in dense areas, KNN is used for trajectory prediction, while in sparse areas, LSTM is adopted. Using publicly available AIS data collected near Xiamen Port, Fujian Province from 2018 to 2019 [5], the results show that this method has better prediction accuracy compared with classical prediction algorithms developed based on Kalman filtering. You et al. [6] proposed a sequence-to-sequence model based on GRU to predict vessel trajectory at a time horizon of 5 min. The model works by first encoding the trajectory as a context vector to maintain the temporal relationship of the trajectory position, and then using GRU as a decoder to output the future trajectory. Numerical results based on the data collected near Chongqing and Wuhan of Yangtze River Channel demonstrate good short-term trajectory predictions. Murray et al. [7] proposed a data-driven trajectory prediction method which first applies Principal Component Analysis (PCA) to each trajectory to generate feature vectors and then uses the extracted feature vectors as inputs for the Gaussian Mixer model to predict multiple possible trajectories at the same time. The time horizon of prediction is about 30 min. All of the above work uses deep learning or data-driven methods and achieves better results.

Compared to the problem of short-term trajectory prediction, the progress on the long-term trajectory prediction problem is much less reported, partly due to the more challenging nature of the long-term prediction problem. Existing short-term prediction algorithms may be used to produce long-term predictions through recursions; however, the accuracies are expected to decrease as the number of prediction steps increases [6]. One relevant work that tackles the problem of long-term trajectory can be found in [22], where DBSCAN was used to cluster historical trajectories, and the trajectory predictions are produced by pre-trained deep learning models. This work was tested using the publicly available AIS data from MarineCadastre in Zones 15 and 16. While [22] also adopted DBSCAN, the overall methodology is completely different from our proposed approach. For instance, in [22], DBSCAN was used to cluster trajectories based on trajectory statistics such as the trimmed mean of the longitude (latitude). In contrast, our work adopts DBSCAN to cluster the key points from many historical trajectories, where key points sharing similar characteristics are grouped into the same cluster. Two other relevant works on long-term vessel trajectory predictions can be found in [17,18]. Both these two works attempt to predict the remaining path of ships to a destination in order to better estimate the remaining traveling distance and to achieve a more accurate estimate of the arrival times. Ogura et al. [17] evaluated the difference between the predicted weather and the weather of historical trajectories, and then adopts the historical trajectory with the smallest difference as the prediction. The AIS data of four cargos traveling to and from Japan was used to test the weather-based algorithm [17]. Alessandrini et al. [18] adopted a graph-based approach by calculating the raster of density and directionality of all the historical data and applying a path-finding algorithm to identify a path from the current location to the destination. This method was tested using AIS data recorded in October 2015 near the port of Trieste, an Italian city on the Northeastern Adriatic Sea, with the assistance of historical LRIT data from 2009 to 2014 [18].

3. Methodology

In this section, we describe the proposed long-term vessel trajectory prediction algorithm. Figure 1 presents a flow chart of the proposed algorithm. As shown in the figure, the algorithm has the following two modules: the key point clustering module and the trajectory finding module. In the remaining part of this section, we present the details of the two modules.

3.1. Key Point Clustering

The key point clustering module aims to extract from pre-processed training trajectories a set of waymark points that are representative of the traveling patterns of the oil tankers in an area of interest. Instead of keeping all the raw AIS data points in the training trajectories, each trajectory is pre-processed to remove redundant data and to reduce the number of points retained in the trajectory. (A detail description of data pre-processing will be presented in Section 4.1 along with the experiments.) The points retained in the training trajectories are called the key points. The key point clustering module treats all key points (from all different training trajectories) as the input and produces a set of waymark points. As a result of clustering, each waymark point not only reflects the Lon and Lat of tankers passing through the vicinity of the corresponding cluster, but also contains information such as the COG to better capture the motion of tankers in the corresponding area.

Each predicted trajectory consists of a number of waymark points. Hence, to achieve accurate predictions, there should be a sufficient number of waymark points that can comprehensively capture all possible routes in the area of interest. Additionally, waymark points that are close in space should have distinguishing features related to the motion of the vessels to make clear selections among the waymark points. To fulfill these purposes, we adopt a COG-based clustering process to generate the waymark points [23].

To facilitate the description of the clustering process, we use ${\lambda }_k$, ${\varphi }_k$ to represent the Lon and Lat of the key point ${\tilde{p}}_k$, respectively, and ${\theta }_k$ the COG of ${\tilde{p}}_k$. Hence, each key point is represented by ${\tilde{p}}_k=\left({\lambda }_k,{\varphi }_k,{\theta }_k\right)$. Let $\mathcal{D}=\left\{{\tilde{p}}_1,{\tilde{p}}_2,{\tilde{p}}_3,\dots ,{\tilde{p}}_K\right\}$ be the set of key points in the training trajectories. The core idea of the clustering algorithm is first to discretize the entire area of interest into a set of grids and then run DBSCAN for each grid to cluster the key points falling into the corresponding area. Specifically, DBSCAN is used to cluster the key points within the grid area based only on the COG to extract representative traveling directions. Each cluster identified by DBSCAN forms a waymark point $q_m$, which contains four-dimensional information, including the average Lon ${\overline{\lambda }}_m$ and average Lat ${\overline{\varphi }}_m$ of all the key points in the corresponding cluster, the median COG of the key points, i.e., ${\overline{\theta }}_m$, and the number of key points ${\overline{\rho }}_m$ in the cluster, i.e., $q_m=\left({\overline{\lambda }}_m,{\overline{\varphi }}_m,{\overline{\theta }}_m,{\overline{\rho }}_m\right)$.

Each of the grid areas can be represented as:

$$\begin{array}{c}\hfill G_{i,j}=\left\{\left(\lambda ,\varphi \right)\right|\lambda \in \left({\overline{\lambda }}_i-{\delta }_{\lambda },{\overline{\lambda }}_i+{\delta }_{\lambda }\right),\varphi \in \left({\overline{\varphi }}_j-{\delta }_{\varphi },{\overline{\varphi }}_j+{\delta }_{\varphi }\right)\},\end{array}$$

(1)

where ${\overline{\lambda }}_i$ and ${\overline{\varphi }}_j$ are the Lon and Lat of the center of the grid, and ${\delta }_{\lambda }$ and ${\delta }_{\varphi }$ are the widths in the longitude and latitude direction. The difference of the center Lon (Lat) between two adjacent grid areas can be made smaller than $2{\delta }_{\varphi }$ ($2{\delta }_{\lambda }$), e.g., $\left|{\overline{\lambda }}_i-{\overline{\lambda }}_{i+1}\right|<2{\delta }_{\varphi }$. In such cases, two adjacent grid areas may partially overlap with each other, and the key points falling into the overlapped areas are used multiple times during clustering. Such a setup makes it easier to obtain a representative set of waymark points from non-uniformly distributed key points by using a relatively large grid width.

It is noted that a customized distance metric is adopted when applying DBSCAN for COG clustering. Specifically, suppose ${\theta }_i\in [0^{°},{360}^{°})$ and ${\theta }_j\in [0^{°},{360}^{°})$ are two COGs. Then, the corresponding DBSCAN distance metric is calculated as $|{\theta }_i-{\theta }_j|$ if $|{\theta }_i-{\theta }_j|<{180}^{°}$ and $360-|{\theta }_i-{\theta }_j|$ otherwise. This customized metric is to avoid the errors caused by COGs close to the north, i.e., COG values close to $0^{°}$ or ${360}^{°}$. Additionally, it is noted that the Lon and Lat are not used in the clustering process. This is due to the fact that the COG clustering is performed for each grid area, where the key points have relatively close values of Lon and Lat.

As an illustrative example, Figure 2 presents the key points obtained from the trajectories of the oil tankers passing through the shown area between 2017 to 2019. It can be seen that there are areas where the key points are sparse, while there are also areas where the key points are much denser. It is known that DBSCAN requires one to specify a value on the minimum number of points in a cluster so the algorithm can operate. If there are less key points than this value, then all points in the grid area will be classified as noise, and no cluster can be formed. Therefore, in a sparse area, the width of the grid area has to be set sufficiently large. However, if overlapping is not permitted, then a large width will lead to a reduced number of waymark points, which will make the representation of the traveling patterns less accurate.

3.2. Trajectory Finding Based on Waymark Points

Without loss of generality, suppose a prediction of the future trajectory is made at time instant 0, where the historical trajectory consisted of the set of previous location points $\{p_{-I}\dots ,p_{-2},p_{-1},p_0\}$, and $p_n=({\lambda }_n,{\varphi }_n)$ contains the Lat and Lon of the n-th position. Denote $\mathcal{Q}=\{q_1,q_2,\dots ,q_M\}$ as the set of all waymark points obtained from the key point clustering. The task of the trajectory prediction is to select from $\mathcal{Q}$ a set of N points to produce the predicted trajectory $T=\{p_1,\dots ,p_n,\dots ,p_N,p_{DES}\}$, where $p_{DES}$ is the destination of interest. This process is implemented in an iterative way by sequentially identifying the N waymark points. Note that the number of point N is determined by the algorithm automatically and is not a preset value.

In the n-th iteration, the algorithm first identifies a candidate waymark point set ${\mathcal{Q}}_n\subseteq \mathcal{Q}$ based on $T_n=\{p_{-I}\dots ,p_0,\dots ,p_n\}$. The candidate set ${\mathcal{Q}}_n$ is then filtered to generate a refined candidate set ${\mathcal{Q}}_n^{\prime }$, from which the predicted point $p_{n+1}$ is chosen based on a weighted distance to be explained later in this section. In this process, the Lat, Lon and the information of direction and point density are all utilized.

For ease of exposition, we introduce the following notation:

•

$dire(p_m,p_n)$: the heading direction from position $p_m$ to position $p_n$;

•

$dist(p_m,p_n)$: the spherical distance between $p_m$ to position $p_n$, calculated using the Python package Geopy;

•

$LDM\{{\overrightarrow{a}}_1,{\overrightarrow{a}}_2,\dots ,{\overrightarrow{a}}_m\}$: the linear directional mean of direction ${\overrightarrow{a}}_1,{\overrightarrow{a}}_2,\dots ,{\overrightarrow{a}}_m$, calculated as:

$$LDM\{{\overrightarrow{a}}_1,{\overrightarrow{a}}_2,\dots ,{\overrightarrow{a}}_m\}=\angle \left(\frac{{\overrightarrow{a}}_1+{\overrightarrow{a}}_2+...+{\overrightarrow{a}}_m}{m}\right),$$

(2)

where ${\overrightarrow{a}}_i$, $i=1,\dots ,m$ is a unit vector from the origin.

•

${\Delta }_{\overrightarrow{a},\overrightarrow{b}}:$ the angular change from direction $\overrightarrow{a}$ to direction $\overrightarrow{b}$, measured in degrees. ${\Delta }_{\overrightarrow{a},\overrightarrow{b}}$ is calculated as:

$${\Delta }_{\overrightarrow{a},\overrightarrow{b}}=\left\{\begin{array}{cc}{360}^{°}-\Delta \hfill & \Delta \ge {180}^{°}\hfill \\ \Delta \hfill & \Delta \le {180}^{°}\hfill \end{array}\right.$$

(3)

where $\Delta =|\angle \overrightarrow{a}-\angle \overrightarrow{b}|$. Here, angle $\angle \overrightarrow{a}$, $\angle \overrightarrow{b}$ both take value in $[0^{°},{360}^{°})$, where $0^{°}$ corresponds to the north.

At the n-th iteration, the distance from $p_n$ to $p_{DES}$ is compared against a pre-determined threshold ${\epsilon }_{DES}$. If 

$$dist(p_n,p_{DES})<{\epsilon }_{DES},$$

i.e, the distance to the destination is sufficiently small, then the iterative process terminates, and $p_{DES}$ is appended to $T_n$ to form the completed trajectory prediction. Otherwise, the algorithm continues to identify $p_{n+1}$. The threshold ${\epsilon }_{DES}$ needs to be chosen with care, as setting it too large can lead to premature terminations of the trajectory finding, while setting it too small may make the prediction susceptible to interference from other routes in the areas of high concentrations of multiple routes towards different directions, e.g., the areas near a port. We recommend to set ${\epsilon }_{DES}$ to be the distance from the destination from which various routes begin to diverge.

3.2.1. Candidate Set ${\mathcal{Q}}_n$ Identification and Navigation Status Update

At iteration $n\ge 1$, the algorithm first identifies the candidate set ${\mathcal{Q}}_n$, which is made up of all the waymark points within a rectangle centered at the latest point in the trajectory, i.e., around point $p_n=({\overline{\lambda }}_n,{\overline{\varphi }}_n)$:

$$\begin{array}{c}\hfill {\mathcal{Q}}_n=\left\{({\overline{\lambda }}_m,{\overline{\varphi }}_m,{\overline{\theta }}_m,{\overline{\rho }}_m)\right|{\lambda }_m\in ({\overline{\lambda }}_n-{\epsilon }_{\lambda },{\overline{\lambda }}_n+{\epsilon }_{\lambda }),{\varphi }_m\in ({\overline{\varphi }}_n-{\epsilon }_{\varphi },{\overline{\varphi }}_n+{\epsilon }_{\varphi })\}\end{array}$$

(4)

where ${\epsilon }_{\lambda }>0$ and ${\epsilon }_{\varphi }>0$ are used to define the size of the rectangular area.

To ensure that the waymark point $p_{n+1}$ locates in the direction that is consistent with the recent trajectory trend, we update the current heading trend ${\overrightarrow{d}}_{{cur}_n}$ as:

$${\overrightarrow{d}}_{{cur}_n}=\left\{\begin{array}{cc}{\overrightarrow{d}}_n\hfill & n=1\hfill \\ LDM\left({\overrightarrow{d}}_{{cur}_{n-1}},{\overrightarrow{d}}_n,{\overrightarrow{d}}_{{\mathcal{Q}}_{n-1}^{\prime }}\right)\hfill & n>1\hfill \end{array}\right.$$

(5)

where ${\overrightarrow{d}}_n=dire(p_{n-1},p_n)$ and ${\overrightarrow{d}}_{{{\mathcal{Q}}^{\prime }}_{n-1}}$ is the mean of the COGs of the waymark points in set ${\mathcal{Q}}_{n-1}^{\prime }$, and ${\mathcal{Q}}_{n-1}^{\prime }$ is the filtered candidate set in the $(n-1)$-th iteration (see Section 3.2.2 for details). From the above equation, it can be seen that when $n>1$, ${\overrightarrow{d}}_{{cur}_n}$ is the linear directional mean of the (predicted) traveling direction of iteration $n-1$, the heading direction ${\overrightarrow{d}}_n$ from point $p_{n-1}$ to $p_n$ (the current position), and the mean COG of the key points in ${\mathcal{Q}}_{n-1}^{\prime }$.

3.2.2. Candidate Set Filtering

In this step, the waymark points in ${\mathcal{Q}}_n$ that deviate significantly from the expected traveling direction are filtered out from further consideration in the current prediction iteration. To do so, we first denote ${\overrightarrow{d}}_{DE{S}_n}\doteq dire(p_n,p_{DES})$ as the direction from the current position to the destination and then introduce the following angle-dependent parameters:

$${\delta }_m^{\prime }={\Delta }_{{\overrightarrow{d}}_{p_n,q_m},{\overrightarrow{d}}_{DE{S}_n}}$$

(6)

$${\delta }_m={\Delta }_{{\theta }_n,{\overline{\theta }}_m}$$

(7)

where ${\overrightarrow{d}}_{p_n,q_m}=dire\left(p_n,q_m\right)$ is the heading direction from $p_n$ to $q_m$. Parameter ${\delta }_m^{\prime }$ measures the degree of difference between ${\overrightarrow{d}}_{p_n,q_m}$ and the direction of the current position towards the destination, while ${\delta }_m$ measures the degree of difference between the COG of $p_n$ and the COG of $q_m$. An illustration of the above notation is given in Figure 3.

Ideally, if a tanker does not make a sharp turn, ${\delta }_m$ is expected to take small values. Additionally, as the tanker is traveling towards the destination, it is expected that the heading from $p_n$ to $q_m$ does not deviate too drastically from ${\overrightarrow{d}}_{DE{S}_n}$. Therefore, it is expected that ${\delta }_m^{\prime }$ does not take extremely large values. Therefore, we identify the refined candidate set by applying filtering based on ${\delta }_m$ and ${\delta }_m^{\prime }$ as follows:

$${\mathcal{Q}}_n^{\prime }=\left\{q_m|q_m\in Q_n,{\delta }_m\le {ϵ}_m,{\delta }_m^{\prime }\le {ϵ}_m^{\prime }\right\},$$

(8)

where ${ϵ}_m$ and ${ϵ}_m^{\prime }$ are two thresholds.

3.2.3. Waymark Point Selection

Once the filtered candidate set of waymark points, ${\mathcal{Q}}_n^{\prime }$, has been identified, the algorithm can proceed to make a selection for the next waymark point. As each waymark point $q_m=\left({\overline{\lambda }}_m,{\overline{\varphi }}_m,{\overline{\theta }}_m,{\overline{\rho }}_m\right)$ represents a cluster of key points, a higher value of ${\overline{\rho }}_m$ means that the corresponding grid area has been visited more frequently. Hence, it is reasonable to allocate more weight to such waymark points in the selection. Additionally, the expected trajectory trend, i.e., ${\overrightarrow{d}}_{cu{r}_n}$, is expected to be close to the heading from $p_n$ to $q_m$. Based on these considerations, we propose a weighted distance to account for both the frequency that a waymark point is visited and the angular change from $p_n$ to $q_m$. Specifically, the weighted distance between the current point $p_n$ and a candidate waymark point $q_m\in {\mathcal{Q}}_n^{\prime }$ can be calculated as:

$$Dis\left(p_n,q_m\right)=\kappa \sqrt{{\rho }_{max}/{\overline{\rho }}_m}+\left(1-\kappa \right){\Delta }_{{\overrightarrow{d}}_{cu{r}_n},{\overrightarrow{d}}_{p_n,q_m}}$$

(9)

where ${\Delta }_{{\overrightarrow{d}}_{cu{r}_n},{\overrightarrow{d}}_{p_n,q_m}}$ is the angular change defined in (3), and $\kappa \in (0,1)$ is the weight. Usually, ${\Delta }_{{\overrightarrow{d}}_{cu{r}_n},{\overrightarrow{d}}_{p_n,q_m}}$ is much larger than $\sqrt{{\rho }_{max}/{\overline{\rho }}_m}$; thus, $\kappa$ can be set close to 1 to better balance the two terms. Parameter ${\rho }_{max}=max\left\{{\overline{\rho }}_m\right|q_m\in {\mathcal{Q}}_n^{\prime }\}$ is the maximum density of all the waymark points in the refined candidate set ${\mathcal{Q}}_n^{\prime }$. Based on the weighted distance, the next waymark point is then selected as the one with the smallest weighted distance, i.e.,

$$p_{n+1}=arg\underset{q_m\in {\mathcal{Q}}_n^{\prime }}{min}Dis(p_n,q_m).$$

(10)

The weighted distance defined in (9) does not consider the geo-distance between two waymark points. This is due to the fact that the refined candidate set ${\mathcal{Q}}_n^{\prime }\subseteq {\mathcal{Q}}_n$ contains waymark points that are sufficiently close to the current reference point, as shown in (4). In this case, it is sufficient to focus on the angular and density differences to identify the next waymark point. The proposed trajectory prediction algorithm is finally summarized in Algorithm 1:

Algorithm 1 Trajectory prediction algorithm

Require: $p_{DES}$, $T_0=\{\dots ,p_{-1},p_0\}$, $\mathcal{Q}=\{q_1,q_2,\dots ,q_M\}$

Parameters: ${\epsilon }_{DES}$, ${\epsilon }_{\lambda }$, ${\epsilon }_{\varphi }$, ${\epsilon }_{{\delta }_k}$, ${\epsilon }_{{\delta }_k^{\prime }}$, $\kappa$, ${ϵ}_m,{ϵ}_m^{\prime }$

Initialization: $T=\left\{p_0\right\}$

while  $dist(p_n,p_{DES})\ge {\epsilon }_{DES}$ do

${\mathcal{Q}}_n=\left\{\left({\overline{\lambda }}_m,{\overline{\varphi }}_m,{\overline{\theta }}_m,{\overline{\rho }}_m\right)\right|{\overline{\lambda }}_m\in \left({\overline{\lambda }}_n-{\epsilon }_{\lambda },{\overline{\lambda }}_n+{\epsilon }_{\lambda }\right),$

${\overline{\varphi }}_m\in \left({\overline{\varphi }}_n-{\epsilon }_{\varphi },{\overline{\varphi }}_n+{\epsilon }_{\varphi }\right)\}$

${\mathcal{Q}}_n^{\prime }=\left\{q_m|q_m\in Q_n,{\delta }_m\le {ϵ}_m,{\delta }_m^{\prime }\le {ϵ}_m^{\prime }\right\}$

for all $q_m\mathbf{in}{\mathcal{Q}}_n^{\prime }$ do

$Dis\left(p_n,q_m\right)=\kappa \sqrt{{\rho }_{max}/{\overline{\rho }}_m}+(1-\kappa ){\Delta }_{{\overrightarrow{d}}_{cu{r}_n},{\overrightarrow{d}}_{p_n,q_m}}$

end for

$p_{n+1}=arg{min}_{q_m\in {\mathcal{Q}}_n^{\prime }}Dis(p_n,q_m)$

Append  $p_{n+1}$ to T

end while

return T

4. Experimental Results

In this section, we present the experimental results that evaluate the performance of the proposed trajectory prediction method. We start this section by introducing the dataset and the pre-processing techniques used to process the raw AIS data. The results obtained from key point clustering are then presented to illustrate the operation of the proposed method, followed by the evaluations of trajectory prediction accuracies.

4.1. Dataset and Pre-Processing

In the experiment, we evaluate the proposed trajectory prediction method using the historical AIS data provided by the Danish Maritime Authority (DMA) [10]. The raw data contain AIS messages from various types of vessels traveling around Danish waters. The AIS data have 26 fields, including latitude, longitude, COG, destination, etc. In these experiments, we will mainly use latitude, longitude and COG fields for trajectory prediction. A comprehensive list and description of other fields can be found in Appendix A of [23]. The experiment uses the AIS data from oil tankers traveling to Port Skagen during the period between 2017 and 2019. This dataset contains the AIS data for 1278 vessels and has over 40 million AIS records. In this dataset, the data between 1 January 2017 and 30 June 2019 are used to extract the training trajectories and produce the waymark points. The AIS data from oil tankers traveling to the same port from 1 July 2019 to 31 December 2019 are used to test the proposed trajectory prediction algorithm.

The raw data are pre-processed before being used to produce the waymark points. Key steps of the pre-processing include trajectory extraction and integrity check and downsampling. We refer to Appendix B of [23] for a detailed description of the pre-processing steps. In the experiments, a total of 1046 training trajectories were obtained after the pre-processing, and there are 232 trajectories in the test period. Each trajectory undergoes a route to the destination, i.e., Port Skagen.

4.2. Results for Grid-Based Key Point Clustering

In the training dataset, there are over 550,000 key points. During the key point clustering, the difference in latitude/longitude between the center of the adjacent grid areas is set to 0.05 degrees, while the width of the grid area is 0.05 degrees for both the latitude and the longitude. The parameters for DBSCAN are listed in

Table 2. Key point clustering: DBSCAN parameter settings.

Minimum Number of Points				$\mathit{ϵ}$
COG clustering	key points n	$100\le n$	40	1
		$10\le n<100$	5
		$n<10$	$\lceil n/3\rceil$

.

Figure 4 plots all the key points from the training trajectories, from which the results show that there are some relatively clear routes in the eastern part of the area. In contrast, the routes are much less clear in the western part of the area. After DBSCAN clustering, a total of 6127 waymark points are obtained, which are plotted in Figure 5. Comparing Figure 4 to Figure 5, the results show that although the number of waymark points is only about 1.1% of the key points, they capture all the main route patterns in both the western and the eastern part of the area.

Figure 6 shows the empirical cumulative distribution (CDF) of the COG standard deviation (STD) before and after clustering. Before applying the clustering algorithm, each COG STD is estimated based on the key points within the same grid area. After the clustering, the COG STD is for the key points within the same cluster. Therefore, as expected and shown in Figure 6, the clustering significantly reduces the STD of the COG. For example, before applying the COG clustering, the 95% STD value is about $67.8^{°}$. This means that for about 5% of the clusters (before applying the COG clustering), the COG spread measured by the STD is ≥67.8${}^{°}$. This result suggests that although the grid areas are relatively small, the traveling directions of the tankers within the same small area can differ significantly. In contrast, after applying the clustering, the 95% STD of the COG reduces to $2.2^{°}$. In other words, after clustering, the traveling direction of the cluster becomes more consistent, which is more indicative for the prediction of future positions of the vessel of interest.

4.3. Results for Trajectory Prediction

We now evaluate the performance of the proposed trajectory prediction. The parameters for waymark point connection are set as follows. The distance threshold ${\epsilon }_{DES}$ = 20 km, the length and the width of the rectangular region used to identify ${\mathcal{Q}}_n$ were set to ${\epsilon }_{\lambda }=0.2^{°}$, ${\epsilon }_{\phi }=0.2^{°}$, and the weight coefficient of the weighted distance is set to $\kappa =10/11$. In the filtering waypoint step, we have chosen two more lenient values to filter out impossible points, which are ${ϵ}_m^{\prime }={180}^{°}$ and ${ϵ}_m={120}^{°}$.

Figure 7 shows the number of test samples versus the prediction horizon measured by the number of hours. It can be seen that the test samples contain a mixed of short-term and long-term trajectory prediction tests. About 53.5% of the tests correspond to trajectory predictions in 10 h or more.

To test the proposed algorithm, we implemented a trajectory prediction algorithm based on graph theory and the Dijkstra algorithm [18] and use it as the baseline. Additionally, we also implemented a straightforward trajectory prediction algorithm that simply selects the training trajectory that most matches the historical segment of the trajectory under prediction. This baseline is termed best-matched. Directed Hausdorff distance is adopted to measure the similarity. Figure 8 presents four examples of the trajectory prediction results of the proposed algorithm (the red curves), the graph-based baseline from [18] (the greed curves), and the baseline of the best matched (the orange curves). The ground truth trajectory is also plotted as a reference in each example (the grey solid curves). In the figures, the grey dotted curves are the historical segment of the trajectory prior to the prediction moment. The results show that while the baseline algorithms can provide trajectory predictions fairly close to the ground truth, the smoothness of the predicted trajectories is not good enough. As a comparison, the proposed algorithm can offer more accurate predictions in the sense that the predicted trajectories are visually closer to the true ones and retain the smoothness of realistic trajectories.

From a quantitative perspective, we further evaluate the similarities between the predicated and the ground truth trajectories for the proposed algorithm and the baseline algorithms. To measure the similarity, we consider the DTW distance [24,25,26,27] and the Hausdorff distance [28], which are commonly used to measure the similarities of two curves/sequences of data. A smaller value of this metric suggests a higher degree of similarity and thus a more accurate trajectory prediction.

As the area of interest contains different routes to Port Skagen from different directions, to better demonstrate the details of the performance evaluation, we classified the trajectories into eight trajectory clusters, as shown in Figure 9, for the tested data.

Table 3. Performance comparison between different trajectory prediction algorithms based on DTW and Hausdorff distance.

Cluster	Sample Size	DTW Distance (Unit: km)			Hausdorff Distance (Unit: km)
		Proposed	Best-Matched	Graph-Based	Proposed	Best-Matched	Graph-Based
0	137	555.68	676.44	772.76	8.93	9.96	11.68
1	1734	759.10	908.50	1181.15	11.68	16.04	19.06
2	903	1872.65	2086.40	2221.03	13.46	15.80	15.16
3	13	201.43	218.75	421.26	5.10	8.18	9.86
4	105	1475.34	1552.61	1668.01	11.80	10.42	12.79
5	162	1075.40	1050.57	1204.27	13.61	14.70	15.07
6	115	522.76	548.78	943.52	7.83	8.85	14.81
7	8	167.82	408.43	359.49	3.28	8.01	14.45
Average	3177	1094.31	1244.72	1462.59	11.98	15.14	17.02
Percentage			−12.08%	−25.18%		−11.05%	−42.04%

shows the number of test instances for each cluster in the experiment and the corresponding average DTW and HD. Each test instance refers to a time instance of a trajectory in the test dataset. Since each track lasts more than 25 h on average, the number of test instances is significantly larger than the number of trajectories.

From Table 3, it can be seen that the proposed algorithm outperforms the baseline algorithms on both the performance metrics and for all the eight trajectory classes. Specifically, there is a 12.08% and 25.18% improvement in the DTW distance compared to the best-history and the graph-based baselines, respectively. For Hausdorff distance, the improvement is 11.05% and 42.04% over the two baselines. It can also be seen that for trajectory clusters 1, 6 and 7, which locate on the eastern side of the sea, the advantages of the proposed algorithm over the graph-based approach are even more significant due to the complexity of the shipping lanes and the extremely narrow passable channels.

To further show the stability of the performance of the algorithms, box plots of the performance metrics of the three algorithms are presented in Figure 10 and Figure 11, where the orange line indicates the median, and the green triangle is the location of the mean. It can be seen that for both metrics, the proposed algorithm has a lower median and 75% quantile of the mean than the baseline algorithm, demonstrating higher robustness.

5. Conclusions and Discussion

In this paper, a long-term trajectory prediction algorithm has been developed for oil tankers traveling to a known destination port. The proposed algorithm utilizes key points from historical training trajectories that are extracted from historical AIS data in an area of interest. The algorithm then applies the DBSCAN clustering algorithm to generate a set of waymark points that are much fewer than the key points but still retain all the traveling patterns of the oil tankers. Based on the waymark points, a novel path-finding algorithm has been developed that sequentially identifies a set of waymark points to form the predicted trajectory. The proposed algorithm was tested on real AIS data for oil tankers in Danish waters with a fixed destination of port Skagen. Experimental results show that compared to existing trajectory prediction algorithms, including the graph-based and the best-matched approach, the proposed method can achieve more accurate trajectory predictions with higher trajectory smoothness. Specifically, measured by DTW distance and Hausdorff distance, the proposed method offers a reduction of 25.18% and 42.04% over the graph-based baseline and 12.08% and 11.05% over the best-matched baseline, respectively.

While the proposed algorithm has been tested using real AIS data, the tested scenario is relatively limited. Further testings with AIS data from different areas, on different types of vessels, and to different destinations are needed to obtain a more comprehensive evaluation of the performance of the algorithm. Additionally, the focus of the current work was only on the prediction of the trajectory path. Further effort is needed to match the positions of the predicted trajectory to timestamps. In this direction, a combination of the current algorithm with machine learning techniques may be needed to better capture the motion dynamics of vessels in different areas.