A New Challenge: Path Planning for Autonomous Truck of Open-Pit Mines in The Last Transport Section

Abstract

During the operation of open-pit mining, the loading position of a haulage truck often changes, bringing a new challenge concerning how to plan an optimal truck transportation path considering the terrain factors. This paper proposes a path planning method based on a high-precision digital map. It contains two parts: (1) constructing a high-precision digital map of the cutting zone and (2) planning the optimal path based on the modified Hybrid A* algorithm. Firstly, we process the high-precision map based on different terrain feature factors to generate the obstacle cost map and surface roughness cost map of the cutting zone. Then, we fuse the two cost maps to generate the final cost map for path planning. Finally, we incorporate the contact cost between tire and ground to improve the node extension and path smoothing part of the Hybrid A* algorithm and further enhance the algorithm’s capability of avoiding the roughness. We use real elevation data with different terrain resolutions to perform random tests and the results show that, compared with the path without considering the terrain factors, the total transportation cost of the optimal path is reduced by 10%–20%. Moreover, the methods demonstrate robustness.

1. Introduction

Cost control has always been the top priority in the mining industry, and how to control the cost is one of the critical points for mining companies to make profits. In the open-pit mines, the transportation cost accounts for about 50% of the total mining cost [1,2,3]. In the absence of mining technology innovation, how to improve transportation efficiency and reduce extra transportation costs are crucial to improve the profits of mining enterprises. The rise of driverless technology makes the problem have a new solution. At present, unmanned ground vehicles (UGVs) have been widely used in port or warehouse automation logistics, disaster relief, military, and other industrial and military fields [4,5,6,7,8].

In the mining field, autonomous navigation of the load-haul-dump (LHD) machine in underground mines has been successfully applied [9,10]. However, it is still a challenge to apply autonomous technology in open-pit mining. The open-pit mining area is an ideal environment for unmanned vehicles due to its closed environment and few external interference factors. Komatsu and other mining equipment manufacturing giants have carried out unmanned transportation experiments in a Chilean copper mine since 2008 [11]. Their results show that the transportation capacity is increased by 40% and the cost has reduced by 15% in the autonomous haulage system (AHS). Another unmanned transportation experiment in Australian mining also shows that unmanned transportation has more considerable economic significance compared with traditional transportation [12]. Nonetheless, Ayres [13] and Juliana Parreira [14] both pointed out that the AHS may have problems such as unable to respond to ground pits or rocks promptly, resulting in extra fuel consumption and tire wear. Also, some studies [3,15] have concluded that, when haulage trucks run on poor road conditions, the vehicle cost will be increased. Kansake [16], through the simulation experiment, concluded that under the general mining road conditions, the dynamic load of the tire would be 60% higher than the static load, which will aggravate the tire wear. For the existing AHS, only the large obstacles on the ground are considered in path planning. In contrast, the tiny obstacles are ignored, like rocks or terrain undulation on the surface, resulting in additional tire wear. Thus, the truck route should be planned on smoother ground.

The path planning of UGV is divided into route planning, trajectory planning, and motion planning [17]. Route planning is point-to-point planning, which regards the vehicle as a particle and does not consider other constraints, so it requires less map accuracy. Trajectory planning is based on route planning, with more constraints of vehicles, such Cering the most constraints, such as vehicle control strategy. Choi et al. [18] used the map from Google Earth to optimize the route of dump trucks in open-pit mines. Google Earth only provides a low-resolution surface elevation model, but the realization of trajectory planning requires more accurate surface topographic maps.

The digital surface model (DSM) is a model that contains the surface elevation information, and it is generally used for the reconstruction of the terrain model as its high precision. There are several methods to build a DSM: traditional digital mapping [19], unmanned aerial vehicle (UAV) lidar mapping [20,21], oblique photogrammetry [22], etc. Among these methods, the oblique photogrammetry has been widely used in urban surveying and mapping [23], disaster monitoring [24], and the open-pit mining industry [25,26,27] because of its high precision and convenient operation. Based on DSM, trajectory planning can be carried out. Jurgen Zoto et al. [28] planned a coverage path for UGV of vineyard based on DSM. Kai Zhang et al. [29] used lidar to measure the surface elevation point cloud and constructed the ground map to carry out the path planning of cross-country vehicles. These DSM-based path planning studies only separated obstacle information but did not consider surface roughness information in the planning map.

The transportation environment of open-pit mines is unstructured, and the Hybrid A* algorithm (HA* algorithm) is a heuristic algorithm widely used in real-time vehicle path planning for unstructured environments [30,31,32,33,34]. This algorithm adds the non-holonomic constraints into the heuristic functions and plans the path in the Voronoi field map. Nevertheless, in the planning process, the regular Hybrid A* algorithm only considers the large obstacles. However, it ignores the little pit or roughness of the surface, which may cause an additional cost for large vehicles.

There are a lot of open spaces in the cutting zone of open-pit mines without any obstacles. In this scenario, the path planner considers the road surface as a passable area without obstacles, thus planning a path directly to the loading position which may result in extra tire cost because of the ground roughness. In order to reduce the extra cost, this paper proposes a high-precision planning map construction method based on DSM and an improved Hybrid A* algorithm. This planning map contains both surface obstacles and roughness information. Firstly, the DSM is divided into an obstacle point set and surface roughness point set with different terrain features. Secondly, the obstacle point set is calculated as obstacle cost map (OCM) according to generalized Voronoi diagram (GVD). Then the sliding window standard deviation method is used to process the surface roughness point set and generate the roughness cost map (RCM) without obstacle. Finally, the two maps are integrated into the cost map (C-MAP) for path planning, which contains the information of obstacles and roughness on the road. With this C-MAP, we modified the Hybrid A* algorithm by adding the contact cost between the tire and ground into the heuristic function, so that the path can be planned on the smooth road surface. The main contributions of this paper are as follows:

• We propose a truck path planning method based on high-precision digital map for the final transportation section of open-pit mines, which can effectively reduce the driving cost of trucks.
• We design a planning map that contains obstacle information and roughness information on the road surface and its construction method.
• We improve the Hybrid A* algorithm to incorporate the roughness information into the cost estimation function of node extension.

The structure of this paper is as follows: Section 2 describes the modified Hybrid A* algorithm in detail, including planning map building method and algorithm improvement. In Section 3, we evaluate the algorithm performance in different real open-pit cutting zone scenes and analyze the algorithm robustness through a thousand random tests with different terrain resolution maps. Conclusions and future work are presented in Section 4 and Section 5.

2. Materials and Methods

To avoid the rough area on the transportation road, it is necessary to make a planning map containing the roughness information and add it to the estimation function of the Hybrid A* algorithm. The method of constructing a planning map with roughness information is as follows:

• DSM data acquisition and generation. The raw elevation data of cutting zone is obtained by some topographic mapping methods (oblique photogrammetry, lidar scanning, etc.), and the corresponding DSM is generated.
• Constructing the obstacle cost map of the cutting zone. Firstly, the obstacle point set in DSM is filtered by setting different terrain factors thresholds. Then the obstacle binary map is constructed by the elevation scanning method, and finally, the obstacle cost map is made from GVD.
• Constructing the roughness cost map of the cutting zone. Firstly, removing the obstacle point set in DSM, filtering the remaining point set with another terrain threshold. Then, using the standard deviation of moving window to traverse calculation, its value is taken as the roughness value of all coordinate points within the window range. Finally, normalizing the roughness matrix to generate the roughness cost Map of the cutting zone.
• Generating the C-MAP. Overlaying and normalizing the value matrix of obstacle cost map and roughness cost map, get the C-MAP, which is used in path planning.

The modified parts of the Hybrid A* algorithm are as follows:

• Modify the g-value estimation function. The contact cost between the tire and the ground surface is added into the g-value estimation function. Consequently, the step length a dynamic value instead of a static value in the part of the node extension.
• Modify the conjugate gradient object function. The original Voronoi term is replaced by a new roughness term to correct the path points passing through rough area.

The primary process of map construction and modified Hybrid A* algorithm is shown in Figure 1.

2.1. Raw Data Collection and DSM Generation

There are many topographic surveys to collect DSM original image data. In this paper, UAV oblique photogrammetry is used for data collection. After that, using terrain processing software to extract the elevation points from the original image and generating the DSM, as shown in Figure 2. Save the elevation points information of DSM in the following matrix:

$$E=\left[\begin{array}{ccc}{e}_{11}& \cdots & e_{1n}\\ ⋮& \ddots & ⋮\\ e_{m1}& \cdots & e_{mn}\end{array}\right],$$

(1)

where $e_{ij}$ is the corresponding coordinate of elevation points.

2.2. Obstacle Cost Map Construction

2.2.1. Elevation Feature Points Extraction

In DSM, the pits or fallen rocks on the ground are expressed in the form of elevation difference. Wang et al. [35] proposed a method to extract the elevation feature points to represent the terrain changes. The feature points were extracted by the four directions (vertical, horizontal, 45° left, 45° right) scanning method based on the section scanning. The key is the selection of terrain thresholds to distinguish different terrains, such as obstacles or little pit. Moreover, for high-precision DSM, Wang pointed out that the equal-interval sampling method can improve the efficiency of subsequent calculations, and the model error is acceptable. Therefore, the original DSM with 0.032 m accuracy is reprocessed to 0.096 m accuracy by the equal-intervals sampling method in this paper. Referring [16], we set $T=0.3\mathrm{m}$ as the terrain threshold of impassable elevation difference, and $\theta ={15}^{°}$ as the slope threshold of the maximum passable slope. The steps of extracting elevation feature points are as follows (take the vertical scanning as an example):

• Reforming the elevation matrix $E$ according to the different scanning directions. The elements of each column in $E$ is put into $E_i$, as shown in Equation (2).

  $$E_i=\left[\begin{array}{c}{e}_{1i}\\ ⋮\\ e_{mi}\end{array}\right],E=\left(E_1,E_2\dots E_n\right),i=1,2\dots n$$

  (2)
• Scanning each data point sequentially in $E_i$: set a searching group with $e_{j-1},e_j,e_{j+1},\left(j=2,3\dots m-1\right)$, the step length is one. If the point is the first or last point of $E_i$, or local extremum of the group, add the point $p$ to the candidate feature points set $F_p$, where $p$ contains the elevation value and coordinate information. The extracted feature points are shown in Figure 3a.
• Filtering the points of $F_p$ by terrain threshold $T$ which is defined by the characteristics of the cutting zone. Then determine whether $p_k$ and $p_{k+1}$ are adjacent in $F_p$. If they are adjacent, judge the absolute value of two-point difference, whether it’s over the threshold $T$. If the value greater than $T$, $p_{k+1}$ will be reserved in $F_p$. Otherwise, $p_{k+1}$ will be removed from $F_p$. The original section scanning line and the final feature point set $F_p$ are shown in Figure 3b.

2.2.2. Binary Map Construction

Binary map construction is the preprocessing work for Obstacle Cost Map construction. Firstly, Compare whether the slope angle $\beta$ of adjacent points $p_k\left(x_{p_k},y_{p_k}\right)$ and $p_{k+1}\left(x_{p_{k+1}},y_{p_{k+1}}\right)$ in $F_p$ greater than $\theta$ sequentially. If so, all the corresponding coordinate points between $p_k$ and $p_{k+1}$ are assigned as 1, otherwise, it is 0, i.e., the value of one is represented obstacle. As shown in Equations (3) and (4). Figure 4 is the schematic diagram of the four-direction scanning binary map.

$$\beta =ta{n}^{-1}\left(\frac{y_{p_{k+1}}-y_{p_k}}{x_{p_k}-x_{p_k}}\right),$$

(3)

$$p_j=\{\begin{array}{l}1,\beta \ge \theta \\ 0,\beta <\theta \end{array},j=1,2\dots m,$$

(4)

where $p_j\in E_i$.

Secondly, fuse four scanning binary maps to form the final binary map in a simple mixture method: Traverse all points of the four maps. If the accumulated value of the same coordinates in four maps is greater than or equal to 2, then the value of the coordinate point is set to 1. Otherwise, it is 0. The final obstacle binary map is shown in Figure 5.

2.2.3. Obstacle Cost Map Generation

For the heuristic algorithm, the map cost value directly affects the accuracy of the estimation function, thus affecting the quality of the final path. Voronoi diagram is a geometry-based route planning method for robot path planning [36,37], which divides the space into several regional units according to the nearest neighbor attributes of the elements in the site collection. Dolgov [30], the inventor of the Hybrid A* algorithm, used the generalized Voronoi diagram to generate the planning map, which made the planning path can pass through the narrow area and without collision. The generation steps for the obstacle cost map are as follows:

• Building the two-dimensional Voronoi diagram from the binary map. The 1-value coordinate points in the binary map of obstacles are regarded as the scatter points. Then, calculate the 2-D Voronoi diagram of all scatter points. The red dot is the Voronoi station, and the blue line is the Voronoi edge, as shown in Figure 6a.
• Cutting off the Voronoi edges which passing through the obstacle polygon, and the remaining edges are the edges of GVD. The black square is the obstacle polygon, and the red line is the GVD edges, as shown in Figure 6b.
• Converting the GVD into grid form by sampling method, as shown in Figure 6c. Then, calculating each grid cost by Equation (5) to construct the obstacle cost map. The grid cost reaches its maximum within obstacles, and the color is close to black. It reaches its minimum on the GVD edges. The color is close to white, as shown in Figure 6d.

$$gridCost\left(x,y\right)=\{\begin{array}{l}(\frac{\alpha }{\alpha +d_o(x,y)})(\frac{d_v(x,y)}{d_o(x,y)+d_v(x,y)})(\frac{{(d_o-d_o^{max})}^2}{{(d_o^{max})}^2}),d_o<d_o^{max}\\ 0,others\end{array}$$

(5)

where $gridCost\left(x,y\right)$ is the cost of each grid. It is continuously distributed in the range of 0 to 1, only takes 1 in the obstacles, and takes 0 on the Voronoi edge. $\alpha$ is the parameter to control the descent rate of cost, $\alpha >0$. $d_o^{max}$ is the maximum range of potential value distribution, $d_o^{max}>0$. For the grid points that exceed the $d_o^{max}$, the grid cost is 0. $d_o$ is the distance from the grid point to the nearest obstacle edge. $d_v$ is the distance from the grid point to the nearest Voronoi edge.

Using this method to process the obstacle binary map of the cutting zone, we get the corresponding obstacle cost map, as shown in Figure 7. However, due to the high grid precision and the maximum range of potential value distribution, only the area at the edge of the obstacle (black area) has an available cost (gray area). For the rest of the open space area (far away from the obstacle), the grid cost is 0 (white area). Therefore, to make the algorithm consider the cost in the open area, we construct the roughness cost map to complement this area.

2.3. Roughness Cost Map Construction

A large number of blank areas exist in the obstacle cost map of the cutting zone, in which there are many possible adverse road conditions. If it is not included in the heuristic function of the planning algorithm, it will result in a deviation of the planning path. Generally, the standard deviation of the sliding window can be used to estimate the ground surface roughness [38]. The smaller the standard deviation value, the smoother the road surface. We use this method to construct the roughness cost map:

• Removing the coordinate of obstacles in the original elevation point set $E$, and obtained a new point set $E_r$, which only keeps the elevation value of the open space area.
• Using the section scanning method (threshold $T=0,\theta =5^{°}$) to process the open space area, deleting the elevation value of a point below the threshold from $E_r$.
• In $E_r$, taking a window and using the sliding standard deviation method to traverse $E_r$, and calculating the standard deviation of the window, then assigning its value back to each point in this window.
• Normalizing all the values of points in $E_r$ with a range of 0 to 1, which is the same as the obstacle cost range.
• Through the above steps, we can get the rough ground cost map of the cutting zone, as shown in Figure 8a.

2.4. C-MAP Generating

After constructing OCM and RCM with the same value range at $cost=\left[0,1\right]$, the C-MAP can be constructed as follows:

• Removing the obstacle points (cost = 1) from OCM, and adding the remaining OCM to RCM.
• Normalizing the obtained matrix.
• Backfilling obstacle points of OCM.

The constructed C-MAP still maintain the value range at $cost=\left[0,1\right]$, and it contains both the obstacle information (black part of Figure 8b) and surface roughness information in the open space area (other color parts of Figure 8b), in which the planning algorithm can take roughness cost into account in the estimation function.

2.5. Hybrid A* Algorithm Improvement

This section introduces the planning process of the modified Hybrid A* algorithm, and how to optimize it to fit the scene described in this article. The Hybrid A* algorithm is an improved A* path planning heuristic algorithm proposed by Dolkov [30], which adds the vehicle non-holonomic constraints to the estimation function. Based on the two-dimensional grid $\left(x,y,\theta \right)$ extension of A* algorithm, it adds a three-dimensional search, considers the three-dimensional state $\left(x,y,\theta ,\lambda \right)$ of the vehicle. Where $x$ and $y$ are coordinate information, $\theta$ is direction information, and $\lambda$ is vehicle forward and backward status information. Nevertheless, the path produced by the Hybrid A* algorithm has local bending, so Dolkov uses the conjugate gradient (CG) method to smooth the path.

We add the cost of tire contact with the road surface into the estimation function of the Hybrid A* algorithm instead of only taking the cost of the center point as the estimation. Also, we replace the Voronoi term with the Roughness term to enhance the path’s avoiding ability for roughness surface.

2.5.1. Node Extension Improvement

The node extension of Hybrid A* is guided by two heuristic methods. One is the 2D heuristic method in the traditional A* algorithm, which considers the information of obstacles in the map, but ignores the physical attributes of the vehicle itself. Its node extension only contains coordinate information and direction information, as shown in Figure 9a. Another is the non-holonomic-without-obstacles heuristic method. This method considers the vehicle non-holonomic constraints and assumes that there are no obstacles in the map, assumes the goal state of $\left(x_g,y_g,{\theta }_g\right)=\left(0,0,0\right)$, computes the shortest path to the goal from every point $\left(x,y,\theta \right)$. The schematic extension diagram is shown in Figure 9b. The second method uses the maximum of the two methods cost as the heuristic cost. Also, in order to save computing time, the algorithm uses Reed-Shepp (RS) curve to judge the current node and the goal node periodically. If there is no obstacle between the two points, an obstacle-free RS path will be generated directly between the two points. We set a few parameters, such as number of motion primitives (NMP), motion length (ML), and extension interval (EI), to control the node extension. NMP determines how many next nodes will be calculated based on the current node, by way of example, in Figure 9b, the NMP is 5, so there are five forward and five backward points. ML is the length of each extension determined by the minimum turning radius. EI determines how many loops to perform RS curve detector.

When the Hybrid A* algorithm estimates the moving cost, it usually takes the grid value of the center point as the estimation cost. We improve the g-value estimation function of the Hybrid A* algorithm, add a new term to estimate the cost of tire contact with the ground when the vehicle is moving. The improved node extension diagram is shown in Figure 10, the yellow area is the coverage area of the tire when extending, and the occupancy grid value within the range is taken as the extension cost.

2.5.2. Estimation Function Improvement

The Hybrid A* algorithm uses the Equation (6) to estimate the cost of each node extension.

$$f(x)=g(x)+h(x),$$

(6)

where the first term $g\left(x\right)$ is the sum of the estimated cost from the start to the current node, the second term $h\left(x\right)$ is the heuristic value from the current node to the goal, and the sum of the two terms is the evaluation function of the whole algorithm. This paper only improves the $g\left(x\right)$, and use the two heuristic functions of the original algorithm to calculate the $h\left(x\right)$. The calculation of $g\left(x\right)$ is realized by Equation (7).

$$g(x)=DirectionCost\times (g(x_{parent})+length(x_{parent},x)+TireCost(x_{parent},x))+\lambda \cdot S\mathrm{witch}Cost$$

(7)

where $DirectionCost$ contains ForwardCost (FC) and ReverseCost (RC), penalizing the motion cost according to a different direction. $g\left(x_{parent}\right)$ represents the sum cost from the start point to the parent node of the current node. $length\left(x_{parent},x\right)$ is the length of the curved path between the parent node and the current node.$S\mathrm{witch}Cost$ is the direction switch cost to penalize the motion state, when the direction of the vehicle is changed $\lambda =1$, otherwise $\lambda =0$. We add a new cost term, $TireCost\left(x_{parent},x\right),$ based on the original g-value estimation function, to estimate the cost of tire contact with the surface. This term represents the sum cost of the occupied grid values by the tire line from the current node to its parent node:

$$TireCost(x_{parent},x)={\displaystyle \sum _{i=x_{parent}}^x(CostMa{p}_i)}+{\displaystyle \sum _{j=x_{parent}}^x(CostMa{p}_j)},$$

(8)

where the first term is the accumulative cost of the left tire, and the second is the right.

2.5.3. Conjugate Gradient Method Improvement

There are still some unnecessary turns in the path planned by the Hybrid A* algorithm. Dolkov used the conjugate gradient (CG) smooth method to optimize the path. The objective function is:

$$P={\omega }_{\rho }{\displaystyle \sum _{i=1}^N\rho v(x_i,y_i)+{\omega }_o{\displaystyle \sum _{i=1}^N{\sigma }_o(\left|x_i-o_i\right|-d_{max})+{\omega }_{\kappa }{\displaystyle \sum _{i=1}^{N-1}{\sigma }_{\kappa }(\frac{\mathsf{\Delta }{\mathsf{\Phi }}_i}{\left|\mathsf{\Delta }{x}_i\right|}-{\kappa }_{max})+{\omega }_s{\displaystyle \sum _{i=1}^{N-1}{(\mathsf{\Delta }{x}_{i+1}-\mathsf{\Delta }{x}_i)}^2}}}},$$

(9)

where the first term is Voronoi term to guide the path away from an obstacle in both narrow and wide passages, the second term is Obstacle term to penalize the collisions, the third term is curvature term to ensure the path curvature satisfies vehicle non-holonomic constraints, and the fourth term is the smoothness term to make the path smoother.

In the cutting zone scene, there is a large open space area instead of narrow passages. Thus, we replace the Voronoi term with roughness term in CG, and the function is:

$$P_r={\omega }_r{\displaystyle \sum _{i=1}^{N-1}{\sigma }_r\left({\displaystyle \sum _{n=1}^4{\omega }_n(\left|x_{in}-r_{in}\right|)}\right)}$$

(10)

The derivative is:

$$\frac{\partial P_r}{\partial x_i}=2{\omega }_r\cdot \left({\displaystyle \sum _{n=1}^4{\omega }_n(\left|x_{in}-r_{in}\right|)}\right)\cdot \left({\displaystyle \sum _{n=1}^4{\omega }_n\frac{x_{in}-r_{in}}{\left|x_{in}-r_{in}\right|}}\right),$$

(11)

$${\omega }_n=1-\left(r_{nmax}/{\displaystyle \sum _{n=1}^4{r}_{nmax}}\right),$$

(12)

where ${\omega }_r$ is the weight of gradient; ${\sigma }_r$ is a simple quadratic penalties function; $x_{in}$ are the locations of four tire contact with the surface which corresponding to each path point; $r_{in}$ are the locations of the maximum value roughness points nearest to the $x_{in}$; $r_{nmax}$ is the distance between $r_{in}$ and $x_{in}$; ${\omega }_n$ is the weight of the gradient.

This term penalizes the contact between the tire and the nearby rough ground, which increases the ability to deviate from rough terrain during path point optimization.

2.5.4. The Modified Hybrid A* Algorithm in Detail

The modified Hybrid A* algorithm considering terrain roughness is improved on the regular Hybrid A* algorithm, which has the same essential parts. Before planning starts, it is necessary to set vehicle parameters, including length, width, tire information, and minimum turning radius, etc. Generally, tire information includes information on wheelbase and tire size. After that, the C-MAP is input, which includes both obstacle and roughness information to the path planner, and then the essential planning parameters are set, e.g., FC, RC, EI, etc. The pseudo code for Modified Hybrid A* algorithm (Algorithm 1) is as follows.

Algorithm 1 Modified Hybrid A*
1:	Set Vehicle properties:
	$VP\leftarrow \left(Length,Width,Tireinformation,Minimumturningradius\right)$
2:	Set planner: $planner\leftarrow \left(CostMap,Planningparameters\left(ML,NMP,FC,RC,SC,EI\right)\right)$
3:	Input: $startNode\left(x_s,y_s,{\varphi }_s\right),goalNode\left(x_g,y_g,{\varphi }_g\right)$
4:	Procedure Hybrid A* planner$\left(startNode,goalNode,VP,planner\right)$
5:	Node extension
6:	Validate the node
7:	RS curve check
8:	return Path Points
9:	Processing Path Points by ConjugateGradientSmoother$\left(CostMap,PathPoints\right)$
10:	return the Path Points after CG optimization
11:	end procedure
12:	Procedure gScore$\left(n_{current},n_{next},planner,VP\right)$
13:	$gScor{e}_{parent}\leftarrow gScore$ of $n_{current}$
14:	$Lengt{h}_{current\to next}\leftarrow ML$
15:	Generate the tire path besides the extension line according to $VP$
16:	Calculate the occupied grid coordinates of the tire path
17:	Delete duplicate grid coordinates
18:	$TireCost$ is the total cost of the occupied gird
19:	if the direction of $n_{next}$ $\ne$ the direction of $n_{current}$ then
20:	$gScore=gScor{e}_{parent}+Lengt{h}_{current\to next}+TireCost+SC$
21:	else
22:	$gScore=gScor{e}_{parent}+Lengt{h}_{current\to next}+TireCost$
23:	end if
24:	end procedure

The procedure gScore is an improved g-value estimation function. In each node extension step, it accumulates the total cost of the occupied grid of the corresponding tire path which is calculated from vehicle parameters. This allows Hybrid A* algorithm can incorporate surface roughness into motion cost estimation. Steps 19 to 23 mean that only when the motion status changes, in other words, when the motion direction changes from forward to backward, the SwitchCost will penalize the motion cost.

3. Performance Analysis of Modified Algorithm

To verify the effectiveness and robustness of the modified algorithm, we experiment with and analyze its performance in two different scenes of the Chengmenshan copper mine in Jiangxi. The original image data of the cutting zone are obtained by quad-rotor UAV oblique photogrammetry (Figure 11), and the corresponding C-MAPs of two scenes are constructed by the method proposed in this paper (Figure 12). In these two scenes, we set three different start positions and a fixed position of truck for path planning and verification. Path length, accumulative path cost, and path total score are used for the path evaluation.

Table 1. Cartographic information.

Camera Pixels of UAV	DSM Original Accuracy	DSM Accuracy after Sampling	C-MAP Unit Size
20 megapixel × 5	0.032 m	0.096 m	0.1 m × 0.1 m

is the cartographic information, and

Table 2. Simulated haulage truck parameters.

Length	Width	Wheelbase	Minimum Turning Radius	Tire Size
8.7 m	4.525 m	3.75 m	7.2 m	18.00-33-32 PR

is the parameters of the simulated haulage truck.

Table 3. Hybrid A* algorithm parameters.

Number of Motion Primitives	Motion Length	Forward Cost	Reverse Cost	Switching Cost	Extension Interval
5	72	1	5	100	30

shows the main parameters of the Hybrid A* algorithm. The parameters of conjugate gradient smoother are set to ${\omega }_r={\omega }_o=0.002,{\omega }_{\kappa }={\omega }_s=0.2$.

3.1. Analysis of Algorithm Capability in Different Scenes

According to the actual mining haulage route, we demarcate three initial starting positions in the green square frame, and a loading position in the red square frame, as shown in Figure 13. Then we verified the algorithm by the following two parts:

• Testing the recognition and avoidance ability of modified Hybrid A* algorithm for rough terrain. Firstly, plan a path in the obstacle cost map without roughness information. Then, after a virtual rough area is artificially added to the map and plan the path again. The result is shown in Figure 14.
• Evaluating the path planned by the modified Hybrid A* algorithm in obstacle cost map and C-MAP, the results are shown in Figure 15.

In the first part of the test, it is evident that after adding the virtual rough area in the obstacle cost map, the path avoids the rough area. It is proved that the modified Hybrid A* algorithm can detect and avoid the large roughness.

In the second part of the test, it is intuitive to see in the comparison diagram that the red line deviates from the rough area better than the black line. In order to quantitatively evaluate the different paths in two scenes, we measure the path length and the accumulative path cost of the tires on both sides.

Since it is challenging to integrate the path length and accumulative cost into actual vehicle operating cost, we use the entropy method [39] to calculate the weight and get the total score of each path to approximate the operating cost with the Equations (13) and (14).

$$E_j=-ln{\left(n\right)}^{-1}{\displaystyle \sum _i^n{p}_{ij}ln\left(p_{ij}\right)},$$

(13)

$${\omega }_i=\frac{1-E_j}{k-{\displaystyle \sum E_j}}(i=1,2,\dots ,k),$$

(14)

where $p_{ij}=Y_{ij}/{{\displaystyle \sum }}_{i=1}^n{Y}_{ij}$, $Y_{ij}=\frac{X_{ij}-min\left(X_i\right)}{max\left(X_i\right)-min\left(X_i\right)}$, $X_i=\{length,cost\}$. $X_i$ is a matrix of path length and accumulative cost information for each path.

In

Table 4. Planning Path Evaluation Table in Two Scene.

		Path Length (m)	Accumulative Cost	Total Score	Optimization Rate (%)
Scene 1	OCM	40	1124.4	1113.8	5.93
Start I	C-MAP	43	1057.7	1047.8
Scene 1	OCM	39	926.9	918.2	5.26
Start II	C-MAP	37	878.1	869.9
Scene 1	OCM	41	316.4	313.7	6.87
Start III	C-MAP	39	294.6	292.1
Scene 2	OCM	57	386.2	382.9	18.23
Start I	C-MAP	61	315.7	313.1
Scene 2	OCM	54	1174.8	1163.6	9.82
Start II	C-MAP	58	1059.4	1049.3
Scene 2	OCM	62	1096.8	1086.4	6.99
Start III	C-MAP	59	1020.1	1010.5

, the optimization rate refers to the optimization degree of the total score of C-MAP compared with the total score of OCM. The results of path evaluation Table 4 show that the planning path with consideration of roughness is better than that without consideration. In these three scenes, the optimization rate is more significant than zero, which shows that the modified algorithm can reduce the potential cost of the planning path.

3.2. Analysis of Algorithm Robustness

To test the optimization of the modified Hybrid A* algorithm under random conditions, we have conducted thousand planning experiments in two scenes with randomly given start position and goal position. We are moreover repeat the experiment at three different terrain resolutions scenes (0.1 m, 0.5 m, 1.0 m). Figure 16 is the schematic diagram of two scenes with lower resolution.

With the decrease of terrain resolution, C-MAP loses a lot of surface roughness information. As a result, in the planning of low-resolution maps, the optimization rate is more concentrated in the 0% to 10% range, as shown in Figure 17c–f. Significantly, when the resolution is from 0.1 m to 0.5 m, the average optimization rate of two scenes decreases rapidly, while the latter decreases slowly from 0.5 m to 1.0 m, as shown in Figure 18. The result shows that the planning algorithm cannot optimize the trajectory well when the terrain resolution (grid accuracy) is reduced from high to low.

Based on the results of the random test, we consider several factors that affect the optimization rate:

• The influence of planning map construction, including the DSM accuracy and filter threshold. The accuracy of DSM directly affects the recognition ability of the construction algorithm for surface roughness. However, the improvement of accuracy will result in more computational work and the program needs to process more data units. Therefore, 0.1 m is an appropriate accuracy value.
• For the filter threshold, it affects the ability of the construction algorithm to distinguish obstacles from the surface. Moreover, it is difficult to construct a surface roughness model if the terrain is flat.
• The influence is caused by the parameter setting of Hybrid A* algorithm. The setting of parameters will directly affect the quality of the final path. To illustrate, if the ReverseCost setting is greater than ForwardCost, the algorithm will consider the backward situation more. If the Extension Interval setting is too large, it will cause extra computation. Switch Cost penalizes the direction change of a vehicle. However, in this paper, since Tire Cost is a dynamic value, the static value of Switch Cost will weaken its effect so that it can be considered as an adaptive value.
• Caused by statistics of path cost. The entropy method may not be able to integrate the path length and the accumulative cost on the surface very well, which makes the calculation error.

4. Conclusions

The trajectory planning of cutting zone for unmanned trucks is essential in the autonomous transportation of open-pit mines. When planning the path, if we only consider the obstacle information, the roughness of ground may lead to extra tire wear. This paper proposes a path planning method considering terrain factors. (1) generate the DSM of the cutting zone using the UAV oblique photogrammetry method. (2) based on the DSM, construct a high-precision cost map (C-MAP) which contains both obstacle and roughness terrain information. (3) On the basis of C-MAP, plan a path by the modified Hybrid A* algorithm that incorporates the surface roughness information into the g-value estimation function. The random tests of actual scenes demonstrate a 10% to 20% reduction of the path cost. Therefore, our proposed method can reduce the transportation cost of the haulage truck and solve the problem of high cost in the last section of open-pit transportation.

5. Future Work

Due to the variability of the cutting zone environment, the method is only applicable to the global static map before transportation. During the mining operation, local elevation will change. Therefore, future work can be divided into two parts:

• We will collect the environmental information of the cutting zone and construct a real-time local DSM map by using detection sensors, such as Lidar, camera, GNSS receiver, etc. This enables continuous path planning by continually updating the global map and maintaining its timeliness and accuracy. Due to the errors in acquisition equipment and mapping procedure, the matching and fusion updating of the real-time local map and the global map is the focus and difficulty of future works.
• We will study the relationship between road roughness and vehicle transportation cost. Establishing a model to predict vehicle transportation costs based on different surface roughness conditions will be the critical feature of future work.