Rapid global path planning algorithm for unmanned surface vehicles in large-scale and multi-island marine environments

Environment map model

In 2D global path planning applications based on the FMM or its improved methods, discrete numerical calculations are based on cartesian grid maps. Therefore, an environment map should first be converted into a binary grid map with a suitable spatial resolution (SR). Free and open-source satellite images (such as Google satellite images) can be used as the data source for the maps in most USV applications, and the corresponding grid maps can be generated using image processing (Shi et al., 2018) combined with manual assistance. The binary cell is set as an obstacle cell when an obstacle exists at the geographic location (0 values); otherwise, it is set as a free cell (1 value).

Fast marching square method

When using FM-based methods to plan a path, the basis is the FMM. The core work of the FMM is calculating solutions of the Eikonal equation (Tsitsiklis, 1995; Adalsteinsson & Sethian, 1995). The Eikonal equation describes a wave front propagation scenario from sources with the speed of a wave front given as F_x at cell x. It can be expressed as $∥\nabla T_x∥F_x=1$, where T_x is the arrival time of the wave from the source to the cell x and ∇ is a vector differential operator. From the perspective of time-cost, the Eikonal equation can be expressed in another form (that of Lin, 2003): (1)${\displaystyle ∥\nabla T_x∥={\tau }_x}$ where τ_x is the time-cost at cell x and is equivalent to 1∕F_x. The solution we want to calculate is T_x, and all of them compose the arrival time map, $\mathcal{T}\left(x\right)$. All time-cost τ_x compose the time-cost function map, $\tau \left(x\right)$.

The FMM was first proposed independently by Tsitsiklis (1995) and Adalsteinsson & Sethian (1995).The solution T_x in cell x can be interpreted as the wave arrival time from the nearest source to cell x. Three cell sets (the accepted set ${\mathcal{S}}_{\mathrm{A}}$, trial set ${\mathcal{S}}_{\mathrm{T}}$, and far set ${\mathcal{S}}_{\mathrm{F}}$) are defined to calculate $\mathcal{T}\left(x\right)$. The set of all cells at which T_x will not change is ${\mathcal{S}}_{\mathrm{A}}.{\mathcal{S}}_{\mathrm{T}}$ is the set of cells to be examined. Every T_x in ${\mathcal{S}}_{\mathrm{T}}$ has been previously computed but may be updated. ${\mathcal{S}}_{\mathrm{F}}$ is the set of all other cells with T_x, which has never been computed (Lin, 2003).

The calculation process includes initialization and loop procedures. During initialization, source cells x₀ with T_x0 = 0 are set to ${\mathcal{S}}_{\mathrm{T}}$, and all other cells are set to ${\mathcal{S}}_{\mathrm{F}}$ with T_x = ∞. After the initialization, the cell x_m with the smallest T value in ${\mathcal{S}}_{\mathrm{T}}$ is selected and moved from ${\mathcal{S}}_{\mathrm{T}}$ to ${\mathcal{S}}_{\mathrm{A}}$. Thereafter, the non-accepted neighbor cells of x_m are updated, including the solutions and cell sets. Considering Tsitsiklis’s deduction (Tsitsiklis, 1995; Lin, 2003), each of the new neighbor solutions is determined by: (2)${\displaystyle T_x=min\left({\tilde{T}}_x,\left(T_{x_m{x}_i},\mathrm{i}=1,2\right)\right)}$ where ${\tilde{T}}_x$ denotes the original solution. T_xmxi, i = 1, 2 are the candidate solutions for the paths that pass through the line segment x_mx_i and propagate to cell x (see Fig. 1), which are calculated by: (3)${\displaystyle T_{x_m{x}_i}=\left\{\begin{array}{cc}\frac{1}{2}\left(T_{x_m}+T_{x_i}+\sqrt{2{\tau }_x^2-{\left(T_{x_m}-T_{x_i}\right)}^2}\right),\hfill & T_{x_m{x}_i}>T_{x_m},\mathrm{and}{T}_{x_m{x}_i}>T_{x_i}\hfill \\ min\left(T_{x_m},T_{x_i}\right)+{\tau }_x,\hfill & \text{others}\hfill \end{array}\right.}$ where T_xm and T_xi are the accepted solutions at cells x_m and x_i, respectively. If cell x_i is in ${\mathcal{S}}_{\mathrm{T}}$, T_xi is ∞. For the cell set, if a non-accepted neighbor cell x is in ${\mathcal{S}}_{\mathrm{F}}$, it moves from ${\mathcal{S}}_{\mathrm{F}}$ to ${\mathcal{S}}_{\mathrm{T}}$. The loop procedure is performed continually until all cells are in ${\mathcal{S}}_{\mathrm{A}}$. The resulting map is the arrival time map, $\mathcal{T}$.

The main disadvantage of the basic FMM is that the computed path is too close to obstacles and forces the vehicle to perform abrupt turns (Garrido, Alvarez & Moreno, 2020). As described in Garrido et al. (2007) and Garrido, Alvarez & Moreno (2020), a smooth path with sufficient safety distances from obstacles can be computed using the FM² method. This method applies the basic FMM twice. The procedure for computing paths is described as follows (Garrido, Alvarez & Moreno, 2020):

1. The environment is modeled as a binary grid map (Fig. 2A).

2. The FM-1st step. All obstacle cells are used as wave sources (T = 0), expanding several waves at the same time at a constant speed. The value of each cell in the resulting map indicates the time required for a wave to reach the closest obstacle (see Fig. 2B). This is proportional to the distance from the obstacles. By reversing the meaning of these values, they can be understood as the maximum admissible speed at each cell. Finally, the speed values are rescaled to fix a maximum cell value of 1.

   Modifications of the speed map with two parameters, α and β, can adjust distances between the computed paths and obstacles. The value of each cell in the speed map F_i,j is adjusted exponentially by α: (4)${\displaystyle new{F}_{i,j}=F_{i,j}^{\alpha }}$ The parameter β is used to saturate the values in the speed map. It is defined within the range of 0 and 1. Every F_i,j with a value greater than β is set to one (see Fig. 2C).

   Comparing with a path without modifications, the path modified by α will be closer to obstacles when α < 1. On the contrary, if α > 1, the modified path will stay further away from obstacles. The parameter β allows the path to move closer to obstacles. When β is smaller, the modified path will be closer to obstacles.

3. The FM-2nd step. The goal point is used as a unique wave source. The wave is expanded over the map until the starting point is reached. The speed at each cell (equivalent to 1∕τ_x) is obtained from the modified speed map computed in the FM-1st step. The resulting arrival time map is shown in Fig. 2D.

4. Finally, a gradient descent is applied over the resulting arrival time map from the starting point to the goal point. An optimal path in terms of the arrival time, smoothness, and safety is obtained.

Gradient descent method

The global path can be extracted by applying the gradient descent method. The path propagates along the gradient descent direction from the starting position P_s to the goal position P_g with a step length d. The value of d can be set to a value equal to the SR of the grid map. The gradient of the cell x_i,j is calculated by: (5)${\displaystyle \nabla T_{i,j}={\left[\begin{array}{cc}\hfill \frac{T_{i+1,j}-T_{i-1,j}}{2}\hfill & \hfill \frac{T_{i,j+1}-T_{i,j-1}}{2}\hfill \end{array}\right]}^{\mathrm{T}}}$ where T_i,j represents the arrival time value at cell x_i,j and the non-italic T in the upper-right corner represents the transpose of vector $\left[xy\right]$.

The path waypoints are not limited to cell centers. Generally, the gradient of a path waypoint $P_n=\left(u,v\right)$ located in the cell x_i,j is approximated by the cell gradient ∇T_i,j. However, a modified gradient ∇T_u,v calculated by the bilinear interpolation can be used to improve the path precision (see Fig. 3).

Mapping of two-level SR grid maps

The two-level SR grid maps are contained in the HSR and LSR grid maps. The HSR grid map is directly obtained from the Google satellite image data. A mapping relationship between the LSR cells and the corresponding L × L HSR cell sub-blocks is established, which can be expressed as: (6)${\displaystyle \left(i_{\mathrm{L}},j_{\mathrm{L}}\right)\sim \left[\begin{array}{ccc}\hfill \left(i_{\mathrm{H}}^{\prime },j_{\mathrm{H}}^{\prime }+L-1\right)\hfill & \hfill \cdots \hfill & \hfill \left(i_{\mathrm{H}}^{\prime }+L-1,j_{\mathrm{H}}^{\prime }+L-1\right)\hfill \\ \hfill ⋮\hfill & \hfill \ddots \hfill & \hfill ⋮\hfill \\ \hfill \left(i_{\mathrm{H}}^{\prime },j_{\mathrm{H}}^{\prime }\right)\hfill & \hfill \cdots \hfill & \hfill \left(i_{\mathrm{H}}^{\prime }+L-1,j_{\mathrm{H}}^{\prime }\right)\hfill \end{array}\right]}$ where $\left(i_{\mathrm{L}},j_{\mathrm{L}}\right)$ are the LSR cell coordinates, and $\left(i_{\mathrm{H}}^{\prime },j_{\mathrm{H}}^{\prime }\right)$ are the original cell coordinates of the mapped HSR cell sub-block, as shown in Fig. 4.

To map the LSR grid map to a sub-block of the HSR grid map, the following is used: (7)${\displaystyle \left\{\begin{array}{c}{i}_{\mathrm{H}}^{\prime }=i_{\mathrm{Ho}}+L{i}_{\mathrm{L}}\hfill \\ j_{\mathrm{H}}^{\prime }=j_{\mathrm{Ho}}+L{j}_{\mathrm{L}}\hfill \end{array}\right.}$

To map from a sub-block of the HSR grid map to a cell of the LSR grid map, we use: (8)${\displaystyle \left\{\begin{array}{c}{i}_{\mathrm{L}}=f_{\text{floor}}\left(\frac{i_{\mathrm{H}}-i_{\mathrm{Ho}}}{L}\right)\hfill \\ j_{\mathrm{L}}=f_{\text{floor}}\left(\frac{j_{\mathrm{H}}-j_{\mathrm{Ho}}}{L}\right)\hfill \end{array}\right.}$ where $\left(i_{\mathrm{H}},j_{\mathrm{H}}\right),i_{\mathrm{H}}^{\prime }\le i_{\mathrm{H}}\le i_{\mathrm{H}}^{\prime }+L-1,j_{\mathrm{H}}^{\prime }\le j_{\mathrm{H}}\le j_{\mathrm{H}}^{\prime }+L-1$ are the HSR cell coordinates. $\left(i_{\mathrm{Ho}},j_{\mathrm{Ho}}\right)$ are the original cell coordinates of the original HSR cell sub-block (see Fig. 4), which are determined by: (9)${\displaystyle \left\{\begin{array}{c}{i}_{\mathrm{Ho}}=\left(i_{\mathrm{Hg}}-f_{\text{floor}}\left(\frac{L}{2}\right)\right)\text{\%}L\hfill \\ j_{\mathrm{Ho}}=\left(j_{\mathrm{Hg}}-f_{\text{floor}}\left(\frac{L}{2}\right)\right)\text{\%}L\hfill \end{array}\right.}$ where $\left(i_{\mathrm{Hg}},j_{\mathrm{Hg}}\right)$ are the HSR cell coordinates of the goal cell x_g, the function $f_{\text{floor}}\left(m\right)$ indicates a rounding down of m, and a%b indicates the modulo operation.

The SR of the LSR grid map is: (10)${\displaystyle D_{\text{LRes}}=L{D}_{\text{HRes}}}$ where D_HRes is the SR of the HSR grid map. In general, L is defined as: (11)${\displaystyle L\le f_{\text{round}}\left(\frac{D_{\mathrm{Th}}}{2D_{\text{HRes}}}\right)}$ where D_Th is the distance threshold of the virtual obstacle influence, and the function $f_{\text{round}}\left(m\right)$ indicates a rounding of m. The cell numbers of the LSR grid map in the X-axis and Y-axis directions are: (12)${\displaystyle M_{\mathrm{L}}=f_{\text{floor}}\left(\frac{M_{\mathrm{H}}-i_{\mathrm{Ho}}}{L}\right)}$ (13)${\displaystyle N_{\mathrm{L}}=f_{\text{floor}}\left(\frac{N_{\mathrm{H}}-j_{\mathrm{Ho}}}{L}\right)}$ where M_H and N_H are the cell numbers of the HSR grid map in the X-axis and Y-axis directions, respectively.

Based on the mapping relationship, an LSR cell type is determined by comparing the setting threshold Γ and obstacle cell proportion γ in the HSR cell sub-block. If γ > Γ, the LSR cell is set as an obstacle cell; otherwise, it is set as a free cell. Γ prefers to select a small value to retain as much obstacle information as possible.

IDC-FM² method

In the basic FM² method, the path is improved using a modified speed map. Two parameters, α and β, are used to modify the speed map. These modifications ensure the flexibility of the computed path. However, these two parameters have no clear physical meaning. The normalization of the speed values in the FM-1st step also results in several shortcomings. For example, all the speed values in the speed map must be calculated. When the map range changes slightly, the entire speed map rescales and differs in the same locations of the original map, and the parameters must be adjusted.

An IDC-FM² method is proposed to address these shortcomings. In this method, a time-cost weighting function, w_x, is introduced to adjust the speed map. Three inshore-distance parameters (the distance threshold D_Th, as well as the virtual strong and weak constraint distance parameters, D_sc and D_wc, respectively) and two weighting function values with respect to D_sc and D_wc (w_sc and w_wc, respectively) are used to determine the function w_x, and D_Th is used to achieve the same function as the speed saturation. Although five parameters are set, only two inshore-distance parameters, D_Th and D_sc, must be considered, while suitable values of w_sc and w_wc can be selected as constant values. Based on this method, the path portion around an obstacle is constrained in a region determined by D_Th and D_sc.

Parameters for the time-cost weighting function

Five parameters, D_Th, D_sc, D_wc, w_sc, and w_wc, are introduced to determine the coefficients a and b, which are determined by: (18)${\displaystyle b=\frac{ln\left(w_{\mathrm{sc}}-1\right)-ln\left(w_{\mathrm{wc}}-1\right)}{ln\left(1-e_{\mathrm{sc}}\right)-ln\left(1-e_{\mathrm{wc}}\right)+ln{e}_{\mathrm{wc}}-ln{e}_{\mathrm{sc}}}}$ (19)${\displaystyle a=\left(w_{\mathrm{sc}}-1\right){\left(\frac{e_{\mathrm{sc}}}{1-e_{\mathrm{sc}}}\right)}^b}$ where e_sc = D_sc∕D_Th, e_wc = D_wc∕D_Th, D_sc < D_wc < D_Th, and w_sc > w_wc > 1.

The parameter D_Th determines the largest virtual influence scope of the obstacles. The suggestion for the selection of D_Th is that there should be a sufficiently safe buffer area for obstacles. The parameter D_sc is very important for the safety of the USV. First, it is suggested that this parameter satisfies: (20)${\displaystyle D_{\mathrm{sc}}\ge v_{\mathrm{U},max}{t}_{\mathrm{r}}-\frac{v_{\mathrm{U},max}^2}{2a_{\mathrm{d}}}}$ where v_U,max is the maximum speed of the USV, t_r is the reaction time of the thruster, and a_d is the negative maximum acceleration of the USV under braking. In an unknown environment, it is better to select a slightly larger value to ensure safety, which is similar in environments with shallow and reef waters. This value may be relatively small in deep and reef-free waters to improve path efficiency. D_wc is a parameter that acts as the adjusted constraint distance. In a non-channel ocean area, the path portions around obstacles will be limited to the regions between D_wc and D_Th. The D_wc value can be set independently or determined jointly by D_Th and D_sc. In this study, it is set as: (21)${\displaystyle D_{\mathrm{wc}}=D_{\mathrm{Th}}-\frac{\sqrt{2}}{2}\left(D_{\mathrm{Th}}-D_{\mathrm{sc}}\right)}$

The area around the obstacles is divided into four parts by the three inshore-distance parameters: D_sc, D_wc, and D_Th. The three parts extending from an obstacle to the outside region are the danger region, S_d; the strong constraint region, S_sc; and the weak constraint region, S_wc (see Fig. 5). The influences of w_wc and w_sc on the path are shown in Fig. 5. As shown in Fig. 5A, the path will be located far from the obstacle until it is near the boundary of S_wc when w_wc is large. In contrast, the path is closer to the obstacle when w_sc is large (see Fig. 5B). The minimum degree of the path close to S_sc is mainly determined by w_wc, which can be inferred by further comparing the paths in Fig. 5B with the closest path to the obstacle (w_wc = 1.1) in Fig. 5A. However, regardless of the change in w_wc and w_sc (w_sc > w_wc > 1), the path portions bypassing obstacles will always be within S_wc or around the outside boundary of S_wc in non-channel ocean areas. In channel areas where there is no S_wc, the path will follow the quasi-centerline of the channel if it passes through the channel. In this study, the values of w_sc = 40.0 and w_wc = 2.0 were selected and fixed. The paths can be flexibly modified by the inshore-distance parameters D_Th and D_sc.

Determination of effective regions within the HSR grid map

The effective region within the HSR grid map is important for improving the computational efficiency of the proposed algorithm. It is acquired based on the mapping between the LSR and HSR grid maps when the corresponding region within the LSR grid map is determined.

Two effective regions that are respectively used in the FM-1st and FM-2nd steps of the IDC-FM² method must be determined. The low-precision initial path, ℓ_ini, is obtained first based on the LSR arrival time map ${\mathcal{T}}_{\mathrm{L}}$. Two effective regions within the LSR grid map, S_{LER_1st} and S_{LER_2nd}, are determined by expanding the cells around the “path-passed” cells of ℓ_ini. As shown in Fig. 6, the nearest neighbor cell among the four neighbor cells of the waypoint $P_i=\left(x_i,y_i\right)$ is defined as the “path-passed” cells. When all “path-passed” cells are determined, the κ_2nd = κ levels of the cells are extended to obtain S_{LER_2nd}. Considering S_{LER_1st}, there are three different situations:

• Situation 1: there is no cell with w_x > 1 in S_{LER_2nd};

• Situation 2: there is at least one cell with w_x > 1 and no obstacle cell in S_{LER_2nd}.

• Situation 3: there is at least one obstacle cell in S_{LER_2nd}.

In Situations 1 and 3, κ_1st = κ. In Situation 2, S_{LER_2nd} is extended continuously with Δκ levels of cells until there is at least one obstacle cell for the first time, and κ_1st = κ + Δκ. Finally, two corresponding effective regions within the HSR grid map, S_{HER_1st} and S_{HER_2nd}, are determined and used in the different FM steps of the IDC-FM² method to obtain the HSR arrival time map ${\mathcal{T}}_{\mathrm{H}}$.

The extended parameter κ affects the paths (recorded as ℓ_κ) based on the proposed algorithm. When κ is not adequately large, these paths may not be consistent with the reference path ℓ_ref, which is acquired by applying the IDC-FM² method in the HSR grid map directly. A navigable global path is used as a case to compare the consistency between ℓ_κ and ℓ_ref. The basic parameters used in the proposed algorithm are listed in Table 1. Using ℓ_ref as the reference, the distances between the corresponding waypoints of ℓ_κ (κ = 3, …, 6) and ℓ_ref are shown in Fig. 7. It is shown that the largest distance decreases when κ increases, and the distance becomes 0 when κ ≥ 7. Other cases were tested and similar results were obtained. These results indicate that there is good consistency when κ is sufficiently large (such as κ = 10). Δκ is a variable parameter that was determined based on three different cases.

Algorithm flow

The proposed algorithm is improved based on the basic FM²-based algorithm. The basic FM²-based algorithm is executed directly on a single grid map. The algorithm flow of this basic algorithm is shown in Fig. 8 in the form of main data stream and used methods. This algorithm flow has two main steps. In Step S1, an arrival time map with obstacle sources and a saturation threshold, ${\mathcal{T}}_{\mathrm{sat}}$, is obtained by first applying the FMM. Then, the approximate time-cost weighting function map, W, is calculated and used to adjust the time-cost map. Based on the adjusted time-cost map, an arrival time map with the goal point source, $\mathcal{T}$, is obtained by applying the FMM again. The entire process is used to obtain $\mathcal{T}$ by applying the IDC-FM² method. Finally, the global path is acquired based on $\mathcal{T}$ by applying the gradient descent method in Step S2.

The main data stream in the proposed algorithm and the methods used to obtain them are shown in Fig. 9. The main flow is as follows:

1. Step T1. Obtain the two-level SR grid maps;

2. Step T2. Obtain the LSR arrival time map ${\mathcal{T}}_{\mathrm{L}}$ in the LSR grid map;

3. Step T3. Determine the effective regions within the HSR grid map (S_{HER_1st} and S_{HER_2nd});

4. Step T4. Obtain the HSR grid arrival time map ${\mathcal{T}}_{\mathrm{H}}$ with the effective region constraint in both FM steps of the IDC-FM² method. There are several adjustments in both FM steps compared to the basic FMM performance. For the FM-1st step, all w_x values within S_{HER_1st} need not be adjusted, and this step can be ignored in Situation 1. In Situation 2 and Situation 3, all free cells out of S_{HER_1st} are set to ${\mathcal{S}}_{\mathrm{A}}$ with T = ∞ in the initialization of this step. In the FM-2nd step, all free cells out of S_{HER_2nd} are set to ${\mathcal{S}}_{\mathrm{A}}$ with T = ∞ in the initialization;

5. Step T5. Acquire the final high precision global path by applying the gradient descent method based on ${\mathcal{T}}_{\mathrm{H}}$.

Simulation environments

Two surrounding spatial areas of Zhucha Island and Changhai County were selected as the simulation environments. The original maps adopted Google satellite images, as shown in Fig. 10. The spatial resolutions in the longitude and latitude directions are approximately 4.8 m and 3.9 m, respectively. Temporary binary grid maps are obtained based on the image processing method (Shi et al., 2018) combined with manual assistance, and then the HSR grid maps with D_HRes = 10m in both the longitude and latitude directions are determined by resampling processing. The corresponding binary grid maps are presented in Fig. 11. Their ranges are 7 km × 7 km and 64 km × 48km, respectively.

Inshore-distance-constraint performances

A path planning from P_s = [4.3 km, 2.8 km] to P_g = [3.6 km, 2.0 km] is used as a typical case to analyze the inshore-distance-constraint performance of the IDC-FM² method using the inshore-distance parameters D_Th and D_sc. As shown in Fig. 12, the path planned based on the basic FMM, ℓ_FMM, is very close to the islands when this path bypasses them. For comparison, all paths planned by the proposed algorithm, ℓ₁ to ℓ₄, are located away from islands by a certain distance. They are therefore significantly better choices from a safety perspective.

Paths ℓ₁ to ℓ₄ are acquired using different inshore-constraint distance parameters (D_Th and D_sc, see Table 2 wherein D_wc is calculated based on D_Th and D_sc) in the IDC-FM² method in which w_sc = 40.0 and w_wc = 2.0. These paths are clearly adjusted by these distance parameters. When the distance constraints (D_Th and D_sc) are small, the path (such as ℓ₁, see Fig. 12) will be a somewhat close to the islands. The estimated quasi-closest distance from ℓ₁ to an island is approximately 38 m, which is shorter than half of the shortest channel width (d_hscw ≈ 101m). In this situation, ℓ₁ is always outside of its S_sc (<28.2 m) value, and the path portions around the islands are located in the region of its S_wc (28.2 m to 60.0 m). When the distance constraints increase, the paths (such as ℓ₂ and ℓ₃, see Fig. 12) in the non-channel areas will be outside of S_sc. The estimated quasi-closest distances (approximately 115 m and 128 m, respectively) are larger than d_hscw. However, ℓ₂ and ℓ₃ will be along the quasi-midline of the channel in the channel area, regardless of whether d_hscw is greater or less than D_wc for these two paths because the weighting time-cost will be smaller than the path around the islands (as in ℓ₄, see Fig. 12). If the distance constraints are further increased, ℓ₄ will be selected for safety, although more time will be required. These results indicate that the path based on the proposed algorithm can be adjusted by setting different inshore-constraint distance parameters to meet safety requirements. Additionally, when the inshore-constraint distance parameters are determined, the path will cross a channel when its width is adequately large; otherwise, the path will bypass around the islands.

One classical approach to plan collision-free paths is applying the FMM within the map with inflated obstacles. For comparison, paths planned by the classical approach and the IDC-FM² method are shown in Fig. 13. The starting and goal positions are P_s = [4.15 km, 2.74 km] to P_g = [4.33 km, 2.40 km] and the island obstacles have been inflated by 94.0 m. As shown in Fig. 13, the path planned by the classical approach, ℓ_FMM+inflated, is obviously influenced by the shape of the inflated obstacle. In some path turns which are marked by ellipses in Fig. 13, they are both affected by the shape of the inflated obstacle and the time-value gradient calculation of cells adjacent to the inflated obstacle. A smooth path, ℓ_{FMM+inflated+smooth}, is also shown in Fig. 13. Although it removes abrupt turn waypoints, the path is still not very smooth. As comparisons, paths planed based on the IDC-FM² method, ℓ₁ and ℓ₂, are obviously smoother than ℓ_{FMM+inflated+smooth}. In addition, when using the proposed rapid path planning algorithm based on two-level SR grid maps to improve the computational efficiency, the channel will be mapped as obstacle cells when mapping the HSR grid map to the LSR grid map in this case (taking the mapping parameter L = 8 as shown in Table 1), because of the influence of the inflated obstacle. This unexpected mapping will result in the loss of the path through the channel.

Global path planning in a large-scale and complex multi-island environment

To verify the path planning ability and the computational efficiency improvement of the proposed algorithm in a large-scale and relatively complex multi-island environment, the simulation environment shown in Fig. 11B is selected, and the long-length path cases of USVs around the islands or across channels are investigated.

Path planning cases

Five typical path cases were selected to demonstrate the performance of the proposed algorithm. The P_s and P_g groups are listed in Table 3. When determining the inshore-distance parameter D_sc, a USV with v_U,max = 6m∕s, t_r = 2s, a_d =  − 1m∕s² is considered as a case. Based on Eq. (20), D_sc is suggested to select a value larger than 30 m. Therefore, a larger value D_sc = 50m is used for the path planning of the selected five path cases. Another inshore-distance parameter D_Th = 200 m is selected empirically. The main parameters used in these path planning cases are shown in Table 1.

Table 3:

Quantified information of the five typical paths.

Path	Ps(km)	Pg(km)	Length (km)
ℓ1	[35.34, 39.25]	[15.31, 11.65]	38.84
ℓ2	[19.42, 41.02]	[17.10, 3.63]	37.90
ℓ3	[42.96, 43.67]	[46.34, 8.24]	38.20
ℓ4	[36.11, 18.77]	[47.44, 41.01]	26.47
ℓ5	[3.95, 26.52]	[50.45, 30.83]	47.82

DOI: 10.7717/peerjcs.612/table-3

As shown in Fig. 14, all paths successfully bypass the islands. The enlarged views (see Fig. 15, every path is displayed with a solid line) provide further details. Every path maintains a relatively safe distance when close to an island, and smooth when turning around the island. For a narrow channel of a certain degree (usually wider than 2D_LRes), the path is planned along the quasi-midline of the channel to ensure that it is as safety as possible. The path lengths range from about 26.52 km to 47.86 km (see Table 3). These lengths can cover the range of most applications of current small-and medium-sized USVs. These results show the effective path planning ability of the proposed algorithm in a large-scale and complex multi-island environment.