4.1. Feature model of vessel trajectory

In order to build the vessel trajectory feature model for visual analysis, a trajectory clustering algorithm is applied. It is based on the tree clusters method proposed by Piciarelli et al. (Reference Piciarelli, Foresti and Snidaro2005). There are two advantages of this method. Firstly, we improve the Euclidean distance measurement (when a match between a trajectory and a cluster is found, the cluster must be updated). This way each cluster is a dynamic approximation of the mean and variance of trajectories that matched it, with an exponentially decreasing weight for the older trajectories. Secondly, we use tree clusters, which eliminates the trajectory division processes without loss of local cluster features. The clusters are dynamically built in real-time as the trajectory data is acquired, saving the trouble of off-line processing.

However, in order to gain better trajectory clustering results from AIS data and increase clustering efficiency, the following improvements are implemented. Firstly, to reduce invalid matches, a method for trajectory pre-classification is proposed: we select the starting points of all trajectories clustered by a Density-Based Spatial Clustering of Applications with a Noise (DBSCAN) algorithm (Ali et al., Reference Ali, Asghar and Sajid2010) and classify the trajectories based on their starting points. The classification can be used for future trajectories' clustering and selection. This method can be regarded as building indices for trajectory clustering, which will significantly reduce the matching failure rate. Secondly, a new method of updating clustering is presented to support the trajectory contribution rate: adding a support ct _j to each feature point j of the cluster. When the cluster C only consists of one trajectory, the ct _j of each point is 1. Each time a feature point is updated, we add 1 to ct _j. This support can also yield the overall arithmetic mean of the whole pool of parameters when the parameters of every point are updated.

The trajectory is defined as follows:

(1)$$T_i = \{t_{i1}, \ldots, t_{in}\}, \hbox{where } t_{ij} = \{x_{ij}, y_{ij}\}$$

Each trajectory T _i is represented by a list of vectors t _ij which gives the spatial position of the vessel i at time j.

The improved cluster is defined as Equation (2).

(2)$$C_i = \{ c_{i1}, \ldots c_{in}\} \,\hbox{where } c_{ij} = \{ x_{ij},y_{ij},\sigma _{ij},ct_{ij}\} $$

where C _i denotes cluster i, and c _ij is the feature point of the cluster. Every feature point comprises the trajectory point (x _ij, y _ij), local approximation σ_ij of the cluster variance (Piciarelli et al., Reference Piciarelli, Foresti and Snidaro2005) and support ct _ij at time j.

In order to check if a trajectory fits a given cluster, a distance measure must be defined. The distance between a trajectory T = {t ₁…t _n} and a cluster C = {c ₁…c _n} is defined as Equation (3), which is as in Piciarelli et al. (Reference Piciarelli, Foresti and Snidaro2005) and is an improved Euclidean distance measure.

(3)$$D(T,C) = \displaystyle{1 \over n}\sum\limits_{i = 1}^n d (t_i,C)$$

where

(4)$$\eqalign{& d(t_i,C) = \mathop {\min }\limits_j \left( {\displaystyle{{dist(t_i,c_j)} \over {\sqrt {\sigma _j} }}} \right) \cr & j \in {\rm \{ \lfloor (1 - \delta )i \rfloor \ldots \lceil(1 + \delta )i \rceil \} }} $$

dist(t _i, c _j) is the usual Euclidean distance. The distance of a trajectory from a cluster is thus the mean of the normalised distances of every trajectory element t _i from the nearest cluster element c _j found inside a temporal window centred in i and with a variable size, in order to permit matching under accumulating temporal differences. The window is clipped in the range [1 … m] and, if it completely falls outside this range, the distance D is set to ∞.

However, a different method is used to update the cluster when a match between a trajectory and a cluster is found. The improved method uses arithmetic mean instead of a fixed parameter α (Piciarelli et al., Reference Piciarelli, Foresti and Snidaro2005). If elements t _i = {x _i, y _i} and c _j = {x _j, y _j, σ _j, ct _j} match (c _j is the cluster element that is cluster feature point, nearest to t _i inside the temporal window), then c _j is updated as follows:

(5)$$\left\{ {\matrix{ {x_j = \left( {\displaystyle{{ct_j} \over {ct_j + 1}}} \right)x_j + \left( {\displaystyle{1 \over {ct_j + 1}}} \right)x_i} \hfill \cr {y_j = \left( {\displaystyle{{ct_j} \over {ct_j + 1}}} \right)y_j + \left( {\displaystyle{1 \over {ct_j + 1}}} \right)y_i} \hfill \cr {\sigma _i = \left( {\displaystyle{{ct_j} \over {ct_j + 1}}} \right)\sigma _j + \left( {\displaystyle{1 \over {ct_j + 1}}} \right){(dist(t_i,c_j))}^2} \hfill \cr {ct_j = ct_j + 1} \hfill \cr } } \right.$$

where dist(t _i, c _j) is the same as in Equation (4).

4.2. Visual analysis of vessel trajectory

Since vessel behaviour is a study on navigation features of vessel groups, not on those of an individual vessel, a visual analysis approach based on clustering vessel trajectories (Clustering Visualisation for short) can be better for showing the vessel trajectory distribution and vessel navigation trends.

Clustering Visualisation is a method for visualising the results of clustering analysis of the vessel trajectories, which is described in Section 4.1. The result is a set SSR = {StR ₁, StR ₂, …, StR _m} which consists of multiple subsets; each subset is a forest of clusters. Visualising the result means plotting every cluster in result set SSR on a map. However, it is necessary to distinguish the trajectory support and the degree of clustering during the plotting. “Cluster” is similar to trajectory and can be represented by a series of point coordinates, but it also contains the support ct _ij and local approximation σ_ij for each point, as shown in Equation (2). Support ct _ij is the number of trajectories in the cluster while σ_ij represents the average variance of the trajectories around the cluster centre. Therefore, support ct _ij reflects the frequency of the trajectories, and local approximation σ_ij the clustering degree of the cluster, can also be regarded as the extent of the trajectory close to the cluster centre. In this paper, clustering degree is distinguished by different trajectory widths, and the frequency is distinguished by colour depths. The description of algorithm of Clustering Visualisation (TrajVisualise for short) is shown in Table 1.

Table 1.

Description of TrajVisualise.

4.3. Interactive visualisation model for analysing vessel trajectory

Visualisation technology cannot provide users with dynamic analysis, so human-computer interaction methods are also added. Based on our algorithm TrajVisualise, an interactive visualisation model for analysing vessel trajectory is defined as follows.

1. (1) Positioning starting point for trajectory cluster subsets. By pre-classifying, trajectories are divided into different subsets. Users can forecast the trajectory direction by fixing the position of the start point of a trajectory cluster subset using this interactive visual analysis approach.

2. (2) Positioning ending point for trajectory cluster subsets. Trajectories which have the same starting point may move towards different directions, which forms a tree structure of trajectory clusters. However, the trajectories which have the same end point usually belong to the same course. Analysis of a trajectory cluster with a particular end point is helpful for the analysis of a particular course.

3. (3) Setting threshold in the clustering algorithm. Setting the threshold for trajectory similarity is the key to the clustering effect. Through human-computer interaction to change the threshold, users can more intuitively and quickly observe the clustering effect under different thresholds and enable the selection of a more suitable threshold.

5.1. Feature model of vessel density distribution

Vessel density refers to the number of vessels per unit water area, such as per square kilometres or per 100 square kilometres. However, vessel density distribution is different from vessel density, which represents changes of vessel density across different locations. Vessel density distribution is very important in analysing waterway condition and maritime accidents.

In this paper, kernel density estimation is used to model vessel density distribution. Vessel density distribution can be treated as the distribution of vessels in 2D. The kernel density estimation function is defined as Equation (6).

(6)$$\hat f_{2D}(x) = \displaystyle{1 \over {nh_{1,1}h_{2,2}}}\sum\limits_{i = 1}^n K \left( {\displaystyle{{(x - x_i)} \over {h_{1,1}}},\displaystyle{{(y - y_i)} \over {h_{2,2}}}} \right)$$

where n is the number of points whose distance from point (x, y) is less than h _{1, 1} along the x-axis and less than h _{2, 2} along the y-axis. h _{1, 1} and h _{1, 2} represent the attenuation threshold along the x-axes and y-axes respectively. K(x, y) is a two-dimensional Gaussian kernel.

Based on the feature of kernel density estimation, this paper focuses on the following problems when modelling vessel density distribution.

1. (1) A kernel density estimation matrix requires too much computation power: to reduce the memory and Graphics Processing Unit (GPU) capacity consumption during plotting, grid partitioning is adopted. Grid partitioning divides the water area into equal-sized grids, and separately calculates a kernel density estimation matrix for each grid.

2. (2) How to obtain instantaneous density information: drawing a heat map for vessel density distribution requires overlaying vessel density information from each time entry. However, the frequency of AIS message transmissions depends on vessel speed, so different vessels' position information is not synchronised. Thus, in the dynamic AIS database, position information for a vessel at time t is queried by combining the data from a range of predefined time zones centred at time t. Then the closest position to t is chosen as the instantaneous position information. Since approximately 5 minutes is the maximum gap between AIS reports, a range of 300 seconds is chosen.

3. (3) Removing the repeated samplings for vessel position information: the data for each vessel position is inserted into a dynamic AIS data table as soon as it is obtained. The number of the grid in which the vessel is located at time t can be acquired by searching the records which have the same Maritime Mobile Service Identity (MMSI) and are closest to time t in the table. Thus, using this method we can avoid the problem that a vessel passes through several grids during the same time quantum.

5.2. Visual analysis of vessel density distribution

A heat map is the most common and effective way to express density distribution. In this paper, a density distribution heat map is plotted using kernel density estimation which has the advantage that users can interact with a heat map in more depth by changing the attenuation threshold of kernel density estimation.

The procedure of the visual algorithm for vessel density distribution analysis based on kernel density estimation is as follows:

1. (1) Mesh the water area, and obtain grids of equal size;

2. (2) Calculate the kernel density estimation matrix for each grid in turn;

3. (3) Take the estimated value of the kernel density estimation matrix at each point and map it with the heat map colour;

4. (4) Mark the points from the kernel density estimation matrix in the map;

5. (5) Plot a heat map for each grid, and then combine them to form a heat map for the overall water area.

5.3. Interactive visualisation model for analysing vessel density

Since kernel density estimation is used to model the features of vessel density distribution, an interactive analysis model is built by means of setting the attenuation threshold. The only parameter h, the attenuation threshold, in kernel density estimation, influences the degree of detail in the heat map. The bigger the attenuation threshold h, the smaller the difference of density estimation from neighbouring points. The smaller the attenuation threshold, the bigger the difference of the density estimation from the neighbouring points. Therefore, changing the attenuation threshold can effectively change view precision and deliver detail information of different vessel densities.

6.1. Feature model of vessel speed distribution

6.1.2. Feature model of vessel speed spatial distribution

First of all, the water area is meshed. Then the model of vessel speed distribution based on this grid is built as follows:

(7)$$\bar v = \displaystyle{{{\bar v}_1 + {\bar v}_2 + {\bar v}_3 + \cdots + {\bar v}_n} \over n}$$

where $\bar{v}_{1}$, $\bar{v}_{2}$, $\bar{v}_{3}\comma \; \ldots\comma \; \bar{v}_{n}$ are the average speeds of single vessels, n represents the total number of vessels in the current water area and $\bar{v}$ represents the average speed.

According to the Equation (7) and the dynamic AIS data table, Equations (8) and (9) can be used to calculate the average vessel speed in the gridded water area.

(8)$$\matrix{{Average \,\,speed \,\,of \,\,a \,\,vessel \,\,in \,\,a \,\,grid = \displaystyle{{\sum {speed \,\,of \,\,the \,\,vesse{l}^{\prime}s \,\,records \,\,in \,\,the \,\,grid} } \over {The \,\,number \,\,of \,\,records}}} \hfill \cr } $$

(9)$$\matrix{{Average \,\,speed \,\,of \,\,vessels \,\,in \,\,a \,\,grid = \displaystyle{{\sum {average \,\,speed \,\,of \,\,every \,\,vessel \,\,in \,\,the \,\,grid} } \over {The \,\,number \,\,of \,\,vessels}}} \hfill \cr } $$

6.2. Visual analysis of vessel speed distribution

6.3. Interactive visualisation model for analysing vessel speed

In this section we aim at visually analysing vessel speed spatial distribution and vessel speed probability distribution. The interactive visualisation model is defined as follows:

1. (1) Setting the grid parameters. In our model of vessel speed spatial distribution which is based on grid-enabled colour coding, grid height hf and grid width wf influence the accuracy of the visualisation result and processing efficiency. When wf and hf are too large, the visualisation result does not reflect the gradient variation of vessel speed in different areas accurately, and when wf and hf are too small, the time consumed in dividing the grid and computing average vessel speed is too high. Therefore, users need to adjust the grid scale to achieve different views of speed distribution.

2. (2) Setting the speed level. Expressing speed gradient variation by different colours is helpful for users in analysing speed spatial and probability distribution. Users can set the colour coding cl, that is, the level of the speed level boundaries, based on their preference so that different visualisations can be shown.

3. (3) Browsing grid details. Each grid has its histogram of vessel speed probability distribution in the view. To ensure clearness of view, the view is only rendered by colour but does not show values by default. The speed probability distribution cannot be displayed in the same view as speed spatial distribution to avoid one data set obscuring the other. The details of vessel speed spatial distribution and vessel speed probability distribution can be shown by interactive display.

7.1. Universal model of interactive analysis

In order to obtain better interaction with the visual analytics system, in addition to the interactive visualisation model for analysing the characteristics of vessel behaviour, a universal interactive model is constructed in the prototype system. The model includes Data presetting, Switch of visual analysis view and Universal interactive operations.

1. (1) Data presetting. Data presetting is to obtain the data sets filtered by the user that need to be analysed. Related operations include filtering the data by time interval, vessel type and geographic location.

2. (2) Switch of visual analysis view. Different views give an analysis of different characteristics of vessel behaviour. A user can switch between multiple views.

3. (3) Universal interactive operations. Basic interactive view operations are supported, including zooming in and out, refreshing, and dragging.

7.2. Architecture of prototype system

The prototype of our visual analysis system for vessel behaviour consists of three subsystems: Visual analysis of vessel trajectory, Visual analysis of vessel density, and Visual analysis of vessel speed distribution. These are built by feature model and interactive visualisation model, which are illustrated in Sections 4, 5 and 6 respectively. The architecture of the visual analysis of vessel trajectory subsystem is shown in Figure 2. This subsystem consists of the following components:

1. (1) Data preprocessing: cleansing data and extracting subset of AIS data interactively;

2. (2) Trajectory feature modelling and visualisation mapping: pre-processing trajectories based on DBSCAN; mining trajectories by means of tree clustering algorithm, and visual mapping;

3. (3) Visualisation: rendering the results of trajectory clustering;

4. (4) Interaction setting: realising interactive analysis, including pre-setting data, setting clustering threshold, and locating trajectory sub-clusters.

Figure 2.

Architecture of Visual analysis for vessel trajectory.

8.1. Experimental environment and data

A prototype system for visual analysis of vessel behaviour was realised according to Figure 1. For the experiment, we used a machine with the following configuration: a 2.80 GHz Intel^® Xeon^® CPU with 8 GB main memory and a 1TB hard disk. Operating System was Windows 7 x64. SQL Server 2008 R2 database was used to store AIS data and other processed data. To realise the related algorithms, we used the Java programming language. To implement the visual analysis of the vessel behaviour features we used the open source library Mapv of Baidu Map (http://huiyan.baidu.com/mapv/). The communication between front-end and back-end was ensured by a Javaweb framework and the application server was built using Apache-Tomcat-6.0.36.

As dataset for our experiment, we used historical AIS data collected from middle Hankou channel to Yangluo channel in the Yangtze River in 2015. The data was then cleaned and organised into three datasets D₁, D₂ and D₃ which are shown in Table 2. It should be noted that the AIS data for our experiment was collected in the framework of the project of Natural Science Foundation of China (No. 51479155) in 2015.

Table 2.

Experimental Data Samples.

8.2. Experimental results and analysis

We present three experiment results, as follows.

8.2.1. Setting cluster thresholds

The subsystem Visual analysis of vessel trajectory requires the user to set two cluster thresholds, σ and ε, which are the key parameters for the vessel trajectory cluster. The definition of σ is given in Equation (2), and ε is the threshold for D(T, C). The illustration of the effect of different values of these two thresholds on vessel trajectory clusters is depicted in Figure 3: (a) $\sqrt{\sigma} = 100\,\, \hbox{m}$, ε = 1.2; (b) $\sqrt{\sigma} = 150\,\, \hbox{m}$, ε = 1.2; (c) $\sqrt{\sigma} = 150\,\, \hbox{m}$, ε = 1.35; (d) $\sqrt{\sigma} = 150\,\, \hbox{m}$, ε = 1.5. It can be seen from Figure 3(a), 3(b) and 3(c) that the number of trajectories absorbed by the trajectory cluster will be more and their supports will be bigger as the cluster threshold $\sqrt{\sigma}$ and ε increase. However, trajectory clusters with different features will be together in a cluster and lose the discrimination between trajectories if the threshold ε is too large, just as the result shown in Figure 3(d). It can also be inferred that σ will become small so that those trajectories which should be in a cluster, will be clustered in different clusters if ε is too small, for example, if it is less than 1. Therefore, by setting cluster thresholds interactively, the user can change and analyse the result of trajectory clustering and choose suitable thresholds.

Figure 3.

Interactive analysis comparison of setting cluster thresholds.

8.2.2. Setting the attenuation threshold of kernel density estimation

In Visual analysis of vessel density distribution, the user can set the attenuation threshold h of kernel density estimation in the heat map.

As shown in Figure 4, the attenuation threshold h in the four figures are respectively 200 m, 400 m, 600 m and 1200 m. Through comparison and analysis, it can be seen that the effect of the heat map is obviously different because of the different attenuation thresholds. The smaller the threshold, the more detailed the density distribution of the heat map. The larger the threshold, the more the heat map can show the trend of density distribution of the whole channel.

Figure 4.

Comparison of attenuation threshold adjusting in a heat map.

8.2.3. Setting the grid parameters

In visual analysis of vessel speed distribution, users can set the grid parameters themselves. As shown in Figure 5, the meshed parameters wf × hf in (a) and (b) are 1000 × 1000 m and 400 × 400 m. From the two figures, we can obtain a view with finer meshing and more detailed speed distribution when the grid parameters are smaller. The view of vessel speed distribution with large grids provides more detailed information of local area speed distribution, while the view with small grids can give information of the trend of speed distribution in the water area in general.

Figure 5.

The result of changing the meshed parameter.

In summary, through the comparison of the results before and after changing various parameters interactively, the user can obtain the information of vessel behaviour more comprehensively, clearly and intuitively.