{'fraction_non_alphanumeric': 0.034662691017517705, 'fraction_numerical': 0.00760342899739098, 'mean_word_length': 4.085481425322214, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
An Alternative Metric for Detecting Anomalous Ship Behavior Using a Variation of the DBSCAN Clustering Algorithm

There is a growing need to quickly and accurately identify anomalous behavior in ships. This paper applies a variation of the Density Based Spatial Clustering Among Noise (DBSCAN) algorithm to identify such anomalous behavior given a ship's Automatic Identification System (AIS) data. This variation of the DBSCAN algorithm has been previously introduced in the literature, and in this study, we elucidate and explore the mathematical details of this algorithm and introduce an alternative anomaly metric which is more statistically informative than the one previously suggested.


Introduction
In this paper we identify anomalous behavior in ships given their Automatic Identification System (AIS) data. AIS data is reported by all ships, and among other things, these data include a ship's position (latitude & longitude), speed, and course (direction) over time. We use this data to identify spatial and behavioral patterns of ships. With these patterns identified from the training data, we can potentially spot anomalous behavior in other ships given their new/incoming AIS data. We do this by applying a variation to the DBSCAN (Density Based Spatial Clustering Among Noise) algorithm ( [6]). The DBSCAN algorithm is used to create spatial clusters of data (see [5], [17], [22]), and because of its simplicity and speed, it is becoming an increasingly popular method used to identify anomalous behavior in ships (see [1], [11], [12], [15]). The variation to the DBSCAN algorithm we consider was introduced by Liu et al. ([9]). In this paper, we review, correct, and expand on the anomaly detection methods given in Liu ([9] & [10]), and we also offer an alternative and more statistically informative way of identifying anomalous behavior. This paper begins by showing the data and briefly discussing the problem at hand. This is done in Section 2. Simply stated, the problem is this: "How does one identify anomalous behavior among the ships that are reporting AIS data?" An easy way to answer such a question would be to create spatial clusters of AIS positional data using the DBSCAN algorithm. If a new ship then reported AIS positional data which was geographically far from any of these clusters, such data may be considered unusual, or anomalous. Methods such as these have been applied in a variety of maritime settings. Lee et al. ([8]) applied it to ships in port, and Wang et al. ( [20]) applied it to ships in waterways. There is so much data to create these positional/spatial clusters that Chen et al. ( [3]) studied how to form such clusters on the Hadoop platform.
But what if one wished to detect anomalies with respect to other variables, such as speed or direction, in addition to the positional/spatial variables? For example, how would it be possible to identify if a ship were traveling in an unusual direction and/or at an unusual speed, yet was (geographically) very close to a cluster of the training data? Capabilities have been created to separately address how to detect anomalies in speed (see [19] and [21]) or anomalies in direction (see [21]), but many of these methods are only suited to identify specific (predefined) departures from normal behavior (see [14]). Spatial clustering algorithms that separately include speed or direction have even been used to identify resting or stationary points (often called "stops") of ships. The clustering algorithm given in Palma ([13]), for example, identifies "stops" as places where ships are close together and have low speed, and the clustering algorithm given in Rocha ([16]) identifies "stops" as places where ships are close together and have multiple directions (if ships are at rest they can have multiple directions and not endanger one another). Liu ([9] & [10]) and Kontopoulos et al. ([7]) jointly cluster on location, speed and direction in an attempt to identify anomalous behavior with respect to any and/or all of these variables. Kontopoulos et al. ([7]) identify anomalous behavior visually; they form convex hulls with the clusters they create, and any ship that has a trajectory which goes beyond the borders of these hulls is said to exhibit anomalous behavior. Liu et al. ([10]) attempt to identify anomalous behavior by reporting a single number. The statistical significance of this number, however, is not evident. It is thus possible for Liu's statistic to report "anomalous" behavior yet (given the training data) not be uncommon. This is the problem we correct for. In this paper, we offer an alternative metric which clearly conveys the statistical significance of the anomaly score, and we prove that the asymptotic distribution of this metric is normal. An additional feature to the asymptotic distribution of our proposed metric is that it is invariant to the distributions of location, speed, or direction in the training data.
To explain how we create this metric, we first review how Liu et al. ( [9]) modify the DBSCAN algorithm to create appropriate clusters. This is done in Section 3. Section 3 also states and clarifies some mathematical assumptions that were made in their methodology. Section 4 then explains how these clusters can be used to test for anomalous behavior. Section 4 specifically discusses (and occasionally corrects for) how Liu ([10]) use these clusters to create gravity vectors and stationary sampling points. These vectors and points are summaries of the training data that are necessary in calculating an anomaly score of new/incoming trajectories. The statistical properties of Liu's anomaly score, and those of the new anomaly score that we propose are also discussed in Section 4.

The Data
In this paper, we look at AIS data that was reported on January 1 of 2017 along the mid-Atlantic coast of the United States. This data is publicly available at MarineCadastre.gov. These data are shown in Figure 1, and this is the data set we will use for training purposes. We will denote the training data as P and assume there are n observations in this data set, making P = z P i : i = 1, 2, . . . , n , where z P i = y P i , x P i , s P i , c P i , y P i is the latitudinal position of the i th observation, x P i is the longitudinal position of the i th observation, s P i is the "s"peed of the i th observation, and c P i is the "c"ourse of the i th observation. 1 Figure 1, of course, only shows the latitude and longitude of the training data.
The goal of this paper is to identify patterns among this training data using as few assumptions as possible. Any other data which exhibits significant departures from these patterns will be considered an anomaly.
A way to identify spatial patterns in the data while making no distributional assumptions is to apply the DBSCAN algorithm. Although this is not the algorithm we ultimately use to create our clusters, we do use a close variation of it. Section 3 reviews the DBSCAN algorithm and then provides details on how Liu et al. ( [10]) altered it to accommodate our needs. Recall that our needs involve clustering the data not just with respect to position, but with respect to position, speed, and course.

Clustering the Data
To cluster spatial data using the DBSCAN algorithm, two parameters need to be specified, and Nmin. The parameter specifies the maximum distance each observation in a cluster is from another observation in that same cluster, and Nmin specifies the minimum number of objects in a cluster. The following definitions (borrowed directly from [6] yet applied to P) will be used to define the DBSCAN algorithm with more mathematical clarity. Definition 1 A point z P i is directly density-reachable from an object z P j with respect to and Nmin in the set of objects P if 1. z P i ∈ N z P j , where N z P j is the subset of P contained in the neighborhood of z P j . 2. |N z P j | ≥ Nmin, where |N z P j | is the cardinality of the set N z P j . To visualize the concept of two points being directly density-reachable, imagine the point z P j and all of the objects in the set P that are within from z P j . If there are at least Nmin objects that are within of z P j , and z P i is one of those elements, then z P i is directly density-reachable from the object z P j . Figure 2 illustrates objects that are directly density-reachable.
z P j z P i Figure 2: z P i is Directly Density-Reachable from z P j with respect to and N min = 5.
Definition 2 An object z P i is density-reachable from an object z P j with respect to and Nmin in the set of objects P if there is a chain of points p1, p2, . . . , pn, with p1 = z P j and pn = z P i such that pi ∈ P ∀ i and pi+1 is directly density-reachable from pi with respect to and Nmin.
Density-reachability is different from two points being directly density-reachable in the sense that it implies that there is a sequence of directly density-reachable points from z P j , and z P i is directly density reachable from one of the points in that sequence. Figure 3 illustrates objects that are density-reachable from one-another. Figure 3: z P i is Density-Reachable from z P j with respect to and N min = 5. In this case, n = 9, and observe that p i+1 is directly density-reachable from p i for all i.

Definition 3 An object z P
i is density-connected to an object z P j with respect to and Nmin in the set of objects P if there is a point z P o ∈ P such that both z P i and z P j are density-reachable from P with respect to and Nmin in P. Figure 4 illustrates the concept of two elements in P being density-connected. In the DBSCAN algorithm, all elements in the same cluster are density-connected. Another (and perhaps more intuitive) way to communicate this is to say that for every element in a DBSCAN cluster, there is at least one other element in that same cluster which is units away. And it must be the case that for at least one of the elements that is units away, there are at least Nmin observations within of it.
Applying the DBSCAN algorithm to the training data shown in Figure 1 with = .02, Nmin = 5, and dist z P we get the results shown in Figure 5. Note that the distance we calculate between two points is the Euclidian distance between the two points' latitude and longitude coordinates. Our methodology thus assumes a flat earth. The more locally we apply our algorithm, the more valid this assumption. The results in Figure 5 cluster points based on their location. We wish to cluster with respect to location, speed, and direction, however. Let us specifically assume that we wish to cluster such that points z P i and z P j are in the same cluster if dist z P i , z P j < Dist, c P i − c P j < Crs, and s P i − s P j < Spd . To do this, an edited version of the DBSCAN algorithm is necessary since DBSCAN is only engineered to consider one distance between two points. Liu et al. ( [9]) edit the DBSCAN algorithm to accommodate these wishes. Their edited version of the DBSCAN algorithm is called DBSCANSD, where the "SD" stands for "S"peed and "D"irection.
The DBSCANSD algorithm applied to the training data set P is given below and requires that the additional thresholds Crs and Spd be specified. Liu et al.'s ( [9]) addition to the standard DBSCAN algorithm is shown in red.
We apply the DBSCANSD algorithm to the "moving" members of P, P mv , and the DBSCAN algorithm to the stationary members of P, P st to get the results shown in Figure 6. It should be noted that in this application, the speed threshold used was 2.5 knts and the direction threshold used was 90 degrees ( Spd = 2.5 knts and Crs = 90 degrees) if cluster.size ≥ N min then Mark z P i as core point j is core point and clusterA contains cluster Q then merge ← TRUE for each point z P l in clusterB do clusterA.add (z P l ) Figure 6: DBSCANSD Results. The different clusters are in different colors.
From this point on, we will denote the set of points put in cluster j as a result of the DBSCANSD algorithm as P mv, cl j , and we will assume that the number of observations in this cluster is n mv j , making P mv, cl j = z mv, cl j 1 , z mv, cl j 2 , . . . , z mv, cl j n mv The set of points put in cluster l as a result of the DBSCAN algorithm will be denoted as P st, cl l , and we will assume that the number of observations in this cluster is n st l , making The next section discusses how we use these results to identify outliers or outlying/anomalous behavior in ships.

Identifying Anomalous Behavior
To identify anomalous behavior in a new ship's trajectory, one first has to separate this new trajectory into a set of stationary points and a set of moving points. Just as Liu et. al ( [10]) did, we identify the set of stationary points as that set such that the speed is less than 0.5 knots. The set of moving points is the complement of that. We will assume that there are a total of m points in this new trajectory, and we will denote the set of points in this trajectory as N (for "N"ew), where The set of points in N that are moving is and the set of stationary points is We will assume there are m mv values in N mv and m st values in N st , making We then see how the points in N st depart from stationary points in the training data set, and how the points in N mv depart from the moving points in the training data set. To do this, Liu et al. ( [10]) first create two sets of points, one set which summarizes the stationary points in the training data, the other set which summarizes the moving points in the training data. These sets are respectively called the stationary sampling points and gravity vectors, and we denote these sets as S and G. The set of new trajectory points, N , are then compared to S and G (the set N st is compared to S, and the set N mv is compared to G), and it is from this comparison that a trajectory is identified as being anomalous or not. Subsection 4.1 describes how the gravity vectors and stationary sampled points are calculated. Subsection 4.2 describes how these set of points, S and G, are compared to the new trajectory, N . The subsection specifically reviews how Liu et al. ( [10]) calculate and assign an anomalous score to a new trajectory and then discusses our alternative anomalous score.

Creating Gravity Vectors and Stationary Sampled Points
Generally speaking, a gravity vector is a point (or vector) that is meant to summarize and describe all of the moving points of the training set around it. We will let G cl j be the set of gravity vectors which summarize cluster j, and we will assume there are n cl j grv gravity vectors associated with cluster j, i.e., To calculate all of the gravity vectors associated with cluster j, one first has to calculate the average course (direction) of the entire cluster. We will call this average direction c mv, cl j and calculate it as After calculating c mv, cl j , one must define a line along this direction and divide this into segments of length δ. Liu et al. ( [10]) set δ = , and we do the same to achieve our results. All of the observations that are then within a particular band (of width δ) are considered and their average latitude, longitude, speed, direction, and median distance from the average position are calculated and reported in one gravity vector. This concept is illustrated in Figure 7, and in this paper we mathematically formalize the calculation of these vectors. This is done in the text box below.  The line along the direction of c mv, cl j is of length L mv, cl j , where L mv, cl j = max i y mv, cl j i / cos c mv, cl j − min i y mv, cl j i / cos c mv, cl j .
We are dividing this line into bands of width δ, making the number of gravity vectors in cluster j, n cl j grv , n cl j grv = L mv, cl j δ.
The observations in cluster j to be considered in the calculation of the l th gravity 1 ≤ l ≤ n cl j grv vector are then where b = mini y mv, cl j i cos c mv, cl j . In words, the observations in P cl j grv l are those in the l th slice of the moving points in cluster j. The direction and orientation of this slice is governed by cos c mv, cl j . The l th gravity vector in cluster j is then the average location, speed, and course over this set. It also includes the median distance of each observation to the location of the gravity vector. We denote the l th gravity vector in cluster j as Stationary sampling points are meant to describe the stationary points in the training set. For each stationary cluster, the number of stationary sampled points meant to summarize it are N/ , and they are randomly selected according to the following algorithm: Algorithm 2: Extracting Stationary Sampling Points from P st, cl j input : P st, cl j , Nmin, Dist output: S cl j S cl j ← empty Lat1, Lat2 ← minimum and maximum of all points' latitude in P st, cl l Lon1, Lon2 ← minimum and maximum of all points' longitude in Dist count ← 0 while count < sample size do Randomly select one point from cluster P st, cl l if Randomly selected point is far from all points in S cl j then Add point to S cl j count = count + 1 Assume there are n cl j ssp stationary sampled points in stationary cluster j. We will call this set of points S cl j , and S cl j = s cl j 1 , s cl j 2 , . . . , s cl j n cl j ssp , where s cl j l = y cl j ssp l , x cl j ssp l , s cl j ssp l , c cl j ssp l . It is with the gravity vectors and stationary points that an anomaly score is calculated. The subsection below explains how Liu ( [10]) calculates this anomaly. We add some mathematical and statistical detail to their calculations and also introduce an alternative and more flexible way to measure anomalous behavior.

Calculating Anomalous Behavior
To assign an anomaly score to N , Liu et al. ( [10]) first split the new track into its stationary and moving parts, N st and N mv . For each point in N st , they calculate the smallest distance between it and the set of stationary sampled points. This distance is called the Absolute Distance Deviation (ADD), and for point i in N st , it is calculated as For each point in N mv , they calculate two distance metrics, the Relative Distance Deviation (RDD) and the Cosine Division Distance (CDD). The RDD is similar to the ADD in that it calculates the smallest distance between a point and the set of gravity vectors, but this metric is different in that it accounts for the variation and geographical spread around the gravity vector. It does this by dividing the distance by the median of the associated gravity vector. For the i th point in N mv , RDD is calculated as The CDD accounts for any difference in heading and/or speed a point in N mv may have from the closest gravity point. For point i in N mv , it is calculated as where α = c * − c N mv i , and c * and s * are the course and speed components of gravity vector g * (the closest gravity vector), where The equation in (1) looks at the ratio of the two different speeds and the magnitude in the difference of their directions. If the speeds and courses are identical, CDD N mv i = 1, and the smaller its value, the greater the anomaly.
It should be noted that our definition of CDD is different than how it is written in Liu et al. ( [10]). They write CDD as but this definition seems unclear. As they have written it, they are looking for a gravity vector with a speed and course which most closely matches the speed and course of point i in N mv . The location of this gravity vector is not considered. Given their definition of CDD, it would be possible for the point i in N mv to be in the proximity of points going in an opposite direction and at a much different speed, yet still have a CDD value that went unnoticed. With this in mind, we modified/rewrote the definition of CDD. For each point i in N mv , we evaluate CDD at the closest gravity vector. The closest gravity vector is that which minimizes RDD (g * in Equation 1). With these three metrics (ADD, RDD and CDD) calculated, Liu et al. ( [10]) score each observation in N depending on whether the calculated ADD, RDD, or CDD are beyond a certain threshold. For points in N st , the score is calculated as For points in N mv , the score is calculated as They then calculate the total anomaly score for the new trajectory, N , as One way they obtain these three thresholds, ADD Threshold , RDD Threshold , and CDD Threshold is by considering an entirely different data set, D (for "D"ifferent), calculating the distribution of ADD, RDD and CDD values in this data set, and then letting ADD Threshold and RDD Threshold be the 95 th percentile of the distribution in the ADD and RDD values, and CDD Threshold be the 5 th percentile of the distribution in the CDD values. This is mathematically formulated in the textbox below.

Calculating Threshold Values
We assume D has r observations, r st which are stationary and r mv which are moving. From this data set, we calculate r st values of ADD, written as and r mv values of RDD and CDD, written as ADD Threshold is the 95 th percentile of ADD D , RDD Threshold is the 95 th percentile of RDD D , and CDD Threshold is the 5 th percentile of CDD D .
A setback to the method that Liu et al. use to measure anomalous behavior is that it fails to highlight the extremity of the anomaly. For example, if N were a set of stationary points, all of which were just barely beneath ADD Threshold , its anomaly score would be 0. Its anomaly score would also be 0 if all these stationary points were significantly below ADD Threshold . Yet another setback to Liu et al.'s method is that the statistical significance of their anomaly is not immediately transparent. For instance, if AnomLiu (N ) = 0.3, it is not obvious from this statistic what the probability is of observing an anomaly as or more extreme than the one observed.
Part of the novelty proposed in this paper is in how we calculate the anomaly of a new trajectory. We calculate the anomaly in such a way that the two extreme cases described above would have considerably different scores. The statistical significance of our anomaly score is also transparent (we ultimately report a z−score).
We begin by scoring each observation not with a 1 or a 0 (as Liu et al. did), but with the fraction of ADD D , RDD D , and CDD D values that are more extreme than the one observed. The smaller this fraction, the more unusual/extreme an observation. These scores are written below in Equations 2 and 3. Theorem 1 in Section A of the Appendix calculates the expected value and variance of AnomLiu (N ). Theorem 2 in Section A of the Appendix justifies the asymptotic normality of AnomBotts (N ). In Section B of the Appendix, some simulation results are provided which illustrate this asymptotic normality and its invariance to the distribution of ADD, CDD, and RDD.
We would expect normal/in-family trajectories to have AnomLiu (N ) and AnomBotts (N ) values close to these expected values. For anomalous trajectories, we would expect AnomLiu (N ) to be considerably larger than its expected value, and AnomBotts (N ) to be considerably less than its expected value. Remember that AnomLiu (N ) counts the fraction of points in N that are beyond a certain threshold. The larger the value of AnomLiu (N ), the more anomalous a trajectory is. AnomBotts (N ) considers the fraction of values in ADD D , RDD D , and CDD D that are more extreme than those observed in N and turns this into a z−score. The smaller the value of AnomBotts (N ), the more anomalous a trajectory. It is thus expected that AnomLiu (N ) and AnomBotts (N ) will jointly go in opposite directions of their respective expectations. Figures 8 and 9 show normal trajectories with values of AnomLiu (N ) and AnomBotts (N ) that are close to their expected values. Figures 10 and 11 show abnormal trajectories. In these cases, AnomLiu (N ) is higher than expected, and AnomBotts (N ) is significantly lower than expected. Assuming the ADD, RDD, and CDD values in N come from the same distribution as those in the training set and D, the probability of observing AnomBotts (N ) scores as or more extreme than the ones observed in In both of these figures, the difference between the two metrics illustrates why AnomBotts (N ) may be a preferred metric to AnomLiu (N ) .
In Figure 12, N is made up entirely of stationary points, and they are all (nearly) at the same location. Each point has an ADD that is approximately .017, which is beneath the ADD threshold of .034. Since they are all beneath the threshold, AnomLiu (N ) = 0. AnomBotts (N ), however, does capture how extreme these values of ADD are, since it does not depend on a threshold and merely counts the number of ADD D values that are greater than it.
The same principle applies in Figure 13. In this case, very few CDD values in the moving part of N are below the CDD threshold of -.746, and no RDD values are above the RDD threshold of 2.255, yet 62% of N 's RDD values are within the 95 th and 70 th percentile of the RDD values in RDD D .

Conclusion
In this paper, we illustrate, correct, and clarify the DBSCANSD clustering algorithm presented in Liu et al. ( [9] ). The DBSCANSD algorithm is one of few algorithms that allows one to cluster AIS data based on location, speed, and course. With these clustering results, we offer an alternative to the anomaly metric presented in Liu et al. ( [10]). The statistical significance of the metric we propose is transparent (unlike the one in Liu et al. ( [10])), and the asymptotic distribution of our proposed statistic is invariant to the distributions of course, speed, or location in the training data. In the future we hope to consider other variables such as heading and/or time in the clustering algorithm. Perhaps they can also be considered when trying to detect anomalous behavior in ships.

A Proof of Theorems
Assume these sets of random variables are independent of one another and also assume that the variables within each set are independent of another. We will also assume that all values in V D and V N follow the common distribution fV (v), all the values in Y D and Y N follow the common distribution fY (y), and all the values in Z D and Z N follow the common distribution fZ (z). With these definitions, consider the statistic T1(α), where: 2 ). With uniform integrability, we can apply expectations to get Proof. For calculating the asymptotic distribution of T2, we first have to remember that from the Dvoretsky-Kiefer-Wolfowitz inequality ([4]), we get Since convergence in probability implies convergence in distribution, we get The quantities S1 and S2 thus, respectively, converge in distribution to the values S * 1 and S * 2 as The expected value and variance of A are 1 3 and 1 18 , respectively, implying (by the central limit theorem)

B Simulations
In this section, we perform simulations which illustrate the asymptotic result proven in Theorem 2.