Unmanned aerial vehicle three-dimensional path planning method and navigation equipment

Abstract

The invention relates to the technical field of path planning, and provides an unmanned aerial vehicle three-dimensional path planning method and navigation equipment. The method comprises the following steps: S100, acquiring a flight path planning model comprising an environment model and a total cost function; s200, solving the flight path planning model by adopting a honeypot algorithm, and initializing population parameters; s300, when the coefficient vector A meets the condition that A is larger than or equal to 1, a search strategy with random individual disturbance is executed, and otherwise, the positions of the badgers are still updated with Xprey as the reference; s400, when the number of iterations t reaches the maximum number of iterations tmax, outputting Xprey and a fitness value thereof; otherwise, increasing 1 to the number t of iterations, and returning to the step S300. Wherein A is linearly reduced along with the number of iterations, and in the earlier stage of iteration, the algorithm executes a random search strategy for multiple times, so that premature population aggregation is avoided, the ergodicity of badger individuals in a search space is enhanced, and the global search ability is improved.

Description

Technical Field

[0001] The present invention relates to the technical field of path planning, and in particular to a three-dimensional path planning method for unmanned aerial vehicles and a navigation device. Background Art

[0002] Three-dimensional path planning is a complex constrained optimization problem, and some algorithms and improvement studies have been proposed in the prior art. Traditional path planning methods include the A* algorithm, the artificial potential field method, and the rapidly-exploring random tree algorithm (RRT), etc. Such algorithms will have problems such as slow convergence speed for path planning in complex environments. Swarm intelligence algorithms have become an important method for solving path planning problems. This type of meta-heuristic algorithm simulates the information sharing and mutual learning among biological individuals in nature, and has stronger self-learning, self-adaptability, and self-organization. Specifically, there are the particle swarm algorithm, the artificial bee colony algorithm (Artificial bee colony algorithm, ABC), the whale optimization algorithm (Whale Optimization Algorithm), the Harris hawks optimization algorithm (Harris Hawks Optimization, HHO), the sparrow search algorithm (Sparrow Search Algorithm, SSA), the dung beetle optimizer (dung beetle optimizer, DBO), the crested porcupine optimizer (Crested Porcupine Optimizer, CPO), etc.

[0003] The honey badger algorithm (Honey Badger Algorithm, HBA) mainly performs optimization by simulating the intelligent foraging behavior of honey badgers. The honey badger algorithm initializes the population by using a random distribution method; in the excavation stage, the honey badger group is guided by the odor intensity of the prey for in-depth development, which can accelerate the convergence of the algorithm and improve the operation efficiency. However, when the food source falls into a local optimum, the update mechanism guided by the odor intensity of the prey will cause the algorithm to fall into a local optimum, and the local development ability of the honey badger algorithm is relatively weak, making it difficult to jump out of the local optimum when optimizing actual complex problems. Summary of the Invention

[0004] The purpose of the present invention is to provide a three-dimensional path planning method for unmanned aerial vehicles and a navigation device, aiming to solve the technical problem that the existing honey badger algorithm is difficult to jump out of the local optimum.

[0005] In a first aspect, the present application provides a three-dimensional path planning method for an unmanned aerial vehicle, and the method includes the following steps:

[0006] S100: Obtain a flight path planning model including an environmental model and a total cost function;

[0007] S200: Solve the flight path planning model using the honey badger algorithm and initialize the population parameters;

[0008] S300: When the coefficient vector A satisfies |A|≥1, the search mathematical expression of the population in the excavation stage is:

[0009] X new =X rand +Dir·β·I·X rand +Dir·r2·α·d i ·|cos(2·π·r3)×(1 - cos(2·π·r4))|;

[0010] Otherwise, the search mathematical expression in the excavation stage is:

[0011] X new =X prey +Dir·β·I·X prey +Dir·r2·α·d i ·|cos(2·π·r3)×(1 - cos(2·π·r4))|;

[0012] Where, X new is the updated position of the honey badger individual, X prey is the global optimal position; α is the density factor; the food acquisition ability factor β≥1; r2, r3, r4, r5 are random numbers between [0,1]; I is the prey control intensity; Dir is the search direction factor; X rand is a randomly selected honey badger individual under the current iteration; A = 2mr5 - m, m is the convergence factor;

[0013] S400: When the iteration number t reaches the maximum iteration number t max output X prey and its fitness value; otherwise, the iteration number t is incremented by 1 and return to step S300.

[0014] In a second aspect, the present application provides a navigation device, including a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the computer program, it implements the three-dimensional path planning method for an unmanned aerial vehicle as described in any one of the above.

[0015] The beneficial effects of the UAV three-dimensional path planning method and navigation device provided by the present invention are as follows: The present application improves the honey badger algorithm by introducing the coefficient vector A. When the coefficient |A|≥1, a search strategy with random individual perturbation is executed to perturb the optimal position of the current population, helping the algorithm to jump out of the local optimal trap; in the early stage of iteration, when |A|≥1, the search strategy with random individual perturbation is executed multiple times to avoid premature aggregation of the population, enhance the traversability of honey badger individuals in the search space, focus on global search, and avoid premature convergence, thereby improving the global search ability and solving the technical problem that the existing honey badger algorithm is difficult to jump out of the local optimum. BRIEF DESCRIPTION OF THE DRAWINGS

[0016] In order to more clearly illustrate the technical solutions in the embodiments of the present invention, the following will briefly introduce the drawings required for use in the embodiments or the description of the prior art.

[0017] Figure 1 It is a schematic flowchart of the UAV three-dimensional path planning method provided by the embodiment of the present invention;

[0018] Figure 2 It is a specific flowchart of step S110 in the embodiment of the present invention;

[0019] Figure 3 It is a specific flowchart of step S112 in the embodiment of the present invention;

[0020] Figure 4 It is a schematic structural diagram of the laser beam splitter in the embodiment;

[0021] Figure 5 It is a partial flowchart of step S100 in the embodiment;

[0022] Figure 6 It is a schematic diagram of the threat cost model;

[0023] Figure 7 It is a schematic diagram of population initialization;

[0024] Figure 8 It is another flowchart of the UAV three-dimensional path planning method provided by the embodiment;

[0025] Figure 9 It is a statistical chart of Friedman test rankings in 30 dimensions and 100 dimensions;

[0026] Figure 10 It is an iterative convergence chart of the CEC2017 function set in 30 dimensions;

[0027] Figure 11 It is an iterative convergence chart of the CEC2017 function set in 100 dimensions;

[0028] Figure 12 It is the box plot of all algorithms on CEC2017;

[0029] Figure 13 It is the exploration and exploitation curve of LRMHBA and HBA;

[0030] Figure 14 (a) is the front view and top view of the optimal path planning in Scenario 1;

[0031] Figure 14 (b) is the front view and top view of the optimal path planning in Scenario 2;

[0032] Figure 14 (c) is the front view and top view of the optimal path planning in Scenario 3;

[0033] Figure 15 It is the iterative curve of the average flight cost in three scenarios;

[0034] Figure 16 It is the structural schematic diagram of the navigation device provided by the embodiment of the present invention. Detailed implementation manners

[0035] The traditional honey badger algorithm has an update mechanism guided by the smell intensity of the prey. When the food source falls into a local optimum, it is easy for the algorithm to fall into the local optimum, and its local development ability is relatively weak, making it difficult to jump out of the local optimum.

[0036] In the first aspect, in combination with Figure 1 , the present application provides a three-dimensional path planning method for an unmanned aerial vehicle. The method includes the following steps:

[0037] S100: Obtain a flight path planning model including an environmental model and a total cost function.

[0038] S200: Solve the flight path planning model using the honey pot algorithm and initialize the population parameters.

[0039] S300: When the coefficient vector A satisfies |A|≥1, update the position of the population in the excavation stage according to Equation (1); otherwise, update the position of the population in the excavation stage according to Equation (2).

[0040] X new = X rand + Dir·β·I·X rand + Dir·r2·α·d i ·|cos(2·π·r3)×(1 - cos(2·π·r4))| (1)

[0041] X new = X prey + Dir·β·I·X prey+Dir·r2·α·d i ·|cos(2·π·r3)×(1 - cos(2·π·r4))| (2)

[0042] Wherein, X new is the updated position of the honey badger individual, and X prey is the global optimal position. Specifically, in the current iteration, the individual with the optimal fitness value is selected as X prey . X rand is a honey badger individual randomly selected in the current iteration. The honey badger individual is a solution in the flight path planning model. The coefficient vector A is calculated by equations (3) and (4), and m is the convergence factor.

[0043] A = 2mr5 - m (3)

[0044]

[0045] α is the density factor, which controls the time-varying randomization to ensure a smooth transition from exploration to exploitation. The density factor α is calculated by equation (6). Wherein, t max is the maximum number of iterations, and the constant C ≥ 1. In some embodiments, C = 2.

[0046]

[0047] The food acquisition ability factor β ≥ 1 represents the ability of the honey badger to acquire food. In some embodiments, β = 6. r1, r2, r3, r4, r5, r6, r7, r8 are random numbers between [0, 1], and the values of the random numbers are independent of each other, and the values can be the same or different.

[0048] I is the prey control intensity, which represents the ability of the honey badger to control the prey and can be obtained according to equations (7), (8), and (9). Wherein, S is the source intensity or concentration intensity; d i represents the distance between the prey and the current honey badger individual. Dir is the search direction factor to utilize a large number of opportunities to let the search individual strictly scan the search space, see equation (10).

[0049]

[0050] S = (X i - X i+1 ) 2 (8)

[0051] d i = X prey - X i (9)

[0052]

[0053] S400: When the iteration number t reaches the maximum iteration number t max , output X prey and its fitness value; otherwise, increase the iteration number t by 1 and return to step S300.

[0054] The improved honey badger algorithm (hereinafter referred to as LRMHBA) of the present application introduces a coefficient vector A. When the coefficient |A| ≥ 1, a search strategy of randomly perturbing individuals is executed to perturb the optimal position of the current population, helping the algorithm to jump out of the local optimal trap, enhancing the robustness of the algorithm, and making it more advantageous in dealing with complex three-dimensional path planning problems. Combining Equation (3) and Equation (4), it can be seen that |A| ≥ 1 is satisfied in the early stage of iteration, and the search strategy of randomly perturbing individuals is executed multiple times to avoid premature aggregation of the population, enhance the traversability of honey badger individuals in the search space, focus on global search, avoid premature convergence, and thus improve the global search ability, solving the technical problem that the existing honey badger algorithm is difficult to jump out of the local optimum.

[0055] In addition, |A| decreases linearly with the iteration number, focusing on local development in the later stage of iteration, improving the convergence accuracy, enabling LRMHBA to achieve a balance between global search and local development, and being able to find a high-quality path planning scheme faster.

[0056] In some embodiments, the flight path planning model in step S100 can be an environmental model drawn from on-site shooting, purchased from a commercial navigation company, or satellite shooting, which is not uniquely limited herein.

[0057] In one of the embodiments, the environmental model in step S100 is constructed by the following steps:

[0058] S110: The reference terrain model in the environmental model is represented by Equation (11):

[0059]

[0060] z1(x,y) is the height value at the terrain model point (x,y). k, a, b, c, d, e, f, g are coefficients used to control the digital terrain ratio and undulation. Different coefficients can simulate different reference landform features. For example, the coefficient k, as a scale factor or height scaling factor, is used to control the height range of the entire terrain; increasing k will raise the terrain as a whole, and decreasing k will lower the terrain as a whole. The coefficients a, b, c are used to control the basic shape and undulation frequency of the terrain. Among them, a controls the main undulation direction of the terrain (such as east-west or north-south). b and c control the secondary undulations or detailed features of the terrain. The coefficients d, e, f are used to introduce the irregularity or complexity of the terrain. d controls the amplitude of random perturbation of the terrain. e and f control the frequency or direction of random perturbation. The coefficient g is usually used to control the overall tilt or offset of the terrain.

[0061] S120: The mountain model in the environmental model is represented by Equation (12):

[0062]

[0063] where z2 is the mountain height value at point (x, y): h i is the height peak of the i-th mountain, x i and y i are the central positions of the i-th mountain; a i and b i are the slope parameters of the i-th mountain peak in the x and y directions. In the natural environment, the most threatening to the UAV is the natural mountain terrain. Using Equation (12) can simulate the complexity and diversity of the real mountain terrain. The UAV can identify the specific position, height, and slope of the mountain, so as to avoid the mountain in path planning and reduce the collision risk.

[0064] In one embodiment, in combination with Figure 2 , the key to the reference terrain model in step S110 specifically includes the following steps:

[0065] S111: Obtain the two-dimensional vector map corresponding to the reference terrain model. Among them, the two-dimensional vector map is represented by a two-dimensional graph composed of a series of lines connected by points. The two-dimensional vector map is unified into a standardized vector format through vectorization processing, reducing data errors and improving the accuracy of the reference terrain model. Optionally, the picture of the two-dimensional map is converted into a digital form through image recognition technology to generate a two-dimensional vector map.

[0066] S112: Obtain the shape features and height features corresponding to the two-dimensional vector map. Among them, the shape features generate information describing its geometric shape, and the geometric shape information covers the coordinates of points, lines, and surfaces, such as the type of building, the width and grade of the road, etc. The height feature is the height value corresponding to the shape feature.

[0067] S113: Calculate the coefficients k, a, b, c, d, e, f, and g according to the shape features and height features. Based on this, by obtaining the two-dimensional vector map and its corresponding shape features and height features, and substituting specific values into Equation (11), the coefficients in Equation (11) can be quickly obtained, which helps to quickly construct an accurate reference terrain model in a complex terrain environment.

[0068] In one of the embodiments, in combination with Figure 3 , the height features are obtained by the following steps:

[0069] S1121: Divide the linearly chirped optical signal into a laser detection signal and a local oscillator signal. Among them, the splitting ratio range of the laser detection signal and the local oscillator signal is in the range of 90:10 to 99:1. The local oscillator signal is used for demodulating the echo signal, and laser detection signals of different frequencies are directly emitted outward to shoot at the obstacle corresponding to the shape feature.

[0070] Specifically, in step S1121, first, the control module sends a modulation signal to the drive module. Second, the drive module combines the modulation signal with a necessary DC bias current, and after appropriate amplification or adjustment, forms a modulation current with the correct bias and modulation amplitude, and delivers the modulation current to the single longitudinal mode laser. Third, the single longitudinal mode laser drives the laser diode inside it to generate a laser signal with a single frequency; under the action of the modulation current, the frequency of the laser signal will periodically change within a preset bandwidth to form a linearly chirped optical signal. Finally, the laser beam splitter divides the linearly chirped optical signal into a laser detection signal and a local oscillator signal with a splitting ratio range of 90:10 to 99:1. The local oscillator signal is delivered to the mixer of the system for subsequent echo signal demodulation, and the laser detection signal is directly emitted at the obstacle corresponding to the shape feature.

[0071] In one example, combined with Figure 4 , the laser beam splitter includes a first coupling device 11, a phase regulator 12, and a second coupling device 13 that are optically connected in sequence. The first coupling device 11 is responsible for splitting the linearly chirped optical signal into two beams of light - the first beam and the second beam. At the second coupling device 13, these two beams of light are further divided into a laser detection signal and a local oscillator signal. The role of the phase regulator 12 is to adjust the phase of the first beam of light, thereby changing the phase difference between the two beams of light, so as to allow the second coupling device 13 to flexibly adjust the power distribution ratio between the laser detection signal and the local oscillator signal, increase the optical power ratio of the laser detection signal, significantly reduce the system loss, improve the utilization efficiency of light energy, and do not need to rely on the configuration of an optical circulator.

[0072] Among them, the splitting ratio can be dynamically adjusted in the range of 90:10 to 99:1, which can not only enhance the intensity of the laser detection signal to prevent the echo signal from being masked by the strong local oscillator signal, but also avoid unnecessary attenuation or discarding of the local oscillator signal, thereby optimizing the utilization of optical power. If there is a large deviation between the actual splitting ratio and the ideal value due to manufacturing process limitations, fine-tuning can be performed through the phase regulator 12 to ensure that the output power of the actual local oscillator signal is consistent or close to the output power in the ideal state. For example, when the optical power ratio of the laser detection signal and the local oscillator signal is set to 99:1, the laser detection signal exhibits excellent detection performance. At the same time, the intensities of the echo signal and the local oscillator signal are similar, which is conducive to the interference between the two to generate a clear interference signal, thereby simplifying the subsequent analysis process of the measurement results.

[0073] S1122: Multiple silicon photomultipliers arranged in an array respectively receive echo signals of different frequencies. Optionally, the multiple silicon photomultipliers are arranged in an M*M detector array, where each silicon photomultiplier is specifically configured to capture a specific, pre-set echo signal frequency. No two silicon photomultipliers will share the same echo signal frequency, ensuring that the echo signal frequencies detected by each silicon photomultiplier are different.

[0074] S1123: Process the echo signals to generate height values at each terrain model point.

[0075] Specifically, the specific process of step S1123 is as follows: First, a phase shifter is used to perform a phase distortion compensation operation on the echo signals of different frequencies, aiming to correct the phase errors that may be introduced during signal propagation. Subsequently, a beam combiner integrates and converges the light waves of the echo signals after phase distortion compensation to accurately calculate the distance to the obstacle. Finally, based on these distance data, the height characteristics of the obstacle at the terrain model points are inferred.

[0076] In one embodiment, after step S113, the following steps are further included:

[0077] S114: Use the data of the lidar to revise the reference terrain model. Since the reference terrain model may not immediately simulate the actual environment. When the data verification pass rate of the lidar is lower than 90%, that is, when the coincidence degree between the data of the lidar and the reference terrain model generated in step S113 is lower than 90%, for the non-coincident areas, use the data of the lidar to modify the above reference terrain model. For example, by adjusting the coefficients k, a, b, c, d, e, f, and g, and / or modifying the height characteristics, so that the data verification pass rate of the lidar ≥ 90%.

[0078] In one embodiment, after step S113, the following steps are further included:

[0079] S115: When there are unknown areas in the reference terrain model, the area of the reference terrain model is divided into a dangerous area, a safe area, and an unknown area. The dangerous area is enlarged by a preset multiple, and the safe area minus the enlarged dangerous area is obtained as the acquisition area. The dangerous area refers to an area with obstacles such as mountains, which makes it difficult for the drone to pass through. The preset multiple is positively correlated with the outer diameter of the lidar, so that the boundary of the enlarged dangerous area moves at least 2 times the outer diameter of the lidar outward compared with the boundary of the dangerous area before enlargement.

[0080] S116: The lidar reaches the acquisition area to map the unknown area, so as to update the reference terrain model and eliminate the unknown area. Among them, the acquisition area has a certain safety distance from the dangerous area before magnification and does not approach the dangerous area, so as to avoid the lidar colliding with the actual dangerous area and improve the safety of mapping.

[0081] In addition, when performing the mapping task in the unknown area, only relying on the lidar to operate independently, the data throughput required for this process is relatively low, and there is no need to configure other auxiliary sensors to participate in the mapping operation, which not only reduces the system complexity, but also significantly improves the environmental perception efficiency of the lidar, enabling it to quickly complete the construction and update of high-precision maps.

[0082] In some embodiments, in combination with Figure 5 , the total cost function is constructed by the following steps:

[0083] S130: Construct a flight distance cost function (see Equation (13)), taking into account flight efficiency and obstacle avoidance requirements.

[0084]

[0085] Among them, (x i , y i , z i ) is the coordinate of the i-th path point; n is the number of path points, and multiple path points are connected in sequence to represent the flight path, and the sum of the distances between every two adjacent path points is the flight distance. Based on this, the flight distance cost function quantifies the flight distance into a specific value, which is conducive to finding a flight path with a smaller flight distance cost in the subsequent algorithm iteration, and controlling the UAV subject to fuel capacity and fuel consumption rate limitations.

[0086] S140: Construct a flight altitude cost function (see Equation (14)).

[0087]

[0088] Among them, z i is the altitude of the i-th path point; is the average altitude of the flight path, Based on this, in order to optimize fuel efficiency and safety at a reasonable flight altitude, the flight altitude cost function calculates the sum of the absolute values of the differences between each path point and the average altitude, quantifies the altitude fluctuation of the UAV, reduces the climbing and descending actions of the UAV, reduces fuel consumption, and the UAV maintains a stable low-altitude flight mode to ensure concealment. By maintaining a relatively stable altitude, the UAV can more easily avoid obstacles and reduce the risk of collision.

[0089] S150: Construct a flight turning cost function (see Equations (15) and (16))

[0090]

[0091]

[0092] In the formula, b i represents the vector of the i-th path segment formed by the i-th path point pointing to the (i + 1)-th path point, and b i and b i+1 are the vectors of the i-th and (i + 1)-th path segments respectively; |b i | and |b i+1 | are the lengths of b i and b i+1 respectively. θ i is the angle between adjacent path segments, that is, the angle between the path segment vectors b i and b i+1 . Ф is the preset maximum turning angle. Based on this, a flight turning cost function C C quantifies the severity of the turn (i.e., the turning angle) into a specific numerical value. The closer the angle θ i is to the maximum turning angle Φ, the more significantly the flight turning cost function C C will increase, thus prompting the algorithm to select a smoother turning path and ensuring flight safety at the same time.

[0093] S160: Construct a terrain constraint cost function (see Eqs. (17), (18), and (19)).

[0094]

[0095] d i = z i - h i (19)

[0096] In the formula, z i is the height of the i-th path point; h i is the height of the terrain at this point; d i is the vertical distance from the i-th path point to the terrain surface, and d i ≤ SD. SD is the set terrain safety distance. Optionally, SD = 50 m to 500 m. For example, SD = 200 m. Based on this, the terrain constraint cost function C M quantifies the height close to the terrain ground into a specific numerical value, prevents terrain collision during task execution, and the planned flight path needs to maintain a terrain height advantage and a set safety redundancy.

[0097] S170: Construct an obstacle threat cost function (see Eqs. (20), (21), and (22)).

[0098]

[0099] Wherein, (x i , y i ) are the coordinates of the i-th path point; (x k , y k ) are the coordinates of the center of the threat area; d ik is the vertical distance from the i-th path point to the coordinates of the k-th threat center; R k is the radius of the k-th threat area; SF is the safety distance coefficient of the threat area. Based on this, when the UAV executes tasks, it needs to avoid colliding with obstacles and avoid threat areas such as radar. Combining Figure 6 , the threat area in the map is set as a cylinder model, and it is set that flight is not allowed at any height within the coordinates of the threat area. O k is the center coordinate (x k , y k ) of the k-th threat source. When d ik ≥ 2R k ·SF, it is regarded as having no threat cost, otherwise the threat cost increases quadratically as the distance decreases. The path planned by the UAV should not only avoid obstacles and threat sources, but also leave an appropriate safety distance from the threat area, just like the terrain constraint cost.

[0100] S180: Construct the total cost function (see Equation (23)).

[0101]

[0102] Wherein, μ L to μ T are respectively the weight coefficients of each cost function, and satisfy μ L + μ H + μ C + μ M + μ T = 1. If the planned path is feasible, the total cost function is the weighted sum of all cost functions, otherwise, the total cost function is a maximum penalty value CMAX. Optionally, CMAX takes the value of 10000.

[0103] It should be noted that the total cost function C total can include only the above five sub-cost functions, or can include only one or more than two of the above five sub-functions. It can be understood that in other embodiments, the total cost function C total can also include other cost functions and / or any combination of the above five sub-cost functions, which is not uniquely limited here.

[0104] In this embodiment, the fitness value is the reciprocal of the total cost function C total , and the total cost function C totalThe smaller it is, the larger and better the fitness value. The purpose of LRMHBA is to obtain a flight path with the largest fitness value.

[0105] In some embodiments, in combination with Figure 7 (a), step S200 includes:

[0106] S201: Implement population initialization using Latin Hypercube Sampling (hereinafter referred to as LHS). Among them, LHS can achieve full coverage of the parameter space through a uniform stratification mechanism. Under a limited sample size, it ensures that the sampling points are more evenly distributed in the search space, improving the search efficiency.

[0107] In some embodiments, in combination with Figure 7 (b), step S200 includes:

[0108] S202: The population is randomly initialized within the designed boundaries according to Equation (24).

[0109] X i = lb i + r1·(ub i - lb i ) (24)

[0110] Where x i is the i-th individual of the population. lb and ub are the lower and upper boundaries of the search. Random initialization can maintain the diversity of the population, enabling the algorithm to consider more possibilities during the search process, helping the algorithm to widely explore the search space, and avoiding falling into local optimal solutions.

[0111] In some embodiments, in combination with Figure 7 (c), step S200 includes:

[0112] S203: The population respectively uses LHS and random sampling to generate two samples. After mixing the two samples, the top 50% of the individuals with high fitness values are selected to achieve population initialization. It can be understood that in other embodiments, step S200 can also use other methods to achieve population initialization, which is not uniquely limited here.

[0113] In some embodiments, when the fitness value has not improved after continuous iteration K times, before returning to step S300, step S400 further includes the following steps:

[0114] S410: If the iteration number t ≤ the maximum iteration number t max*2 / 3, individuals with the first preset ratio perform differential mutation according to formula (25), and the remaining individuals perform differential mutation according to formula (26). Optionally, the first preset ratio is 30% - 80%. For example, if the first preset ratio is 50%, then 50% of the individuals perform differential mutation according to formula (25), and the remaining 50% of the individuals perform differential mutation according to formula (26). Optionally, sorted by fitness value, the individuals with the first preset ratio in the front perform differential mutation according to formula (25), and the remaining individuals perform differential mutation according to formula (26), which is not uniquely defined here.

[0115] V i,G = X i,G + F·(X best,G - X ri,G ) + F·(X r9,G - X r10,G ) (25)

[0116] V i,G = X r9,G + F·(X r10,G - X r11,G, ) (26)

[0117] According to the evolutionary strategy of formula (25), using the information of the current optimal individual to guide the search helps to accelerate the convergence speed. However, excessive dependence on the optimal individual may lead to a decrease in population diversity. According to the evolutionary strategy of formula (26), randomly selecting individuals for differential evolution, without relying on the optimal individual information, has a strong exploration ability, helps to maintain population diversity, and avoids the algorithm from converging to a local optimum prematurely. Based on this, the combination of the two evolutionary strategies can accelerate convergence by using the optimal individual information while maintaining population diversity through random exploration and avoiding falling into a local optimum; especially in the initial stage of the algorithm operation (t ≤ maximum iteration number * 2 / 3), the algorithm focuses more on exploring the search space to find potential optimal solution regions.

[0118] S420: If the iteration number t > maximum iteration number t max *2 / 3, individuals with the second preset ratio perform differential mutation according to formula (27), and the remaining individuals perform differential mutation according to formula (28). Optionally, the second preset ratio is 30% - 80%. For example, if the second preset ratio is 50%, then 50% of the individuals perform differential mutation according to formula (27), and the remaining 50% of the individuals perform differential mutation according to formula (28). Optionally, randomly select individuals with the second preset ratio to perform differential mutation according to formula (27), and the remaining individuals perform differential mutation according to formula (28); or sorted by fitness value, the individuals with the second preset ratio in the front perform differential mutation according to formula (27), and the remaining individuals perform differential mutation according to formula (28), which is not uniquely defined here.

[0119] V i,G = Xc1,G +F·(X c1,G -X i,G )+F·(X c2,G -X i,G ) (27)

[0120] V i,G =X i,G +F·(X best,G -X ri,G )+F·(X r9,G -X r10,G ) (28)

[0121] In the above embodiments, r9, r10, r11, and r12 are individual indices randomly selected from the G-th generation; ri represents the current individual index; F is the scaling factor of the difference vector; X best,G is the optimal individual of the G-th generation; X i,G is the i-th individual of the G-th generation; V i,G is the i-th mutated individual of the G-th generation; K is a positive integer. Generally, K ≥ 50. Optionally, the value range of K is 100 to 300. For example, in one embodiment, K is 150, that is, when the fitness value has not improved after 150 consecutive iterations, before returning to step S300, steps S410 and S420 are executed.

[0122] Based on this, in the later stage of the algorithm operation (t > 2 / 3 of the maximum number of iterations), the algorithm focuses more on using the information of the optimal solution already found for fine search to improve the solution accuracy. According to the strategy of formula (27), X c1,G and X c2,G are respectively constructed by using random individuals and the current optimal individual, which not only considers the overall distribution of the population, helps to maintain the diversity of the population, but also helps to accelerate the convergence speed. According to the strategy of formula (28), the current optimal individual is mainly used to accelerate the convergence speed for fine search, and at the same time, the random individuals X r9,G and X r10,G are partially combined to maintain the exploration ability to a certain extent and avoid the algorithm missing better solutions during the fine search process.

[0123] In some other embodiments, in combination with Figure 8 , when the fitness value has not improved after K consecutive iterations, before returning to step S300, step S400 includes the following steps:

[0124] S430: Sort all individuals according to the fitness value and divide them into the first population and the second population.

[0125] S440: The first population evolves according to the first differential evolution strategy.

[0126] S450: The second population evolves according to the second differential evolution strategy, and the first differential evolution strategy is different from the second differential evolution strategy.

[0127] S460: The individuals of the next generation are updated according to Equation (29). Among them, f(V i,G ) and f(X i,G ) represent the fitness values of the mutant crossover individual and the original individual respectively.

[0128]

[0129] Based on this, LRMHBA adopts a hybrid initialization strategy that combines LHS and random sampling. The hybrid strategy not only ensures the uniformity of the space but also increases the diversity of the population, enabling the initial population to achieve a good balance between quality and diversity. At the same time, sorting the first population and the second population according to fitness helps the algorithm adopt a differential strategy in the subsequent search process and converge to the optimal solution or approximate optimal solution faster.

[0130] In the following embodiments, r9, r10, r11, and r12 are individual indices randomly selected from the G-th generation; ri represents the current individual index; F is the scaling factor of the differential vector; X best,G is the optimal individual of the G-th generation; X i,G is the i-th individual of the G-th generation; V i,G is the i-th mutant individual of the G-th generation; X elite is the individual (see the individuals in the red circles in Figure 7 (c)) located in the top third preset proportion sorted by fitness value in the hybrid pool composed of the first sample and the second sample, that is, the elite individuals. The phased double-population co-evolution strategy that combines multiple differential evolution methods under the guidance of elite individuals selects different mutation strategies according to the evolution process and population adaptability, better balancing the exploration and exploitation of the algorithm and enhancing the global optimization ability of the algorithm. Optionally, the third preset proportion is 5% - 50%. For example, the third preset proportion is 5%, 10%, 15%, 20%.

[0131] In one of the embodiments, if the number of iterations t ≤ 2 / 3 of the maximum number of iterations t max *, the first differential evolution strategy is Equation (30), and the second differential evolution strategy is Equation (31).

[0132] V i,G = X c1,G + F·(X c1,G - X i,G ) + F·(X c2,G - X i,G ) (30)

[0133] V i,G = X r9,G + F·(X r10,G - X r11,G, ) (31)

[0134] In the first two-thirds of the iteration process, the overall population tends to conduct global exploration in a wide area to discover potential global optimal solutions, with emphasis on exploration rather than exploitation. Based on this, the first population has better adaptability, so the strategy of Equation (30) with relatively balanced exploration and exploitation capabilities is used for the first population. The second population has relatively weak adaptability. Therefore, at this stage, it can conduct global exploration in a wide area to discover potential global optimal solutions without relying on the positions of any other individuals in the group, helping the algorithm to jump out of the local optimum. So, Equation (31) with excellent global exploration capabilities is used for the second population.

[0135] In one embodiment, if the number of iterations t > the maximum number of iterations t max *2 / 3, the first differential evolution strategy is Equation (32), and the second differential evolution strategy is Equation (33).

[0136] V i = X i,G + F·(X best,G - X ri,G ) + F·(X r9,G - X r10,G ) (31)

[0137] V i,G = X c1,G + F·(X c1,G - X i,G ) + F·(X c2,G - X i,G ) (32)

[0138] In the last one-third of the iteration process, the population gradually begins to gather around the optimal solution. Therefore, emphasis should be placed on exploitation rather than exploration. So, Equation (31) with excellent local exploitation performance is applied to the first population, and Equation (32) is applied to the second population to balance a certain degree of exploration.

[0139] In some embodiments, in combination with Figure 8 , step S300 further includes the following steps:

[0140] S301: When the adaptive coefficient A satisfies |A| ≥ 1, search according to Equation (34) during the honey collection stage; otherwise, search according to Equation (35) during the honey collection stage. d i represents the distance between the prey and the current honey badger individual. In Equation (34), d i = X rand - X i ; in Equation (35), di = X prey -X i .

[0141] X new = X rand + Dir·r7·α·d i (34)

[0142] X new = X prey + Dir·r7·α·d i (35)

[0143] |A| decreases linearly with the number of iterations. In the early stage of iteration, |A| ≥ 1. The algorithm executes the random search strategy multiple times. Equation (34) is used to avoid premature aggregation of the population, enhance the traversal of honey badger individuals in the search space, and improve the global search ability.

[0144] Next, the CEC2017 test set is used to verify Figure 8 the effectiveness of LRMHBA in. CEC2017 is a test set containing a total of 29 functions, including 2 unimodel functions, 7 simple multimodel functions, 10 hybrid functions, and 10 composition functions, which can comprehensively evaluate the performance of the algorithm. Among them, four types of algorithms are selected for comparative analysis, specifically including: (1) classical highly cited algorithms and their variants: PSO, DE, WOA, AOA, GQPSO; (2) new algorithms and their variants in the past two years: PO, DBO, QHDBO; (3) the champion algorithm: LSHADE; (4) HBA algorithms and their variants: HBA, SaCHBA_PDN. The analysis experiments include test analysis on test functions in different dimensions, ablation experiments, and exploration and exploitation experiments. The population size of all algorithms is set to 100, and the maximum number of iterations is set to 100,000. Each algorithm runs independently 30 times. When the test dimensions are relatively low-dimensional (30 dimensions) and high-dimensional (100 dimensions), the optimal value, average value, and standard deviation of the test functions are calculated.

[0145] The data results when the test dimensions are 30 and 100 respectively are as Figure 9 shown. At 30 dimensions, LRMHBA ranked first in 18 functions, and at 100 dimensions, it was 24. Among all algorithms, it ranked first in both dimensions with average ranking values of 1.58 and 1.27 respectively. All the above data show that the LRMHBA algorithm performs better than other algorithms on the CEC2017 test set, especially in high-dimensional problems, and this advantage is more significant.

[0146] The convergence curve graphs of all the test algorithms are as Figure 10 and Figure 11 shown. In the case of 30 dimensions, only the convergence curve graphs of 20 functions in the test set are listed. From Figure 10 and Figure 11 it can be seen that LRMHBA found better average fitness values on most of the benchmark test functions. Compared with the three most competitive algorithms, LSHADE, HBA, and SaCHBA_PDN, LRMHBA generally has a faster convergence speed than LSHADE and HBA. Although the convergence speed of LRMHBA is slower than that of SaCHBA_PDN in the early stage, SaCHBA_PDN is more likely to fall into local optimal solutions. On the contrary, LRMHBA can gradually converge to a better value as the number of iterations increases, indicating that LRMHBA has stronger global optimization ability.

[0147] The data in Table 1 are the p-values of the significance statistics of the Wilcoxon Rank-Sum Test at 100 dimensions. If p < 0.05, it means there is no similarity between the two algorithms. We can use the symbols "+", "=", and "-" to represent that the performance of the LRMHBA algorithm is better than, equivalent to, and worse than other algorithms respectively. The statistical results show that except for two p-values greater than 0.05 in the comparison with LSHADE, the rest are less than 0.05, indicating that there are significant differences between the LRMHBA algorithm and other algorithms.

[0148]

[0149]

[0150] Table 1 Significance statistical results of Wilcoxon Rank-Sum Test at 100 dimensions

[0151] Figure 12 shows the box plots of all algorithms on the 30-dimensional CEC2017 benchmark test set. From this, it can be seen that although there are a few outliers for LRMHBA on some test functions, its box length is very short, and its median and average rank are better than those of other comparison algorithms on most test functions, indicating that the LRMHBA algorithm has very good stability.

[0152] Exploration means that the algorithm tries to access new regions to discover potentially better solutions, while exploitation means that the algorithm concentrates resources to conduct a detailed search in a known promising region in order to find the optimal solution within that region. Excessive exploration may lead to waste of computing resources, while excessive exploitation may cause the algorithm to converge prematurely to a sub-optimal solution and be unable to jump out of the local optimum. Therefore, an effective heuristic algorithm needs to find an appropriate balance between the two. Combining equations (36) to (38), this section adopts a dimensional diversity measurement method to better evaluate the ability of LRMHBA to balance exploration and exploitation.

[0153]

[0154] where Div represents the diversity value of the population; Div max represents the maximum diversity value in the iteration; N is the population size; D is the variable dimension; is the j-th dimensional position of the i-th individual; median(x j ) is the median of the j-th dimensional variables of all individuals in the population.

[0155] One unimodal function, six simple multimodal functions, seven hybrid functions, and six composite functions in the CEC2017 test set were selected to conduct experimental comparison tests on the two algorithms of LRMHBA and HBA. Each algorithm was run 30 times repeatedly, and the number of iterations was set to 2,000. The test results are as Figure 13 shown. It can be seen from Figure 13 that LRMHBA has a strong exploration ability in the initial stage of iteration, which is beneficial for it to find potential optimal solutions in different regions of the entire solution space. As the iteration progresses, LRMHBA gradually shifts from exploration to exploitation, carefully optimizing the found high-quality solutions, and finally concentrating near the global optimal solution. In addition, it can be seen that compared with the HBA algorithm, LRMHBA has a faster convergence speed, and the population quickly approaches the optimal solution, indicating that the Latin hypercube sampling method homogenizes the population, improves the optimization efficiency of the LRMHBA algorithm, and the random perturbation strategy and the phased double-population co-evolution strategy that combines multiple differential mutation methods better balance exploration and exploitation, improving the global optimization ability of the algorithm.

[0156] Furthermore, three map scenarios with a size of 100*100*3 km and increasing complexity were constructed, and algorithms that performed well in this problem, such as PSO, SSA, HHO, DBO, and CPO, as well as the original HBA algorithm and the improved version SaCHBA_PDN of HBA, were selected for comparative analysis of the three-dimensional path planning of UAVs. The population size of each algorithm was set to 100, and the maximum number of iterations was 500. The starting point of the map was (5, 5, 0.3), the target point was (90, 90, 0.8), and 8 intermediate path points were set. The mountain model and threat area parameters are shown in Table 2. To reduce the randomness of the algorithm, each algorithm was run independently 30 times, and the optimal value, average value, variance, and Friedman test ranking were introduced as evaluation indicators.

[0157]

[0158] Table 2 Mountain Model Parameters and Threat Model Parameters

[0159] The optimal path planning diagrams and average cost convergence curves in the three scenarios are shown respectively as Figure 14 and Figure 15 shown. From Figure 14 it can be seen that the LRMHBA algorithm found the shortest and smoothest optimal path in the simplest scenario 1. Although the shortest path was not found in scenarios 2 and 3, Figure 15 it shows that the minimum average cost value of path planning was obtained in all three scenarios, indicating that the LRMHBA has good stability. In addition, in scenarios 1 and the most complex scenario 3, the LRMHBA had a faster convergence speed than the original HBA algorithm and found a feasible path through fewer iterations. In scenario 2, the convergence speed of HBA was slightly better than that of LRMHBA, but it fell into a local optimum.

[0160] Table 3 summarizes the planning results under three scenarios, where the best value of each metric is shown in bold. LRMHBA found the optimal shortest path in Scenario 1, achieved the best average flight cost among the three scenarios, and most metrics outperformed the original HBA algorithm, which is consistent with the results of the optimal path planning diagram and the iteration curve diagram. The sparrow search algorithm has the best stability because it found a feasible path every time, while on the contrary, the HHO, PSO, and DBO algorithms have large averages and variances in all three scenarios for HHO, and large averages and variances in the relatively complex Scenarios 2 and 3 for PSO and DBO, indicating that these algorithms encountered more search failures during the search process in the corresponding scenarios, especially in complex maps. Generally speaking, LRMHBA performed the best, and its Friedman test ranking was also the first. Therefore, the proposed LRMHBA algorithm demonstrated excellent adaptability in both simple and complex three-dimensional obstacle environments and could plan an efficient and stable path.

[0161]

[0162] Table 3 Path Planning Results

[0163] Second, in combination with Figure 16 , this application provides a navigation device 20, including a memory 22, a processor 21, and a computer program 23 stored in the memory 22 and executable on the processor 21. When the processor 21 executes the computer program 23, it implements the three-dimensional path planning method for drones as described in any one of the above.

[0164] Among them, the navigation device 20 can be a mobile phone, a tablet computer, a wearable device, an augmented reality (AR) / virtual reality (VR) device, a notebook computer, an ultra-mobile personal computer (UMPC), a netbook, a personal digital assistant (PDA), etc., and is not uniquely limited here.

Claims

1. A three-dimensional path planning method for an unmanned aerial vehicle, characterized in that, The method includes the following steps: S100: Obtain a flight path planning model including an environmental model and a total cost function; S200: Use a honeycomb algorithm to solve the flight path planning model and initialize population parameters; S300: When the coefficient vector A satisfies |A|≥1, the search mathematical expression of the population in the mining stage is: X new = X rand + Dir·β·I·X rand + Dir·r2·α·d i ·|cos(2·π·r3)×(1 - cos(2·π·r4))|; Otherwise, the search mathematical expression in the mining stage is: X new = X prey + Dir·β·I·X prey + Dir·r2·α·d i ·|cos(2·π·r3)×(1 - cos(2·π·r4))|; Among them, X new is the updated position of the honey badger individual, and X prey is the global optimal position; α is the density factor; the food acquisition ability factor β ≥ 1; r2, r3, r4, r5 are random numbers between [0, 1]; I is the prey control intensity; Dir is the search direction factor; X rand is the honey badger individual randomly selected under the current iteration; A = 2mr5 - m, where m is the convergence factor; S400: When the iteration number t reaches the maximum iteration number t max , output X prey and its fitness value; otherwise, increment the iteration number t by 1 and return to step S300.

2. The method for three-dimensional path planning of an unmanned aerial vehicle according to claim 1, wherein When continuously iterating K times and the fitness value has not improved, before returning to step S300, step S400 further includes the following steps: S410: If the number of iterations t ≤ the maximum number of iterations t max *2 / 3, individuals with the first preset ratio mutate according to the following differential mutation strategy: V i,G = X i,G + F·(X best,G - X ri,G ) + F·(X r9,G - X r10,G ); The remaining individuals mutate according to the following differential mutation strategy: V i,G = X r9,G + F·(X r10,G - X r11,G, ); S420: If the number of iterations t > the maximum number of iterations t max *2 / 3, individuals with the second preset ratio mutate according to the following differential mutation strategy: V i,G = X c1,G + F·(X c1,G - X i,G ) + F·(X c2,G - X i,G ); The remaining individuals mutate according to the following differential mutation strategy: V i,G = X i,G + F·(X best,G - X ri,G ) + F·(X r9,G - X r10,G ); Among them, r9, r10, r11, and r12 are individual indices randomly selected from the G-th generation; ri represents the current individual index; F is the scaling factor of the difference vector; X best,G is the optimal individual of the G-th generation; X i,G is the i-th individual of the G-th generation; V i,G is the i-th mutated individual of the G-th generation; 3. The drone three-dimensional path planning method according to claim 1, characterized in that When continuously iterating K times and the fitness value has not improved, before returning to step S300, step S400 further includes the following steps S430: Sort all individuals according to the fitness value and divide them into a first population and a second population; S440: The first population evolves according to a first differential evolution strategy; S450: The second population evolves according to a second differential evolution strategy, and the first differential evolution strategy is different from the second differential evolution strategy; S460: Update the next-generation individuals according to the following formula: Among them, f(V i,G ) and f(X i,G ) represent the fitness values of the mutated and crossed individual and the original individual, respectively.

4. The UAV three-dimensional path planning method according to claim 3, characterized in that: If the number of iterations t ≤ the maximum number of iterations t max *2 / 3, the first differential evolution strategy is as follows: V i,G = X c1,G + F·(X c1,G - X i,G ) + F·(X c2,G - X i,G ); The second differential evolution strategy is: V i,G = X r9,G + F·(X r10,G - X r11,G, ); If the number of iterations t > the maximum number of iterations t max *2 / 3, the first differential evolution strategy is as follows: V i = X i,G + F·(X best,G - X ri,G ) + F·(X r9,G - X r10,G ); The second differential evolution strategy is: V i,G = X c1,G + F·(X c1,G - X i,G ) + F·(X c2,G - X i,G ); Among them, r9, r10, r11, and r12 are individual indices randomly selected from the G-th generation; ri represents the current individual index; F is the scaling factor of the difference vector; X best,G is the optimal individual of the G-th generation; X i,G is the i-th individual of the G-th generation; V i,G is the i-th mutated individual of the G-th generation; X elite is the individual ranked in the top third preset proportion according to the fitness value in the mixed pool composed of the first population and the second population.

5. The method for three-dimensional path planning of an unmanned aerial vehicle according to claim 1, wherein Step S300 further includes the following steps: S301: When the adaptive coefficient A satisfies |A|≥1, the search mathematical expression in the honey collection stage is: X new = X rand + Dir·r7·α·d i ; At this time, d i represents the distance between the prey and the current honey badger individual, d i = X rand - X i ; Otherwise, the search mathematical expression in the honey collection stage is: X new = X prey + Dir·r7·α·d i ; At this time, r7 is a random number between [0, 1]; d i = X prey - X i .

6. The method for three-dimensional path planning of an unmanned aerial vehicle according to claim 1, wherein, The environmental model is constructed by the following steps: S110: The reference terrain model in the environmental model is represented by the following formula: Where z1 is the height value at the terrain model point (x, y); k, a, b, c, d, e, f, g are coefficients used to control the digital terrain ratio and undulation; S120: The mountain model in the environmental model is represented by the following formula: where z2 is the mountain height value at the point (x, y): h i is the height peak of the i-th mountain, x i and y i are the central positions of the i-th mountain; a i and b i are the slope parameters of the i-th mountain peak in the x and y directions.

7. The method for three-dimensional path planning of an unmanned aerial vehicle according to claim 6, wherein Step S110 specifically includes the following steps: S111: Obtain the two-dimensional vector map corresponding to the reference terrain model; S112: Obtain the shape feature and height feature of the two-dimensional vector map; S113: Calculate the coefficients k, a, b, c, d, e, f, and g according to the shape feature and the height feature.

8. The method for three-dimensional path planning of an unmanned aerial vehicle according to claim 7, wherein: The acquisition of the height feature adopts the following steps: S1121: Divide the linear chirped optical signal into a laser detection signal and a local oscillator signal, where the splitting ratio range of the laser detection signal and the local oscillator signal is in the range of 90:10 to 99:1, the local oscillator signal is used for echo signal demodulation, and the laser detection signals of different frequencies are directly emitted outward; S1122: Receive echo signals of different frequencies through a plurality of silicon photomultiplier tubes arranged in an array; S1123: Process the echo signals to generate height features at each terrain model point.

9. The method for three-dimensional path planning of an unmanned aerial vehicle according to any one of claims 1 to 8, characterized in that, The total cost function is constructed by the following steps: S130: Construct a flight distance cost function: Among them, (x i , y i , z i ) is the coordinate of the i-th path point; n is the number of path points; S140: Construct a flight height cost function: where z i is the height of the i-th point; is the average height of the path, S150: Construct a flight turning cost function: Among them, θ i is the included angle between adjacent path segments; Ф is the maximum turning angle; b i and b i+1 are the vector of the i-th and (i + 1)-th path segments; |b i | and |b i+1 | are the lengths of b i and b i+1 respectively. S160: Construct a terrain constraint cost function: d i = z i - h i ; z i is the height of the i-th path point; h i is the height of the terrain at this point; d i is the vertical distance from the i-th path point to the terrain surface; SD is the set terrain safety distance; S170: Construct an obstacle threat cost function: where (x i , y i ) are the coordinates of the i-th path point; (x k , y k ) are the coordinates of the center of the threat area; d ik is the perpendicular distance from the i-th path point to the k-th threat center area; R k is the radius of the k-th threat area; SF is the safety distance coefficient of the threat area S180: Construct the total cost function: where, μ L to μ T are the weight coefficients of the respective cost functions, and satisfy μ L + μ H + μ C + μ M + μ T = 1.

10. A navigation device, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that: When the processor executes the computer program, it implements the three-dimensional path planning method for an unmanned aerial vehicle according to any one of claims 1 to 9.

Images