1. Introduction

Drones come in different shapes, sizes and functions, ranging from small civilian drones to large military drones, and have been widely used over the past few years due to their low cost and low risk nature. For example, military drones in the current Russian-Ukrainian war are able to drop supplies into war zones and perform combat missions instead of single soldiers, greatly reducing the risk of casualties and transportation costs. At the same time, small civilian drones can be used for drone cluster light show performances, which have gained general recognition from the public. And no matter what kind of mission the UAV performs, it is inevitable to realize that it has to reach from one target area to another designated target area, and it is always important to keep the trajectory of the UAV away from any kind of threats and restricted areas in the meantime. Algorithm design is the core aspect of solving the complex route planning for UAVs. Due to the complex and variable nature of the UAV flight environment, traditional methods struggle to meet these conditions. Therefore, it is challenging to efficiently and accurately plan a safe and effective trajectory for the UAV during its mission.

Intelligent optimization algorithms for bionic populations have emerged in recent years and are widely used in various industrial and production fields. Intelligent algorithms are especially important in solving optimization problems in real engineering, such as particle swarm optimization algorithm , genetic algorithm , ant colony optimization algorithm, cuckoo search algorithm , firefly algorithm , artificial bee colony optimization algorithm and so on[1-6]. However all these algorithms fall into local optimal paths to some extent due to randomness as the number of iterations increases.

In this paper, a UAV trajectory planning algorithm based on modeling bird migration improvement is proposed. Birds migrate in organized flocks and their main behaviors include foraging and avoiding predators, which is similar to the situation of UAVs in complex terrains, where the UAV flies towards a set target while avoiding trap areas. The size and habitat of the birds exhibit a strong resemblance to the dimensions and terrain of UAV movement. The algorithm is proposed to simulate the migration process of birds using the simulation method and the algorithm is simulated and evaluated using matlab.

The widespread use of UAV technology has led to the development of path planning algorithms to meet diverse application requirements. Machine learning and deep learning techniques have become key drivers in the path planning field, enabling UAVs to learn planning strategies from sensor data. In addition, multi-intelligence collaborative planning, high-precision maps and sensor fusion, and dynamic environment adaptation techniques have all driven advances in UAV path planning. In the future, with continuous technological innovations, we can expect the development of smarter, adaptive, and safe path planning algorithms to cater to the requirements of UAVs in a variety of application areas. This development will further promote the extensive utilization of UAV technology in military, commercial, and civil applications.

2. Methodology

2.1. Terrain modeling

Since the UAV's working location is variable, it is located in different three-dimensional spatial environments, with differences in altitude, terrain complexity, etc., which can be specifically subdivided into plain areas, hilly areas, and mountainous areas [7]. Plain areas have small terrain ups and downs, and the vertical climbing and descending requirements for UAVs are not high. Hilly and mountainous regions have complex terrain, similar to the intricate environment in which UAVs perform their tasks. Therefore, the hilly terrain is simulated as the environment for UAV mission execution. The function simulation method is used to construct the 3D terrain features with the following formulas:

\( z(x,y)=sin⁡(y+a)+bsin(x)+ccos(d\sqrt[]{{y^{2}}+{x^{2}}})+ecos(y)+fsin(f\sqrt[]{{y^{2}}+{x^{2}}})+gcos(y) \) (1)

Where (x, y) is the point coordinate of a point on the terrain projected on the horizontal plane, and z is the height of the corresponding point coordinate. Where a , b , c , d , e , f , g are constant coefficients, which can be changed by altering the size of the constant coefficients in order to get different geomorphological features. At the same time, in order to simulate the effect of higher randomness and more closer to the actual conditions, the peak model is constructed on the basic terrain, using the same function simulation method. The specific formula is as follows:

\( h(x,y)={h_{i}}exp⁡[-\frac{(x-{x_{oi}}{)^{2}}}{a_{i}^{2}}-\frac{(y-{y_{oi}}{)^{2}}}{b_{i}^{2}}]+{h_{o}} \) (2)

h₀ and h_i denote the base terrain and the height of the i-th peak, respectively, (x_oi,y_oi) denotes the center coordinate position of the i-th peak, and ai and bi are the slopes of the i-th peak along the x-axis and y-axis, respectively.

2.2. Threat area modeling

While avoiding obstacles, drones often encounter areas that threaten flight safety, which we call threat zones. These threat zones can be the detection threat zones of enemy radar and anti-aircraft missile systems or some other threats, and once the UAV enters these zones, it is likely to be shot down or crashed. In order to simplify the model, this paper adopts a number of cylindrical regions with radius r to represent the threat areas, which can change the size of the radius r and the coordinates of the projected position of the center region of the cylinder, and the size of its radius determines the coverage of the threat areas[8]. The center of each cylinder is the place that poses the greatest threat to the UAV, and decreases outward sequentially. The code snippet is as follows:

%% Mapping of the threat areas

%Coordinates of the threat area center

ThreatAreaPostion = [80,50;50,100];

%Radius of threat area

ThreatAreaRadius = [40;20];

%0verlaying threat areas onto the map

figure

mesh(Map);

hold on;

for i= 1:size( ThreatAreaRadius)

[X,Y,Z] = cylinder(ThreatAreaRadius(i),50);

X = X + ThreatAreaPostion(i,1);

Y = Y + ThreatAreaPostion(i,2);

Z(2,:) = Z(2,:) + 50;%Height of threat area

mesh(X,Y,Z)

end

2.3. Algorithmic mathematical model construction

The main behaviors of birds include foraging for food and avoiding predators. This is similar to the flight of UAVs in complex terrain, where UAVs fly towards a set target while avoiding trap areas. The size and living area of the birds show a high similarity to the size and mission terrain of the UAV. I will use the artificial swarm algorithm as a reference and look forward to better simulation results and practical applications. Thus, I had the idea of simulating bird migration to plan the path of the UAV.

The basic principle of the algorithm is as follows. There are generally three types of individuals in a bird flock: discoverers, joiners, and early warning. Discoverers are responsible for finding food and communicating information about the location of the foraging area to other individuals in the flock, and joiners use the information provided by discoverers to obtain food.In the natural state, individuals engage in mutual monitoring, and members of a population typically compete for the food resources of those with high intake in order to enhance their own foraging rate.While foraging, all individuals are alert to their surroundings to prevent the arrival of natural predators[9], which is referred to as early warning behavior.

The main rules are as follows:

Discoverers usually have high energy reserves and are responsible for searching for food-rich areas in the entire population, providing foraging areas and directions for all joiners. The level of energy reserves in the modeling depends on the Fitness Value of the individual.

Once a predator is detected, the individual begins to chirp as an alarm signal. When the alarm value is greater than the safety value, the finder will take the joiner to other safe areas for foraging.

The identities of discoverers and joiners are dynamic. The roles of discoverers and joiners are constantly changing. Any individual has the potential to become a discoverer if they find a better food source, but the ratio of discoverers and joiners to the total population remains constant. In other words, when one individual becomes a discoverer, another must take on the role of a joiner. The lower the energy of the accessions, the worse their foraging position in the overall population. Some hungry joiners are more likely to fly elsewhere to forage for more energy.

During foraging, accessions were always able to search for the finder that offered the best food, and then either feed from the best food or forage around that finder. At the same time, some joiners may constantly monitor the finder in order to enhance their foraging efficiency and thus compete for food resources.

When aware of the danger, individuals at the edge of the group move quickly towards a safe area to gain a better position, and individuals in the middle of the population walk randomly to get closer to other individuals.

Figure 1. Technical route of the algorithm

First, initialize the coefficients of the population and set the number of population, the maximum number of iterations, the ratio of the discoverer and the follower, and the warning value. Then initialize the position of the population, calculate the fitness value function, compare and get the information of the optimal and worst fitness value; use the position update formula to update the position of the discoverer, follower and early warning; finally, get the global optimal position during the current iteration period, and if a better optimal position occurs, determine whether the replacement operation is needed until the preset maximum number of iterations occurs to terminate the loop, or until the optimal solution of the fitness function is searched.

Initialize the position of individuals in the population, assuming there are n individuals[10].

\( X=[\begin{matrix}{x_{1,1}} & {x_{1,2}} & ⋯ & {x_{1,d}} \\ {x_{2,1}} & {x_{2,2}} & ⋯ & {x_{2,d}} \\ ⋮ & ⋮ & ⋮ & ⋮ \\ {x_{n,1}} & {x_{n,2}} & ⋯ & {x_{n,d}} \\ \end{matrix}] \) (3)

Initialize fitness values for each individual in the population[10].

\( {F_{x}}=[\begin{matrix}f[({x_{1,1}} & {x_{1,2}} & ⋯ & {x_{1,d}})] \\ f[({x_{2,1}} & {x_{2,2}} & ⋯ & {x_{2,d}})] \\ ⋮ & & & \\ f[({x_{n,1}} & {x_{n,2}} & ⋯ & {x_{n,d}})] \\ \end{matrix}] \) (4)

An equation for updating the position of a discoverer in a population[10].

\( X_{i,j}^{t+1}=\begin{cases}X_{i,j}^{t}\cdot exp⁡(\frac{-i}{α\cdot ite{r_{max}}})    {R_{2}} \lt ST \\ X_{i,j}^{t}+Q\cdot L    {R_{2}}≥ST\end{cases} \) (5)

An equation for updating the position of a joiner in a population[10].

\( X_{i,j}^{t+1}=\begin{cases}Q\cdot exp⁡(\frac{X_{worst}^{t}-X_{i,j}^{t}}{{i^{2}}})    i \gt \frac{n}{2} \\ X_{p}^{t+1}+|X_{i,j}^{t}-X_{p}^{t+1}|\cdot {A^{+}}\cdot L    i≤\frac{n}{2}\end{cases} \) (6)

An equation for updating the position of an early warning in a population[10].

\( X_{i,j}^{t+1}=\begin{cases}X_{best}^{t}+β\cdot |X_{i,j}^{t}-X_{best}^{t+1}|    {f_{i}} \gt {f_{g}} \\ X_{i,j}^{t}+K\cdot (\frac{|X_{i,j}^{t}-X_{worst}^{t}|}{({f_{i}}-{f_{w}})+ε})    {f_{i}}={f_{g}}\end{cases} \) (7)

Table 1. Meaning of the mathematical symbols in the formulas

3. Simulation and results

The UAV mission environment is simulated using matlab, including terrain simulation and threat area simulation.Subsequently, the simulation algorithm is executed with the number of iterations set to 50 in order to obtain the optimal path trajectory and fitness value. The convergence relationship curve between the optimal fitness value and the number of iterations is then plotted.. The terrain and threat area simulation parameters are changed to compare the optimal path selection of UAVs under different mission environments, and at the same time, the similarities and differences of the iteration curves are compared to evaluate the performance of the algorithm. Part of the code for terrain simulation is as follows:

%topographic constant

a=1;b=1;c=1;d=1;e=1;f=1;g=1;

z= sin (y+ a) + b*sin (x) + c*cos ( d*( y^2+ x^2) ) + e*cos (y) + … f*sin(f*(y^2+x^2) ) + g*cos(y);

%\% Peak modeling

A = [60,20;100,80;150,150;]; % Peak Coordinate Points

h0 = 0; %Lowest position of the peak

hi =[50;60;80]; %Height of the three peaks

ai =[20;30;20]; %X-direction slopes of the three peaks

bi =[20;30;20]; %Y-direction slopes of the three peaks

h = sum( hi. * exp( - ( x- A( : , 1) ) . ^2. / ai. ^2 - ( y- A( : , 2) ) . ^2. / bi. ^2) ) + h0;

Value = max(z,h);

end

Part of the code for the threat area simulation is as follows:

%% Mapping of threat areas

%Coordinates of the threat area center

ThreatAreaPostion =[100,50;50,150];

%Radius of threat area

ThreatAreaRadius =[30;20];

%Overlaying threat areas onto the map

figure

mesh(Map);

hold on;

for i= 1:size(ThreatAreaRadius)

[ X , Y , Z ] = cylinder( ThreatAreaRadius(i) , 50) ;

X=X + ThreatAreaPostion(i,1);

Y=Y + ThreatAreaPostion(i,2);

Z( 2, : ) = Z( 2, : ) + 50; %Height of threat areal

mesh (X,Y,Z)

end

Set the coordinates of the start point and the end point, set the number of populations to 30 and the maximum number of iterations to 50, and also set the percentage of discoverers to 0.7, the percentage of early warning to 0.2, and the alert threshold to 0.6, and refer to the parameter settings of Xue [11]. The optimal path curves, optimal fitness values, and convergence relation curves obtained before and after changing the terrain and threat area simulation parameters are as follows:

4. Conclusion

According to the convergence relation curve, the following conclusions can be drawn. The high skewness of the convergence curve indicates that the algorithm is more capable of finding the optimum and has high generality, which is more suitable for solving some complex optimization problems than other simulation algorithms.The shortcoming lies in its inability to escape from local optimal solutions.The convergence speed is faster and the approximation to the optimal solution is faster. The disadvantage is that it is easy to fall into the local optimal solution prematurely. The algorithm is averaged over multiple runs and results in a stable approximation to the optimal solution. The algorithm maintains the correlation between the test data points better. Compared with other algorithms it has better global stability and greater development capability.

It is verified through a series of simulations and the experimental results are analyzed. It can be observed that the algorithm is able to give optimization results with superiority, stronger convergence, better global optimization performance and good stability. The comprehensive comparison results are highly satisfactory and can serve as an effective tool for addressing complex practical problems. Compared with some other algorithms, the advantages are obvious, the number of iterations is small, easy to program and implement, while avoiding the disadvantage of slower convergence. So the approximate solution of the problem can be found quickly and the global optimal point can be obtained very quickly without any a priori knowledge.

This study still has many shortcomings that need to be improved. This study only simulates the static terrain and the known shape, size and location of the threat area, i.e., the environment is known and the complexity is low. However, the actual UAV mission environment is complex and variable, with sudden environmental changes and unexpected threats that are difficult to predict in advance. Therefore, it is hoped that an adaptive algorithm will be designed subsequently to make corresponding strategy adjustments for the actual situation. To address the issue of easily falling into a local optimal solution, it is suggested that a randomization strategy be implemented to mitigate the impact of the local optimal solution.