2.1.1 Graph Theory-Based Path Search

Particle Swarm Optimization (PSO), Genetic Algorithm (GA) and ACO are classical methods of path search. By comparing the performance of path planning based on different algorithms in raster map, the algorithm more suitable for graph theory is selected. Specifically, the static simulation environment (20 x 20) is designed, and the starting point and ending point of the target are known in the environment. PSO, GA and ACO are used for shortest path planning. The parameters of ACO are set as follows: the number of ants $M=10$, the maximum number of iterations $N_{\max }=50$, information heuristic operator $\alpha =1$, expectation heuristic factor $\beta =5$, pheromone intensity $Q=5$, pheromone residue operator $\rho =0.7$. The parameters of PSO are set as follows: particle size $S=20$, number of iterations $G_{\max }=50$, inertia weight operators $w=0.5$, acceleration operators $c_1=c_1=1$, $r_1=r_2=0.3$. The parameters of GA are set as follows: population size $NP=20$, iteration times $G_{\max }=50$, crossover rate $p_c=0.7$, mutation rate $P_m=0.05$.

Figure 1 Path planning results of different algorithms.

The path planning results based on PSO, GA and ACO after 50 iterations are shown in Figure 1. The path obtained by PSO occupies 38 grids and is turned 16 times. The path obtained by GA algorithm occupies 42 grids and is turned 18 times. The path obtained by ACO algorithm occupies 36 grids and is turned 7 times. The planning results show that ACO is better than PSO and GA under the same number of iterations, so it is more suitable for solving the path planning problem.

Figure 2 The framework of the two-layer path planning.

The selection operator in GA can ensure that all the selected individuals are excellent, but the mutation operator and crossover operators only introduce new individuals, and their operations cannot guarantee whether the new individuals are excellent. If the newly generated individuals are not sufficiently excellent, the introduction of new individuals will become an interference factor, which will slow down the search speed of GA. For PSO, the biggest problem in the application of graph theory path planning is that there is no unified prior knowledge of parameter setting (such as inertia weight, acceleration coefficient) and application background, resulting in unstable algorithm performance. ACO is also a type of feedback random search algorithm, but it relies on the pheromone (based on fitness function) left by ants for optimization. The initialization setting of the pheromone and the selection probability are closely related to the length of the path. This means that a shorter path attracts more ants. At the same time, the amount of residual pheromone is greater, which speeds up the optimization process. Therefore, ACO is chosen as the basic method to solve the problem of searching paths on graphs.

3.2.1 Adaptive spatial decomposition strategy

1) The coordinate system is built with the target as the center, and the distance between the initial position of the UAV and the target is the far boundary. Extraction of terrain, wind speed, and visibility within the airspace. Due to the wide range of airspace, an adaptive spatial distribution decomposition of airspace based on obstructive factors is proposed. Among them, the obstacle factor is the environmental factor that causes the UAV's flight direction to be impossible and the heading direction must be changed [28]. In this paper, elevation greater than 5200 m, visibility less than 0.8 km, and wind speed greater than 25 m/s, are considered as obstruction factors. As shown in Figure 3(a), there are two obstacle factors in the airspace, which is divided into three-ring layers.

2) To reduce the evaluation errors among sub-regions, each ring layer is divided into 32, 16, and 8 sub-regions from outer to inner layers. As illustrated in Figure 3(b), the first layer (o-R1) of the ring region is divided into eight sub-regions. The second layer of the ring region (R2-R1) is divided into 16 sub-regions. The third layer (R3-R2) is divided into 32 sub-regions.

Figure 3 Spatial decomposition based on obstacle ring.

Figure 4 Connectivity relationship of sub regions.

3) Connectivity of subregions. The $region(i,j)$ represents the relative position between sub-regions, indicating that the current sub-region is located in the j-th discrete region of the i-th ring layer. There are only three choices from the current sub-region to the next sub-region. The connectivity of sub-regions is illustrated in Figure 4. If the numbers of sub-regions in two adjacent layers are the same, the next possible sub-regions are shown in Figure 4(a). If the numbers of sub-regions in two adjacent layers are different, the next sub-regions are demonstrated in Figure 4(b). Therefore, the mathematical description is shown in Eq.(11). There is no connectivity between the innermost sub-regions. The ring layers are unidirectionally connected, which can only be connected from the outer ring to the inner ring.

where $num\left(L_{i+1}\right)$ are the total number of neutron regions in $L_i$ and $L_{i+1}$ layer respectively, $n-1$, $n+1$ are $num\left(L_i\right)$-ary operation, $\lceil \frac{n}{2}\rceil$ is a ceiling operation.

3.2.2 Calculation of heading selection probability

During UAV flight operations, this study primarily considers the impacts of elevation, visibility, and wind speed on reconnaissance mission performance. Through information fusion, we subsequently evaluate the flight selection probability within each sub-region.

(1) Most UAVs adopt low-altitude penetration tactics, so the terrain elevation has a decisive effect on the heading. The lower the ground elevation, the stronger the shielding ability, and the greater the selection possibility $P_e$.

(2) When the visibility is below 1.22 km, the clarity of vision is poor, which seriously threatens the flight safety of the UAV, and its flight direction selection is limited. When the visibility is between 1.22 km and 6 km, it significantly affects the clarity of vision and the degree of radar attenuation. When the visibility exceeds 6 km, the possibility $P_v$ of choosing the direction is higher.

(3) With the increase in wind speed, the control performance and maneuvering stability of the UAV deteriorate. Therefore, when selecting the flight direction, it often avoids regions with high wind speeds. When the wind speed is less than 1.5 m/s, it does not affect the UAV's flight, and the selection possibility is high. Wind speed starts to impact operational performance in the range of 1.5 m/s to 20 m/s, and the selection possibility decreases as wind speed increases. When the wind speed exceeds 30 m/s, it becomes unsafe to fly, and the selection possibility $P_w$ is low.

Combined with the grading standards of terrain and meteorology, this paper makes a fuzzy representation of environmental factors, as shown in Table 1.

Then the calculation formula based on multi-source information fusion of the heading selection probability is shown in Eq.(12), where $S_{(i,j)}^{\prime }$, $S_{(i,j)}$ represent the barrier area and total area of the subregion $(i,j)$ respectively.

3.3.1 Heading planning objective function

The course of UAV is described by selecting the sub-region sequence set, $course=\{regio{n}_1,\cdots ,regio{n}_N\}$.

The quality of feasible solution is evaluated from sequence length, possibility, and tortuosity. The length of the sequence is not only limited by the flight range of UAV, but also affects the flight time; a shorter sequence length results in a shorter flight time, which is crucial for completing the flight task efficiently. A higher mean probability of the sequence reduces the risk to the UAV. The sequence's tortuosity should be as smooth as possible due to flight maneuverability constraints. We propose that the direction sequence quality be evaluated according to Eq.(13).

where $P_L\left(course\right)$, $P_{cs}\left(course\right)$, $P_{\theta }\left(course\right)\in (0,1)$ represent the length ratio, mean probability, tortuosity of sequence respectively. The calculation formulas are as shown in (14). And $w_{P_L}$, $w_{P_{cs}}$, $w_{P_{\theta }}$ are the corresponding weights.

where $L_{course}$ represents total length of sequence sets, $L_{min}$ is the straight-line distance from the flight start point to the end point, $L_{max}=L_{min}+2\pi L_{min}$. $C{S}_{course}$ represents total probability of sequence sets, $N$ is the number of sub-regions in the course sequence set. ${\theta }_{max}$ is a constant, representing the maximum allowable direction tortuosity of UAV.

3.3.2 IACO

When the flight environment is complex, ACO is challenged due to the large search space when selecting the next region, leading to extended processing time and subpar real-time performance. To address this issue, an improved ant colony algorithm (IACO) is proposed.

Figure 5 Schematic of path safety.

(1) The new heuristic function is critical in inspiring ants to select nodes, directly affecting the solution quality. This new heuristic information considers the heading selection possibility of sub-regions and the distance between the sub-region and the flight destination. This allows for selecting sub-regions beneficial to the flight task, layer by layer.

where $P_{c{s}_{ij}}$ is the selection possibility of sub-region $(i,j)$, and $P_{d_{ij}}$ is the distance between the center of the sub-region $(i,j)$ and the flight end point, $P_{c{s}_{ij}}$, $P_{d_{ij}}\in (0,1)$, $w_1$ and $w_2$ are weights of two heuristic factors, and $w_1+w_2=1$. The calculation of $P_{d_{ij}}$ is shown in (16).

When the flight start point and end point are known, $L_{min}$ is a constant, and $d_{ij}$ represents the distance between selected sub-region and flight end point. According to Eq.(16), the closer the last selected sub-region is to the flight endpoint, the greater the probability of selecting the sub-region, which reduces the search range and improves convergence speed.

(2) To solve the problem that the global pheromone update rule based on positive feedback cannot guide ants to search for the optimal solution in time, a new pheromone update strategy based on positive and negative feedback is constructed.

1) For ants searching for feasible solutions from the start to the end, the pheromones are updated based on positive feedback, as shown in Eqs. (17) and (18).

where $J_{best}$ is the quality evaluation value of the optimal sequence of current iteration. Compared with the traditional ACO, the search scope of ants in IACO is more concentrated in the neighborhood of the optimal heading sequence to achieve global optimization.

2) For ants that do not find a feasible solution, pheromone concentrations in corresponding subregions are updated based on negative feedback, as shown in Eq (19) .

where $P_{d_{ij}}$ represents the distance between the last selected sub-region traversed by the deadlock ant and the flight end point. The attenuation of pheromone concentration is related to $P_{d_{ij}}$, which can not only avoid local optimization, but also avoid blind updating of global pheromone.

3.4.1 Objective function of 3D path planning

(1) The objective function of 3D path planning is similar to that of heading planning, as both have specific requirements for path length and smoothness. The difference lies in the calculation method. Moreover, flight path planning also takes into account the safety of UAVs.

where, $F_l$, $F_e$, $F_s$ are the three cost functions, ${\omega }_1$, ${\omega }_2$, ${\omega }_3$ are the corresponding weight coefficients, ${\lambda }_1$, ${\lambda }_2$ are the scaling factors. The smaller the value of the objective function, the better the quality of the path.

1) Path length

where $n$ represents the number of track points, and $(x_i,y_i,z_i)$ represents the position of track point $P_i$.

2) Path smoothness The 3D path design must meet its maximum climb angle and climb rate limits to ensure the curvature of the path and avoid sharp changes in direction and height fluctuations. $l_i$ represents the vector between two track points. Eqs.(22) and (23) respectively represent the yaw angle ${\phi }_i$ and the pitch angle ${\varphi }_i$. The cost function of path smoothness can be defined as Eq.(24).

3) Path safety In order to ensure the safe flight of UAVs, a safety-oriented mechanism must be incorporated into the path planning process. As shown in Figure 5, there is an obstacle K in the airspace, whose central coordinate is $O_k$ and radius is $R_k$. The vertical distance $d_k$ between the flight node of the UAV and the obstacle should be greater than the safety distance threshold $S=800$ m.

Figure 6 The influence of three main parameters.

(2) constraint condition In order to achieve the feasibility of 3D path planning, constraints such as the flight performance of UAV must be considered, and the objective function exceeding the constraints should be punished with a penalty coefficient $P_f$.

1) Maximum pitch angle. The pitch angle represents the angle of rotation around the X-axis. Constraining the angle change in the vertical direction when the UAV moves between adjacent nodes is crucial to prevent potential damage or crashes.

2) Maximum yaw angle. The yaw angle of the UAV involves its lateral rotation around the vertical axis (y-axis). The change in yaw angle must be strictly controlled to ensure that it stays within a safe range and prevents flight instability.

3) Flight altitude constraints. In order to ensure safe flight, the UAV should have a certain safe distance from the ground.

3.4.2 AMDBO

(1)Initialization strategy based on Tent Chaos When dealing with complex optimization problems, the original DBO randomly generates a population for initialization. However, this approach may lead to low population diversity, uneven population distribution, and rapid convergence to local optimal solutions. Tent Chaotic mapping offers better uniformity and correlation. Therefore, we utilize it as a new initialization method to enhance the population diversity of DBO. The formula is presented in Eq.(29).

(2)Reverse learning strategy

Since the global search ability of reproduction and stealing behavior decreases with the increase of iterations, a random reverse learning strategy is introduced to enhance the global search ability of these two behaviors.

where, $x_r^t$ is the reverse solution, $lb$ and $ub$ are D-dimensional vectors representing the lower and upper bounds of each dimension, and $rd$ is the D-dimensional random vector of $(0,1)$, $x_i^t$ is the current feasible solution.

Figure 7 Heading planning results in M1.

Figure 8 Heading planning results in M2.

(3) Adaptive population variation strategies In the original DBO, the ratio of the four behaviors in the population is not specified, and each individual can only perform one behavior. This limitation may lead to an inadequate search in the solution space or slow convergence. To address this issue, this paper proposes an adaptive behavior variation strategy.

1) Behavior variation strategy is based on the number of iterations. In this paper, the iteration number is used to simulate time, and the behavior variation of dung beetles is carried out every M iterations. This process changes the current behavior into the following behavior, allowing each dung beetle to perform four action behaviors. This approach maximizes the advantages of different behaviors and enhances the optimization ability.

2) Behavior variation strategy based on population similarity. Cosine similarity is introduced to measure population similarity, as shown in Eq.(31).

When the value is less than or greater than 0.5, individuals engaging in reproduction and foraging behaviors will transition into rolling ball behaviors. This transformation enhances the algorithm's global search capability, promoting species diversity. After finding the new best individual or reaching the iteration threshold of the variant individual, the variant individual will return to the original behavior and continue searching.

4.1.1 Parameter setting

The setting of hyperparameters significantly impacts on the performance of evolutionary algorithms. Since there are no theoretical methods to determine the best parameter combination, hyperparameters must be fine-tuned and optimized through experiments to ensure the algorithm can find high-quality solutions efficiently [19]. In the ACO, pheromone and heuristic information play crucial roles in ant selection. Therefore, we focus on three key parameters: pheromone important parameter $\alpha$, heuristic factor important parameter $\beta$, pheromone evaporation coefficient $\rho$. The other initial parameters are: the total number of ants is 50, the maximum number of iterations is set to 100, and pheromone increase intensity coefficient $Q=20$.

This paper analyzes the influence of different values of $\alpha \in \{1,1.5,2,2.5,3\}$, $\rho \in \{0.1,0.2,0.3,0.5,0.8\}$, and $\beta \in \{1,3,5,7,9\}$ on the convergence results. Set a set of default values: $\alpha =2$, $\beta =3$, $\rho =0.2$. When analyzing the values of multiple parameters using the control variable method, one parameter's value is changed at a time while keeping the other parameters at their default values. To enhance the algorithm's robustness, the selected parameter combinations from each group are simulated in 20 random two-dimensional probability maps, and the average convergence values are compared and analyzed. The experimental results are shown in Figure 6.

It can be seen from the results in Figure 5 that when $\alpha =2$, $\beta =5$, $\rho =0.3$, the average convergence value reaches the maximum value, and the influence of $\alpha$ on the average convergence result fluctuates the most. Therefore, the combination of $\alpha =2$, $\beta =5$, $\rho =0.3$ is adopted in the subsequent experiments in this paper.

4.1.2 Contrastive analysis

In order to further verify the effectiveness of heading planning based on IACO, the heading sequence is planned in the same 2D probability map (M1 and M2) based on D*, ACO, and GA. Figure 6 and Figure 7 display the planning results in M1 and M2, respectively.

Figure 9 Convergence curves of three algorithms.

In Figure 7, the heading sequence sets of ACO and IACO exhibit high similarity. The results of D* and the other three methods show the most significant difference. In Figure 8, the difference between the heading sequence sets of different methods is more significant than in M1, with the overall distribution ranging from 0° to 180° and 180° to 360°. The similarity between ACO and IACO is higher in the 8 ring layers near the target. The heading sequences obtained by GA and D* also align consistently within the 8 ring layers.

In this paper, the quality of the heading sequence is evaluated based on the length (HL), tortuous degree (HT), and selection possibility (HP/) of the heading sequence. The evaluation results of the heading sequence in Figure 7 and Figure 8 are presented in Table 2. It can be seen that the heading quality obtained by the four different algorithms in M2 is lower than that obtained in M1. The sequence length and sequence tortuous degree based on D* are the maximum values in the two maps, and the choice possibilities in the set of sub-regions vary widely. This result indicates that the heading sequence based on D* is trapped in a local optimum. IACO achieves the shortest length, the smallest tortuous degree, the largest average choice possibility in the two probability maps, and the highest comprehensive quality.

In addition, the convergence of GA, ACO, and IACO in the two maps is compared, as shown in Figure 9. The convergence quality and convergence speed based on the three algorithms in M2 are not as good as in M1. The convergence rate of GA is slow, and the convergence value of GA in M2 remains stable at a low level, indicating that GA also gets trapped in local optima. The main reason is that the mutation and crossover operators in GA cannot guarantee the quality of new individuals, thereby slowing down the optimization process. In comparison, ACO and IACO exhibit faster convergence speeds and higher quality. This result shows that ACO and IACO are more suitable for solving graph search problems.

Figure 10 Comparison of iterative convergence curves.

The convergence performance based on IACO in the two maps is the best, converging to 0.815 and 0.755 in the 8th and 11th iterations, respectively. The results show that the heading planning based on the IACO can provide reliable decision support for the subsequent 3D path planning.

Figure 11 Environmental information.

Figure 12  Multi-source environmental information of subregions.

4.2.1 Test analysis based on benchmark functions

Compared with heading prediction in a 2D map, the algorithm's global search capability and search efficiency are more critical for the precise planning of a 3D path. DBO, Improved Grey Wolf Optimizer (IGWO), and Sparrow Search Algorithm (SSA) are comparisons to verify the optimization performance of AMDBO. Six functions with different optimal characteristics are compared, and their specific details are presented in Table 3. F1 and F2 are single-peak reference functions that are used to evaluate the algorithm's single-objective solving capability. F3-F4 represent multi-modal reference functions, while F5-F6 are mixed reference functions that assess the algorithm's ability to escape local optima.

(1) Performance analysis In order to enhance the reliability of the experimental results, this paper adopts a consistent experimental setting by fixing the population size of each algorithm at 30 and the number of iterations at 500. Thirty independent experiments are performed on each function to obtain the target solution under each run. The average value (AVE), best value (BEST), and standard deviation (STD) are then calculated. These values serve as indicators to evaluate the reliability and stability of the optimization algorithm. Table 4 displays the test results of the four optimization algorithms on the six test functions.

Figure 13 Heading sequence planning results derived from the multi-source data-fused two-dimensional map.

Figure 14 Path planning based on AMDBO.

The AVE, BEST, and STD in F1 based on AMDBO are superior to the other three algorithms. When compared with the test results of F2, it is evident that the BEST of AMDBO is 57.76% lower than that of DBO, but the AVE and STD of AMDBO are 73.80% and 109.52% higher than those of DBO, respectively. Regarding the overall performance of the two single-peak functions, the overall performance of AMDBO is better than that of other three algorithms. The STD of AMDBO is slightly lower than that of SSA on F4. However, the AVE and the BEST of AMDBO are the highest, and the accuracy is greater than that of DBO and IGWO in F3 and F4. In F5 and F6, the best values of the four algorithms can achieve the theoretical optimal solution. Although the AVE and STD of AMDBO are slightly lower than IGWO on F5, they are superior to SSA and DBO. Additionally, the STD of AMDBO in F6 outperforms the other three algorithms.

Figure 10 illustrates the average convergence curves of 4 algorithms in 30 independent tests across 6 test functions. It can be seen that AMDBO achieves a higher objective function value within the same number of iterations during single-objective optimization. It exhibits the fastest optimization speed and the highest convergence accuracy. AMDBO can swiftly escape local optimal solutions in multi-objective optimization, leading to faster convergence with fewer iterations.

Figure 15 Path planning in 3D map based on other meathods.

To sum up, the AVE and BEST tested by AMDBO on the test functions are mostly better than the other three algorithms. With the same number of iterations, the optimization accuracy of AMDBO is higher. Most test results of the DBO algorithm are better or even superior to the SSA and IGWO, but worse than AMDBO.

4.2.2 Ablation Study

AMDBO introduces three improvement strategies to enhance the DBO algorithm. By implementing these strategies, three different algorithms derive through individual enhancements to DBO. The algorithms are as follows: 1) DBO1, which incorporates the Tent Chaos Initialization strategy. 2) DBO2, which integrates the reverse learning strategy of breeding and stealing. 3) DBO3, which includes the adaptive mutation strategy. These four algorithms are evaluated on test functions , and the AVE, BEST and STD of each algorithm are calculated after 30 independent runs. The experimental results are presented in Table 5.

In summary, the average optimization accuracy of a single improved strategy is higher than that of DBO. However, its variance is worse than that of DBO in test functions, indicating that the stability of a single improved strategy is insufficient. Compared with DBO, the performance indicators of AMDBO are improved, indicating that the comprehensive improvement of the three strategies is effective.

4.3.1 Simulation scene test

Taking the endpoint as the center of the circle and the distance between the starting point (-674,156) and the endpoint as the radius, the environmental information such as elevation, visibility, wind speed, and obstruction factors in this vast airspace is simulated, as shown in Figure 11. Then, the UAV trajectory is planned based on the fused multi-source environmental information.

According to the spatial distribution of obstacle factors in Figure 11(d), the airspace is adaptively decomposed into 80 km, 200 km, 310 km, 460 km, 540 km, and 700 km. The environmental situation information in each sub-area is shown in Figure 12. Comparing it with the actual scene grid information in Figure 11(a-c), it can be observed that the decomposition method accurately reflects the changing trend of environmental information while reducing the number of subregions. The results indicate that the adaptive spatial decomposition method effectively addresses the issue of sparse environmental information.

Based on Eq.(12), the 2D probability map is constructed, as shown in Figure 13(a). The IACO is used to search the heading sequence, and the result is shown in Figure 13(b). The flight direction gradually changed from 292.5° to 225°-270°. The average choice possibility of the route is 0.46, and the minimum value is 0.32, indicating that flight safety is relatively reliable. The tortuous degree of the heading sequence is 0.52, and the length is 1075.1 km, which is relatively smooth and short.

The AMBDO is used to plan a 3D path in regions where the heading sequence is located. The results are shown in Figure 14. The length of the 3D path is 4536 km, the average height is 4674 m, and the average tortuosity is 7.9. The 3D path is relatively smooth, and obstructions can be avoided.

4.3.2 Contrastive analysis

In order to verify the advantages of the two-layer path planning model proposed in this paper, the AMDBO is used to plan the path directly in 3D space. The planning results are shown in Figure 15. It is evident from Figure 13 that feasible paths are successfully planned using the three methods. The paths generated by DBO are distributed in regions with elevated terrain. To ensure a safe distance from the ground, the resulting paths are at a higher altitude and exhibit a greater tortuous degree. The other two methods exhibit a similar path trajectory, opting to first navigate towards lower-lying areas.

Figure 16 Comparison of iterative convergence curves.

In this paper, path quality is evaluated based on path length (L), average flight height (H), tortuosity (T), and minimum distance from obstacles (Df). The evaluation results for the three methods are presented in Table 6. The path generated by DBO is shorter than that of AMDBO, but it exhibits high tortuosity. The path produced by AMDBO has a lower average height, enhancing its ability for concealed flight, but it has a longer path length. The proposed method demonstrates the shortest path length and tortuosity, with a flight altitude comparable to other methods, resulting in the highest overall quality.

In addition, the convergence curves of the two models are compared, as shown in Figure 16. The DBO jumps out of the local solution twice and finally converges to 2091 in the 574th iteration. Although the objective function value of AMDBO is higher than that of DBO in the 1-73 iterations, its convergence rate is faster, and the convergence value is lower than that of DBO. The initial value of the proposed method is close to the convergence value of AMDBO and converges to 1007 at the 81st iteration. Compared with AMDBO and DBO, the convergence speed of the proposed method is increased by 75% and 85.89%, respectively, and the convergence values decreased by 31.87% and 52.05%, respectively.

Based on the above analysis, the two-layer path planning framework proposed in this paper can not only fully consider the impact of the environment on the UAV mission but also carry out fine path planning combined with the maneuvering performance of the UAV. Simultaneously, the environmental information of large airspace is analyzed layer by layer from fine-grained to coarse-grained to fine-grained. This effectively reduces the computational complexity and the execution efficiency of the optimization algorithm, thus improving the efficiency and quality of the flight path planning.