*3.6. Occupancy Grid Map Construction*

In order to apply the environmental information provided by the electronic chart for path planning, it is necessary to process the chart into a binary occupancy grid map as shown in Figure 2. In this process, we first set an appropriate binarization threshold that converts an RGB image into a binary image. Such a threshold [38] is vital to constructing an accurate grid map that ensures the feasibility of path planning. If the threshold is inappropriate, as shown in Figure 2c, then a shoal is identified as a passable area, greatly increasing the navigation risk. Second, the grid size determines the resolution of path planning. Too large a grid cannot capture the subtle details of environments, e.g., small obstacles, resulting in unsafe path searching. On the other hand, too small a grid greatly increases the search space, reducing the computation speed.

**Figure 2.** Converting an electronic chart into an occupancy grid map, where the black areas in the grid map represent obstacles and the white areas represent passable areas.

Considering the maneuverability of a ship, this paper chooses the minimum turning radius as the criterion to measure the grid size. The ship's minimum turning radius can be measured through the ship's maneuverability experiment. The turning radius depends on the sailing speed and water flow velocity. According to [39], in our simulation, we set the minimum turning radius of the ship as three times the length of the ship. Finally, we map the geographic chart represented in terms of longitude and latitude to the grid map using the following equation:

$$\begin{cases} l\alpha\_{(i,j)} = tl\,L\alpha + \frac{|br\,L\alpha - tl\,L\alpha|}{\frac{\upsilon}{\upsilon}} \cdot (i - 1) \\\ l\alpha\_{(i,j)} = br\,L\alpha + \frac{|tl\,L\alpha - br\,L\alpha|}{h} \cdot (j - 1) \end{cases} \tag{20}$$

where *lon*(*i*,*j*) and *lat*(*i*,*j*) denote the longitude and latitude of position (*i*, *j*). *tlLon* and *tlLat* denote the longitude and latitude at the top-left corner of the selected area. *brLon* and

*brLat* denote the longitude and latitude at the bottom-right corner of the selected area. *w* and *h* denote the width and height of the electronic chart.

In our simulated forward exploration, we use eight discretized directions in the grid map to search for paths as shown in Figure 3.

**Figure 3.** 8 discretized directions in path searching for autonomous ships.

#### *3.7. Quantification of Turning Cost*

The A\* algorithm uses Equation (18) as the evaluation function to obtain a path with the shortest distance.

$$f(N\_i) = \lg(N\_i) + h(N\_i) \tag{21}$$

where *f*(*Ni*) denotes the estimated cost from the starting point to the target point, *g*(*Ni*) the actual cost from the initial node to node *Ni*, and *h*(*Ni*) the estimated cost of the best path from node *Ni* to the target node.

However, in path searching for ships, the turning is much more difficult than for cars or aerial vehicles. Thus, we have to add the cost of turning in order to better evaluate the path. Here we use the diagonal distance to compute the turning cost as shown in Figure 4.

The yaw *φ<sup>i</sup>* between waypoint *pi* and *pi*+<sup>1</sup> can be computed as follows.

$$\phi\_i = \arctan(\left|\frac{p\_{i+1}(y) - p\_i(y)}{p\_{i+1}(x) - p\_i(x)}\right|)\tag{22}$$

where *pi*(*x*), *pi*(*y*) denote the coordinates of waypoint *pi*, *pi*+1(*x*), *pi*+1(*y*) the coordinates of waypoint *pi*+1, *pi*−1(*x*), *pi*−1(*y*) the coordinates of waypoint *pi*−1. Preventing collisions is still of the highest priority in path planning. Thus, it is not reasonable to simply pursue the minimum turning cost in planning. We add a penalty to the turning cost:

$$\mathcal{L}(N\_i) = \varepsilon \cdot \max(0, \Delta \phi\_i - \phi) \tag{23}$$

where *φ* is the penalty threshold of the yaw and *ε* the penalty coefficient. Empirically, we set *φ* = 30◦ and *ε* = 0.8. Δ*φ<sup>i</sup>* is computed as follows.

$$\Delta \phi\_i = \left| \arctan(\frac{p\_{i+1}(y) - p\_i(y)}{p\_{i+1}(x) - p\_i(x)}) - \arctan(\frac{p\_i(y) - p\_{i-1}(y)}{p\_i(x) - p\_{i-1}(x)}) \right| \tag{24}$$

Finally, the new evaluation function of our path searching algorithm is defined as:

$$f(N\_i) = g(N\_i) + h(N\_i) + c(N\_i) \tag{25}$$

**Figure 4.** An illustration of how to compute the turning cost.

The improved A\* algorithm is shown in Algorithm 1.


#### *3.8. Global Trajectory Optimization*

The previous path searching gives us a discrete path with the minimum cost. However, the path does not consider the dynamic feasibility with respect to time, velocity, and acceleration. This part requires a further step to optimize the searched path into a smooth trajectory. The trajectory can be defined as a -order polynomial.

$$p(t) = p\_0 + p\_1 t + p\_2 t^2 + \dots + p\_n t^n = \sum\_{i=0}^n p\_i t^i \tag{26}$$

where *p*0, *p*1, ... , *pn* are the coefficients of this trajectory. We denote *P* = [*p*0, *p*1,..., *pn*] *T*, and then Equation (23) can be rewritten as

$$p = \left[1, t, t^2, \dots, t^n\right] \cdot P \tag{27}$$

Then, Equation (14) can be expressed as

$$\int\_0^T \left(p^{(4)}(t)\right)^2 dt$$

$$= \sum\_{i=1}^k \int\_{t\_{i-1}}^{t\_i} \left(p^{(4)}(t)\right)^2 dt$$

$$= \sum\_{i=1}^k \int\_{t\_{i-1}}^{t\_i} \left(\left[0,0,0,0,24,\ldots,\frac{n!}{(n-4)!}t^{n-4}\right] \cdot p\right)^T \left[0,0,0,0,24,\ldots,\frac{n!}{(n-4)!}t^{n-4}\right] \cdot p dt$$

$$= \sum\_{i=1}^k p^T \int\_{t\_{i-1}}^{t\_i} \left[0,0,0,0,24,\ldots,\frac{n!}{(n-4)!}t^{n-4}\right]^T \left[0,0,0,0,24,\ldots,\frac{n!}{(n-4)!}t^{n-4}\right] dt \cdot p \quad \text{(28)}$$
Let

$$Q\_i = \int\_{t\_{i-1}}^{t\_i} \left[ 0, 0, 0, 0, 24, \dots, \frac{n!}{(n-4)!} t^{n-4} \right]^T \left[ 0, 0, 0, 0, 24, \dots, \frac{n!}{(n-4)!} t^{n-4} \right] dt \tag{29}$$

We have

$$\int\_{0}^{T} \left( p^{(4)}(t) \right)^{2} dt = \sum\_{i=1}^{k} p^{T} \mathbb{Q}\_{i} p \tag{30}$$

However, the polynomial expression cannot explicitly control the shape of the trajectory. To gain better control, we choose the Bezier curve using Bernstein polynomials. The *k*-th segment of the trajectory can be expressed as

$$B\_k(t) = \sum\_{i=0}^n c\_k^i b\_n^k(t) \tag{31}$$

where *b<sup>k</sup> <sup>n</sup>*(*t*) = ' *n k* ( ·*t i* ·(1 − *t*) *n*−*i* , <sup>t</sup> <sup>∈</sup> [0,1] *<sup>c</sup><sup>i</sup> <sup>k</sup>* denotes the control point of the *k*-th segment of the Bezier curve.

Since the trajectory must pass through the first and last control points, it can satisfy the positional constraints of the initial and final states. In addition, based on the hodograph of the Bezier curve, we impose constraints on the velocity and acceleration of the trajectory, ensuring the multi-order continuity of the trajectory.

#### *3.9. Real-Time Obstacle Avoidance*

In a static chart, a ship can navigate safely along the aforementioned global trajectory. However, the marine environment is complex and changeable. Ships need to deal with unexpected risks when sailing, e.g., avoiding islands and reefs. As shown in Figure 5, if the ship maintains the planned global trajectory, it will collide with a temporary obstacle. To address this issue, we perform local path planning based on the B-spline curve. The advantage of the B-spline trajectory is that it can change the curve locally by adjusting few control points, while any control point of a Bezier curve will change the shape of the whole trajectory. Moreover, it guarantees that the locally replanned trajectory still satisfies the ship's kinematic and dynamic constraints. This not only achieves the goal of real-time obstacle avoidance, but also satisfies all optimization objectives.

**Figure 5.** A ship can hit an expected obstacle by sailing along the generated global trajectory.

A B-spline can be expressed as

$$\mathcal{C}\_{\mathcal{P}}(\boldsymbol{\mu}) = \sum\_{i=0}^{n} \mathcal{N}\_{i,\mathcal{P}}(\boldsymbol{\mu}) P\_i \tag{32}$$

where *Pi* denotes the *i*-th control point and *Ni*,*p*(*u*) is the B-spline basis function of degree *p*. *u* = [*u*0, *u*1,..., *u*m] is the knot vector. Typically, a three-degree B-spline can ensure the smoothness of accelerations. Thus, we have

$$P\_{(0,3)} = \frac{1}{6} \begin{bmatrix} 1 & t & t^2 & t^3 \end{bmatrix} \begin{bmatrix} 1 & 4 & 1 & 0 \\ -3 & 0 & 3 & 0 \\ 3 & -6 & 3 & 0 \\ -1 & 3 & -3 & 1 \end{bmatrix} \begin{bmatrix} P\_0 \\ P\_1 \\ P\_2 \\ P\_3 \end{bmatrix} \tag{33}$$

Figure 6 illustrates the collision avoidance algorithm. In this figure, *P*<sup>0</sup> = [*x*0, *y*0] and *P*<sup>3</sup> = [*x*3, *y*3] are the start and the end of the local planning. The gray area *ABCD* represents an obstacle. *Ψ<sup>i</sup>* denotes the yaw at position *Pi*. *di* denotes the distance between *Pi*−<sup>1</sup> and *Pi*. The geometrical relationship among these positions can be expressed as

$$\begin{cases} \mathbf{x}\_1 = \mathbf{x}\_0 + d\_1 \cos \Psi\_0\\ y\_1 = y\_0 + d\_1 \sin \Psi\_0 \end{cases} \tag{34}$$

$$\begin{cases} \mathbf{x}\_2 = \mathbf{x}\_3 - d\_3 \cos \Psi\_2\\ \mathbf{y}\_2 = \mathbf{y}\_3 - d\_3 \sin \Psi\_2 \end{cases} \tag{35}$$

First, based on random sampling [40] and collision detection [41], we obtain *d*1, *d*<sup>3</sup> and *Ψ*2(*Ψ*<sup>0</sup> = 0), and then use the geometric relations in Equations (34) and (35) to solve for the position of *P*<sup>1</sup> = [*x*1, *y*1] and *P*<sup>2</sup> = [*x*2, *y*2]. At last, the locally replanned B-spline can be generated by the control points *P*0, *P*1, *P*2, and *P*3.

**Figure 6.** An illustration of the collision avoidance algorithm.

#### **4. Simulation and Results**

We conduct the evaluation of ESP by implementing a data-driven simulator in MAT-LAB. All simulation experiments are run on a quad-core 2.40 GHz Intel i5-1135G7 processor and 16 GB RAM. We input the data from the database of our developed GIS. The data include image and vector maps, longitude and latitude coordinates, and ship route information. During the simulation, we use environmental data such as shoals and whirlpools that are not currently marked in the GIS database to evaluate the effectiveness of ESP. The simulated settings are listed in Table 1. Specifically, the ship length is 30 m, the maximum turning radius is 3 times the ship length. The maximum velocity is 7.7 m/s. The maximum acceleration and jerk are 1 m/s2 and 10 m/s3 [42], respectively.

**Table 1.** Simulation parameters.


We simulate ESP in two scenes using the GCJ-02 coordinate system. In both scenes, the velocities and accelerations at the start and the end are 0. The simulated occupancy grid map is of size 50 × 50. We compare the path searching results of ESP with the A\* algorithm as shown in Figure 7. The results show that ESP effectively reduces the number of turning points and the planned path is safe.

**Figure 7.** In (**a**), the black line represents the searched path by A\*. The green dot denotes the start and the red star denotes the goal. In (**b**), the red line represents the searched path by ESP, considering the turning cost. (**c**) shows the searched path by ESP in the chart.

Figure 8 shows the performance of RRT. Due to the randomness of RRT, the planned paths have many unnecessary turning points, which is detrimental to the safe navigation of ships.

**Figure 8.** The results of the RRT algorithm. (**a**–**c**) are the three stochastic planning results of the RRT algorithm. The solid black line is the planned route, the solid red line records the sampling process of the algorithm, and the solid blue line shows the planning results.

Figure 9 shows the result of RRT\*. RRT\* needs to rewire parents to find asymptotically optimal paths. The result will be close to the optimal solution with more iterations.

Figure 10 shows the optimized trajectory of our proposed ESP. It can be seen that the optimized trajectory (the green curve in Figure 10a) meets the requirements of safety, feasibility, and smoothness.

Table 2 shows the numeric comparison in Scene 1 among ESP, A\* [12], RRT [43] and RRT\* [44] in terms of computation time, number of turns, and fuel consumption. It can be seen that due to the need to measure the turning cost of the path, the searching time of ESP is 0.105 s longer than that of the A\* algorithm, and the algorithm's operating efficiency is reduced by 36.71%. Nevertheless, it is still 1.771 s shorter than that of the RRT algorithm, and 6.021 s shorter than that of the RRT\* algorithm. The efficiency of the algorithm is improved by 4.35 times and 15.42 times, respectively. In addition, the number of turns of ESP is obviously less than that of the A\* algorithm, which reduces the number of largeangle steers to 8 and improves the safety of ship navigation. The fuel consumption of ESP is 164.6008 kg less than that of the A\* algorithm, 387.1543 kg less than that of the RRT algorithm, and 24.3311 kg less than that of the RRT\* algorithm.

**Figure 9.** Three different searched paths from the RRT\* algorithm. (**a**–**c**) are the three stochastic planning results of the RRT\* algorithm. The solid black line is the planned route, the solid red line records the sampling process of the algorithm, and the solid blue line shows the planning results.

**Figure 10.** (**a**) The blue rectangles represent the sailing corridor. (**b**) The black line is the searched path. The optimized trajectory is shown as the red curve.

To evaluate the smoothness, Figure 11 shows the generated positions, velocities, accelerations, and jerks. All these curves are continuous and satisfy the ship's dynamic constraints.


**Table 2.** Performance comparison in Scene 1 with 1000 trials.

**Figure 11.** The optimized smooth trajectory in terms of position, velocity, acceleration, and jerk.

To highlight the effectiveness of ESP, Figure 12 compares the result of ESP with that which did not consider the ship's dynamic constraints. The blue line denotes the searched path. The green line shows that Bezier curve without considering the dynamic constraints. The red line is the optimized trajectory from ESP.

From the enlarged part, we can see that there are knots in the result without considering the dynamic constraints, making sailing control very difficult. In contrast, ESP generates a smooth and continuous trajectory with the lowest number of turns, which is safer.

Figures 13 and 14 show the collision avoidance in Scene 1 and 2. In Scene 1, ESP generates a feasible and smooth local trajectory (the orange line) that avoids the unexpected circular obstacle. In Scene 2, the original planned trajectory is very close to the shoal, increasing the risk of the ship running aground. The local replanned trajectory effectively solves the problem. In both scenarios, the goal of safely avoiding temporarily appearing static obstacles is achieved.

The computation times of the local replanning in both scenes are listed in Table 3. Over 1000 trials, the best calculation time is 48 ms, and the maximum computation time is 265 ms. The mean computation time is 192 ms in Scene 1 and 194 ms in Scene 2, ensuring real-time processing.

**Figure 12.** ESP vs. *w*/*o* considering the dynamic constraints.

**Figure 13.** Avoid the unexpected circular obstacle in Scene 1.

**Figure 14.** Avoid the shoal to reduce the risk of ship running aground in Scene 2.

**Table 3.** Computation time of the local replanning in both scenes with 1000 trials.


Avoiding multiple obstacles requires multiple iterations of the above process, and the response time will be multiplied by the number of iterations.

#### **5. Conclusions**

This paper proposes a GIS-data-driven method for the efficient and safe path planning of autonomous ships in maritime transportation, which makes up for the shortcomings of existing methods that ignore the motion dynamic limitations of ships in order to achieve the shortest path, leading to sudden changes in the planned route and thus lacking practical applicability. To this end, we propose ESP, a new path planner that provides comfortable sailing while saving fuel. The key intuition of our proposal is to reduce the expensive turning in path searching. The expensiveness comes from the inertia exhibited by the huge weight of the ship. To realize the above intuitive idea, we design three components for ESP. First, we quantify the ship's turning cost based on its kinematic and dynamic model and develop a modified A\* path search algorithm. Second, we formulate an optimization problem subject to dynamic ship constraints and environment constraints to produce a safe and smooth trajectory that consumes minimal fuel. Finally, we use the B-spline representation to perform real-time local replanning, enabling autonomous ships to quickly respond to unexpected risks while maintaining the previous optimization objectives. The data-driven experiments demonstrate the effectiveness of ESP. However, we currently only consider the effect of ocean currents on the dynamic ship model. In future, we will consider the influence of other environmental factors, e.g., sea wind, to build a more robust model for autonomous ships. The avoidance of dynamic obstacles and the real-time avoidance of multiple static obstacles will also be investigated on this basis.

**Author Contributions:** Conceptualization, X.H. and K.H.; methodology, K.H., X.H. and D.T.; software, D.T., X.H. and K.H.; validation, X.H., K.H. and D.T.; formal analysis, X.H., K.H. and D.T.; investigation, D.T., Y.Z. and Y.H.; resources, X.H., K.H., D.T. and Y.H.; data curation, X.H., K.H. and D.T.; writing—original draft preparation, X.H. and D.T.; writing—review and editing, D.T. and X.H.; visualization, D.T. and K.H.; supervision, X.H., K.H. and Y.Z.; project administration, X.H. and Y.Z.; funding acquisition, Y.Z. and Y.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the Research Project of Wuhan University of Technology Chongqing Research Institute (No. YF2021-06). This work was also supported by a grant from the National Natural Science Foundation of China (Grant No. 61801341).

**Data Availability Statement:** The simulation data used to support the findings of this study are available from the corresponding author upon request.

**Conflicts of Interest:** The authors declare no conflict of interest.
