Extended linear regression model for vessel trajectory prediction with a-priori AIS information

ABSTRACT

As maritime activities increase globally, there is a greater dependency on technology in monitoring, control, and surveillance of vessel activity. One of the most prominent systems for monitoring vessel activity is the Automatic Identification System (AIS). An increase in both vessels fitted with AIS transponders and satellite and terrestrial AIS receivers has resulted in a significant increase in AIS messages received globally. This resultant rich spatial and temporal data source related to vessel activity provides analysts with the ability to perform enhanced vessel movement analytics, of which a pertinent example is the improvement of vessel location predictions. In this paper, we propose a novel strategy for predicting future locations of vessels making use of historic AIS data. The proposed method uses a Linear Regression Model (LRM) and utilizes historic AIS movement data in the form of a-priori generated spatial maps of the course over ground (LRMAC). The LRMAC is an accurate low complexity first-order method that is easy to implement operationally and shows promising results in areas where there is a consistency in the directionality of historic vessel movement. In areas where the historic directionality of vessel movement is diverse, such as areas close to harbors and ports, the LRMAC defaults to the LRM. The proposed LRMAC method is compared to the Single-Point Neighbor Search (SPNS), which is also a first-order method and has a similar level of computational complexity, and for the use case of predicting tanker and cargo vessel trajectories up to 8 hours into the future, the LRMAC showed improved results both in terms of prediction accuracy and execution time.

KEYWORDS:

• Automatic Identification System (AIS) data
• Linear Regression Model (LRM)
• trajectory mining
• spatial map
• historic data
• trajectory prediction

1. Introduction

The world’s oceans are of critical importance to humanity as they are key to fisheries, shipping, and the environment. From an economic perspective, it is estimated that 90% of all global goods and energy transportation are done by sea, with millions of people being dependent on maritime-related activities for their livelihood. As maritime activities increase globally, there is a greater dependency on technology in monitoring, control, and surveillance of vessel activities. One of the most prominent systems for monitoring vessel activity is the Automatic Identification System (AIS). AIS operates in the VHF band and transmits messages from vessels to other vessels, terrestrial shore stations, and satellites. Due to the global increase in vessels fitted with AIS transmitters and the proliferation of satellite and terrestrial receiving stations, there has been a significant increase in AIS messages received globally. This increased data volume makes it possible to track the real-time movement of vessels and opens the door for improving vessel location predictions via historic vessel movement patterns. Many algorithms have been developed in recent years to aid in improved vessel coordinate prediction. Methods range from simplistic models as done by Burger, Kleynhans, and Lups Grobler (2020) to machine learning (ML) models as done by Xin (2020).

Pallotta et al. (2014) presented a vessel prediction method based on Ornstein-Uhlenbeck stochastic processes, where the parameters of these processes are estimated from historical patterns in historic AIS data. The data is clustered into three types: vessels, waypoints, and routes. Route extraction is done using Traffic Route Extraction and Anomaly Detection from AIS Data (THREAD). The implementation details of THREAD can be found in Pallotta, Vespe, and Bryan (2013a, 2013b). The three types of clustered data aid in vessel prediction and empirical calculations. The maximum prediction time window of a vessel depends on the mean duration of the historically observed route.

Lee, Han, and Whang (2007) presented a trajectory clustering method, where similar trajectories are clustered together. The method works by employing a partition-and-group framework for clustering trajectories. It works by first grouping trajectories into a set of line segments, and similar segments are then grouped into clusters. The method consists of two phases, partitioning and grouping. In the first phase, the minimum description length principle is used (Grünwald, Jae Myung, and Pitt 2005). In the second phase, a density-based line-segment clustering algorithm is used (Ester et al. 1996)based on a Density-Based Algorithm for Discovering Clusters in Large Spatial Databases with Noise (DBSCAN). The model’s sensitivity on its parameters is improved in Jiashun (2012).

Rong, Teixeira, and Guedes Soares (2019) presented a probabilistic trajectory prediction model that describes the uncertainty in future spatial locations along a vessel’s trajectory by using continuous probability distributions. A non-parametric Bayesian model based on a Gaussian Process (GP) is presented. The GP is used to describe the uncertainty of the vessel’s motion, while the longitudinal uncertainty is derived from the dubiety of the vessel’s acceleration. The parameters of the model are derived from historic AIS data. Cholesky decomposition is used to reduce the computational complexity of the algorithm. Prior probabilities are used to predict vessel coordinates in real-time.

Jaskolski (2017) implemented a Discrete Kalman Filter (DKF) (Kalman 1960) to predict future locations of vessels. The DKF, in the context of vessel coordinate prediction, constantly adjusts itself for an improved prediction as new observations are observed. It is assumed that a vessel fitted with an AIS sensor will not constantly send updates. The DKF consists of two sets of equations: prediction and measurement update equations used in prediction and parameter updates, respectively. Burger, Kleynhans, and Lups Grobler (2020), showed that there is no significant improvement in the prediction capability that can be achieved by using a DKF instead of a Linear Regression Model (LRM) to predict linear trajectories.

Xiao et al. (2020) conducted an extensive review of maritime knowledge mining and traffic forecasting technologies. The LRM is compared to several non-linear approaches. Three broad categories of non-linear algorithms are considered: machine learning approaches, knowledge-based approaches, and control theory assisted methods. The predictions range from long to short-term predictions. It is also shown that more complex methods are more accurate than the LRM but have a higher computational cost than the LRM.

Forti et al. (2020) made use of a deep learning (DL) neural network (NN) approach to predict trajectories of vessels. A sequence-to-sequence model that utilizes a Long Short-Term Memory (LSTM) encoded-decoder recurrent neural network (RNN) is proposed. Historic AIS data is used to train the LSTM model. The method aims to learn the predictive distribution of maritime traffic patterns using historic AIS data. Learning the predictive patterns enables the model to predict more accurately. It was shown that the model could predict more accurately than the Ornstein-Uhlenbeck process, given a time window of 20 observations.

Murray and Prasad Perera (2020) presented a novel dual linear autoencoder approach to predict a vessel’s trajectory. The method predicts a future trajectory using historic AIS data. The method implements unsupervised learning for trajectory clustering and classification.

Xin (2020) presented a context-based trajectory prediction algorithm utilizing LSTM networks. Real-valued target trajectories are converted into discrete path sets. Distinctive patterns are clustered hierarchically using historical AIS data. Two models are compared, an RNN consisting of one LSTM and another RNN consisting of $\mathit{k}$ LSTMs. In the RNN with $\mathit{k}$ LSTMS, one LSTM is created for each distinct path. Yaun et al. (2019) used an LSTM model to reconstruct vessel trajectories. This complex method requires pre-processing and data clustering.

Dimitrios, Xidias, and Lekkas (2016) aim to accurately predict future geo-coordinates of a vessel by using artificial NNs (ANNs). Their model learns in real-time whilst predicting. It has a prediction time horizon of up to 15 minutes. Different types of model pre-processing and construction are implemented. Historic AIS data was used to train the ANN. Alizadeh, Asghar Alesheikh, and Sharif (2021) proposed three novel prediction methods based on historic AIS data. The first method proposed is a Point-based Similarity Search Prediction (PSSP), which was inspired by Wijaya and Nakamura (2013). The historical points are measured in terms of their spatial location, SOG, and COG. The second method proposed is called Trajectory-based Similarity Search Prediction (TSSP), where each recorded AIS trip is regarded as a trajectory. The PSSP is a point-based method, whereas the TSSP is a trajectory-based method. Finally, a trajectory-based similarity search prediction is proposed using an RNN LSTM (TSSPL). Alizadeh, Asghar Alesheikh, and Sharif (2021) point out that vessel movement is affected by external movements such as wind, waves, and sea currents. The PSSP and the TSSP are not able to account for these external factors. Another RNN LSTM model was, therefore, built to take this into account (i.e. TSSPL). The TSSPL has an additional input, a measure of similarity between trajectories (similar to what was done in Tang, Yin, and Shen 2019).

The Safety of Life at Sea (SOLAS) regulation V/19 states that passenger vessels and vessels larger than 300 Gross Tonnage must have AIS transmitters fitted. All cargo and tanker vessels should adhere to the SOLAS regulations. The regulations were adopted and are now regulated by the International Maritime Organization (IMO).¹ There are, however, some drawbacks when working with AIS data. Vessels may deactivate their AIS transponders, and some vessels may not be fitted with transponders (Chen et al. 2019).

Other vessel tracking information systems include Long-Range Identification and Tracking (LRIT), Coastal radar, Satellite-borne Synthetic Aperture Radar (SAR), and optical satellites. Moreover, methods such as video tracking at ports exist where maritime video surveillance is present, as done in Chen et al. (2019) and Chen et al. (2020).

Finally, the method presented by Hexeberg, Flåten, and Brekke (2017) uses historic AIS data to predict future locations of vessels. The method is called Single-Point Neighbor Search (SPNS). The method does a close neighbor (CN) search by extracting historic observations within a certain radius of the current vessel’s spatial location. Vessels in the CN set that do not adhere to prespecified SOG and COG range values are removed. Using the CN set, the median COG and SOG values are calculated. Using the median COG and SOG, the predicted longitude and latitude are calculated. The method predicts at constant distance intervals, where the SOG is used to calculate the time passed between two observations. The method can be confidently predicted with a time horizon of up to 15 minutes.

The aforementioned methods are complex in nature, being programmatically challenging to implement from first principles. The setup, initialization, and parameter estimation of these methods require fine-tuning to obtain optimal and accurate results. Burger, Kleynhans, and Lups Grobler (2020) showed that the complexity trade-off of using a DKF over an LRM resulted in no significant performance improvement in prediction accuracy, where the DKF was much more complex than the LRM. The DKF parameters are sensitive, and getting the initial set of parameters to be a reasonable estimate is crucial for the model’s performance. Xiao et al. (2020) showed that the increase in performance gains comes with a significant increase in computational cost.

Therefore, this paper proposes a novel vessel location prediction method that extends the Linear Regression Model (LRM) proposed by Burger, Kleynhans, and Lups Grobler (2020). Moreover, this method uses historic AIS vessel information to improve prediction accuracy. It is easy to implement, and its initial parameters require little to no fine-tuning. The proposed algorithm also enables the LRM to predict non-linear trajectories. The a-priori information used includes a Spatial Map (SM) of historic cell counts, COG values, and COG standard deviations (SDs). The SMs the method uses are easy to generate and update, which is essential as data and movement patterns of vessels change over time. We refer to this method as the LRM with added COG (LRMAC) throughout the paper.

Due to the fact that the SPNS is also a simplistic first-order method (proposed by Hexeberg, Flåten, and Brekke 2017), the LRMAC will be compared to it. When comparing the models, an in-depth comparison of the model accuracies and time complexities will be made. Our tests will be implemented on cargo and tanker vessel data only. This paper is structured as follows: first, we discuss the dataset and the required processing. Second, we summarize the a-priori information the LRMAC requires. Third, we present the LRMAC and the comparison method, the SPNS. Fourth, we discuss the experimental design. Finally, we end the paper with our results and a conclusion.

2. Data

In this section, we will introduce the dataset that was used to test the performance of the proposed and benchmark algorithms. We will also be discussing the pre-processing steps that are necessary to ensure the proposed algorithm will run as intended.

An open-source dataset by Ray et al. (2019) was used. It consists of AIS messages recorded in the Celtic Sea, the North Atlantic Ocean, the English Channel, and the Bay of Biscay (France). In Figure 1, the spatial range of the data is depicted visually (also see the details in Table 1). The specific data characteristics are given in Table 2, and the observational period was over 6 months starting on 1 October 2015.

Figure 1. Data Spatial range – courtesy of Google Earth 2021.

Table 1. Spatial range details for Google Earth extract.

Name	Abbreviation	Longitude	Latitude
Upper Left-Hand Corner	ULHC	51	10
Upper Right-Hand Corner	URHC	51	0
Lower Right-Hand Corner	LRHC	45	0
Lower Left-Hand Corner	LLHC	45	10

Table 2. Dataset attributes.

Description	Measurement Unit	Attribute Range
MMSI		9-digit values
Latitude	DD.dddd	[−10.00  , 0.00]
Longitude	DD.dddd	[45.00  , 51.00]
Course Over Ground	Degrees	0–360
Speed Over Ground	Knots	0–110
Timestamp	UTC	[2015–10-01 00:00:00, 2016–03-31 23:59:59]
Ship type		[10, 99]

In Table 2, the dataset attributes and characteristics are introduced. The mathematical symbols that we use in this paper are presented in Table 3.

Table 3. Variables and Symbol Descriptions.

Symbol	Meaning	Type	Dimensions
$\mathit{t}$	Time (DD-MM-YYYY HH:MM:SS)	Scalar
$\mathrm{\Delta }{k}_t$	Prediction time interval size (s)	Scalar
$\mathit{\varphi }$	Latitude (DD.dddd)	Scalar
$\mathit{\lambda }$	Longitude (DD.dddd)	Scalar
$\mathit{\psi }$	Course Over Ground (°)	Scalar
$\mathit{V}$	Speed Over Ground (kt)	Scalar
$\mathit{V}{ }^{\mathrm{\prime }}$	Speed Over Ground in Degrees	Scalar
$\mathit{\omega }$	Window Size	Scalar
$\mathit{\kappa }$	SM cell length	Scalar
$\mathit{\xi }$	Single SM cell	Scalar
$\mathit{\eta }$	Neighborhood search size	Scalar
$\mathit{H}$	Neighborhood grid index values	Matrix	$\left(2\mathit{\eta }+1\right)\times \left(2\mathit{\eta }+1\right)$
$\mathit{K}$	Vessel counts per cell SDM Matrix	Matrix	$1250\times 1250$
$\mathbf{\Psi }$	Course Over Ground SM Matrix (°)	Matrix	$1250\times 1250$
$\mathrm{\Sigma }$	Course Over Ground SD SM Matrix (°)	Matrix	$1250\times 1250$
$i_{\mathit{\lambda }}$	Denotes the Longitude index of a SM matrix	Scalar
$i_{\varphi }$	Denotes the Latitude index of a SM matrix	Scalar

2.1. Pre-processing steps

Observations were removed from the dataset that did not adhere to the following criteria:

• Ship type within [70,89], where cargo vessel’s type is in [70,79] and tankers in [80,89].

• SOG > 0.5kn, removing stationary observations, including stationary vessels experiencing drift due to currents and other natural phenomena due to being anchored.

• SOG < 60kn, observations with high speeds are likely outliers as cargo and tanker vessels move at relatively low speeds. If more than 60 kn is observed for cargo or tanker vessels, it is most likely due to technical errors in the recorded AIS data.

The remainder of the data were grouped according to vessel MMSI and sorted in ascending order according to the timestamp recorded by each observation. All trajectories with less than 20 observations or those that span less than 5 min in total were removed. The data spans over a period of 6 months, as denoted in Table 2, which implies that there will be more than one trajectory for a given vessel at different time periods and spatial locations. Pre-processing was done to remove any non-sensical data that may lead to incorrect SMs.

2.2. Vessel Interpolation for Spatial Maps (SMs)

After pre-processing, an augmented dataset was constructed by interpolating between observations for each unique vessel MMSI, creating less sparse trajectory observations per vessel. This augmented dataset is then used to construct physically realistic and usable spatial maps (SMs) of the dataset in question. The original dataset is too sparse (as is) for a meaningful and sensical SM to be created from it. The sparsity of the original data is due to several reasons. Firstly, vessels closer to receivers can share their position more often than vessels further away. Furthermore, vessels close to the shoreline travel at lower speeds than vessels further away. Naturally, more AIS messages will be recorded for slower vessels than faster vessels over the same distance traveled.

Observational interpolation was only performed in the following cases:

• The time difference between the two observations is no longer than six hours,

• The distance between the two observations is within 15 km, and

• The distance between observations is no smaller than the size of one grid cell. The grid cell size used was $0.88\mathrm{km}\times 0.88\mathrm{km}$. This specific constraint prevents the over-representation of a grid cell, only adding one observation to a cell if the interpolated trajectory passed through it.

We made use of linear interpolation models by using scikit-learn (Pedregosa et al. 2011). Longitude, latitude, and SOG were interpolated if the cases above were met. Gaussian filtering was also used to have smoother versions of the SMs.

The COG of the interpolated and recorded observations was calculated via

(1) ${\psi }_t=\mathit{arctan}\left(\frac{{\mathit{\varphi }}_{\mathit{t}-1}-{\varphi }_t}{{\mathit{\lambda }}_{\mathit{t}-1}-{\mathit{\lambda }}_t}\right)$(1)

We did not make use of the COG values present in the dataset. These values yielded inaccurate results, more accurate results were obtained by calculating the actual COG based on the longitude and latitude. Having a new augmented dataset with a better representation of the historical locations and trajectories, we can create representative SMs, essential to the algorithm we propose.

3. Spatial Maps (SMs)

Any two-dimensional grid that spans the earth’s surface will be referred to as a spatial map. Each cell in the two-dimensional grid is associated with a specific range of longitudinal and latitudinal coordinates. In this paper, our SMs will be set up from the observations recorded in the dataset, which spans a latitude and longitude of $\mathit{\varphi }\in \left[-{10}^{\circ },0^{\circ }\right]$ and $\mathit{\lambda }\in \left[{45}^{\circ },{55}^{\circ }\right]$, respectively.

In this paper, we use a square SM, meaning the SM’s width and length have the same number of cells. The dimensions of the SMs depend on the range of the $\mathit{\varphi }$ and $\mathit{\lambda }$. The dimensions of the SMs we constructed were $1250\times 1250$ cells as shown in Table 3. The size of the SMs is ${10}^{\circ }\times {10}^{\circ }$ square degrees. Each cell’s resolution, therefore, is equal to $0.008\times 0.008$ square degrees. An extract of the upper left-hand corner of an artificially created SM is depicted in Figure 2. This figure denotes the number of recorded observations in a cell, within a longitude and latitude range based on the historic AIS data.

Figure 2. Spatial distribution map matrix extract example.

Haversine’s formula for distance on a sphere is given by Equation (2)(2) $\begin{array}{rl}& \mathit{d}=2\mathit{r}\cdot \mathit{arcsin}\\ & \sqrt{\mathit{si}{n}^2\left(\frac{{\varphi }_2-{\varphi }_1}{2}\right)+\mathit{cos}\left({\varphi }_1\right)\mathit{cos}\left({\varphi }_2\right){\mathit{sin}}^2\left(\frac{{\lambda }_2-{\lambda }_1}{2}\right)}\end{array}$(2) . This formula is used to calculate the distance between two coordinates on a sphere.

(2) $\begin{array}{rl}& \mathit{d}=2\mathit{r}\cdot \mathit{arcsin}\\ & \sqrt{\mathit{si}{n}^2\left(\frac{{\varphi }_2-{\varphi }_1}{2}\right)+\mathit{cos}\left({\varphi }_1\right)\mathit{cos}\left({\varphi }_2\right){\mathit{sin}}^2\left(\frac{{\lambda }_2-{\lambda }_1}{2}\right)}\end{array}$(2)

where,

• $\mathit{r}$ represents the mean radius of the earth ($\approx 6371\mathrm{km}$).

• ${\varphi }_1$ and ${\varphi }_2$ represent the latitudinal coordinates of two observations, the latitudes of points one and two, respectively.

• ${\mathit{\lambda }}_1$ and ${\mathit{\lambda }}_2$ represent the longitudinal coordinates of two observations, the longitudes of points one and two, respectively.

It now follows from the Haversines formula that the area associated with each cell is roughly $0.89\times 0.89$ km², where $0.89$ km $\approx 0.48$ nautical miles.

3.1. Vessel counts SDM ($\mathit{K}$)

The first type of SM that we use represents the number of observations recorded within each grid cell $\mathit{\xi }$. Each $\mathit{\xi }$ in $\mathit{K}$ records the number of observations that were historically within the longitude and latitude range of that specific cell. Cell counts will be higher in spatial areas where vessels are slower moving compared to faster-moving areas. Figure 3 shows the logarithmically scaled SDM of the vessel counts for each cell. The cell counts of the current vessel location and its neighbors are used in the prediction equations of the LRMAC, scaling the contribution of historic COG information.

Figure 3. Log scale vessel counts SDM.

3.2. The COG $\left(\mathbf{\Psi }\right)\mathit{\hspace{0.17em}}$ and COGSD $\left(\mathbf{\Sigma }\right)$ SM

The second SM that we will be using will be the course over ground SM (represented by $\mathbf{\Psi }$). Each cell of the $\mathbf{\Psi }$ represent the mean COG value that was recorded in that cell. COG is measured in degrees and will always be positive.

The third SM we used is the so-called COG standard deviation SM (represented by $\mathbf{\Sigma }$). The value of each of its cells was computed as follows:

(3) ${\mathbf{\Sigma }}_{{\mathrm{i}}_{\varphi },{\mathrm{i}}_{\mathit{\lambda }}}=\sqrt{\frac{1}{n}\sum _{\mathit{j}=1}^n{\left({\psi }_{i_{\varphi },{{\mathit{i}}_{\mathit{\lambda }}}_{\mathit{j}}}-{\mathbf{\Psi }}_{i_{\varphi },i_{\mathit{\lambda }}}\right)}^2}$(3)

where

• $\mathit{n}$, represents the number of observations observed in a cell as determined by $\mathit{K}$.

• $i_{\varphi },i_{\mathit{\lambda }}$, represents the index values associated with $\mathit{\varphi }$ and $\mathit{\lambda }$ respectively on the SM grid.

• ${\mathbf{\Psi }}_{i_{\varphi },i_{\mathit{\lambda }}}$, represents the mean COG at a specific index value on the SM grid.

• ${{\psi }_{(i_{\varphi },i_{\mathit{\lambda }})}}_{\mathit{j}}$, represents the $j^{\mathit{th}}$ COG value in the cell with index values $i_{\varphi },i_{\mathit{\lambda }}$

• ${\mathbf{\Sigma }}_{i_{\varphi },i_{\lambda }}$, refers to the COG standard deviation at index value $i_{\varphi },i_{\mathit{\lambda }}$.

Loosely speaking, the entries in $\mathbf{\Sigma }$ can be interpreted as how “confident” we ought to be in the corresponding entry in $\mathbf{\Psi }$. Higher cell values in $\mathbf{\Sigma }$ mean that historically there were many vessels traveling in different directions as the SD is higher.

In Figure 4, we can see a visual representation of $\mathbf{\Psi }$. Looking at the figure, we can see that cargo and tanker vessels move in a specific direction in certain geographic locations. Two distinct highways are clearly visible close to the center of Figure 4 (Grobler and Kleynhans 2019). A highway is a route that many vessels traverse. The aqua green highway is used by cargo and tankers to travel upward (North), while the blue highway is used to travel downward (South).

Figure 4. Course over ground SM.

In Figure 5, we see a visual representation of $\mathbf{\Sigma }$, as calculated by Equation (3)(3) ${\mathbf{\Sigma }}_{{\mathrm{i}}_{\varphi },{\mathrm{i}}_{\mathit{\lambda }}}=\sqrt{\frac{1}{n}\sum _{\mathit{j}=1}^n{\left({\psi }_{i_{\varphi },{{\mathit{i}}_{\mathit{\lambda }}}_{\mathit{j}}}-{\mathbf{\Psi }}_{i_{\varphi },i_{\mathit{\lambda }}}\right)}^2}$(3) . The standard deviation of the areas that contain more traffic in different directions is larger than those containing less traffic (this is especially true for the areas surrounding harbors). Yellowish colors represent higher SD values. The highways mentioned, however, are not associated with high SD values. This implies that these highways are highly directional as shown in Figure 5.

Figure 5. Course over ground standard deviation SM.

4. The proposed algorithm

In this section, we present the proposed LRMAC algorithm. As mentioned before, this algorithm uses historic (a-priori) AIS data to predict the trajectories of vessels. The algorithm can predict non-linear vessel trajectories. To be more specific, the LRMAC uses three SMs discussed in the previous section $\mathit{K},\mathit{\hspace{0.17em}}\mathbf{\Psi },$ and $\mathbf{\Sigma }$. We start this section by discussing the different unit conversions that the algorithm requires.

4.1. Speed Over Ground (SOG)

We converted SOG from knots to meters per second via Equation (4)(4) $V_t^{\prime }=0.51\dot{4}\cdot {\mathrm{V}}_{\mathrm{t}}^{{ }^{\prime \prime }}$(4) .

(4) $V_t^{\prime }=0.51\dot{4}\cdot {\mathrm{V}}_{\mathrm{t}}^{{ }^{\prime \prime }}$(4)

We then converted the SOG into degrees/s (°/s) using

(5) ${\mathrm{V}}_{\mathrm{t}}=\frac{V_t^{\prime }}{\mathit{l}}$(5)

with

(6) $\mathit{l}=\frac{2\mathit{\pi }}{360}\times 6378000=111137\mathrm{m}$(6)

The constant $\mathit{l}$, can be interpreted as the average number of meters that one degree of latitude and longitude span on Earth.

4.2. The Linear Regression Model (LRM)

The method presented in this paper (LRMAC) is an extension of the LRM presented by Burger, Kleynhans, and Lups Grobler (2020). We will summarize it below for the benefit of the reader. The LRM can only be predicted linearly with regular time intervals. The prediction interval that we used was 1 s, i.e. $\mathrm{\Delta }{k}_t=1\mathit{s}$. The LRM-based approach works by estimating $V_t$ using a rolling window linear model.

One of the parameters that must be set is the prediction duration, i.e. the total number of hours one would like to predict into the future. The other parameter will be the step size of the predictions, measured in seconds between any two observations. The prediction interval can then be computed.

The LRM consists of two sets of equations: the predictor and measurement update equations. Let the vector ${\mathit{x}}_t$ be defined as:

(7) ${\mathit{x}}_t={\left[{\mathit{\lambda }}_t,{\varphi }_t\right]}^{\mathit{T}}$(7)

The symbols $\mathit{\varphi }$ and $\mathit{\lambda }$ are defined in Table 3. Subscript $\mathit{t}$ refers to the observations at timestep $\mathit{t}.$ If there is no observation at timestep $\mathit{t},{\mathit{x}}_t$ is assumed to be equal to the all-zero vector.

The prediction equations for the LRM are

(8) ${\hat{\mathit{x}}}_{\mathit{t}}^{-}={\hat{\mathit{x}}}_{\mathit{t}-1}+{\hat{V}}_{\mathit{\omega },\mathit{t}}\cdot {\mathbf{\Lambda }}_{\mathrm{t}}$(8)

(9) ${\hat{V}}_{\mathit{\omega },\mathit{t}}=\mathrm{\Delta }{\hat{V}}_{\mathit{\omega },\mathit{t}-1}\cdot k_t+{\hat{V}}_{\mathit{\omega },c_{\mathit{t}-1}}$(9)

(10) ${\mathbf{\Lambda }}_{\mathrm{t}}={\left[cos\left({\psi }_t\right)sin\left({\psi }_t\right)\right]}^{\mathit{T}}$(10)

The variables in the equations above are defined below,

• ${\hat{\mathit{x}}}_{\mathit{t}}^{}-$ denotes the vessel’s predicted position vector using all observations up until timestep $\mathit{t}-1.$

• ${\hat{\mathit{x}}}_{\mathit{t}-1}$ denotes the vessel’s updated estimated position vector using all observations up until timestep $\mathit{t}-1$.

• ${\hat{V}}_{\mathit{\omega },\mathit{t}}$ denotes the vessels predicted SOG using all observations up until timestep $\mathit{t}-1$.

• $\mathrm{\nabla }{\hat{V}}_{\mathit{\omega },\mathit{t}-1}$ denotes the vessel’s updated estimated gradient of our LRM using all observations up until timestep $\mathit{t}-1$.

• ${\hat{V}}_{\mathit{\omega },c_{\mathit{t}-1}}$ denotes the updated estimated $\mathit{y}$-intercept of our LRM using all observations up until timestep $\mathit{t}-1$.

• ${\psi }_t$, denotes the COG of a vessel at timestep $\mathit{t}$, the COG remains the same until we observe a new COG value, i.e. getting an update from the vessel.

• $\mathit{\omega }$, denotes the window size of our LRM.

• $k_t$, denotes the elapsed time in seconds at timestep $\mathit{t}$.

The measurement update equations are

(11) ${\hat{\mathit{x}}}_{\mathit{t}}={\hat{\mathit{x}}}_{\mathit{t}}^{-}+\left({\mathit{x}}_t-{\hat{\mathit{x}}}_{\mathit{t}}^{-}\right)={\mathit{x}}_t$(11)

(12) $\mathrm{\nabla }{\hat{V}}_{\mathit{\omega },\mathit{t}}=\frac{{\sum }_{\mathit{i}={\text{b}1}_{n_t>\mathit{\omega }}\left(n_t-\mathit{\omega }\right)}^{n_t}\left(V_i-{\bar{V}}^{\mathit{\omega }}\right)\left(k_i-{\bar{k}}^{\mathit{\omega }}\right)}{{\sum }_{i=1_{n_t>\mathit{\omega }}\left(n_t-\mathit{\omega }\right)}^{n_t}\left(k_i-{\bar{k}}^{\mathit{\omega }}\right)}$(12)

(13) ${\hat{V}}_{\mathit{\omega },c_t}={\bar{V}}^{\mathit{\omega }}-\mathrm{\nabla }{\hat{V}}_{\mathit{\omega },\mathit{t}}\cdot {\bar{k}}^{\mathit{\omega }}$(13)

Moreover

(14) $1_{n_t>\mathit{\omega }}=\left\{\begin{array}{c}1,\mathrm{if}{n}_t>\mathit{\omega }\\ 0,\mathrm{otherwise}\end{array}\right.$(14)

(15) ${\bar{V}}^{\mathit{\omega }}=\frac{1}{{\mathit{n}}_{\mathit{\omega }}}\sum _{i=1_{n_t>\mathit{\omega }}\left(n_t-\mathit{\omega }\right)}^{n_t}{V}_i$(15)

(16) ${\bar{k}}^{\mathit{\omega }}=\frac{1}{{\mathit{n}}_{\mathit{\omega }}}\sum _{i=1_{n_t>\mathit{\omega }}\left(n_t-\mathit{\omega }\right)}^{n_t}{k}_i$(16)

(17) $n_{\mathit{\omega }}=\left\{\begin{array}{c}{n}_t,\mathrm{if}{n}_t<\mathit{\omega }\\ \mathit{\omega },\mathrm{otherwise}\end{array}\right.$(17)

The variables in the equations above are defined below:

• $V_i$ denotes the $i^{\mathit{th}}$ true SOG observation that was recorded for a particular vessel. We assumed $V_0=4\mathrm{m}/\mathrm{s}$.

• $k_i$ denotes the total time that has elapsed after having the $i^{\mathit{th}}$ true observation.

• $n_t$ denotes the total number of true observations recorded after $\mathrm{\Delta }{k}_t\cdot \mathit{t}$ seconds, where in this paper $\mathrm{\Delta }{k}_t=1\mathrm{s}$. $\mathrm{\Delta }{k}_t$ refers to the prediction time interval between two subsequent observations.

• Note that subscript $\mathit{t}$ is used as a time step index, while subscript $\mathit{i}$ is used as an observational index.

The measurement update equation values are only updated once a new observation is recorded. If there are no new observations recorded, the algorithm will evaluate the prediction equations. During the period of no observations, it is assumed that the COG remains unchanged. Moreover, the predicted position vector at $\mathit{t}-1$ is used to obtain a new estimate of the position vector at $\mathit{t}$.

4.3. LRM with a-priori COG information (LRMAC)

In this section, we will present our novel prediction method that is an extension of the LRM. The proposed method uses spatial maps as they can easily be loaded into memory, and their sizes are relatively small. The proposed method will allow for non-linear trajectory predictions of vessels. Note that the symbols used in this section should not be confused with those used in the previous section for the LRM.

The LRMAC algorithm no longer assumes a constant COG. The value of the COG is dynamically updated using a-priori information, using the three SM matrices $\mathit{K}$, $\mathbf{\Psi }$ and $\mathbf{\Sigma }$. As the COG value is dynamically updated, the LRMAC can predict non-linear trajectories (i.e. the SOG in the longitudinal and latitudinal directions will change). It is assumed that the SOG will remain constant over the prediction period, where the constant SOG used during prediction is derived from the last $\mathit{\omega }$ observed SOG values of the vessel under consideration.

The predictor equations are modified for the LRMAC compared to the LRM. Whilst in prediction mode, the COG value is updated via a-priori information. The COG value computed in step $\mathit{t}$ is only applied during step $\mathit{t}+1$.

Let $\widehat{x}\stackrel{-}{{\mathit{ }}_t}$ be the predicted position vector as in Equation (8)(8) ${\hat{\mathit{x}}}_{\mathit{t}}^{-}={\hat{\mathit{x}}}_{\mathit{t}-1}+{\hat{V}}_{\mathit{\omega },\mathit{t}}\cdot {\mathbf{\Lambda }}_{\mathrm{t}}$(8) . Let,

(18) ${\mathit{n}}_{\hat{\varphi },\hat{\lambda }}=\left[n_{\hat{\varphi }},n_{\hat{\lambda }}\right]$(18)

denote the SM index positions associated with ($\mathit{\varphi }$,$\mathit{\lambda }$). In other words, with $\left[n_{\hat{\varphi }},n_{\hat{\lambda }}\right]$ we can extract the elements in $\mathit{K}$, $\mathbf{\Psi }$ and $\mathbf{\Sigma }$ associated with ($\mathit{\varphi }$, $\mathit{\lambda }$) by making use of matrix subscripting. Let us now construct the following index matrix H:

(19) $\mathit{H}=\left[\begin{array}{c}{\mathit{n}}_{\hat{\varphi }-\mathit{\eta }\cdot \mathit{\kappa },\hat{\lambda }-\mathit{\eta }\cdot \mathit{\kappa }}\dots {\mathit{n}}_{\hat{\varphi }-\mathit{\eta }\cdot \mathit{\kappa },\hat{\lambda }}\dots {\mathit{n}}_{\hat{\varphi }-\mathit{\eta }\cdot \mathit{\kappa },\hat{\lambda }+\mathit{\eta }\cdot \mathit{\kappa }}\\ ⋮\ddots ⋮\text{.}\text{.}\text{.}⋮\\ {\mathit{n}}_{\hat{\varphi }\cdot \mathit{\kappa },\hat{\lambda }-\mathit{\eta }\cdot \mathit{\kappa }}\dots {\mathit{n}}_{\hat{\varphi },\hat{\lambda }}\dots {\mathit{n}}_{\hat{\varphi },\hat{\lambda }+\mathit{\eta }\cdot \mathit{\kappa }}\\ ⋮\text{.}\text{.}\text{.}⋮\ddots ⋮\\ {\mathit{n}}_{\hat{\varphi }+\mathit{\eta }\cdot \mathit{\kappa },\hat{\lambda }-\mathit{\eta }\cdot \mathit{\kappa }\dots {{\displaystyle \mathit{n}}}_{\hat{\varphi }-\mathit{\eta }\cdot \mathit{\kappa },\hat{\lambda }}\dots {{\displaystyle \mathit{n}}}_{\hat{\varphi }+\mathit{\eta }\cdot \mathit{\kappa },\hat{\lambda }+\mathit{\eta }\cdot \mathit{\kappa }}}\end{array}\right]$(19)

where

• $\mathit{\eta }$ in Equation (19)(19) $\mathit{H}=\left[\begin{array}{c}{\mathit{n}}_{\hat{\varphi }-\mathit{\eta }\cdot \mathit{\kappa },\hat{\lambda }-\mathit{\eta }\cdot \mathit{\kappa }}\dots {\mathit{n}}_{\hat{\varphi }-\mathit{\eta }\cdot \mathit{\kappa },\hat{\lambda }}\dots {\mathit{n}}_{\hat{\varphi }-\mathit{\eta }\cdot \mathit{\kappa },\hat{\lambda }+\mathit{\eta }\cdot \mathit{\kappa }}\\ ⋮\ddots ⋮\text{.}\text{.}\text{.}⋮\\ {\mathit{n}}_{\hat{\varphi }\cdot \mathit{\kappa },\hat{\lambda }-\mathit{\eta }\cdot \mathit{\kappa }}\dots {\mathit{n}}_{\hat{\varphi },\hat{\lambda }}\dots {\mathit{n}}_{\hat{\varphi },\hat{\lambda }+\mathit{\eta }\cdot \mathit{\kappa }}\\ ⋮\text{.}\text{.}\text{.}⋮\ddots ⋮\\ {\mathit{n}}_{\hat{\varphi }+\mathit{\eta }\cdot \mathit{\kappa },\hat{\lambda }-\mathit{\eta }\cdot \mathit{\kappa }\dots {{\displaystyle \mathit{n}}}_{\hat{\varphi }-\mathit{\eta }\cdot \mathit{\kappa },\hat{\lambda }}\dots {{\displaystyle \mathit{n}}}_{\hat{\varphi }+\mathit{\eta }\cdot \mathit{\kappa },\hat{\lambda }+\mathit{\eta }\cdot \mathit{\kappa }}}\end{array}\right]$(19) denotes the neighborhood parameter, and

• $\mathit{\kappa }$ represents the width and length of an SM cell (see Table 3 and Figure 2).

The index matrix $\mathit{H}$ is used to select a specific sub-grid/matrix from $\mathit{K}$ and $\text{b}\mathrm{\Sigma }$.

Let the a-priori cell counts, of the area surrounding ${\hat{\mathit{x}}}_{\mathit{t}}^{-}$, therefore, be denoted by

${\mathit{K}}_H,\mathrm{where}{\mathit{K}}_{\mathit{H}}\subset \mathit{K}$

The size of the aforementioned area is determined by $\mathit{\eta }$. Moreover, let the a-priori average COG value associated with ${\hat{\mathit{x}}}_{\mathit{t}}^{-}$, be denoted by

${\mathbf{\Psi }}_{{\mathit{n}}_{\hat{\varphi },\hat{\lambda }}},\mathrm{where}{\mathbf{\Psi }}_{{\mathit{n}}_{\hat{\varphi },\hat{\lambda }}}\in \mathbf{\Psi }$

Finally, let the COG SD associated with ${\hat{\mathit{x}}}_{\mathit{t}}^{-}$, be defined as,

${\mathbf{\Sigma }}_{{\mathit{n}}_{\hat{\varphi },\hat{\lambda }}},\mathrm{where}{\mathbf{\Sigma }}_{{\mathit{n}}_{\hat{\varphi },\hat{\lambda }}}\in \mathbf{\Sigma }$

4.3.1. Updating the COG using a-priori information

We first need to calculate the confidence we have in the a-priori COG value ${\mathbf{\Psi }}_{{\mathit{n}}_{\hat{\varphi },\hat{\lambda }}}$. This confidence measure $\mathit{\rho }$ allows ranges between $\left[0,1\right]$, and it is used to scale the contribution of the a-priori COG information. It is determined by the number of observations in the neighborhood, as shown by Equation (20)(20) $\mathit{\rho }=1_{\mathrm{\Psi }}\cdot \frac{{\mathit{K}}_{{\mathit{n}}_{\hat{\varphi },\hat{\lambda }}}}{\max ({\mathit{K}}_{\mathit{H}})}$(20) . If the a-priori grid count for our current cell is high compared to the other cell counts in the neighborhood, we can more confidently say that the predicted location is in an area where historically many vessels have traveled before. If the observed cell count for the current position is lower than in the surrounding cells, we will give less confidence in the a-priori COG update and rather assign more confidence in the LRM with the scaling factor. We can achieve this by calculating the confidence factor $\mathit{\rho }$:

(20) $\mathit{\rho }=1_{\mathrm{\Psi }}\cdot \frac{{\mathit{K}}_{{\mathit{n}}_{\hat{\varphi },\hat{\lambda }}}}{\max ({\mathit{K}}_{\mathit{H}})}$(20)

where

• $\mathit{\rho }$ denotes the confidence (scaling factor) that we have in our prediction as determined by the current position of the vessel.

• ${\mathit{K}}_{{\mathit{n}}_{\hat{\varphi },\hat{\mathit{\lambda }}}}$ denotes the cell count value associated with $\left(\hat{\varphi },\hat{\lambda }\right)$. The index ${\mathit{n}}_{\hat{\varphi },\hat{\lambda }}$ is used to extract it from $\mathit{K}$.

• $1_{\mathrm{\Psi }}$ denotes the indicator function that sets $\mathit{\rho }$ to zero. The indicator function enforces the restrictions that we impose on whether the COG should be updated or not.

$\mathit{\rho }$ determines the total weight the a-priori COG information should have in the LRMAC. If $\mathit{\rho }$ is close to one, we can be confident in the a-priori value ${\mathbf{\Psi }}_{{\mathit{n}}_{\hat{\varphi },\hat{\lambda }}}$, and if it is close to zero, we are less confident, and it should have less of a contribution in our predictions, as seen in Equation (22)(22) ${\hat{\psi }}_{\mathit{t}+1}=\left(1-\mathit{\rho }\right){\psi }_t+\mathit{\rho }{\mathbf{\Psi }}_{{\mathit{n}}_{\hat{\varphi },\hat{\lambda }}}$(22) . Two factors influence $\mathit{\rho }$. First, if a few vessels have traversed the cell associated with ${\mathit{x}}_t$ relative to its neighboring cells, $\frac{{\mathit{K}}_{{\mathit{n}}_{\hat{\varphi },\hat{\lambda }}}}{\max ({\mathit{K}}_{\mathit{H}})}$ would be relatively close to zero, indicating that the a-priori COG information should have less impact in updating ${\mathbf{\Psi }}_{\mathit{t}+1}$. The same is true for the opposite: values closer to one would have a more significant contribution to updating ${\mathbf{\Psi }}_{\mathit{t}+1}$. The second factor is determined by the value of the indicator function in Equation (21)(21) $1_{\Psi }=\{\begin{array}{c}0,{\Psi }_{n_{\widehat{\varphi },\widehat{\lambda }}}\notin \mathit{\Re }\\ 0,{\Sigma }_{n_{\widehat{\varphi },\widehat{\lambda }}}\notin \mathit{\Re }\\ 0,{\Sigma }_{n_{\widehat{\varphi },\widehat{\lambda }}}>{10}^{°}\\ 0,\mathit{max}\left(K_H\right)=01,\mathit{otherwise}\end{array}$(21) .

(21) $1_{\Psi }=\{\begin{array}{c}0,{\Psi }_{n_{\widehat{\varphi },\widehat{\lambda }}}\notin \mathit{\Re }\\ 0,{\Sigma }_{n_{\widehat{\varphi },\widehat{\lambda }}}\notin \mathit{\Re }\\ 0,{\Sigma }_{n_{\widehat{\varphi },\widehat{\lambda }}}>{10}^{°}\\ 0,\mathit{max}\left(K_H\right)=01,\mathit{otherwise}\end{array}$(21)

This equation will evaluate to zero if,

• $\mathbf{\Psi }$ or $\mathbf{\Sigma }$ contains no information at index ${\mathit{n}}_{\hat{\varphi },\hat{\lambda }}$. This implies that there is no a-priori information available for us to make use of.

• The square root of the COG SD at ${\mathit{n}}_{\hat{\varphi },\hat{\lambda }}$ is larger than 10° i.e. standard deviation. This implies that many vessels have traversed through the cell associated with ${\hat{\mathit{x}}}_{\mathit{t}}$, all going in different directions. This implies that the a-priori COG value is unreliable.

The COG value can now be updated as follows:

(22) ${\hat{\psi }}_{\mathit{t}+1}=\left(1-\mathit{\rho }\right){\psi }_t+\mathit{\rho }{\mathbf{\Psi }}_{{\mathit{n}}_{\hat{\varphi },\hat{\lambda }}}$(22)

where

• $\left(1-\mathit{\rho }\right)$ indicates the role that the previous observed COG should have in the COG update.

• ${\mathit{\psi }}_t$ denotes the previously observed or predicted COG.

• $\mathit{\rho }$ denotes the scaling factor that the historic a-priori COG information should have.

• ${\mathbf{\Psi }}_{{\mathit{n}}_{\hat{\varphi },\hat{\lambda }}}$, denotes the a-priori COG value at index ${\mathit{n}}_{\hat{\varphi },\hat{\mathit{\lambda }}}$.

If our confidence in the a-priori COG value is high, we weigh it accordingly. If it is low, we rather put more trust in the COG value at time step $\mathit{t}$ (weighted by $\mathit{\rho }$). Note the estimated COG will replace the COG value that is used in Equation (10)(10) ${\mathbf{\Lambda }}_{\mathrm{t}}={\left[cos\left({\psi }_t\right)sin\left({\psi }_t\right)\right]}^{\mathit{T}}$(10) .

4.3.2. A flow diagram representation of the LRMAC

In Figure 6, a flowchart of the LRMAC methodology is depicted. The flowchart depicts the LRM and how a-priori course information is added to extend the LRM into the LRMAC. Initializations are in green, functions are in blue, the predicted location is in gray, and parameter extracts are denoted in orange. It is assumed that all pre-processing has already been applied to the dataset.

Figure 6. LRMAC flowchart.

The first step is to construct the SMs from the dataset containing all historic AIS data and initialize all parameters. The last $\mathit{\omega }$ recorded observations are used as additional input parameters in the LRM, specifically for the measurement update equations. The LRM consists of two sets of equations, measurement update and predictor equations. The measurement equations are used to update the predictor parameters and the predictor equations to predict the next set of coordinates.

The algorithm starts at the indicated red dot in the flowchart. The LRM is used to estimate the longitude and latitude at the next time step. All predictions are made at regular spaced time intervals $\mathrm{\Delta }{k}_t$. The LRM assumes a constant SOG and COG, based on the last $\mathit{\omega }$ observations.

Given the predicted longitude and latitude, we extend the LRM into the LRMAC by dropping the constant COG assumption. The COG will now be updated based on a-priori COG information at the location predicted by the LRM. The extension of the LRM allows the predictions to follow historic movement trends in the SMs, allowing for non-linear trajectory prediction.

Given the predicted location ${\hat{\mathit{x}}}_{\mathit{t}}$, the corresponding index ${\mathit{n}}_{\hat{\varphi },\hat{\lambda }}$ in the SMs can be calculated, and the neighboring indexes matrix $\mathit{H}$. Using ${\mathit{n}}_{\hat{\varphi },\hat{\lambda }}$ and $\mathit{H}$, the a-priori values: ${\mathit{K}}_{\mathit{H}}$, ${\mathbf{\Psi }}_{{\mathit{n}}_{\hat{\varphi },\hat{\lambda }}}$, ${\mathbf{\Sigma }}_{{\mathit{n}}_{\hat{\varphi },\hat{\lambda }}}$, and ${\mathit{K}}_{{\mathit{n}}_{\hat{\varphi },\hat{\lambda }}}$ is extracted and used in Equations (20)(20) $\mathit{\rho }=1_{\mathrm{\Psi }}\cdot \frac{{\mathit{K}}_{{\mathit{n}}_{\hat{\varphi },\hat{\lambda }}}}{\max ({\mathit{K}}_{\mathit{H}})}$(20) and (21(21) $1_{\Psi }=\{\begin{array}{c}0,{\Psi }_{n_{\widehat{\varphi },\widehat{\lambda }}}\notin \mathit{\Re }\\ 0,{\Sigma }_{n_{\widehat{\varphi },\widehat{\lambda }}}\notin \mathit{\Re }\\ 0,{\Sigma }_{n_{\widehat{\varphi },\widehat{\lambda }}}>{10}^{°}\\ 0,\mathit{max}\left(K_H\right)=01,\mathit{otherwise}\end{array}$(21) ) to calculate the updated COG value ${\hat{\psi }}_{\mathit{t}+1}$. The updated COG is used in the next iteration by updating the measurement equations of the LRM, which in turn updates the predictor equations. This allows for a COG that updates dynamically based on historic AIS data.

All other parameters, including the COG, will be updated once new observations are received from the vessel, updating all the measurement equations.

4.4. Single-Point Neighbor Search (SPNS)

We will now be summarizing the method presented by Hexeberg, Flåten, and Brekke (2017) called Single-Point Neighbor Search (SPNS). We will be comparing the LRMAC to the SPNS, as it was one of the most similar methods we found. We will compare the prediction accuracy and the time complexity of each. Note, do not confuse the variables defined in this section with the variables used in the rest of the paper.

Let

(23) $\mathit{X}={\left[{\mathit{X}}_1\mathit{\hspace{0.17em}}{\mathit{X}}_2\dots \mathit{\hspace{0.17em}}{\mathit{X}}_M\right]}^{\mathit{T}}$(23)

be defined as a matrix with all the historic AIS data observations, where $\mathit{M}$ indicates the number of AIS messages recorded.

Let

(24) ${\mathit{X}}_i=\left[\mathit{MMS}{I}_i{t}_i{\mathit{p}}_{\mathit{i}}^T\mathit{\hspace{0.17em}}{\mathit{\chi }}_{\mathit{i}}{v}_i\right]$(24)

be defined as a vector, where $\mathit{i}\in \left\{1,2,\dots ,\mathit{M}\right\}$. The elements in ${\mathit{X}}_i$ refer to the MMSI, timestamp, position vector, COG, and SOG of timestep i. The position vector ${\mathit{p}}_i^T$ denotes the spatial location of a vessel, where ${\mathit{p}}_i={\left[{\mathit{\lambda }}_i,{\varphi }_i\right]}^{\mathit{T}}$ denotes the longitude and latitude at message $\mathit{i}$.

A predicted trajectory consists out of $K^s$ predicted positions, which are instants of time. At every iteration $\mathit{k}$, a prediction is made, where $\mathit{k}\in 1,\dots ,K^s$. The predicted state is divided into an a-priori state ${\hat{\mathit{X}}}_{\mathit{i}}^{\mathit{k}-}$ and a-posteriori state ${\hat{\mathit{X}}}_{\mathit{i}}^{\mathit{k}+}$. The states are denoted as:

(25) ${\hat{\mathit{X}}}_{\mathit{i}}^{\mathit{k}-}=\left[\mathit{MMS}{I}_i{\hat{t}}^{\mathit{k}}{\hat{\mathit{p}}}^{\mathit{k}}{\hat{\chi }}^{\mathit{k}-}{\hat{v}}^{\mathit{k}-}\right]$(25)

and

(26) ${\hat{\mathit{X}}}_{\mathit{i}}^{\mathit{k}+}=\left[\mathit{MMS}{I}_i{\hat{t}}^{\mathit{k}}{\hat{\mathit{p}}}^{\mathit{k}}{\hat{\chi }}^{\mathit{k}+}{\hat{v}}^{\mathit{k}+}\right]$(26)

The only difference between the two equations above, is the COG and SOG. The predicted ${\hat{\chi }}^{\mathit{k}-}$ and ${\hat{v}}^{\mathit{k}-}$ at iteration $\mathit{k}$ in the a-priori state, represent the predicted COG and SOG between the previous position ${\hat{\mathit{p}}}^{\mathit{k}-1}$ and the current position ${\hat{\mathit{p}}}^{\mathit{k}}$. The a-posteriori predicted COG and SOG at iteration $\mathit{k}$ is the difference between the current position ${\hat{\mathit{p}}}^{\mathit{k}}$ and the next position ${\hat{\mathit{p}}}^{\mathit{k}+1}$.

The SPNS makes use of a close neighbor search, where a set radius is defined as a parameter and observations within the radius is queried, given that the observations adhere to a set of predefined constraints.

Let the close neighbors (CNs) at prediction step $\mathit{k}$ be defined as:

(27) $C^k=\{X_i{\mathit{}}_{|}\mathit{d}\left({\widehat{p}}^{\mathit{k}},p_i\right)\le r_c,{\chi }_i\in \mathit{S},X_i\in \mathit{X}\}$(27)

where

• $\mathit{d}\left({\hat{\mathit{p}}}^{\mathit{k}},{\mathit{p}}_i\right)$ is defined as the Haversine distance between the longitude and latitude in ${\hat{\mathit{p}}}^{\mathit{k}}$ and ${\mathit{p}}_{\mathit{i}}$. The Haversine distance is defined in Equation (2)(2) $\begin{array}{rl}& \mathit{d}=2\mathit{r}\cdot \mathit{arcsin}\\ & \sqrt{\mathit{si}{n}^2\left(\frac{{\varphi }_2-{\varphi }_1}{2}\right)+\mathit{cos}\left({\varphi }_1\right)\mathit{cos}\left({\varphi }_2\right){\mathit{sin}}^2\left(\frac{{\lambda }_2-{\lambda }_1}{2}\right)}\end{array}$(2) .

• $r_c$, is defined as the search radius in meters, to search for all the CNs within $r_c$ of the current position. The value of $r_c$ is a predefined parameter.

• $\mathit{S}$ is defined as the interval of course angles. Only observations within the course angle interval will be included in the CN set.

Let the course angles $\mathit{S}$ be defined by:

(28) $\mathit{S}=\left[{\hat{\chi }}^{\mathit{k}-}-\mathrm{\Delta \chi },\mathit{\hspace{0.17em}}{\hat{\chi }}^{\mathit{k}-}+\mathrm{\Delta \chi }\right]$(28)

where $\mathrm{\Delta \chi }>0$, is the maximum course angle deviation, which is a predefined parameter.

The above pre-processing steps of the distance and COG deviation will filter out all the observations that we do not need for the remaining steps. All the CNs are extracted from historic AIS data. Let every state that belongs to the set of CNs at the prediction step $\mathit{k}$ be denoted as ${\mathit{X}}_i\in {\mathit{C}}^k$, where ${\mathit{X}}_c^k=\left[\mathit{MMS}{I}_c^k{t}_c^k{\chi }_c^k{v}_c^k\right]$ and $\mathit{c}\in \left\{1,\dots ,\mathit{C}\right\}$. $\mathit{C}$ is the number of CNs at $\mathit{k}$.

The predicted trajectory at a state ${\mathit{X}}_i$ is defined as,

(29) ${\hat{\mathit{T}}}_{\mathit{i}}=\left\{\left[{\hat{\mathit{p}}}^1{\hat{t}}^1\right],\dots ,\left[{\hat{\mathit{p}}}^{\mathit{K}}{\hat{t}}^{{\mathit{K}}^s}\right]\right\}$(29)

Let the true trajectory, ${\mathit{T}}_i$, given state ${\mathit{X}}_i$, be defined as

(30) ${\mathit{T}}_i=\left\{\left[{\mathit{p}}^1{t}^1\right],\dots ,\left[{\mathit{p}}^L{t}^L\right]\right\}$(30)

$\mathit{L}$ denotes the number of AIS states recorded. $K^s$ and $\mathit{L}$ are not necessarily equal as several prediction steps can be made between two subsequent AIS messages. The first elements in both ${\hat{\mathit{T}}}_{\mathit{i}}$ and ${\mathit{T}}_i$ are equal, as it is the starting point given by state ${\mathit{X}}_i$.

4.4.1 SPNS prediction

A new parameter is introduced, $\mathrm{\Delta }\mathit{l}$, which denotes the step length from the current observation to the next predicted observation in meters. $\mathrm{\Delta }\mathit{l}$ decides how far the next position should be propagated. Let the predicted position be denoted by

(31) ${\hat{\mathit{p}}}^{\mathit{k}+1}={\mathit{p}}^k+\mathrm{\Delta }\mathit{l}\cdot {\left[sin({\hat{\chi }}^{\mathit{k}+})\cdot \mathit{f}\left({\hat{\varphi }}^{\mathit{k}}\right)cos({\hat{\chi }}^{\mathit{k}+})\cdot \mathit{g}\left({\hat{\varphi }}^{\mathit{k}}\right)\right]}^{\mathit{T}}$(31)

$\mathit{f}\left({\hat{\varphi }}^{\mathit{k}}\right)$ and $\mathit{g}\left({\hat{\varphi }}^{\mathit{k}}\right)$ are functions of the current latitude ${\hat{\varphi }}^{\mathit{k}}$, which transforms from meters to degrees longitude and latitude, respectively. The step length, $\mathrm{\Delta }\mathit{l}$, reflects the curvature of the sea lanes ahead.

4.4.1.1 Course prediction

The COG value ${\hat{\chi }}^{\mathit{k}+}$ is used when calculating the predicted position ${\hat{\mathit{p}}}^{\mathit{k}+1}$. ${\hat{\chi }}^{\mathit{k}+}$ is the a-priori course calculated from the extracted CNs at position ${\mathit{p}}^k$. Note that the course is periodic in $\left[0^{\circ },{360}^{\circ }\right]$; special care must be taken when calculating the CN set’s mean COG. The mean COG is calculated as follows:

(32) ${\bar{\mathit{\chi }}}_{\mathit{c}}=\left\{\begin{array}{c}\mathrm{ta}{\mathrm{n}}^{-1}\left(\frac{\bar{s}}{\bar{c}}\right)\mathit{if}\bar{s}>0,\bar{c}>0\\ \mathrm{ta}{\mathrm{n}}^{-1}\left(\frac{\bar{s}}{\bar{c}}\right)+{180}^{\circ }\mathit{if}\bar{c}<0\\ {tan}^{-1}\left(\frac{\bar{s}}{\bar{c}}\right)+{360}^{\circ }\mathit{if}\bar{s}⟨0,\bar{c}⟩0\end{array}\right.$(32)

where

(33) $\bar{s}=\frac{1}{C}\sum _{\mathit{c}=1}^{C}sin\left({\mathit{\chi }}_c\right)$(33)

(34) $\bar{c}=\frac{1}{C}\sum _{\mathit{c}=1}^{C}cos\left({\mathit{\chi }}_c\right)$(34)

A constant velocity model is used whenever $C^k\notin \mathbb{R}$ as done by Hexeberg, Flåten, and Brekke (2017). The median course ${\tilde{\chi }}_{\mathit{c}}$ can be calculated by calculating $\tilde{s}$ and $\tilde{c}.$ Hexeberg, Flåten, and Brekke (2017) recommended to use the median course for non-linear trajectories.

4.4.1.2 Speed prediction

The median ${\tilde{v}}_{\mathit{c}}$ of the CNs is used to calculate the predicted speed. Note that the speed prediction is only used in the time update equation shown in Algorithm 1. The predicted speed ${\tilde{v}}_{\mathit{c}}$ is used to calculate the time passed between the current observation and the predicted observation, where the time passed is denoted by $\frac{\mathrm{\Delta }l}{{\hat{v}}^{\mathit{k}+}}$. The time update equation is defined as:

(35) ${\hat{t}}^{\mathit{k}+1}={\hat{t}}^{\mathit{k}}+\frac{\mathrm{\Delta }l}{{\hat{v}}^{\mathit{k}+}}$(35)

where

• ${\hat{t}}^{\mathit{k}}$ denotes the current time,

• $\mathrm{\Delta }\mathit{l}$ denotes the distance between the current and predicted observation, and

• ${\hat{v}}^{\mathit{k}+}={\tilde{v}}_{\mathit{c}}$, given set $C^k$ at $\mathit{k}$.

The calculation of the predicted time allows the SPNS to have regular spaced distance intervals of length $\mathrm{\Delta }\mathit{l}$. The algorithm for the SPNS is presented below.

Table

Algorithm 1 Single Point Neighbor Search Prediction
1:	${\mathit{X}}_i$ given	The state we predict from
2:	Set decision parameters
	(a) $\mathrm{\Delta }\mathit{l}$ (b) $r_c$ (c) $\mathrm{\Delta }\text{Chi}$ (d) $K^s$	Step length [m] Search radius [m] Maximum course angle deviation [deg] Number of prediction steps
3:	Set ${\hat{\mathit{X}}}_{\mathit{i}}^{\mathit{k}-}={\mathit{X}}_i$
4:	for $\mathit{k}=1$ to $\mathit{K}$ do
5:	Find all CNs $X_c^k$ around ${\hat{\mathit{X}}}_{\mathit{i}}^{\mathit{k}-}$
6:	Calculate ${\hat{\mathit{X}}}_{\mathit{i}}^{\mathit{k}+}\mathit{by}:$ Calculating ${\hat{\chi }}^{\mathit{k}+}$ based on $X_c^k$ Calculating ${\hat{v}}^{\mathit{k}+}$ based on $X_c^k$
7:	Calculate the next predicted position at its predicted point in time: Calculate ${\hat{\mathit{p}}}^{\mathit{k}+1}$ according to Equation (31)(31) ${\hat{\mathit{p}}}^{\mathit{k}+1}={\mathit{p}}^k+\mathrm{\Delta }\mathit{l}\cdot {\left[sin({\hat{\chi }}^{\mathit{k}+})\cdot \mathit{f}\left({\hat{\varphi }}^{\mathit{k}}\right)cos({\hat{\chi }}^{\mathit{k}+})\cdot \mathit{g}\left({\hat{\varphi }}^{\mathit{k}}\right)\right]}^{\mathit{T}}$(31) Calculate ${\hat{t}}^{\mathit{k}+1}={\hat{t}}^{\mathit{k}}+\frac{\mathrm{\Delta }l}{{\hat{v}}^{\mathit{k}+}}$
8:	Set ${\widehat{X}}_{\mathit{i}}^{\left(\mathit{k}+1\right)-}=\mathit{MMS}{I}_i{\widehat{t}}^{\mathit{k}+1}{\widehat{p}}^{\mathit{k}+1}{\widehat{\chi }}^{\mathit{k}+}{\widehat{v}}^{\mathit{k}+}$
9:	end for

In Table 4, the set of hyperparameters that are used for curved trajectories is shown. These decision parameters are identical to the ones presented by Hexeberg, Flåten, and Brekke (2017) and were used to compare the SPNS to the LRMAC.

Table 4. Curved trajectory decision parameters.

Decision Parameter	Value	Explanation
$r_c$	50 m	Search radius of for the CNs
$\mathrm{\Delta }\mathit{l}$	$2r_c$	Prediction step length [m]
$\mathrm{\Delta \chi }$	25°	Maximum course deviation
${\hat{\chi }}_{\mathit{i}}^{\mathit{k}+}$	${\widetilde{\text{Chi}}}_{\mathit{c}}$	Course prediction used at every iteration $\mathit{k}$
${\hat{v}}_{\mathit{i}}^{\mathit{k}+}$	${\tilde{v}}_{\mathit{c}}$	Speed prediction used at every iteration $\mathit{k}$

5. Experimental design

This section will discuss how we set up our experiments to get the optimal set of parameters for the proposed LRMAC model. We will also be discussing how we will compare the LRMAC to the SPNS. We will be comparing the LRM and LRMAC to show the increase in prediction accuracy of the LRMAC over the LRM and then compare the LRMAC to the SPNS model comparing prediction accuracy and time complexity. We will be discussing the vessels used, how we extracted their trajectories, the subsampling approach that we made use of, and how we determined the performance of the two approaches.

To compare the LRMAC and SPNS, 40 unique vessel trajectories were extracted at random from the dataset published by Ray et al. (2019) after having been pre-processed. The pre-processing removed stationary observations and vessels with too few observations, as mentioned earlier in the paper. The extracted vessels are a mix of cargo and tanker vessels. The extracted trajectories of the vessels were all from the high-density areas depicted in Figure 3. We restricted our use case to this example to show the performance of the LRM, LRMAC, and SPNS on curved trajectories. The LRMAC algorithm can be applied to various use cases and will perform well on any SM containing clear regions exhibiting path directionality.

5.1 SPNS vessel query setup

Since the SPNS algorithm requires querying historic observations within a specified radius, all the data was loaded into a PostgreSQL database.² An extension was added to PostgreSQL called PostGIS,³ which allows for improved spatial queries with a datatype called geometry. The PostGIS plugin uses a unique kind of indexing. Querying observations in PostgreSQL with PostGIS allow that vessels within a given radius from the search point can be extracted into $C^k$.

5.2 Trajectory subsampling method

In order to have a large test set of trajectories to identify the ideal parameters for the LRMAC and compare it to the SPNS, we created a method that would subsample trajectories. We are effectively creating multiple trajectories from one observed trajectory. A given trajectory ${\mathit{T}}_i$ will be subsampled into different time subsets given a prediction length $\mathit{h}$. Each time a subset’s starting observation will differ by 1 hour compared to the subset that precedes or follows it. We refer to this hour difference as a stride length $\mathit{s}$, where $\mathit{s}$ defines the starting point of the stride hour from ${\mathit{T}}_i$. Stride values are measured in hours. This method allows us to extract multiple sub trajectories ${\mathit{T}}_{\mathit{s},\mathit{h}}$ from ${\mathit{T}}_i$. Let the number of sub trajectories that can be created with a prediction length of $\mathit{h}$ from ${\mathit{T}}_i$ be denoted by

(36) $\mathrm{\#}{T}_h={⌊\mathit{max}({\mathit{t}}_{{\mathit{T}}_i})⌋}_{\mathit{hour}}-\mathit{h}+1$(36)

and let the total number of sub trajectories from ${\mathit{T}}_i$ with different starting positions $\mathit{s}$ be denoted by,

(37) $\mathrm{\#}{T}_{\mathit{s},\mathit{h}}={\sum }_{\mathit{h}=1}^{{⌊\mathit{max}({\mathit{t}}_{{\mathit{T}}_i})⌋}_{\mathit{hour}}}\mathrm{\#}{T}_h$(37)

where

• ${⌊\mathit{max}({\mathit{t}}_{{\mathit{T}}_i})⌋}_{\mathit{hour}}$ denotes the closest floored hour to the maximum observed time in ${\mathit{T}}_i$.

• $\mathit{h}$ refers to the prediction length, where $\mathit{h}\in \{1,2,\dots ,{\lfloor max(t_{T_i})\rfloor }_{\mathit{hour}}$, and

• $\mathit{s}$ denotes the stride starting position measured in hours,

$\mathit{s}\in \left\{0,1,\dots ,{⌊\mathit{max}({\mathit{t}}_{{\mathit{T}}_i})⌋}_{\mathit{hour}}-\mathit{h}\right\}$

In Figure 7, an example is shown of the trajectory subsampling and the number of subsets we can create given an observed trajectory between 6 and 7 hours. Using Equation (37)(37) $\mathrm{\#}{T}_{\mathit{s},\mathit{h}}={\sum }_{\mathit{h}=1}^{{⌊\mathit{max}({\mathit{t}}_{{\mathit{T}}_i})⌋}_{\mathit{hour}}}\mathrm{\#}{T}_h$(37) , the total number of sub trajectories that can be created from this trajectory is $21$.

Figure 7. Subsample visualization.

This method of sub-trajectory sampling allows us to have multiple trajectories to compare the prediction and time performance between the LRMAC and the SPNS model. The prediction performance is determined by the Haversine distance between the expected spatial location at time $\mathit{t}$ compared to the predicted location at time $\hat{t}$. The time performance is computed by the time the methods take from execution until the predicted trajectory of length $\mathit{h}$ has been calculated.

First, we will run the LRMAC method on all sub-trajectories with different sets of parameters $\mathit{\eta }$ and $\mathit{\omega }$ to determine the optimal set of parameters that results in the smallest median Haversine distance on average. The LRMAC with its optimal parameters and SM densities will be compared to the LRM, and SPNS with the parameters defined by Hexeberg, Flåten, and Brekke (2017). The SM densities were not optimized for these tests, the cell sizes remained constant, and the number of cells fixed at $1250\times 1250$. We will also be testing the speed performance between the two algorithms, measuring the execution times.

5.3 Parameter estimation of LRMAC

To obtain the optimal set of hyperparameters for the LRMAC, the LRMAC was tested extensively on different combinations of $\mathit{\eta }$ and $\mathit{\omega }$. The combination of $\mathit{\eta }$ and $\mathit{\omega }$ that gave the lowest overall Haversine distance error will be chosen as the ideal set of hyperparameters. The considered parameter values were: $\mathit{\eta }=\left[1,2,3,4\right]$ and $\mathit{\omega }=\left[1,3,5,7,9\right]$, resulting in 20 combinations.

The parameters that yielded the best results were: $\mathit{\omega }=3$ and $\mathit{\eta }=1$. The LRMAC was tested on all sub trajectories. The combination of $\mathit{\omega }$ and $\mathit{\eta }$ with the smallest average median Haversine distance over all sub-trajectories and timeframes was chosen as the ideal set of hyperparameters. Figure 8 compares the average median Haversine distances, and we can see that there is no significant performance increase after the combination $\mathit{\omega }=3$ and $\mathit{\eta }=1$. The chosen hyperparameters have a lower complexity as the $\mathit{\omega }$ and $\mathit{\eta }$ are small, compared to the other parameters with similar performance. Smaller values of $\mathit{\omega }$ means that our function is more flexible to incorporate new information, as there are fewer elements to consider. Large values of $\mathit{\omega }$ will mean that our algorithm is less flexible, using more observed observations from the current trajectory to make a prediction. Larger values of $\mathit{\eta }$ means more historic SM cells to include when making a prediction, as the neighborhood size is larger.

Figure 8. Average median haversine error for each hyperparameter combination.

5.4 SPNS prediction adjustment

Since the SPNS predicts in constant distance intervals $\mathrm{\Delta }\mathit{l}$ instead of constant time intervals as the LRMAC does, we had to adjust the final prediction of the SPNS to allow for an exact comparison at $\mathit{h}$ (prediction horizon). To get the predicted spatial location of the SPNS after $\mathit{h}$ has passed, we let the SPNS predict until ${\hat{t}}^{\mathit{k}+1}>\mathit{h}$. We calculate the time that has passed between the last two predicted observations ${\hat{\mathit{X}}}_{\mathit{i}}^{\left(\mathit{k}\right)-}$ and ${\hat{\mathit{X}}}_{\mathit{i}}^{\left(\mathit{k}+1\right)-}$, where time is indicated by ${\hat{t}}^{\mathit{k}}$ and ${\hat{t}}^{\mathit{k}+1}$, and where, ${\hat{t}}^{\mathit{k}}<\mathit{h}<{\hat{t}}^{\mathit{k}+1}.$, let he total time passed be denoted by $\mathrm{\Delta }{t}_{|\hat{y}|}$. The Haversine distance between ${\hat{\mathit{X}}}_{\mathit{i}}^{\left(\mathit{k}\right)-}$ and ${\hat{\mathit{X}}}_{\mathit{i}}^{\left(\mathit{k}+1\right)-}$ is calculated $\mathrm{\Delta }{l}_{\left|{\hat{\mathit{X}}}_{\mathit{i}}^{\left(\mathit{k}\right)-}-{\hat{\mathit{X}}}_{\mathit{i}}^{\left(\mathit{k}+1\right)-}\right|}$. The SPNS will be rerun from ${\hat{\mathit{X}}}_{\mathit{i}}^{\left(\mathit{k}\right)-}$ with the same set of parameters and only changing the prediction step length $\mathrm{\Delta }\mathit{l}=\mathrm{\Delta }{l}_{\left|{\hat{\mathit{X}}}_{\mathit{i}}^{\left(\mathit{k}\right)-}-{\hat{\mathit{X}}}_{\mathit{i}}^{\left(\mathit{k}+1\right)-}\right|}$. This will give us the predicted spatial location at $\mathit{h}$, as ${\hat{t}}^{\mathit{k}+1}=\mathit{h}$ determined by Equation (35)(35) ${\hat{t}}^{\mathit{k}+1}={\hat{t}}^{\mathit{k}}+\frac{\mathrm{\Delta }l}{{\hat{v}}^{\mathit{k}+}}$(35) .

5.5 Prediction performance comparison algorithm

The algorithm below shows how the errors were calculated for both the LRMAC and SPNS on each sub trajectory. We used the median error results for each predicted timeframe $\mathit{h}$ to compare the performance between the two algorithms.

Table

Algorithm 2 Prediction Performance Comparison
1:	Set different parameters: (a) $\mathit{\omega }=3$ (b) $\mathit{\eta }=1$ (c) ${\in }_{\mathrm{LRMAC}}=\left[{\in }_{\mathrm{LRMA}{\mathrm{C}}_1},\dots ,{\in }_{\mathit{LRMA}{\mathit{C}}_{\mathit{max}\left(\mathit{h}\right)}}\right]$ (d) ${\in }_{\mathrm{SPNS}}=\mathit{\hspace{0.17em}}\left[{\in }_{SPN{S}_1},\dots ,{\in }_{\mathit{SPN}{\mathit{S}}_{\mathit{max}\left(\mathit{h}\right)}}\right]\mathit{\hspace{0.17em}}$	LRMAC Jagged array for each $\mathit{h}$ SPNS Jagged array for each $\mathit{h}$
2:	for ${\mathit{T}}_i\mathit{in}\mathit{T}$ do
3:	Get ${\mathit{T}}_{\mathit{s},\mathit{h}}$ which is a set of sub trajectories from ${\mathit{T}}_i$
4:	for $\mathit{h}$ in $1,2,\dots ,{⌊\mathit{max}({\mathit{t}}_{{\mathit{T}}_i})⌋}_{\mathit{hour}}$ do
5:	for $\mathit{s}$ in {$0,1,\dots ,{⌊\mathit{max}({\mathit{t}}_{{\mathit{T}}_i})⌋}_{\mathit{hour}}-\mathit{h}$} do
6:	$T_{\mathit{test}}=T_{\mathit{s},\mathit{h}}\left[0:3\right]$	Initial observations for the algorithms
7:	${\hat{\lambda }}_{\mathit{LRMAC}},$ ${\hat{\varphi }}_{\mathit{LRMAC}}$ = $\mathit{LRMAC}\left(T_{\mathit{test}},\mathit{\eta },\mathit{\omega },\mathit{h}\right)$	LRMAC, given the parameters and h
8:	${\hat{\lambda }}_{\mathit{SPNS}},$ ${\hat{\varphi }}_{\mathit{SPNS}}$ = $\mathit{SPNS}\left(T_{\mathit{test}},\mathit{h}\right)$	SPNS, given the parameters and h
9:	Calculate the haversine distance error: ${\in }_{\mathrm{LRMAC}}\mathit{\hspace{0.17em}}$ append havdist(${\mathit{T}}_{\mathit{s},\mathit{h}}$,${\hat{\lambda }}_{\mathit{LRMAC}},$ ${\hat{\varphi }}_{\mathit{LRMAC}}$)	Error to corresponding $\mathit{h}$ array in $\mathit{\hspace{0.17em}}{\in }_{\mathrm{LRMAC}}\mathit{\hspace{0.17em}}$
10:	${\in }_{\mathrm{SPNS}}$ append havdist(${\mathit{T}}_{\mathit{s},\mathit{h}}$,${\hat{\lambda }}_{\mathit{SPNS}},$ ${\hat{\varphi }}_{\mathit{SPNS}}$)	Error to corresponding $\mathit{h}$ array in $\mathit{\hspace{0.17em}}{\in }_{\mathrm{SPNS}}$ havdist() denotes the haversine distance between the observed and predicted location.

5.6 System specifications

All the tests will be carried out on a system with the specifications listed in Table 5.

Table 5. Testing system specifications.

Description	Name	Specification
CPU	Intel i7-10700K	8 core 16 threads @ 5.1 GHz
Memory	Corsair Vengeance LPX	DDR4 @ 3600 MHz
Storage	HP EX920	NVMe M.2 SSD @ 512 GB
Software	Python	3.8.5
Database	PostgreSQL	12.6
Database plugin	PostGIS	3.1.1
Operating System	Ubuntu Desktop	20.04 LTS

6. Results

In this section, we will compare the LRMAC with the SPNS.

The LRMAC is an improved version of the LRM, where the LRMAC allows for the prediction of non-linear trajectories. In Figure 9 below, the median Haversine error is shown for both the LRM (red) and LRMAC (blue), together with the standard deviation (SD) from the median. The LRMAC has a reduced SD and median error compared to the LRM. In terms of short-term prediction accuracy, we see that the LRM and LRMAC are not significantly different, but for longer prediction windows (>120 min), the proposed method shows an improvement. To further support this observation, when looking at Figure 10, we can see that the LRMAC reproduced the trajectory, where the LRM went off course. The LRMAC assumes a constant speed calculated which was estimated from the first $\mathit{\omega }$ observed observations. The LRMAC predicted where the vessel would be after 8 hours, should it have traveled at a constant speed. The prediction error comparison does not accurately represent the performance of the LRMAC. When looking at the predicted trajectory, we can see that the predicted LRMAC trajectory is more similar to the actual trajectory when compared to the LRM.

Figure 9. LRM vs LRMAC.

Figure 10. LRM vs LRMAC on SDM (MMSI 419689000).

The error is measured as the Haversine distance between the actual observation at time $\mathit{t}$ compared to the predicted observation at time $\mathit{t}$ for each method. The LRMAC has a lower SD for longer range predictions than the LRM. This is expected as our test set contains vessels that are in the highways depicted in Figure 3. Shorter prediction periods $\mathit{h}$ will mean the subsets will be near-linear as cargo and tanker vessels have a slow rate of turn. We do not expect a vessel trajectory to stay linear for more extended prediction periods, as there are obstacles like landmasses and other vessels. The LRM will keep predicting linearly, whereas the LRMAC will adapt to predict non-linear trajectories using historical COG information. We can see from $\mathit{h}=240$ onwards that the median difference between the LRM and LRMAC increases. We expect that the LRMAC will be able to predict with any SM that has path directionality.

In Figure 10, the LRM (lime) and LRMAC (aqua) predictions are shown of MMSI 419689000, together with the actual trajectory (magenta) that the vessel had. Looking at the error lines of the LRMAC (yellow) and the LRM (red), we see that their distances to the expected location differ by $3.47\mathrm{km}$. The error distances are not significantly different; however, the LRMAC was still able to predict the correct path, given the constant speed assumption. The LRMAC only used $\mathit{\omega }=3$ historic observations to predict the 8-hour trajectory. The background of Figure 10 depicts the vessel counts SM $\mathit{K}$, shown in Figure 3. Given a use case where the historic data have paths that are longer and highly non-linear (not just a piecewise linear historic path), we believe that the LRMAC would do significantly better than the LRM as it will follow the trend of the historic data.

In Figure 11, the median prediction error of the LRMAC (blue) is compared to the SPNS (green). Similar to the LRM and LRMAC comparison, the error was calculated over all the sub trajectories with a prediction length of $\mathit{h}$. The median error for both methods was extracted at each $\mathit{h}$ together with the standard deviation (SD). Note that the SD is plotted around the median, instead of the mean absolute deviation, and this was done to reflect the algorithmic stability.

Figure 11. SPNS vs LRMAC prediction results.

Figure 11, which compares the LRMAC and the SPNS, is a better comparison as both the LRMAC and SPNS follow the historic route. Looking at Figure 11, we see no significant prediction performance difference between the LRMAC and the SPNS for short prediction periods. However, when considering longer prediction periods (>120 min), the median error and SD of the LRMAC are lower than that of the SPNS. The SPNS’s SD starts to increase dramatically, compared to the LRMAC’s SD, which increases at a lower rate. When looking at the prediction interval $\left[5,15\right]$ min, the recommended time prediction horizon for the SPNS, the SPNS has a smaller SD than the LRMAC.

One assumption that the LRMAC makes is that a vessel keeps a constant speed, and the prediction is based on this constant speed assumption. The SOG of a vessel is derived from the last $\mathit{\omega }$ observed observations. The SPNS calculates the SOG based on the historic observations, but the speed is not used in the coordinate prediction step but instead in the time predictions.

In Table 6 below, we see the time complexity of both methods for each prediction length $\mathit{h}$. The time was measured over all the sub-trajectories, and the median time is shown. The SPNS computation time is significantly longer to predict than the LRMAC. The speed of the SPNS is limited by the time it takes to query the CNs from a database. Larger database table sizes will lead to more significant query times as there are more observations to search through. The LRMAC uses SMs whose sizes stay fixed, given more observations for the same longitude and latitude range. The SM sizes will only increase if the longitude and latitude range of the observed data increase.

Table 6. LRMAC vs SPNS median prediction times.

$\mathit{h}$	LRMAC time		SPNS time		Time difference
min	s	min	s	min	min
5	0.11	0.0018	30.92	0.5154	−0.5135
15	0.33	0.0055	95.68	1.2614	−1.5892
30	0.66	0.0110	193.07	3.2178	−3.2068
60	1.31	0.0218	371.87	6.1978	−6.1760
120	2.64	0.0440	712.29	11.8715	−11.8275
180	3.98	0.0663	1026.38	17.1063	−17.0399
300	6.56	0.1093	1724.46	28.7410	−28.6317

Looking at Figure 11 and Table 6, we see that even though the two methods are not significantly different in their prediction accuracies for short prediction periods, the LRMAC had a shorter execution time. In Figure 11, we see that the SPNS has a smaller SD than the LRMAC for the $\left[5,15\right]$ min timeframe.

We believe that in higher density areas closer to harbors, the SPNS will have more accurate prediction results as the set of CNs will only contain vessels moving in the same direction, where the LRMAC would default to the LRM as the a-priori COG SD $\mathbf{\Sigma }$, will be high. In areas closer to harbors, vessels tend to move slower, and prediction timeframes are usually smaller with increased AIS coverage. We think that the SPNS can be used together with the LRMAC. If the LRMAC encounters a cell with a high SD, the SPNS can be deployed until a cell with a lower SD is encountered. A hybrid approach between the SPNS and LRMAC may improve prediction times and improve prediction accuracy in areas with a significant amount of traffic in different directions.

It should be noted that the LRMAC is limited by the spatial resolution of the SMs it employs. In the case of sparse historic data, the SMs generated from this data would be inaccurate, resulting in a decrease in performance of the LRMAC. Also, when the prediction time interval between two consecutive predictions is too large, the LRMAC will skip past important a-priori information and have inaccurate prediction results. Furthermore, the LRMAC was only tested on cargo and tanker vessels in this study, as the movement of other vessel types, such as fishing vessels are of an erratic nature and the SMs that will reflect this behavior, with high COG standard deviation.

The LRMAC can be used to predict trajectories of vessels or to impute historic trajectories. The LRMAC would be more accurate as observations are recorded, and the observations in the window size $\mathit{\omega }$ is updated. Currently, the predictions of the LRMAC and SPNS were done on the assumption that only the first $\mathit{\omega }$ observations’ information will be used, simulating AIS transponders that are switched off for extended periods of time.

7. Conclusions

In this paper, we presented an extension of the LRM method proposed by Burger, Kleynhans, and Lups Grobler (2020), which uses historic AIS vessel information (particularly the historic COG information), which is robust and easy to implement. The LRMAC implemented in this paper had a smaller incurred prediction error and associated standard deviation than the LRM implemented by Burger, Kleynhans, and Lups Grobler (2020). The LRMAC can predict a trajectory similar to the actual trajectory where the LRM goes off course. The LRMAC makes use of SMs to improve its prediction capability. The LRMAC was compared to Single-Point Neighbor Search (SPNS), which has a similar level of computational complexity and, for the use case of predicting tanker and cargo vessel trajectories up to 8 hours into the future, showed improved results both in terms of the prediction accuracy and execution time. The LRMAC can be used to predict trajectories of vessels or impute vessel trajectories. The LRMAC will have reduced prediction accuracy in high-density areas where historically vessels traveled in different directions. We think a hybrid approach between the LRMAC and the SPNS may lead to improved prediction accuracy and execution time. Future work includes exploring the possibility of a hybrid approach between the LRMAC and SPNS.

Acknowledgement

We would like to thank MUNUS International for all the financial support given *(munus.ai).

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The data that support the findings of this study are available at https://doi.org/10.5281/zenodo.1167595 (Ray et al. 2019).

Additional information

Notes on contributors

Christiaan Neil Burger

Christiaan Neil Burger an MSc Computer Science student at Stellenbosch University, South Africa, in 2021. The contents of this paper form part of his Masters’ thesis. He is in the field of Data Science with a focus on data mining, prediction, and machine learning.

Waldo Kleynhans

Waldo Kleynhans received a Ph.D. (Electronic Engineering) from the University of Pretoria, South Africa as well as an MBA from Heriot-Watt University, Scotland in 2012 and 2011 respectively. He is a statistical signal processing and machine learning expert focusing on various application areas including maritime domain awareness and telecommunications. He holds adjunct faculty positions with the University of Pretoria, South Africa as well as San Diego State University, USA.

Trienko Lups Grobler

Trienko Lups Grobler received his PhD in Engineering at the University of Pretoria in 2013. He was a postdoctoral fellow at Rhodes University from 2013 to 2016. He joined the Computer Science Division of Stellenbosch University in 2017. His research interests include remote sensing and interferometry.

Notes

1. For more on the International Maritime Organization (IMO), see: www.imo.org . For additional information on AIS transponders see: https://www.imo.org/en/OurWork/Safety/Pages/AIS.aspx.

2. For more on PostgreSQL see: www.postgresql.org

3. For more on PostGIS see: www.postgis.net