1. Introduction

2. Literature Review

3. Methodology：Dynamic Calculation of Navigation Risk

3.1. Probabilistic Conflict Detection

The CD approach is first developed from a probabilistic risk viewpoint, in which the influence of uncertainty related to various sources on the potential conflict is considered. It is designed to quantify the conflict criticality between ship pairs and provide essential functionality for extract ship group. The probabilistic CD approach is characterized by a conflict criticality measure and conflict probability estimation.

3.1.1. Conflict Criticality Measure

A conflict occurs when the predicted trajectories of two ships violate the specified separation distance. The paper defines ship conflict based on the MSD model [39] The method of conflict identification is shown in Figure 2. If the sum of the MSD of the two ships and the distance between the two ships within the prediction time satisfy equation (1), then ship A and ship B are considered to be in a conflict state:

$$D\left(t\right)\le MS{D}_A\left(t\right)+MS{D}_B\left(t\right),$$

(1)

Where D(t) is the distance between the two ships, $MS{D}_A\left(t\right)$and $MS{D}_B\left(t\right)$represent the distances between the ship centers and their dynamic domain boundaries, respectively. Because we focus on identifying ship conflicts with ship motion uncertainty, inequality (1) is a probabilistic issue.

Given that the probability density function (PDF) for the unsafe separation (i.e., $\left(D\left(t\right)\le MS{D}_A\left(t\right)+MS{D}_B\left(t\right)\right)$) between the two ships at time t is denoted as $f_{\gamma L\left(t\right)}$, the instantaneous occurrence probability of a conflict is expressed as follows:

$$PC\left(t\right)=\mathrm{Pr}\left(L\left(t\right)\le 0\right)={\displaystyle {\int }_{-\infty }^0{f}_{L\left(t\right)}}\left(\rho \right)d\rho ,$$

(2)

where L(t) equals $\left(D\left(t\right)-MS{D}_A\left(t\right)-MS{D}_B\left(t\right)\right)$.

The maximum conflict probability within the predicted range is used to measure the severity of the conflict as follows:

$$P_{con}\left(\tau \right):=\underset{t\in \left[0,T\right]}{\max }PC\left(t\right),$$

(3)

3.1.2. Conflict Probability Estimation

Probabilistic CD is based on the prediction of ship trajectories with uncertainty. Therefore, considering the ship motion model with uncertain mode of ship heading and ship position, we use the Monte Carlo simulation algorithm which can accurately calculate the conflict probability to extend CD to the probability framework.

The initial position of the two ships is determined, but the ship is affected by many uncertain factors, such as pilot’s maneuver and intention, wind and waves, mechanical failure and so on. Considering the fluctuation of the ship’s trajectory, the equation of ship’s motion is expressed as follows:

$${\dot{X}}_O=\left(x_O,y_O\right)+W_O=\left(V_O\sin {\phi }_O,V_O\cos {\phi }_O\right)+W_O,$$

(4)

$${\dot{X}}_T=\left(x_T,y_T\right)+W_T=\left(V_T\sin {\phi }_T,V_T\cos {\phi }_T\right)+W_T,$$

(5)

Among them, ${\phi }_O$,${\phi }_T$is the heading angle of $X_O$ and $X_T$ respectively. $W_O~\Psi (0,{\mathcal{Q}}_O)$、$W_T~\Psi (0,{\mathcal{Q}}_T)$ is the error of ship’s predicted position.${\mathcal{Q}}_O$and${\mathcal{Q}}_T$ is an error covariance matrix related to trajectory uncertainty.

The probability of ship conflict depends on the relative motion between ships. The distance between ships at any time, that is, the relative position, can be expressed as follow:

$$D(t)=\sqrt{A\cdot t^2+B\cdot t+C},$$

(6)

Among the top formula：

$$A={\left(V_{xO}-V_{xT}\right)}^2+{\left(V_{yO}-V_{yT}\right)}^2,$$

(7)

$$B=2\left[\left(x_O-x_T\right)\left(V_{xO}-V_{xT}\right)+\left(y_O-y_T\right)\left(V_{yO}-V_{yT}\right)\right],$$

(8)

$$C={\left(x_O-x_T\right)}^2+{\left(y_O-y_T\right)}^2,$$

(9)

The dynamic domain boundary between ships has a significant influence on the severity of the conflict. The following formula is used to calculate the dynamic domain boundary distance of the ship in all directions at a certain time.

$$MS{D}_{\mathrm{A}}\left(t\right)=\left(1.19539\cdot L_A\cdot V_A{\left(t\right)}^{0.33618}\cdot \left(1.866-0.61874\cdot \sin \left({\displaystyle \frac{2.32601-Q_{BA}\left(t\right)}{1.59485}}\cdot \pi \right)\right)\right),$$

(10)

$$MS{D}_{\mathrm{A}}\left(t\right)=\left(1.19539\cdot L_A\cdot V_A{\left(t\right)}^{0.33618}\cdot \left(1.866-0.61874\cdot \sin \left({\displaystyle \frac{2.32601-Q_{BA}\left(t\right)}{1.59485}}\cdot \pi \right)\right)\right),$$

(11)

Where $V_A\left(t\right)、V_B\left(t\right)$ are the predicted speeds of ships A and B at t time respectively, and are assumed to be zero-mean Gaussian distribution.$Q_{BA}\left(t\right)$、$Q_{AB}\left(t\right)$ Represents the predicted relative position of the target ship relative to the own ship at time t.

Under the influence of various uncertain factors, the ship’s course will change (affecting the relative orientation between ships), thereby affecting the MSD. The uncertainty component of the course prediction is taken into account and is described as follows:

$${\phi }^T\left(t\right)=\phi \left(t\right)+\epsilon \left(t\right),$$

(12)

$\epsilon \left(t\right)$represents the course prediction error component of the ship at time t.

In [40], the uncertainty distribution PDF of the ship’s course and position extracted based on AIS data drive is used as the assumption of the uncertainty component of this article. It uses Kernel Density Estimation (KDE), a non-parametric estimation method, to identify the uncertain components of the PDF of a data set via the following equation:

$$g\left(x\right)={\displaystyle \frac{1}{k^{\prime }}}{\displaystyle \sum _{i=1}^{k^{\prime }}{\varphi }_h}\left(x-x_i\right)={\displaystyle \frac{1}{k^{\prime }h}}{\displaystyle \sum _{i=1}^{k^{\prime }}{\varphi }_h}\left({\displaystyle \frac{x-x_i}{h}}\right),$$

(13)

Among them, ${\varphi }_h$is the kernel function with bandwidth h that satisfies ${\varphi }_h\left(x\right)＞0$ and ${\displaystyle {\int }_R}{\varphi }_h\left(x\right)dx=1$, $K^{\prime }$represents the number of elements in the data set within the bandwidth h.

Note that$D\left(t\right)\le MS{D}_A\left(t\right)+MS{D}_B\left(t\right)$in the inequality$D\left(t\right)$and$MSD$are functions of the position and heading prediction error components. Therefore, whether equation is established as a probabilistic event, its occurrence probability needs to be determined by a probability-based calculation method.

Considering the advantages of the Monte Carlo (MC) method as a quantitative risk assessment method in dealing with uncertainty in the ship navigation process, this paper calculates the ship conflict probability distribution based on the MC method. Generate random sampling points based on the PDF of the course and ship position, insert each random sample point into the ship’s uncertainty motion equation model, convert the random problem into a deterministic problem. When the number of simulations is large enough, the statistical information of the uncertain relative positions between ships can be calculated by using a set of deterministic solutions. According to the literature [40], 15000 iterations are needed to achieve an accuracy of 0.01, Algorithm 1 describes in detail the MC simulation process for calculating ship conflict probability, as shown in the following table:

Table 1. Calculation of ship conflict probability.

Table 1. Calculation of ship conflict probability.

Algorithm 1: MC simulation algorithm

1: Function$f\left({\stackrel{·}{V}}_O,{\stackrel{·}{V}}_T,{\stackrel{·}{\phi }}_O,{\stackrel{·}{\phi }}_T,L,lo{g}_O,la{t}_O,lo{g}_T,la{t}_T\right)$Return $L\left(t\right)$ 2: $hits\leftarrow 0$ 3: $N\leftarrow 15000$ 4: For $t=1,2,\cdot \cdot \cdot ,T$ do 5: Generate a random sample vector of position/course for each ship 6: For $i\leftarrow 1$to $N$ do 7:   ${\stackrel{·}{V}}_O\sim V_O+N\left(0,{\sigma }_O\right)$ 8:   ${\stackrel{·}{V}}_T\sim V_T+N\left(0,{\sigma }_T\right)$ 9: Call function $f$ to calculate based on random sample vector$L\left(t\right)$ 10: if $L\left(t\right)\le 0$ then 11:    $hits=hits+1$ 12:  End if 13: $PC\left(t\right)=hits/N$ 14: End for 15: $P_{con}=\max \left\{PC\left(t\right)\right\},t\in T$ Output: $P_{con}$

3.2. Ship Group Extraction

The framework for grouping regional ships is shown in Figure 3. First, the method proposed in Section 3.1 is used to detect conflicts between ships in the core port area of Ningbo-Zhoushan Port. Secondly, a ship topology network is generated based on the conflict relationship between ships. Then, coupled Fast Modularity Optimization and Spectral Clustering (FMOASC) algorithms is used to group the regional ships.

The implementation of FMOASC includes the following three steps: solving the optimal modularity, constructing the new dimension information graph of node and edge relationship, extracting the features of Laplacian matrix and clustering, as shown in Figure 4:

3.2.1. Preliminary Group Partition based on FMO Algorithm

One of the principles of regional ship grouping is to find the ship combination with strong conflict connection within the group and relatively weak connection between different groups of ships. Based on the conflict connection relationship between ships, the FMO algorithm is used to divide the regional ships into sub-communities. The preliminary division of the ship group is more concerned with the conflict connection between ships, so grouping across spatial distances may occur.

Modularity is a strength index to measure the division of a network into sub-communities, which can be used to evaluate the division of communities. The calculation method of modularity is as follows:

$$Modularty={\displaystyle \frac{1}{2m}}{\displaystyle \sum _{i,j}\left[w_{ij}-\frac{k_i{k}_j}{2m}\right]}\delta \left(c_i,c_j\right),$$

(14)

Where$w_{ij}$represents the edge weight between node$V_i$and$V_j$,$k_i$is the sum of edge weights attached to node$V_i$,$C_i$represents the community to which node$V_i$is divided, $2m$is the sum of edge weights in the network, and$\delta \left(u,v\right)$is defined as a simple incremental function as follows:

$$\delta \left(u,v\right)=\{\begin{array}{c}1,u=v\\ 0,u\ne v\end{array},$$

(15)

In the regional ship undirected weighted graph$G=\left(V,E\right)$, taking the ship as the node, $V=\left\{v_1,v_2,\cdot \cdot \cdot ,v_{\mathrm{n}}\right\}$as the node set, $v_i\in V,i=1,2,\cdot \cdot \cdot ,n$as the ship node, the ship with conflict relationship is connected, $E=\left\{e_1,e_2,\cdot \cdot \cdot ,e_{\mathrm{m}}\right\}$is the edge set, $e_{ij}=\left(v_{j1},v_{j2}\right)$,$e_i\in E,j=1,2,\cdot \cdot \cdot ,m$is the conflict relation edge, and the weight is the conflict probability value.

FMO can find the highly modular division of the network in a short time and show a complete hierarchical community structure. It is a bottom-up greedy optimization algorithm based on modularity maximization. Starting from each node as a community, the algorithm merges the neighboring communities in turn, calculates the modularity change $\Delta M$, and continuously merges the communities along the direction of the maximum increase in modularity $M$. Until the result of community division corresponding to the maximum modularity $M$is obtained. The iterative process is shown in the following figure:

Figure 5. Visual diagram of FMO algorithm iteration.

Figure 5. Visual diagram of FMO algorithm iteration.

The modularity gain $\Delta M$gained by moving orphaned nodes to the community can be calculated by the following formula:

$$\Delta M=\left[{\displaystyle \frac{{\Sigma }_{in}+k_{i,in}}{2m}}-{\left({\displaystyle \frac{{\Sigma }_{tot}+k_i}{2m}}\right)}^2\right]-\left[{\displaystyle \frac{{\Sigma }_{in}}{2m}}-{\left({\displaystyle \frac{{\Sigma }_{tot}}{2m}}\right)}^2-{\left({\displaystyle \frac{k_i}{2m}}\right)}^2\right],$$

(16)

Among them, ${\Sigma }_{in}$is the sum of the weights of the connected edges in the community C, ${\Sigma }_{tot}$is the sum of the weights of the connections associated to the node in C, $k_i$is the sum of the weights of the connections associated to the node $i$, $k_{i,in}$is the sum of the weights of the connected edges from the node $i$to the nodes in C, and $m$is the sum of the weights of all the connection relations in the network. When a node $i$is removed from its community, expressions are used to evaluate the change in modularity.

Through the preliminary division of regional ships through the FMO algorithm, sub-communities based on the conflict connection relationships between ship nodes can be obtained. It is the basis for subsequent optimal division using the SC algorithm, and the basis for constructing new dimensional information on the connection relationships between ships.

3.2.2. Optimal Group Partition based on Spectral Clustering

Based on the preliminary grouping results of FMO, combined with the goal of conflict connectivity and spatial compactness, a new dimension of information is constructed to represent the edge weight between ship nodes. The function representation of the edge weight between ships is constructed according to the probability of conflict between ships and the spatial distance. The new dimension information graph is represented as $G^{\prime }=\left(V^{\prime },E^{\prime }\right)$, where $V^{\prime }=\left\{v_1,v_2,\cdot \cdot \cdot ,v_n\right\}$is the node set, $v_x\in V^{\prime },x=1,2,\cdot \cdot \cdot ,n$is the node,$E^{\prime }=\left\{{e^{\prime }}_1,{e^{\prime }}_2,\cdot \cdot \cdot ,{e^{\prime }}_z\right\}$is the edge set, and ${e^{\prime }}_y=\left(v_{y1},v_{y2}\right),{e^{\prime }}_y\in E^{\prime },v_{y1}\ne v_{y2},y=1,2,\cdot \cdot \cdot ,z$is the edge of the graph.

In the second stage of group division, the grouping results obtained by FMO algorithm need to be fine-tuned. The more important goal of this stage is to consider the influence of ship spatial distance on the basis of conflict relations. Therefore, it is necessary to weaken the influence of conflict connections between ships and get spatially compact grouping results. Based on the influence of weakening the conflict connection relationship, the distance between ships $dis{t}^{\prime }$, the collision probability $p_{ij}$between ships and the preliminary division group label $labe{l}_{vi}$are obtained, and the similarity feature vector $D_{ev}=(dis{t}^{\prime },p_{ij},labe{l}_{vi})$is proposed, where $labe{l}_{vi}$is the group tag obtained based on FMO algorithm. The calculation of $dis{t}^{\prime }$is as follows:

$$dis{t^{\prime }}_{ij}=f\left(x\right)=\{\begin{array}{c}0,p_{ij}\ge 0.5\\ dis{t}_{ij}\ast \left(1-p_{ij}\right),0.1\le p_{ij}＜0.5\\ dis{t}_{ij}，p_{ij}＜0.1\end{array},$$

(17)

Next, standardize each feature vector, that is, $z_i=\left(t_i-\mu \right)/\sigma$.

Among them, $z_i$represents the standardized eigenvalues corresponding to each eigenvector; $t_i$and $\mu$are the eigenvalues and the average values of eigenvalues corresponding to each eigenvector, respectively, $\sigma$is the standard deviation of the corresponding eigenvalues of all eigenvectors. The Euclidean distance corresponding to the normalized feature vector between all nodes is $d_{ij}=\Vert D_{evi}-D_{evj}\Vert$, the edge weight is $w_{ij}=\exp \left(-d_{ij}/std\left(d_{ij}\right)\right)$, where$std\left(d_{ij}\right)$is the standard deviation of $d_{ij}$.

Through the feature extraction and decomposition of Laplacian Matrix L (the calculation formula is L=D-W. W is the adjacency matrix represents the adjacency relationship between nodes. D is the degree matrix, and the concept of degree indicates how many edges of a node are connected to it.) and K-means clustering, the optimal division result of ship group is obtained.

3.2.3. Group Quantity Adjustment Indicator

Using the FMOASC algorithm to divide the ship groups can obtain schemes with different numbers of groups. In order to select the ship group division scheme with excellent performance in terms of modularity and the number of noise points, an index I for comprehensive evaluation of these two important considerations is introduced. The modularity and the number of noise points of each grouping scheme are used as evaluation parameters, and the weighted normalization formula is used for calculation, as shown below:

$$I=w_m{r^{\prime }}_m-w_c{r^{\prime }}_c,$$

(18)

$${r^{\prime }}_m={\displaystyle \frac{r_m-r_m^{\min }}{r_m^{\max }-r_m^{\min }}},$$

(19)

$${r^{\prime }}_c={\displaystyle \frac{r_c^{\max }-r_c}{r_c^{\max }-r_c^{\min }}},$$

(20)

Among them, $r_m$and$r_c$are the modularity and the number of noise points corresponding to the ship grouping scheme respectively; ${r^{\prime }}_m$and${r^{\prime }}_c$are the modularity and the number of noise points after normalization respectively; $w_m$and$w_c$are the weights corresponding to the modularity and the number of noise points respectively.

By comparing the indicators under different grouping schemes, the number of ship groups is selected, and a balanced result in terms of modularity and the number of noise points is obtained.

3.3. Dynamic Calculation of Ship Navigation Risk

In order to identify the navigation risks of ships in different local waters and regions in complex navigable waters, explore the evolution law of ship navigation risks, assist maritime regulatory authorities in tracking potential high-risk ships and grasping local high-risk waters, this section carries out the following work based on the ship conflict probability detection and the extraction of ship groups from the perspective of water navigation management: First, the problem of ship navigation risk assessment is described, and the framework of global risk assessment is explained; secondly, the Shapley value method based on cooperative game is introduced to identify the contribution of each ship in a multi-ship encounter situation; finally, the navigation risks of single ships, different local waters and regions are quantified, and the navigation risks of the research waters are visualized using data visualization technology. This will help maritime regulatory authorities to have a more comprehensive understanding of real-time collision risks, effectively grasp the risks and traffic complexity of different local waters and regions, and thus have a deeper understanding of ship navigation risks, provide guidance for the formulation of risk mitigation measures.

3.3.1. Risk Assessment Framework

The ship navigation risk system is composed of the set of all ships sailing in the research waters at a certain time, as shown in Figure 6. For example: if there are n ships in the system, the system can be represented by a graph$G\left(t\right)=\left\{S,E\left(t\right)\right\}$, where S is a set of ships, represented as$\left\{S:s_{ 1},s_{ 2},s_{ 3}\cdot \cdot \cdot ,s_{ n}\right\}$, and the elements in the set are each ship, where the edge set is$E\left(t\right)\subset S\times S$, represents the conflict relationship between ships, and the edge weight is the conflict probability value between ships. The conflict probability of ship$s_i$and ship$s_j$at time t is represented by$p_{ij}\left(t\right)$. Only when there is a conflict between$s_i$and$s_j$at time t, $\left(s_i,s_j\right)\in E\left(t\right)$, that is, when and only when$p_{ij}\left(t\right)\ne 0$,$\left(s_i,s_j\right)\in E\left(t\right)$.

Figure 7 (a) (b) are the spatial position distribution and conflict connection relationship diagrams of ships at a certain moment. Combining Figure 7 (a) (b), we can see that at a certain moment, there is no connection between ships that are far away or even close but have a tendency to move away from each other. Therefore, the ship system graph is not connected, but is divided into several subgraphs. Assume$G_i\left(t\right)=\left\{S_i,E_i\left(t\right)\right\}\left(i=1,2,3,\cdot \cdot \cdot ,m\right)$ is a connected subgraph of$G\left(t\right)$, a subgraph (that is, a ship group divided according to the FMOASC algorithm) is a subset (subsystem) of the ship system, and m is the number of subsystems. The relationship between the ship navigation risk system and subsystems in the study waters is$G\left(t\right)=G_1\left(t\right)\cup G_2\left(t\right)\cup \cdot \cdot \cdot \cup G_m\left(t\right)$.For any two subgraphs $G_i\left(t\right)$and$G_j\left(t\right)$, there exists$S_i\cap S_j=\varnothing$.

The collision risk of the kth subsystem$G_i\left(t\right)$ is represented by$p_k\left(t\right)$ . For any pair of ships with a non-zero collision probability, the risk of the subsystem in which they are located will be affected by the collision probability between them, which in turn affects the system risk of the regional ship composition. This paper will model$p_k\left(t\right)$as the function of$p_{ij}\left(t\right)\ne 0$, that is$p_k\left(t\right)=f_k\left(p_{ij}\left(t\right),\cdot \cdot \cdot ,p_{ul}\left(t\right)\right)$, where$\left(s_i,s_j\right)\in E_k\left(t\right)$and$p_{ij}\left(t\right)\ne 0$.

3.3.2. Global risk Assessment

The process of evaluating global collision risk is shown in the figure below:

Figure 8. Global collision risk assessment process.

Figure 8. Global collision risk assessment process.

It is assumed that the global collision risk of a multi-ship encounter is composed of the collision risk of each ship in the situation. However, due to the different encounter types and collision avoidance maneuvers between different ships, the contribution of each ship to the global collision risk will also be different. Based on the above assumptions, the Global Collision Risk (GCR) of a multi-ship encounter can be expressed by the following equation:

$$GCR={\displaystyle \sum _{i=1}^{n}C{R}_i\times C_i},$$

(21)

Among them, $C{R}_i$and$C_i$refer to the collision risk and contribution of ship $i$ in the multi-ship encounter scenario, respectively.

In order to evaluate the global collision risk in multi-ship encounter, it is necessary to first determine the collision risk of each ship in the encounter and its contribution. For complex navigable waters with high ship density, a large number of ships will be in the same waters at the same time. The calculation of the global collision risk needs to traverse all ships in the waters at the same time. Therefore, the calculation complexity is very large. In order to simplify the calculation, this paper adopts the ship group division method from the perspective of water navigation management proposed in Section 3.2, divides the ships in the study waters into different ship groups, and then identifies the collision risk of each ship based on the ship conflict probability detection method proposed in Section 3.1. Subsequently, the Shapley value method in the cooperative game (to be introduced in the next subsection) is used to determine the contribution of each ship. The collision risk within the ship group is calculated according to formula (21). On this basis, all collision risks of different ship groups are integrated through formula (22) to obtain the Regional Collision Risk (RCR).

$$RCR={\displaystyle \sum _{j=1}^{m}GC{R}_j\times C_j},$$

(22)

Among them, $GC{R}_j$is the Collision Risk of each Group, $C_j$is the contribution of each ship group.

3.3.3. Local/Regional Navigation Risk Quantification

The scope of ship navigation risk events discussed is aimed at two or more ships, that is, at least two ships are the system objects of ship navigation risk events. The navigation status of the ship can be obtained based on AIS data, and the probability of conflict between ships can be obtained using the ship conflict probability detection method proposed in Section 3.1. This degree of conflict is a measure of the severity of the conflict from the perspective of each ship. The forms of local and regional system risks are characterized by$p_k\left(t\right)=f_k\left(p_{ij}\left(t\right),\cdot \cdot \cdot ,p_{ul}\left(t\right)\right)$and$F\left(t\right)=f\left(p_k\left(t\right),\cdot \cdot \cdot ,p_m\left(t\right)\right)$ respectively.

3.3.4. Contribution Determination

The Shapley value method was proposed by Shapley L.S in 1953 to solve the problem of conflicts caused by the distribution of benefits among multiple participants in the process of cooperation. It is a solution method in the field of cooperative games. The Shapley value method represents the contribution of each participant to the entire group. A major advantage of applying the Shapley value is that the benefits are distributed according to the marginal contribution rate of the members to the combination. This solves the problem of benefit distribution because the method distributes the cooperation amount by estimating the contribution of each participant. This method is used to calculate the contribution of each ship to the navigation risk in multi-ships encounter. The general form of the Shapley value is as follows:

$$S_i\left[A\right]={\displaystyle \sum _{T\subseteq N,i\in T}\frac{\left(t-1\right)!\left(n-t\right)!}{n!}}\left[A\left(T\right)-A\left(T-\left\{i\right\}\right)\right],$$

(23)

Among them, $i$refers to the participants in the game, $T$refers to the combination consisting of participants$i$, $t$refers to the number of participants in the combination$T$, $n$refers to the number of participants in the entire group, $A\left(T\right)$refers to the amount generated by the combination$T$,$A\left(T-\left\{i\right\}\right)$refers to the amount generated by the combination before the participants$i$ join$T$, and $S_i\left[A\right]$refers to the Shapley value corresponding to the participants$i$.

When applying the Shapley value method, the multi-ships encounter is regarded as a cooperative game, each ship in the encounter scenario is regarded as a game subject, and its collision risk is regarded as the game amount. Ships can form different combinations through permutations and combinations. According to formula (23), when calculating the Shapley value of each ship, it is first necessary to find the risk $A\left(T\right)$of each combination. The number of combinations is the sum of all the permutations of the encountering ships, expressed as follows:

$$A\left\{1,2,\cdot \cdot \cdot n\right\}={\displaystyle \sum _{i=1}^{n}C{R}_{i|1,2,\cdot \cdot \cdot ,n}},$$

(24)

Among them, $C{R}_{i|1,2,\cdot \cdot \cdot ,n}$refers to the collision risk of the ship$i$ in n-ships encounter.

Each ship’s contribution is calculated as follows:

• Find out all ship combinations by arranging them.
• For each combination, the combined risk value is calculated by summing the risk of each ship according to Formula (24).
• Calculate the Shapley value of each ship according to equation (23).
• Normalized the Shapley value to obtain the contribution degree of each ship to the global collision risk.

In addition to vessels with a high risk of collision, VTS supervisors also need to focus on vessels with a high Shapley value, which is essential for the control and management of global collision risk. Because the reduction of collision risk of ships with high Shapley value is conducive to the mitigation of global collision risk.

4. Applications and Case Study Results

4.4. Quantifying Risks in Different Local Waters

Using 7,680,496 AIS data from the study area for one month, conflict data was calculated. The conflict distribution visualization diagram for each hour is shown in Figure 16. The numbers under each sub-figure in the figure represent 24 time periods, such as “1” represents the conflict distribution during the 00:00-01:00 time period. As shown in sub-figures 16 (1)-(4), the number of ship conflicts during the 00:00-04:00 time period is relatively small; in Figures 16 (5)-(15), the conflict distribution area is wider and is roughly distributed in the Fudu Waterway, Luotou Waterway, and Xiazhimen Channel. Among them, the frequency of ship conflicts in Luotou Waterway and Xiazhimen Channel is higher. The reason for the higher frequency of conflicts in Luotou Waterway is that ship traffic from other directions converges in this water area, increasing the frequency of multiple ship encounters. In addition, ship turning due to terrain geometric constraints in this water area may also be one of the reasons for this result. In fact, Luotou Waterway is an official warning area issued by the Ningbo-Zhoushan VTS Center, which to some extent verifies the effectiveness of the conflict probability detection method proposed in this paper in identifying high collision risk areas. Xiazhimen Channel is the main channel for large ships to enter and leave the port, but the navigation width is narrow, so it becomes one of the areas with high collision risk. One possible reason for the high frequency of ship conflicts in Fudu Waterway is that the traffic width in this water area is narrow, resulting in a reduction in the minimum passing distance between encountering ships, so, many high-severity conflicts occur. The degree and distribution range of ship conflicts are expanding in the time period of 20:00-24:00, which may be due to the increase in the number of ship conflicts caused by reduced visibility at night.

The distribution of ship conflicts during the day and at night is explored, that is, the ship conflicts in the time periods of 06:00-18:00 and 18:00-06:00 are counted, and two conflict distribution images are generated every 10 days, as shown in Figure 17, where “Daytime” in the sub-image represents the conflict distribution in the time period of 06:00-18:00, and “Night” represents the conflict distribution image in the time period of 18:00-06:00. As can be seen from the figure, the degree of conflict is more severe during the day than at night, the number of conflicts is larger, and the high-risk areas are concentrated in the waters of the intersection of Luotou Waterway and Zhitouyang, Fodu Waterway, Xiazhimen Channel, and the southern part of Damao Island. Among them, the waters south of Damao Island are the main roads connecting the port hub with offshore ship traffic. It is the area with the highest ship traffic density and the highest frequency of ship encounters in the port; the intersection of Luotou Waterway and Zhitouyang waters, and the intersection of Xiazhimen Channel and Zhitouyang waters are both navigation warning areas officially established by the maritime department, which verifies the accuracy and effectiveness of this article’s conflict research at the micro-time and space scale.

Figure 18 shows the characteristics of the distribution of the hourly conflict probability values in the study waters. It can be observed that the median conflict probability fluctuates slightly over time. In the 11:00-15:00 time period, the median of conflicts between ships is larger than that in other time periods. The possible reason is that the ship traffic volume in this time period is large and the possibility of dangerous encounters between ships is high. The maximum conflict probability is basically maintained at 1, which means that ship encounters with high conflict severity will occur in every time period. VTS centers can use indicators such as median and maximum values to more comprehensively understand the real-time collision risks and traffic complexity in the region, thereby improving their ability to deal with dangerous encounters caused by high traffic intensity.

Figure 18. Time distribution of conflict.

Figure 18. Time distribution of conflict.

Figure 19. Time distribution of the number of conflicts.

Figure 19. Time distribution of the number of conflicts.

Figure 19 shows the number of conflicts per hour in the study waters. It can be seen that ship conflicts occur more frequently in the time periods of 05:00, 07:00-11:00, and 13:00-17:00, which is consistent with the actual traffic conditions in the study waters, because the ship traffic density is higher during the day. Based on the identification results of the conflict distribution in time and space, marine navigators and maritime departments can have a deeper understanding of when and where to enhance situational awareness during the navigation process of ships, providing valuable reference for traffic regulatory departments to formulate and implement appropriate regulatory strategies.

The distribution of conflict severity and the corresponding cumulative probability distribution are shown in Figure 20. It can be seen from the figure that the severity of most conflicts is below 0.2, which indicates that the probability of a ship being involved in a higher conflict is much lower than the probability of being involved in a minor conflict. In order to further demonstrate the characteristics of the detected conflicts, Figure 21 shows the number of different conflict types among different ship types. It can be seen that the number of crossing encounter conflict types is the largest, followed by overtaking, and then head-on. Compared with overtaking and head-on conflicts, the proportion of conflicts in crossing encounter scenarios is higher, mainly because the relatively large angle range is classified as crossing encounters. In addition, the number of conflicts between cargo ships in the three encounter scenarios is much higher than that of other ships, which may be due to the large number of cargo ships in the study waters and the frequent encounters.

In order to verify the effectiveness and feasibility of the global risk assessment method proposed in this paper, the complexity of maritime traffic is introduced as a collision risk indicator for evaluation. The complexity of maritime traffic is used to indicate the degree of traffic conditions, including the degree of congestion and risk. Although the traffic complexity model is not a traditional risk assessment model, it can reflect the instantaneous collision risk to a certain extent. Generally speaking, the higher the traffic complexity, the higher the collision risk, because when the traffic situation becomes complicated, the difficulty of avoiding collision increases. In addition, compared with simple indicators such as traffic density and ship distance, it can more fully reflect the risk of ship collision. Therefore, this paper uses traffic complexity as an indicator of the effectiveness of the global risk assessment framework.

The definition of the maritime traffic complexity model includes two parts: traffic density factor and traffic conflict factor. The traffic density factor is described by the relative distance between ships, and the traffic conflict factor is described by the length, speed and relative orientation of the ship. The traffic complexity is expressed as follows:

$$Com\left(D_{ij},V_{ij},{\theta }_{ij}\right)=De{n}_{ij}\left(D_{ij}\right)+Pcon{f}_{ij}\left(D_{ij},V_{ij},{\theta }_{ij}\right),$$

(26)

$$Com\left(i\right)={\displaystyle \frac{1}{N}}{\displaystyle \sum _{i\ne j,j=1}^{N}Co{m}_{ij}\left(D_{ij},V_{ij},{\theta }_{ij}\right)},$$

(27)

Among them, $Co{m}_{ij}$、$De{n}_{ij}$、$Pcon{f}_{ij}$are the traffic complexity, traffic density complexity, and traffic conflict complexity of the ship $i$and ship$j$ , respectively. $D_{ij},V_{ij},{\theta }_{ij}$ are the relative distance, relative speed, and relative orientation between the ship $i$and ship$j$, respectively; $Com\left(i\right)$ is the traffic complexity of the ship$i$, is the average of the sum of all complexities involving ship$i$.

According to the ship group extraction algorithm proposed in this paper, the ships in the study waters at a certain time are divided into 10 ship groups, as shown in Figure 22. The contribution value of each ship in each ship group obtained by the Shapley value method in cooperative game is shown in Figure 23. The collision risk and traffic complexity of each ship group are shown in Figure 24.

It can be seen from Figure 24 that the group with the highest collision risk is Group 1, which reaches 0.5058. It can be seen that: the ship traffic density in Group 1 is large, there are many ship conflicts, and the narrow channel width leads to a reduced passing distance between ships, resulting in a high risk within the group; Group 7 has the lowest collision risk of 0.09, and the distance between ships in the group is relatively large, resulting in a conflict of lower severity; and the collision risk in Group 3 is significantly greater than that in Group 7. This is because the distance between ships in Group 3 is much smaller than that in Group 7. When approaching each other, the distance is negatively correlated with the collision risk. Therefore, the collision risk calculated according to this framework is consistent with the actual traffic conditions in the study waters. In addition, it can be seen from Figure 24 that the traffic complexity curve and collision risk curve of each group have the same change trend. In order to quantitatively analyze the correlation between the change trend of collision risk and traffic complexity of the 10 ship groups, a Pearson correlation analysis was performed. The analysis results show that at a confidence level of 99%, the change trend of collision risk of the 10 ship groups is strongly correlated with the change trend of traffic complexity, with a correlation coefficient of 0.923 and a P value of less than 0.001. Therefore, it can be proved that the collision risk of different ship groups in the waters calculated according to the proposed framework is effective.

5. Conclusions

References

1. Goerlandt, F.; Montewka, J. Maritime Transportation Risk Analysis: Review and Analysis in Light of Some Foundational Issues. Reliab Eng Syst Saf 2015, 138, 115–134. [CrossRef]
2. Liu, Z.; Li, Y.; Zhang, Z.; Yu, W.; Han, B.; Du, Y. Spatial Modeling and Analysis Approach for Ship-Encounters Dynamic Spatial–Temporal Domain. Ocean Engineering 2023, 278, 114373. [CrossRef]
3. Yim, J.-B.; Kim, D.-S.; Park, D.-J. Modeling Perceived Collision Risk in Vessel Encounter Situations. Ocean Engineering 2018, 166, 64–75. [CrossRef]
4. Chen, P.; Huang, Y.; Mou, J.; van Gelder, P.H.A.J.M. Probabilistic Risk Analysis for Ship-Ship Collision: State-of-the-Art. Saf Sci 2019, 117, 108–122. [CrossRef]
5. Zhao, Y.; Suo, Y.; Xian, B. Ship Encounter Scenario and Maneuvering Behavior Mining Based on AIS Data. In; 2022; pp. 58–65.
6. Liu, J.; Zhang, J.; Yang, Z.; Wan, C.; Zhang, M. A Novel Data-Driven Method of Ship Collision Risk Evolution Evaluation during Real Encounter Situations. Reliab Eng Syst Saf 2024, 249, 110228. [CrossRef]
7. Lei, P.-R. Mining Maritime Traffic Conflict Trajectories from a Massive AIS Data. Knowl Inf Syst 2020, 62, 259–285. [CrossRef]
8. Chen, P.; Huang, Y.; Mou, J.; van Gelder, P.H.A.J.M. Ship Collision Candidate Detection Method: A Velocity Obstacle Approach. Ocean Engineering 2018, 170, 186–198. [CrossRef]
9. Zhang, W.; Goerlandt, F.; Kujala, P.; Wang, Y. An Advanced Method for Detecting Possible near Miss Ship Collisions from AIS Data. Ocean Engineering 2016, 124, 141–156. [CrossRef]
10. Kundakçı, B.; Nas, S.; Gucma, L. Prediction of Ship Domain on Coastal Waters by Using AIS Data. Ocean Engineering 2023, 273, 113921. [CrossRef]
11. Fujii, Y.; Tanaka, K. Traffic Capacity. Journal of Navigation 1971, 24, 543–552. [CrossRef]
12. Wang, Y.; Chin, H.-C. An Empirically-Calibrated Ship Domain as a Safety Criterion for Navigation in Confined Waters. Journal of Navigation 2016, 69, 257–276. [CrossRef]
13. Wang, N. A Novel Analytical Framework for Dynamic Quaternion Ship Domains. Journal of Navigation 2013, 66, 265–281. [CrossRef]
14. Zhang, L.; Meng, Q. Probabilistic Ship Domain with Applications to Ship Collision Risk Assessment. Ocean Engineering 2019, 186, 106130. [CrossRef]
15. Szłapczyński, R.; Niksa-Rynkiewicz, T. A Framework of A Ship Domain-Based Near-Miss Detection Method Using Mamdani Neuro-Fuzzy Classification. Polish Maritime Research 2018, 25, 14–21. [CrossRef]
16. Dinh, G.H.; Im, N. The Combination of Analytical and Statistical Method to Define Polygonal Ship Domain and Reflect Human Experiences in Estimating Dangerous Area. International Journal of e-Navigation and Maritime Economy 2016, 4, 97–108. [CrossRef]
17. Zhen, R.; Riveiro, M.; Jin, Y. A Novel Analytic Framework of Real-Time Multi-Vessel Collision Risk Assessment for Maritime Traffic Surveillance. Ocean Engineering 2017, 145, 492–501. [CrossRef]
18. Chin, H.C.; Debnath, A.K. Modeling Perceived Collision Risk in Port Water Navigation. Saf Sci 2009, 47, 1410–1416. [CrossRef]
19. Wang, N.; Yuen, K.F.; Yuan, J.; Li, D. Ship Collision Risk Assessment: A Multi-Criteria Decision-Making Framework Based on Dempster–Shafer Evidence Theory. Appl Soft Comput 2024, 162, 111823. [CrossRef]
20. Zhang, M.; Montewka, J.; Manderbacka, T.; Kujala, P.; Hirdaris, S. A Big Data Analytics Method for the Evaluation of Ship - Ship Collision Risk Reflecting Hydrometeorological Conditions. Reliab Eng Syst Saf 2021, 213, 107674. [CrossRef]
21. Gil, M.; Montewka, J.; Krata, P.; Hinz, T.; Hirdaris, S. Determination of the Dynamic Critical Maneuvering Area in an Encounter between Two Vessels: Operation with Negligible Environmental Disruption. Ocean Engineering 2020, 213, 107709. [CrossRef]
22. Gang, L.; Wang, Y.; Sun, Y.; Zhou, L.; Zhang, M. Estimation of Vessel Collision Risk Index Based on Support Vector Machine. Advances in Mechanical Engineering 2016, 8, 168781401667125. [CrossRef]
23. Li, B.; Pang, F.-W. An Approach of Vessel Collision Risk Assessment Based on the D–S Evidence Theory. Ocean Engineering 2013, 74, 16–21. [CrossRef]
24. Lenart, A.S. Analysis of Collision Threat Parameters and Criteria. Journal of Navigation 2015, 68, 887–896. [CrossRef]
25. Larson, J.; Bruch, M.; Ebken, J. Autonomous Navigation and Obstacle Avoidance for Unmanned Surface Vehicles.; Gerhart, G.R., Shoemaker, C.M., Gage, D.W., Eds.; May 5 2006; p. 623007.
26. Huang, Y.; van Gelder, P.H.A.J.M.; Wen, Y. Velocity Obstacle Algorithms for Collision Prevention at Sea. Ocean Engineering 2018, 151, 308–321. [CrossRef]
27. Szlapczynski, R.; Szlapczynska, J. A Target Information Display for Visualizing Collision Avoidance Maneuvers in Various Visibility Conditions. Journal of Navigation 2015, 68, 1041–1055. [CrossRef]
28. Huang, Y.; Chen, L.; van Gelder, P.H.A.J.M. Generalized Velocity Obstacle Algorithm for Preventing Ship Collisions at Sea. Ocean Engineering 2019, 173, 142–156. [CrossRef]
29. Huang, Y.; van Gelder, P.H.A.J.M. Time-Varying Risk Measurement for Ship Collision Prevention. Risk Analysis 2020, 40, 24–42. [CrossRef]
30. Idiri, B.; Napoli, A. The Automatic Identification System of Maritime Accident Risk Using Rule-Based Reasoning. In Proceedings of the 2012 7th International Conference on System of Systems Engineering (SoSE); IEEE, July 2012; pp. 125–130.
31. Hänninen, M. Bayesian Networks for Maritime Traffic Accident Prevention: Benefits and Challenges. Accid Anal Prev 2014, 73, 305–312. [CrossRef]
32. Kujala, P.; Hänninen, M.; Arola, T.; Ylitalo, J. Analysis of the Marine Traffic Safety in the Gulf of Finland. Reliab Eng Syst Saf 2009, 94, 1349–1357. [CrossRef]
33. Bye, R.J.; Aalberg, A.L. Maritime Navigation Accidents and Risk Indicators: An Exploratory Statistical Analysis Using AIS Data and Accident Reports. Reliab Eng Syst Saf 2018, 176, 174–186. [CrossRef]
34. Du, L.; Goerlandt, F.; Kujala, P. Review and Analysis of Methods for Assessing Maritime Waterway Risk Based on Non-Accident Critical Events Detected from AIS Data. Reliab Eng Syst Saf 2020, 200, 106933. [CrossRef]
35. Terndrup Pedersen, P. Collision and Grounding Mechanics;
36. Martins, M.R.; Maturana, M.C. Application of Bayesian Belief Networks to the Human Reliability Analysis of an Oil Tanker Operation Focusing on Collision Accidents. Reliab Eng Syst Saf 2013, 110, 89–109. [CrossRef]
37. Fan, S.; Blanco-Davis, E.; Yang, Z.; Zhang, J.; Yan, X. Incorporation of Human Factors into Maritime Accident Analysis Using a Data-Driven Bayesian Network. Reliab Eng Syst Saf 2020, 203, 107070. [CrossRef]
38. Weng, J.; Li, G.; Chai, T.; Yang, D. Evaluation of Two-Ship Collision Severity Using Ordered Probit Approaches. Journal of Navigation 2018, 71, 822–836. [CrossRef]
39. Chen Y.; Liu K.; Xin X.; Miao P.; Zhang M. Ship collision probability detection method based on AIS data. Chinese Nav. 2022,45(04):15-22.
40. Xin, X.; Liu, K.; Yang, Z.; Zhang, J.; Wu, X. A Probabilistic Risk Approach for the Collision Detection of Multi-Ships under Spatiotemporal Movement Uncertainty. Reliab Eng Syst Saf 2021, 215, 107772. [CrossRef]