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

: 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 modiﬁed Hybrid A* algorithm. Firstly, we process the high-precision map based on di ﬀ erent 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 ﬁnal 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 di ﬀ erent 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.


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 1.
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.

2.
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.

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: . . . . . . . . .
where e ij is the corresponding coordinate of elevation points.

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: where is the corresponding coordinate of elevation points.

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: where is the corresponding coordinate of elevation points.

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 m as the terrain threshold of impassable elevation difference, and θ= 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): 1.
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). 2.
Scanning each data point sequentially in E i : set a searching group with e j−1 , e j , e j+1 , ( j = 2, 3 . . . m − 1), 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.

3.
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.
Appl. Sci. 2020, 10, x 6 of 25  , the red marks are the feature points filtered according to the threshold , and only the points with large elevation change are retained. The x-axis, which depends on the size of the point set represents the length of the scanline; the y-axis represents the elevation value of the point on the scanline.

Binary Map Construction
Binary map construction is the preprocessing work for Obstacle Cost Map construction. Firstly, Compare whether the slope angle of adjacent points , and , in greater than sequentially. If so, all the corresponding coordinate points between and 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.

Binary Map Construction
Binary map construction is the preprocessing work for Obstacle Cost Map construction. Firstly, Compare whether the slope angle β of adjacent points p k x p k , y p k and p k+1 x p k+1 , y p k+1 in F p greater Appl. Sci. 2020, 10, 6622 6 of 22 than θ 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.
where p j ∈ E i .
Appl. Sci. 2020, 10, 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. 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. 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.

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: 1.
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.

2.
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.

3.
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.
where gridCost(x, y) 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. α is the parameter to control the descent rate of cost, α > 0. 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. Appl. Sci. 2020, 10, x 9 of 25 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.

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:

1.
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.

2.
Using the section scanning method (threshold T = 0, θ = 5 • ) to process the open space area, deleting the elevation value of a point below the threshold from E r . 3.
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.

4.
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. 5.
Through the above steps, we can get the rough ground cost map of the cutting zone, as shown in Figure 8a.

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: 1. Removing the coordinate of obstacles in the original elevation point set , and obtained a new point set , which only keeps the elevation value of the open space area. 2. Using the section scanning method (threshold = 0, = 5°) to process the open space area, deleting the elevation value of a point below the threshold from . 3. In , taking a window and using the sliding standard deviation method to traverse , and calculating the standard deviation of the window, then assigning its value back to each point in this window. 4. Normalizing all the values of points in with a range of 0 to 1, which is the same as the obstacle cost range. 5. Through the above steps, we can get the rough ground cost map of the cutting zone, as shown in Figure 8a.

C-MAP Generating
After constructing OCM and RCM with the same value range at = 0,1 , the C-MAP can be constructed as follows: 1. Removing the obstacle points (cost = 1) from OCM, and adding the remaining OCM to RCM.

C-MAP Generating
After constructing OCM and RCM with the same value range at cost = [0, 1], the C-MAP can be constructed as follows: 1.
Removing the obstacle points (cost = 1) from OCM, and adding the remaining OCM to RCM.
The constructed C-MAP still maintain the value range at cost = [0, 1], 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.

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 (x, y, θ) extension of A* algorithm, it adds a three-dimensional search, considers the three-dimensional state (x, y, θ, λ) of the vehicle. Where x and y are coordinate information, θ is direction information, and λ 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.

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 x g , y g , θ g = (0, 0, 0), computes the shortest path to the goal from every point (x, y, θ). 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.
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.

Estimation Function Improvement
The Hybrid A* algorithm uses the Equation (6) to estimate the cost of each node extension.
where the first term ( ) is the sum of the estimated cost from the start to the current node, the second term ℎ( ) 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 ( ), and use the two heuristic functions of the original algorithm to calculate the ℎ( ). The calculation of ( ) is realized by Equation (7).

Estimation Function Improvement
The Hybrid A* algorithm uses the Equation (6) to estimate the cost of each node extension.
where the first term g(x) is the sum of the estimated cost from the start to the current node, the second term h(x) 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(x), and use the two heuristic functions of the original algorithm to calculate the h(x). The calculation of g(x) is realized by Equation (7).
g(x) = DirectionCost × (g(x parent ) + length(x parent , x) + TireCost(x parent , x)) + λ · SwitchCost (7) where DirectionCost contains ForwardCost (FC) and ReverseCost (RC), penalizing the motion cost according to a different direction. g x parent represents the sum cost from the start point to the parent node of the current node. length x parent , x is the length of the curved path between the parent node and the current node. SwitchCost is the direction switch cost to penalize the motion state, when the direction of the vehicle is changed λ = 1, otherwise λ= 0. We add a new cost term, TireCost x parent , x , 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: where the first term is the accumulative cost of the left tire, and the second is the right.

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: 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: The derivative is: where ω r is the weight of gradient; σ 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 ; ω 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.

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. 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.

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 is the cartographic information, and Table 2 is the parameters of the simulated haulage truck. Table 3 shows the main parameters of the Hybrid A* algorithm. The parameters of conjugate gradient smoother are set to = = 0.002, = = 0.2.   Table 1 is the cartographic information, and Table 2 is the parameters of the simulated haulage truck. Table 3 shows the main parameters of the Hybrid A* algorithm. The parameters of conjugate gradient smoother are set to = = 0.002, = = 0.2. Table 1. Cartographic information.   Table 1 is the cartographic information, and Table 2 is the parameters of the simulated haulage truck. Table 3 shows the main parameters of the Hybrid A* algorithm. The parameters of conjugate gradient smoother are set to ω r = ω o = 0.002, ω κ = ω s = 0.2.

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: 1.
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.

2.
Evaluating the path planned by the modified Hybrid A* algorithm in obstacle cost map and C-MAP, the results are shown in Figure 15.

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: 1. 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. 2. Evaluating the path planned by the modified Hybrid A* algorithm in obstacle cost map and C-MAP, the results are shown in Figure 15.

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: 1. 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. 2. 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). where cost . X i is a matrix of path length and accumulative cost information for each path.
In Table 4, 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.

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: 1.
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.

2.
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.

3.
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.

4.
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. 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:   1. 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. 2. 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. 3. 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. 4. 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.

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.

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.

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: 1.
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. 2.
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.
Author Contributions: Z.Z. and L.B. conceived, designed, and performed the experiments. Z.Z. analyzed the data and wrote the paper. All authors have read and agreed to the published version of the manuscript.
Funding: This research was supported by the Fundamental Research Funds for the Central Universities of Central South University, project number 2020zzts709.