A Global Path Planning Method for Unmanned Ground Vehicles in Off-Road Environments Based on Mobility Prediction

Abstract

In a complex off-road environment, due to the low bearing capacity of the soil and the uneven features of the terrain, generating a safe and effective global route for unmanned ground vehicles (UGVs) is critical for the success of their motion and mission. Most traditional global path planning methods simply take the shortest path length as the optimization objective, which makes it difficult to plan a feasible and safe route in complex off-road environments. To address this problem, this research proposes a global path planning method, which considers the influence of terrain factors and soil mechanics on UGV mobility. First, we established a high-resolution 3D terrain model with remote sensing elevation terrain data, land use and soil type distribution data, based on a geostatistical method. Second, we analyzed the vehicle mobility by the terramechanical method (i.e., vehicle cone index and Bakker’s theory), and then calculated the mobility cost based on a fuzzy inference method. Finally, based on the calculated mobility cost, the probabilistic roadmap method was used to establish the connected matrix and the multi-dimensional traffic cost evaluation matrix among the sampling nodes, and then an improved A* algorithm was proposed to generate the global route.

1. Introduction

Due to their unnecessary requirement of human interference and their ability to replace human operations in risky environments, unmanned ground vehicles (UGVs) have several applications, such as disaster rescue, nuclear inspection, planet exploration and military operations, etc. In order to enable a UGV to traverse complex off-road environments quickly and safely and achieve its missions, an effective global path planning method is required, of which the key objective is to generate a safe and feasible route from the initial position to the specified targets based on given off-road environmental map information [1,2].

In recent decades, path planning methods developed rapidly, resulting in numerous new such methods [3]. Path planning methods can be mainly divided into two categories, i.e., graph search-based and sampling-based [3]. Typical graph search-based methods include Dijkstra, A*, D*, field D* and several other variants [4,5,6,7]. The graph search-based path planning methods can find the optimal path accurately, but such algorithms should traverse the surrounding of the current node and select the minimum cost node, which leads to large computational efforts. To reduce the calculation time the algorithm takes to search each interval, an efficient A* algorithm was proposed [8]. Guruji et al. [9] proposed an online path planning method for a mobile robot in a grid-map environment, using a modified real time A* algorithm. In obstacle-avoidance path planning, Duan et al. introduced a safety–cost matrix in the A* algorithm and the heuristic function was modified to ensure the optimal safety path in different environments [10].

Compared with graph search-based methods, sampling-based path planning methods have fast search efficiency, and classical sample-based approaches include rapidly exploring random trees (RRT) [11], probabilistic road maps (PRM) [12] and several other variants [13]. However, the random nature of sampling may lead to the tendency of falling into the local minimum in complex environments [14]. To address this problem, a potentially guided intelligent bi-directional RRT* algorithm was proposed [15], and the proposed algorithm greatly improved the convergence rate and had more efficient memory utilization. Also, Park et al. proposed a hierarchical roadmap (HRM) approach that enabled the indoor mobile robot to cover navigable areas using fewer nodes than a probabilistic roadmap [16].

However, most existing developed planning methods focus on structured environments, such as indoor and urban scenarios. Due to the low bearing capacity of the soil and the uneven features of the terrain, existing methods are very difficult to apply to complex unstructured off-road environments. To adapt to the off-road environments, researchers have proposed several improved planning methods based on traditional methods. For example, to evaluate the comprehensive influence of terrain slope and surface attributes, an improved A* algorithm with a window moving method was proposed; this method can plan a gentle global route for the mobile robot, while the robot can effectively avoid moving obstacles [17]. Similarly, a method was proposed to systematically deal with the effects of kinematic vehicle models, off-road terrains and obstacles to UGV path planning by constructing artificial potential fields [18]. Also, there is an increased risk of vehicle obstruction due to soil deformation. Gonzalez et al. [19] proposed a path planning method based on the classification of terrain into GO/NOGO areas by soil strength value (CI); it is worth mentioning that the method was integrated into the next-generation NATO reference mobility model (NG-NRMM). With the purpose of generating a high confidence global route in the off-road, Jiang et al. [20] employed a global path planning method based on the RRT* algorithm with the cost function of the uncertainty in terrain elevation. In terms of global path smoothing, the geometric curve method transforms the path planning problem into a two-point boundary value problem, such as the B-spline curve [21] and fifth degree polynomial curve [22].

In spite of progress, most existing global path planning methods ignore the effects of soil deformation on vehicle mobility. Although Gonzalez et al. proposed classifying off-road terrain into GO and NOGO based on the soil strength value (CI), it is still difficult to plan an optimal global route for UGV within the GO areas because of the lack of a well-designed cost function to quantify the effects of vehicle–terrain interaction. To address this problem, the contributions of this study are as follows: (1) In order to obtain 3D models of unknown off-road environments effectively, this paper proposes a Kriging terrain interpolation method based on open-source remote sensing DEM data with multi-source data fusion, which is general and easy to implement. (2) To address the problem that the current global path planning methods for unmanned vehicles do not take into account the impact of soil deformation on mobility, this paper proposes an analysis method based on terramechanical theory to calculate the impact of soil deformation on the speed of a UGV. Then, a method to quantify the cost of mobility based on fuzzy inference fusion of multiple influencing factors is proposed. (3) In terms of a global path generation algorithm, an improved A* algorithm based on mobility cost is proposed to generate a topological map by sampling terrain points, and to introduce mobility costs, in both extended cost and heuristic cost, to generate an efficient, feasible and safe global route. Finally, the effectiveness of the proposed method is validated by simulation and comparison experiments.

The remainder of this paper is organized as follows. Section 2 describes the details of the proposed method. The simulation experiments and discussions are given in Section 3. Section 4 concludes the paper and proposes future work.

2. Methods

The proposed global planning method include three phases, i.e., terrain reconstruction, mobility prediction and path generation. In the terrain reconstruction phase, the original sparse digital elevation model (DEM) is firstly converted into a high-resolution 3D point cloud using a geostatistics interpolation method, and then the land use and soil type distribution data are integrated to generate a multisource terrain model. In the mobility prediction phase, the vehicle cone index and Bakker’s theory-based terramechanical model is proposed to calculate the maximum attainable speed of the UGV on different soil areas, and a fuzzy inference method is utilized to generate the global mobility cost map. In the path generation phase, an improved A* algorithm with extended cost and heuristic cost based on the mobility cost map is proposed to generate an optimal global route. Figure 1 shows the method framework.

2.1. Terrain Reconstruction

Terrain geometry information typically comes from remote sensing sources (e.g., radar, imagery camera, etc.). Since the resolution of the original DEM is too sparse for accurate prediction of vehicle mobility, an interpolation process is required to obtain a refined terrain model with a higher resolution. In our study, open-source remote sensing topographic data with a resolution of 30 m was used as the raw data; the required resolution for mobility prediction is usually about 10 m [23]. To further characterize other features of the selected terrain, the obtained land use data and soil type distribution data were projected onto the resulting refined DEM data. The reconstruction result was saved in the form of a 3D point cloud for the following processes.

In this research, we selected a target area of about 40 km² within the range of 97.2°–97.3° E and 31.4°–31.5° N. After determining the target area, the satellite map and the DEM were downloaded on the geospatial data cloud website. In addition to elevation, the file contains a layer with different soil regions classified according to the National Earth System Science Data Center. Specifically, there are four soil types include sandy loam, clayey soil, lean clay and upland sandy. Also, land use data can be obtained using the method described above, and the data include water area, forest land, grass land and cultivated land. All the data are shown in Figure 2.

In the process of reconstructing the terrain, we needed to obtain the coordinates of the DEM data from the geographic coordinates (Longitude–Latitude). This preprocessing process generates the new DEM data in the Cartesian format. Each location was then formatted in terms of planar coordinates (X-Y) and elevation (Z), all of which are in metric units. The resulting 3D terrain model is displayed in Figure 3; the soil type distribution data and the land use data were also projected onto the corresponding terrain points (see Appendix A.1 the Projection of each Data Type.

The process of natural terrain modeling begins with a sparse set of measurements obtained from remote sensors on the terrain area of interest. An important step in simulating the performance of a UGV on such terrain is to generate a continuous surface; this process can be achieved by interpolating the elevations of unknown terrain points from the known terrain points. The most widely used terrain interpolation method in geostatistics is the Kriging [24], which is based on the principle of estimating the elevations of interpolated points from the variance of the known points. The variogram represents the statistical relationship between terrain points. This shows how the heterogeneity between $z\left(s\right)$ and $z\left(s+h\right)$ evolves with the lag distance $h$; its calculation formula is

$$\gamma \left(h\right)=\frac{1}{2N}{\displaystyle \sum }_{i=1}^N{\left(z\left(s_i\right)-z\left(s_i+h\right)\right)}^2$$

(1)

where N is the number of terrain points, $z\left(s_i\right)$ is the elevation at point $s_i$ and $h$ represents the separation distance, which is the vector separating the locations of two variables. For example, $h=1\times d, \dots ,n\times d$ d is the separation distance between points. If points are sampled at 10 m, then d = 10 and the semi-variogram is calculated for the separation distance $h$.

When the variogram model has been built, the Kriging method can be used for the terrain interpolation. There are three main forms of Kriging, depending on whether the mean of the data is known (simple Kriging) or unknown (ordinary Kriging and universal Kriging). In geostatistics, the most common method used is the ordinary Kriging [24], because the mean is not a priori and only spatial variation is meaningful. According to the ordinary Kriging formula, every terrain point with an unknown elevation is interpolated based on a weighted combination of the known terrain points, the expression being

$$\widehat{Z}\left(s_0\right)={\displaystyle \sum }_{i=1}^N{\omega }_{i}Z\left(s_i\right)$$

(2)

where $\widehat{Z}\left(s_0\right)$ is the predicted elevation of the sample located at $s_0$(interpolated terrain point) and ${\omega }_i$ are the weights for the combination of the known terrain points $Z\left(s_i\right).$ The parameter i represents the samples available in the dataset. Ordinary Kriging minimizes the mean squared error of prediction, and its expression is

$${\sigma }_e^2=2{\displaystyle \sum }_{i=1}^N{\omega }_i\gamma \left(s_0-s_i\right)-{\displaystyle \sum }_{i=1}^N{\displaystyle \sum }_{j=1}^N{\omega }_i{\omega }_j\gamma \left(s_i-s_j\right)$$

(3)

After differentiating Equation (3) and setting the derivative to zero, it follows that

$${\displaystyle \sum }_{j=1}^N{\omega }_j\gamma \left(s_i-s_j\right)+\lambda =\gamma \left(s_0-s_i\right)$$

(4)

where $\lambda$ is the Lagrange parameter used for ensuring ${\sum }_{i=1}^N{\omega }_i=1$. By using matrix expressions, the previous set of equations can be written as

$$\begin{gathered} W=C^{-1}D, \\ W=\left[{\omega }_1,{\omega }_2,\dots ,{\omega }_N,\lambda \right], \\ C=\gamma \left(s_i-s_j\right), i=1,\dots ,N; j=1,\dots ,N, \\ D={\left[\gamma \left(s_0-s_1\right), \gamma \left(s_0-s_2\right),\dots ,\gamma \left(s_0-s_N\right)\right]}^{\mathrm{T}} \end{gathered}$$

(5)

Once the weights are calculated, the predicted elevations at the unknown locations can be estimated using Formula (2). The reconstructed DEM is shown in Figure 4a. The slope of the terrain is also an important factor in limiting the mobility of the vehicle, as large slope angles can cause the vehicle to roll over. Therefore, after obtaining the refined 3D terrain model, the slope value of each terrain point is calculated. Specifically, a $3\times 3$ window template is used to traverse all points of the terrain; the slope map is shown in Figure 4b (see Appendix A.2 Slope Calculation Method).

2.2. Terramechanic Analysis

The off-road vehicle mobility prediction model is usually employed in the mobility-analysis phase to map environmental variables, such as soil properties, for example, into the vehicle mobility (vehicle maximum attainable speed). Various theoretical models were developed; for instance, the U.S. Army Engineer Research and Development Center (ERDC) defines vehicle mobility in two different ways. One is the single vehicle cone index (VCI₁) and the other is the fifty-pass vehicle cone index (VCI₅₀) [25]. VCI₁ and VCI₅₀ represent the minimum soil strength required for a vehicle to consistently make a specific number of passes. In the NATO Reference Mobility Model (NRMM) [26], VCI is used as a metric to evaluate the off-road mobility of military ground vehicles on different soils. The mobility performance is judged by the following:

·

If CI $\ge$ VCI₅₀, pass easily

·

If VCI₅₀ > CI $\ge$ VCI₁, pass difficultly

·

If CI < VCI₁, not pass

Where CI (cone index) refers to the ratio of the force required to vertically press the cone penetrometer into the soil at a certain depth at a speed of 0.5 cm/s to the area of the cone. The value of VCI is related to the vehicle mobility index (MI), which is related to the ground pressure coefficient, axle load coefficient, tire coefficient, anti-skid coefficient, engine coefficient and transmission coefficient of the vehicle, reflecting the overall mobility of the vehicle. The expression of MI is

$$\mathrm{MI}=\left(\frac{P_{FG}W}{TG}+L-H\right)EX$$

(6)

where $P_{FG}$ is grounding pressure coefficient, $E$ is the engine coefficient, $G$ is the track tooth coefficient or wheel spur coefficient, $T$ is the tire coefficient, $L$ is the load factor, $H$ is the ground clearance and $X$ is the transmission coefficient.

According to the MI, the empirical formula can be used to calculate the VCI values of a single pass and fifty passes [25], which are VCI₁ and VCI_50, respectively.

For single-pass wheeled vehicles with MI ≤ 115:

$${\mathrm{VCI}}_1=11.48+0.2\mathrm{MI}-\left(\frac{39.2}{\mathrm{MI}+3.74}\right)$$

(7)

For single-pass wheeled vehicles with MI > 115:

$${ \mathrm{VCI}}_1=4.1{\mathrm{MI}}^{0.446}$$

(8)

For wheeled vehicles with 50 passes:

$${\mathrm{VCI}}_{50}=28.23+0.43\mathrm{MI}-\left(\frac{92.67}{\mathrm{MI}+3.67}\right)$$

(9)

To express the interaction between the vehicle and the ground, Bakker measured the mechanical properties of foundation soil experimentally [27], then Bakker’s derived terramechanic model (BDTM) was created. This model was applied to the NATO Reference Mobility Model (NRMM) to obtain the mobility prediction of military vehicles. In particular, when the structural parameters of the vehicle are known, the maximum attainable speed of the vehicle on a soft soil foundation mainly depends on soil characteristics. Considering these, BDTM is used to calculate the maximum attainable speed.

This paper takes wheeled unmanned vehicles as the object of analysis, when the pneumatic tire traverses off-road, it rolls in different states due to the different strengths of the soil. When the soil is soft and the sum of the tire inflation pressure generated by the tire stiffness is greater than the support force of the soil to the lowest point of the tire, the pneumatic tire rolls like a rigid tire; this phenomenon is referred to as the rigid mode of the tire, shown in Figure 5a. In contrast, on relatively solid soil, the tread-grounding parts of tires will be pressed into a plane, also known as the elastic mode, as shown in Figure 5b.

To calculate the maximum attainable speed of the vehicle, it is necessary to predict rolling resistance ($F_{rc}$). When the tire rolls on the soft soil, it forms a rut with a depth of $z_0$, and the contact surface between the tire and the soil is divided into two areas, the curved surface and the plane. According to the pressure and sinkage relationship of BDTM, the calculated unit pressure in the curved surface area is

$$p=\left(\frac{k_c}{b}+k_{\phi }\right)z_0^n$$

(10)

The unit pressure in the plane area is

$$p=p_f+p_c$$

(11)

where $k_c$, $k_{\phi }$ and n are Bakker’s soil parameters, b is the tire width and $z_0$ is the sinkage.

Assuming that the diameter of the rigid wheel is $D$, if the reaction force of the soil on the rigid wheel is only radial, the soil resistance ($F_{rc}$) of the wheel can be obtained as

$$F_{rc}=b{{\displaystyle \int }}_0^{z_0}\left(\frac{k_c}{b}+k_{\phi }\right)z^{n}dz=b\left(\frac{k_c}{b}+k_{\phi }\right)\frac{z_0^{n+1}}{n+1}$$

(12)

In the rigid mode, the soil is compressed only in the vertical direction, so the value of $F_{rc}$ calculated from Equation (12) should be equal to the work done when a flat plate of width $b$ is pressed vertically into the soil to a depth of $z_0$. The vertical load $W$ on the wheel is

$$W=-b{{\displaystyle \int }}_0^{Z_0}pdx=-b{{\displaystyle \int }}_0^{Z_0}\left(\frac{k_c}{b}+k_{\phi }\right)z^{n}dz$$

(13)

This can be obtained from the geometric relationship of each variable:

$$W=\frac{1}{3}b\left(\frac{k_c}{b}+k_{\phi }\right)\sqrt{z_{0}D}{z}_0^n\left(3-n\right)$$

(14)

Combine Equations (13) and (14):

$$Z_0={\left[\frac{3W}{b\left(\frac{k_c}{b}+k_{\phi }\right)\sqrt{D}\left(3-n\right)}\right]}^{\frac{2}{2n+1}}$$

(15)

$$F_{rc}=\frac{{\left(3W\right)}^{\frac{2n+2}{2n+1}}}{{\left(3-n\right)}^{\frac{2n+2}{2n+1}}\left(n+1\right){\left(k_c+b{k}_{\phi }\right)}^{\frac{1}{2n+1}}{D}^{\frac{n+1}{2n+1}}}$$

(16)

Since the speed of a UGV traversing soft soil is very low, the air resistance factor is usually ignored. According to the principle of dynamic balance, the maximum speed that a UGV may achieve on soft soil is

$$v=\frac{p_e}{F_{rc}+W\sin \beta }$$

(17)

where $p_e$ is the rated power of the engine and $\beta$ is the terrain slope.

In the elastic mode, the resistance effect caused by the elastic deformation of the tire must be considered. If a part of the tire is flattened and its projection length is L, the contact stress of this part is

$$p_f+p_c=\frac{W}{bL}$$

(18)

According to Equation (11), the sinkage $z_0$ of the tire at this time is

$$Z_0={\left(\frac{p_i+p_c}{\frac{k_c}{b}+k_{\phi }}\right)}^{\frac{1}{n}}={\left(\frac{W}{L\left(k_c+b{k}_{\phi }\right)}\right)}^{\frac{1}{n}}$$

(19)

Substitute Equation (19) into Equation (12) to obtain:

$$F_{rc}=b\left(\frac{z_0^{n+1}}{n+1}\right)\left(\frac{k_c}{b}+k_{\phi }\right)={\left(\frac{W}{L}\right)}^{\frac{n+1}{n}}\left(n+1\right){\left(k_c+b{k}_{\phi }\right)}^{-\frac{1}{n}}$$

(20)

Similarly, after obtaining the resistance of the soil, the maximum speed that the vehicle may reach can be calculated according to Equation (17).

2.3. Mobility Cost Quantification

As described above, mobility prediction methods usually divide the terrain into GO and NOGO [19], where GO represents the area where the vehicle can traverse and NOGO is the opposite. The terrain points with slopes greater than the kinematic limits of the UGV are marked as NOGO, while soil areas with a CI less than VCI are marked as NOGO. Nevertheless, it is difficult for decision-makers to determine an optimal route for the UGV to travel from the GO areas. Consequently, it is necessary to quantify the mobility of the UGV in the GO areas. In the process of quantifying mobility, the influence of terramechanics and terrain parameters should also be considered.

However, the multi-factor fusion process cannot be simply weighted, so a fuzzy logic-based multi-factor fusion method is used. Fuzzy logic is a technique which uses the degree of truth instead of discrete values such as 0 or 1; it also uses “linguistic” variables such as “low”, “medium”, and “high”, as well as numerical variables for the calculations. The relationship between inputs and outputs is given by some simple statements rather than complex mathematical equations [28]. Let $A$ be a fuzzy subset. We can represent $A$ as

$$A=\frac{\mu }{y}$$

(21)

where $\mu$ is the degree of membership of $y$ in fuzzy subset $A$. The membership degree has a value between 0 and 1. For a fuzzy subset $A$, if $x\notin A$, then $\mu \left(x\right)=0$, if $x\in A$, then $\mu \left(x\right)=1$, if $x$ possibly in $A$, but not sure, then $\mu \left(x\right)\in \left(0, 1\right)$. The block diagram of the proposed fuzzy mobility cost quantification process is shown in Figure 6.

Combining the terramechanics and the fuzzy inference method, the terrain is then divided into NOGO, safe (high), low-risk (medium) and high-risk areas (low), and the mobility cost of each terrain point is converted into a distribution between [0, 1]: the closer the value is to 1, the worse the mobility. The rules need to be further analyzed in the context of specific vehicle and terrain characteristics, and this process is presented in Section 3.2 of this paper. For the determination of the membership function, in general, the steeper its shape, the higher the resolution of the system, the slower the change and the better the stability of the system [29]. Therefore, the membership functions corresponding to “low”, “medium” and “high” are Z-type, triangular type and S-type, respectively, as shown in Figure 7, and the corresponding expressions are listed in Equation (22).

$$\mu \left(x\right)=\left\{\begin{array}{cc}1,& x<x_0\\ {\left(\frac{x_2-x}{x_2-x_0}\right)}^k,& x_0\le x\le x_2\\ 0,& x_2<x\end{array} \Rightarrow \right.\mu \left(x\right)=\left\{\begin{array}{c}\begin{array}{cc}0,& x<x_1\\ \frac{x-x_1}{x_3-x_1},& x_1\le x\le x_3\end{array}\\ \begin{array}{cc}\frac{x_2-x}{x_5-x_3},& x_3<x\le x_5\\ 0,& x_5<x\end{array}\end{array}\right.\Rightarrow \mu \left(x\right)=\left\{\begin{array}{cc}0,& x<x_4\\ {\left(\frac{x_6-x}{x_6-x_4}\right)}^k,& x_4\le x\le x_6\\ 1,& x_6<x\end{array} \right.$$

(22)

For the fuzzy inference method, no consistent theory has been established to guide the construction of fuzzy relations, implying that the construction of fuzzy relations is diverse, the related fuzzy operations are different and the final inference results are not the same. A relatively simple and effective algorithm based on “Max–min” criterion is widely adopted in the practical use of fuzzy inference [30]. The “Max–min” criterion, also known as the pessimistic criterion, is one of the decision criteria for uncertainty-based decision making, and the basic attitude of this approach is conservative. It is implemented by finding the worst possible outcome of each option, i.e., the minimum value, and then selecting the option that provides the maximum payoff, i.e., the maximum value, among its worst outcomes. This guides one to maximize the minimum possible outcome. In order to avoid the risk of damage to the UGV during the off-road mission as much as possible, the “Max–min” criterion is, therefore, used in this case for defuzzification to obtain relatively conservative mobility prediction results.

2.4. Path Generation

In the global path planning method, it is necessary to build the topological relationships of the grids and plan the optimal route according to the topological relationships. In the traditional method, the topological relationships are usually established based on the distance between grids and the distribution of obstacles. Nevertheless, in large off-road terrains (>5 × 5 km) considering multiple influencing factors, traditional topological relationship-building methods cannot meet the requirement of generating a feasible, safe and effective global route. In this section, we propose a method to construct topological relations based on the above quantified results of mobility. Because of the large number of terrain points, the terrain points are sampled by the Latin hypercube sampling (LHS) method, and the sampled terrain points are used as the nodes for establishing topological relations. The topological relations of each node are then constructed based on the above quantified results of the mobility cost. The specific topological relationships are built based on the extended cost evaluation matrix and the heuristic cost evaluation matrix. With the topological relations, an improved A* algorithm is used for global path planning.

According to the results of the fuzzy process, the terrain is divided into NOGO area ($U_B$), low-risk area ($U_L$), high-risk area ($U_H$) and safe area ($U_S$). All the areas are coupled into the total environmental area $U$.

To build the topological relationship between each node, the sampling terrain points are first generated by Latin hypercube sampling, and the topological map $M=\left(P,E\right)$ is formed: P is the set of sampling nodes and E is the set of road sections between sampling points. The connectivity matrix ${\mathit{A}}_p$ is then used for connectivity evaluation, as shown in Figure 8a and the expression is

$${\mathit{A}}_p=\left[\begin{array}{ccc}{a}_{ii}& \cdots & a_{ji}\\ ⋮& \ddots & ⋮\\ a_{ij}& \cdots & a_{jj}\end{array}\right]$$

(23)

where ${\mathit{A}}_p$ is the connectivity matrix between the sampling nodes. If the straight line between two certain nodes $p_i$ and $p_j$ passes through $U_B$, then its connectivity cost $a_{ij}$ = 0, otherwise $a_{ij}$ = 1.

The extended cost evaluation matrix includes the extended low-risk cost matrix ${\mathit{L}}_p$, extended high-risk cost matrix ${\mathit{H}}_p$, extended distance cost matrix ${\mathit{D}}_p$ and extended safe-cost matrix ${\mathit{S}}_p$. For instance, the extended low-risk cost matrix ${\mathit{L}}_p$ consists of the low-risk cost $l_{ij}$ between the nodes in set P. The evaluation method is to take n uniformly distributed sampling nodes on the straight line between two nodes $p_i$ and $p_j$; the low-risk is evaluated by the cost accumulation of each child node in the low-risk area and the evaluation results are indicated as the low-risk cost matrix of the path between multiple nodes, as shown in Figure 8b; its expression is

$$l_{ij}=\left\{\begin{array}{cc}{\displaystyle {\displaystyle \sum }_{k=1}^n}{U}_L\left(k\right),& a_{ij}=1\\ \infty ,& a_{ij}=0\end{array}\right. \stackrel{ }{\Rightarrow } {\mathit{L}}_p=\left[\begin{array}{ccc}{l}_{ii}& \cdots & l_{ji}\\ ⋮& \ddots & ⋮\\ l_{ij}& \cdots & l_{jj}\end{array}\right]$$

(24)

Similarly, the extended high-risk cost matrix ${\mathit{H}}_p$ and the extended safe-cost matrix ${\mathit{S}}_p$ can be obtained, as shown in formulas:

$$h_{ij}=\left\{\begin{array}{cc}{\displaystyle {\displaystyle \sum }_{k=1}^n}{U}_H\left(k\right),& a_{ij}=1\\ \infty ,& a_{ij}=0\end{array}\right. \stackrel{ }{\Rightarrow } {\mathit{H}}_p=\left[\begin{array}{ccc}{l}_{ii}& \cdots & l_{ji}\\ ⋮& \ddots & ⋮\\ l_{ij}& \cdots & l_{jj}\end{array}\right]$$

(25)

$$s_{ij}=\left\{\begin{array}{cc}{\displaystyle {\displaystyle \sum }_{k=1}^n}{U}_S\left(k\right),& a_{ij}=1\\ \infty ,& a_{ij}=0\end{array}\right. \stackrel{ }{\Rightarrow } {\mathit{S}}_p=\left[\begin{array}{ccc}{s}_{ii}& \cdots & s_{ji}\\ ⋮& \ddots & ⋮\\ s_{ij}& \cdots & s_{jj}\end{array}\right]$$

(26)

The distances between the nodes also need to be considered in node expansion, so the extension distance cost matrix ${\mathit{D}}_p$ is composed of Euclidean distance $d_{ij}$ between the nodes in set P, as shown in Figure 8c, and its expression is

$$d_{ij}=\left\{\begin{array}{cc}\Vert p_i-p_j,\Vert & a_{ij}=1\\ \infty ,& a_{ij}=0\end{array}\right. \stackrel{ }{\Rightarrow } {\mathit{D}}_p=\left[\begin{array}{ccc}{d}_{ii}& \cdots & d_{ji}\\ ⋮& \ddots & ⋮\\ d_{ij}& \cdots & d_{jj}\end{array}\right]$$

(27)

To improve the efficiency of path planning and to optimize the expansion direction of nodes, it is necessary to consider the heuristic cost evaluation between the new expansion node $p_i$ and the goal node $p_G$, and search along the path direction with the least heuristic cost. Therefore, the heuristic cost evaluation matrix includes the heuristic barrier cost matrix $\dot{{\mathit{B}}_p}$, the heuristic distance cost matrix $\dot{{\mathit{D}}_p}$ and the heuristic safe cost matrix $\dot{{\mathit{S}}_p}$. $\dot{{\mathit{B}}_p}$ is composed of the heuristic barrier cost $b_{iG}$ between the node $p_i$ and the goal node $p_G$. It is evaluated using the cumulative sum of the cost of the NOGO areas and the cost of the risk area where each child node is located in the path between the sampling nodes and the goal node, as shown in Figure 8d, so the expression is $b_{iG}=\sum _{k=1}^n{U}_B\left(k\right)+\sum _{k=1}^n{U}_L\left(k\right)+\sum _{k=1}^n{U}_H\left(k\right)$, and $\dot{{\mathit{B}}_p}={\left[\begin{array}{ccc}{b}_{iG}& \dots & b_{jG}\end{array}\right]}^{\mathrm{T}}$.

The heuristic distance cost matrix $\dot{{\mathit{D}}_p}$ consists of the Euclidean distance $d_{iG}$ between the node $p_i$ and the goal node $p_G$; its expression is $d_{iG}=\Vert p_i-p_G\Vert$, and $\dot{{\mathit{D}}_p}={\left[\begin{array}{ccc}{d}_{iG}& \dots & d_{jG}\end{array}\right]}^{\mathrm{T}}$.

The heuristic safe cost matrix $\dot{ {\mathit{S}}_p}$ consists of the heuristic safe cost $s_{iG}$ between the node $p_i$ and the goal node $p_G$. It is evaluated by accumulating the cost of the safe area, where each child node is located in the path between the sampling nodes and the goal node, and its expression is $s_{iG}=\sum _{k=1}^n{U}_S\left(k\right)$ and $\dot{{\mathit{S}}_p}={\left[\begin{array}{ccc}{s}_{iG}& \dots & s_{jG}\end{array}\right]}^{\mathrm{T}}$.

The global path planning method based on the extended cost evaluation function and the heuristic cost evaluation function designed in this paper is

$$F\left(p_i\right)=G\left(p_i\right)+H\left(p_i\right)$$

(28)

where $F\left(p_i\right)$ is the cost evaluation function from the start node $p_S$ to the current node $p_i$ and $G\left(p_i\right)$ is the extended cost evaluation function from the start node $p_S$ to the current node $p_i$; $H\left(p_i\right)$ is the heuristic cost evaluation function from the current node $p_i$ to the goal node $p_G$.

The node extension method connects the current node $p_i$ to its surrounding nodes in the topology graph M. The node extension method is shown in Figure 9. The start node $p_S$ is extended to a feasible connected node $p_i$ according to the connectivity matrix. If the connection path $a_{Si}$ between $p_S$ and $p_i$ is feasible, a candidate node set Q is generated. The optimized node $p_{\min }$ with the lowest cost is searched in the Q set, expanded to the connected nodes around it, and then the new extended node is added to set Q, until it extends to the goal node $p_G$.

In the process of node extension, for the new extension nodes, the extension cost evaluation function is used to evaluate the extension cost from the parent node to the next node using weighted calculation, as shown in the formula

$$G\left(p_i\right)=G\left(p_{min}\right)+g_i$$

(29)

$$G_i=k_b·l_{ij}+k_h·h_{ij}+k_s·s_{ij}+k_d·d_{ij}$$

(30)

where $G\left(p_{min}\right)$ is the extension cost from the parent node to node $p_{min}$, and $G_i$ is the extension cost from node $p_{min}$ to $p_i$. $k_b$, $k_h$, $k_s$ and $k_d$ are low-risk weight coefficient, high-risk weight coefficient, safe-weight coefficient and distance weight coefficient, respectively. The heuristic function is shown in the formula:

$$H\left(p_i\right)=k_b·w_{ij}+k_d·d_{ij}+k_s·s_{ij}$$

(31)

where $H\left(p_i\right)$ is the heuristic cost evaluation function between the new extension node $p_i$ and the goal node $p_G$, $k_b$ is the barrier weight coefficient, and $w_{ij}$, $d_{ij}$ and $s_{ij}$ are the heuristic barrier cost, heuristic distance cost and heuristic safe cost of node $p_i$ in the heuristic cost evaluation matrix, respectively.

After the optimized path nodes are connected sequentially, the vehicle motion trajectory contains multi-segment polylines. Because the trajectory has hard inflection nodes, it cannot meet the requirements of vehicle kinematic constraints. In this paper, the clamped cubic B-spline smoothing method is used to optimize the trajectory of path planning; as shown in Figure 10, the black dotted line is the initial path from $p_{i-1}$ to $p_i$, while the green curve is the result of smoothing.

3. Simulation Experiments and Discussion

3.1. Experimental Setting

The UGV model used for the simulation experiment is shown in Figure 11. Its main parameters are given in

Table 1. Parameters of the UGV model.

Parameters	Gross Mass/kg	Clearance/m	Number of Axles	Number of Tires	Tire Load/N	Tire Width/m	Tire Radius/m	Engine Power/kW
Value	1242	0.45	2	4	3200	0.40	0.51	262

.

3.2. Mobility Prediction Results

According to Equations (8) to (10), and the UGV parameters given in Table 1, the MI index of the UGV is 67.9, so VCI₁ is 389.4 kN/m² and VCI₅₀ is 534.1 kN/m². The soil types in this terrain are mainly divided into sandy loam, clayey soil, lean clay and upland sandy, as illustrated in

Table 2. Soil characteristic parameters.

Soil Type	Moisture Content (%)	Deformation Index (n)	Cohesive Force Modulus ${\mathit{k}}_{\mathit{c}}$ (kN/mn+1)	Internal Friction Modulus ${\mathit{k}}_{\mathit{\phi }}$ (kN/mn+2)	Sinkage (m)	Rolling Resistance (kN)	Cone Index (kN/m2)
Sandy loam	26	0.30	2.79	141.11	0.51	18.84	223.40
Clayey soil	38	0.50	13.19	692.15	0.13	8.86	732.50
Lean clay	22	0.20	16.43	1724.69	0.02	3.11	1070.50
Upland sandy	51	1.10	74.60	2080.00	0.16	9.32	528.50

.

Recall that the 3D terrain model we created is a multi-dimensional model containing elevation data, soil type data (soil strength values, CI), slope values and land use data. Therefore, according to the comparison between the calculated VCI and the CI of each area, the area with a CI less than VCI₁ is defined as the NOGO area. According to the obtained slope map in Figure 5, combined with the limitations of the vehicle kinematic model, if the slope value exceeds the maximum pitch angle or roll angle of the UGV, it is then determined as a NOGO area. Additionally, based on land use data, the areas covered by water and forests are also defined as NOGO areas, while the rest are GO areas. Thus, a GO/NOGO map is finally generated, as shown in Figure 12a. In this terrain, the GO areas account for 71% and the NOGO areas account for 29%. Combined with the slope value and the soil parameters in Table 1, the maximum attainable speed of the UGV in each area can be obtained using Equation (18). As a result, the maximum attainable speed distribution map of the UGV can be obtained, as shown in Figure 12b. The maximum speed that the UGV can reach is 32.4 km/h, and the average speed is 9.5 km/h.

The quantification method of mobility results adopts the fuzzy inference described above; three states are defined for the three independent evaluation indexes, namely “low”, “medium” and “high”, to make the final result more accurate. These are combined with the above determination for NOGO areas, with four states being defined for the mobility cost, namely, “NOGO area”, “high-risk area”, “low-risk area” and “safe area”. Reviewing the previous section, the factors affecting mobility are summarized as the maximum attainable speed of the UGV, CI, terrain slope and height. When it comes to the fuzzy rule setting, this study references the expert rules for off-road vehicles in off-road terrain evaluation studies, such as the slope being defined as high, if in the range of 0°–8°, medium 5°–15°, and low 10°–30° [31]. Similarly, terrain height is defined as high if in the range of 0–300 m, medium 200–600 m, and low if greater than 500 m [32]. Under off-road traversing conditions, the vehicle speed is defined as low, if in the range of 0–15 km/h, medium 10–30 km/h, and high 20–40 km/h [28]. By comparing the calculated VCI₁ and VCI₅₀ to the CI, CI is defined as low if in the range of 0–500 kN/m², medium 400–700 kN/m², and high if greater than 700 kN/m² (See Appendix B for details of the specific rules setting). When determining the membership function, we normalize the input data, then use the S-type membership function, triangular membership function and Z-type membership function to express the fuzzy rules, as shown in Figure 13.

The defuzzification process uses the “Max–min” criterion proposed in Section 2.3. After defuzzification, the global terrain is divided into four categories, as shown in Figure 14a. The mobility cost evaluation result is also expressed as the distribution of [0, 1], and the closer to 1, the worse the mobility, as shown in Figure 14b.

3.3. Path Generation Results

After obtaining the mobility cost results, the improved A* algorithm described above is used for global path planning (min-mobility cost route). We also compare two other algorithms, the traditional A* algorithm for searching the shortest route (min-distance route), and the A* algorithm with the consideration of the terrain profile (min-slope route) [33], as shown in Figure 15.

The final experimental result shows that the three path planning strategies generate three global routes, namely, MDR (min-distance route), MSR (min-slope route) and MMCR (min-mobility cost route), as given in Figure 16. To intuitively evaluate the routes generated at different costs, each route node is stored and statistically analyzed. All routes are normalized and segmented in

Table 3. Summary of the three cost functions considered in this research.

Route	Distance (m)	Average Height (m)	Average Slope (°)	Average Speed (km/h)	Average Mobility Cost	Traversal Time (h)
MDR	6724	628.8	20.9	5.3	0.48	1.47
MSR	11,670	593.2	8.5	11.5	0.40	1.01
MMCR	8306	551.1	12.4	13.8	0.32	0.62

to show the average indicators of the UGV in different route segments, and the detailed comparison of the data is shown in Figure 17. After ignoring the NOGO areas, MDR adopts the planning strategy of the shortest route. The comparison of experimental results also proves that this route is the shortest. However, only the distance factor is considered, and the influence of slope and soil strength on vehicle mobility are ignored. The route has the highest average slope (20.9°), the highest height (628.8 m), the lowest average speed (5.3 km/h) and the longest traversing time (1.47 h). Combined with the abovementioned fuzzy rules, the MDR is classified as “high risk”. MSR adopts the planning strategy of the slowest route. It is obvious that the route has the smallest average slope (8.5°); because only the influence of slope is considered in the cost function, the slowest route leads to the longest path length (11,670 m). Simultaneously, when searching for the slowest route, it inevitably selects areas with low soil strength, resulting in an increase in traversing resistance and reduced average speed (11.5 km/h). Combined with fuzzy rules, MSR is classified as “low risk”. Conversely, the minimum mobility cost planning strategy proposed in this paper comprehensively considers the influence of terrain. On this basis, the terrain points with higher speed are searched as the route nodes. Therefore, this route has the highest average speed (13.8 km/h) and the least traversing time (0.62 h), and its fuzzy result is “low risk”.

In addition, the computer execution time for the three methods in Matlab shows the longest to be MSR, which takes 75 min, followed by MDR, which takes 44 min, whereas MCCR combines the characteristics of the sampling algorithm and, therefore, takes the least time, less than 4 min; these experiments were run on a general-purpose computer (Intel Core i7 3 GHz,16 MB RAM).

In order to further prove the effectiveness of the algorithm in this paper, the generated routes are extracted and compared with regions A and B, with a large difference in fitting degree, as shown in Figure 18.

Region A has two types of soil, lean clay (CI = 1070.5 kN/m², 77%) and upland sandy (CI = 528.5 kN/m², 23%), while the average slope value in this area is 15.6° with a standard deviation of 4.2. The comparison results of the three routes are shown in

Table 4. Comparison of the parameters of the three routes in region A.

Route	Distance (m)	Average Height (m)	Average Slope (°)	Average Speed (km/h)	Traversal Time (h)
MDR	1465	518.7	16.1	4.8	0.31
MSR	2734	563.4	8.8	8.9	0.30
MMCR	2267	500.2	9.1	13.3	0.17

.

In region A, MDR only considers distance cost, resulting in traversing areas with low soil strengths and high slopes, so it has the lowest average speed (4.8 km/h) and the longest traversal time. In the MSR search process, only the slope cost is considered, and the result satisfies the requirement of searching for the slowest route; however, this process also produces the longest route and a route with low soil strength. MMCR proposed in this paper, and which considers the comprehensive factors (mobility cost), avoids the soft soil areas well, so that it can search along the lean clay area; this route has the best mobility, with the highest average speed (13.3 km/h) and the lowest time (0.17 h). It also has a smoother route by using the B-spline curve.

In the case of the large terrain relief in region B (mean slope value of 21.1°, standard deviation of 6.7), the soil contains lean clay (CI = 1070.5 kN/m², 56%), upland sandy (CI = 528.5 kN/m², 27%) and clayey soil (CI = 732.5 kN/m², 17%); the results of the three planning strategies are shown in

Table 5. Comparison of the parameters of the three routes in region B.

Route	Distance (m)	Average Height (m)	Average Slope (°)	Average Speed (km/h)	Traversal Time (h)
MDR	2890	708.7	18.3	3.9	0.74
MSR	4518	686.4	10.8	7.2	0.63
MMCR	3321	644.2	13.3	12.1	0.28

.

In region B, MDR still searches for the shortest route within the upland sandy; its average speed is further limited (3.9 km/h) due to the increase in terrain undulation; although its route distance is the shortest, it takes the longest time (0.74 h). Due to the large change of slope, MSR only pursues a gentler route and also increases traversing distance. Meanwhile, part of the route enters the area of upland sandy and clayey soil with relatively low soil strength, which leads to a lower average speed. MMCR can avoid this situation well on the basis of the mobility cost quantification method, while being able to find a faster route on the high strength soils (average speed 12.1 km/h, traversing time 0.28 h).

4. Conclusions and Future Work

In this paper, a mobility prediction-based global path planning method for UGVs is proposed. The experimental results show that the proposed method, MMCR, can effectively generate a more feasible, efficient and safer route when considering the mobility cost. It is worth noting that the current path-planning technologies usually rely on terrain modeling by the UGV’s own sensory devices, which can build a more detailed terrain model, but the data storage volume required is huge and the data parsing method is relatively complex. Sensing the surrounding environment in the actual unknown off-road terrain is likely to cause the risk of failure of the unmanned equipment. This paper proposed an effective method to reconstruct 3D off-road terrain by fusing multi-source remote sensing data based on the Kriging interpolation method. Also based on terramechanics and fuzzy inference, the effect of soil deformation on the UGV is considered to reduce the risk of obstruction. In summary, this paper provides an effective and safe reference method for UGVs to perform off-road tasks. The method is not only applicable to the terrain selected in this paper, but also applicable to other terrains, so it is generalized, and has a certain value in engineering applications.

Future work will focus on several limitations of this study. (1) The limitation of this study is that the attitude change of the vehicle during the driving process is not taken into account; therefore, to obtain more accurate and more informative mobility prediction results, discrete element or finite element techniques will be considered for modeling the soil and coupling the soil model with the vehicle model for simulation. Although several techniques for coupled simulation have been proposed [34,35], these approaches are, nevertheless, time-consuming when predicting mobility over large areas. To obtain mobility prediction results efficiently, it is necessary to combine these methods with machine learning algorithms. (2) On the other hand, there are perceptual errors in satellite remote sensing data and interpolation errors generated using terrain interpolation methods, as well as spatial variability in physical soil properties (water content, season, etc.); all of these may lead to uncertain results of mobility prediction. Therefore, how to reduce or quantify the uncertainty of mobility prediction will be another future research direction.