CN103913172B - A kind of it is applicable to the paths planning method of aircraft under complicated low latitude - Google Patents

Abstract

The invention discloses and a kind of be applicable to the paths planning method of aircraft under complicated low latitude, belong to flight management method and system technical field.Described method initially sets up aircraft path planning model；Then use adaptive differential target evolution algorithm that single rack aircraft is carried out path planning.Based on described single rack aircraft paths planning method, the present invention also provides for a kind of multi rack aircraft paths planning method.The present invention is a kind of conflict Resolution method that multi-machine collaborative is evolved, and operation efficiency is high, it is possible to meet the demand in pre-planning path；Disclosure satisfy that the unit under the conditions of low altitude airspace free flight and multimachine path planning demand；There is certain feedback arrange, it is possible to while ensureing the safety of path planning, also obtain higher flight efficiency.

Description

A kind of it is applicable to the paths planning method of aircraft under complicated low latitude

Technical field

The present invention relates to a kind of flight management method and system, be a kind of all purpose aircraft being applicable to complicated low altitude airspace The method of aircraft path planning.

Background technology

Low altitude airspace refers to the spatial domain scope of the highest less than 1000 meters in principle, and in this region, activity is the most general Aviation aircraft.Low altitude airspace circumstance complication, common aero vehicle can free flight (be to be navigated by the United States Federal in this region Empty office (FAA, Federal Aviation Administration) proposes, and it is a kind of brand-new aviation theory, gives flight Personnel select the right of the state such as air route and the flight speed of a ship or plane the most in real time, i.e. carry out autonomous management), autonomy is strong.At certain A little scene (such as emergency disaster relief) density of getting off the plane are high, and safety and the high efficiency of flight are difficult to ensure that.Along with China progressively opens Putting low altitude airspace, General Aviation certainly will be flourish, for ensureing the safety and efficiency that flight runs, closes the planning of Flight device The path of reason economy is particularly significant.

Path planning refers in given planning space, finds movable body and arrives impact point from starting point and meet constraint Optimum or the most viable track of condition and certain performance indications.Flight device is carried out path planning and is to ensure that aircraft flight Safety and put forward high efficiency important means, the at present both at home and abroad research to this field mainly includes traditional planning algorithm and intelligence rule The method of calculating.Traditional planning algorithm includes dynamic programming algorithm, method in optimal control and derivative method of correlation etc., and this kind of method to be carried out greatly Iterating of scale, lack of wisdom guidance capability, amount of calculation is the biggest, and it is long that algorithm calculates the time, is not suitable for low absolutely empty The free flight in territory, and robustness is poor；Intelligent planning algorithm includes heuristic search search, ant group algorithm, neutral net Algorithm and genetic algorithm.Owing to the essence of path planning is np problem.Both along with the growth of problem scale, computation complexity is drastically Strengthen.Genetic algorithm has proved to be the effective ways of solution np problem because of its heuristic search operator and good optimizing ability. The essence of the path planning problem of common aero vehicle is again a series of conflicting objective optimizations, uses multiple-objection optimization Algorithm is more excellent than to the single object optimization weighting of multiple targets traded off.

Summary of the invention

The present invention is to solve problems of the prior art, it is provided that a kind of be applicable to the road of aircraft under complicated low latitude Footpath planing method, Multipurpose Optimal Method worked in coordination with in term, in described paths planning method, it is assumed that (1) all aircraft all exist Sustained height layer flies with friction speed, and in each comfortable flight course, flight speed does not changes；(2) every frame aircraft is all Determine respective beginning and end；(3) low altitude airspace flight environment of vehicle is it is known that include prohibited flight area position and scope, prestige Side of body point position and radiation radius；(4) general purpose vehicle physical property limits and includes maximum deflection angle and maximum flying distance；Described Paths planning method specifically includes following steps:

The first step, sets up aircraft path planning model；

Second step, uses adaptive differential target evolution algorithm that single rack aircraft is carried out path planning, and concrete steps are such as Under:

S1: determine planning space, starting point S and impact point E, utilizes heuristic initialization algorithm to initialize population, and complete Become the individual population pretreatment converted to absolute value coordinate individuality of contrast pole coordinate；

S2: according to headway v and sampling time step delta t, population at individual is completed track points discretization；

S3: according to the aircraft path planning model set up in the first step, utilize the individual discretization of initial population Track points carries out flight path constraint checking and desired value calculates, and calculates constrained binding occurrence sum err_ of each individuality ind；

S4: carry out parent population carrying out genetic cross and mutation operation according to adaptive differential operator, form sub-population；

S5: complete sub-population generation of neutrons individuality track points discretization, bring the track points of discretization into problem model, complete With individuality constraint, target function value and the calculating of constrained binding occurrence sum err_ind；

S6: parent is formed with filial generation and mixes group, and mixing group implements the quick non-dominated ranking improved, and determines that mixing group is each The individual Pareto layer being distributed；

S7: the crowding distance that improves of individuality and the Crowing Mechanism to each layer of group of mixing, according to dominance relation and individual Body crowding is chosen the individuality of Population Size number from mixing group and is formed new parent population；

S8: judge whether evolutionary process terminates, if reaching the algebraically of greatest iteration, proceeds to S9；Otherwise proceed to S4；

S9: select individuality in non-dominant layer from population, its flight path represented is required flight path；

Initialization population described in step S1, uses two kinds of initial population individuality expression waies, i.e. contrast pole coordinate to compile Code mode and absolute value codes co-ordinates mode, comprise contrast pole Coordinate generation initial chromosome and initial population individuality Coordinate Conversion Two steps:

(S101) by formula (26) acquisition contrast pole coordinate, (r, θ, value z) represents the three-dimensional of its correspondence to flight path node initializing Flight path node, wherein rand () represents and is uniformly distributed value, the lower limit Minseg of r and the upper limitRepresent between adjacent track points The shortest permission distance, the bound of turn_ θ is relevant with aircraft mobility, the θ obtained_iPositive and negative represent next flight path section The yawing moment of more previous flight path section, is just counterclockwise, is negative clockwise；

$\left\{\begin{array}{c}{r}_i=rand(Minseg,\frac{MaxL}{N_w})\\ {\mathrm{\θ}}_i=rand(-turn\_\mathrm{\θ},turn\_\mathrm{\θ})\\ z_i=rand\left(\min \right(start.z,end.z),\max (start.z,end.z\left)\right)\end{array}\right.---\left(26\right)$

Wherein, (start.z, end.z) represents the height coordinate value of starting point and impact point.

(S102) initial population individuality coordinate system conversion:

After utilizing contrast pole coordinate method to generate initial population, utilize contrast pole coordinate and the corresponding relation of absolute value coordinate, By formula (27) and formula (28), the individual two-dimentional flight path of the contrast pole coordinate representation to generating carries out coordinate system conversion, will generate Contrast pole coordinate initial population change into the initial population of absolute value coordinate representation；First flight path node W₁Absolute value sit Mark (x₁,y₁,z₁) by its contrast pole coordinate (r₁,θ₁,z₁) obtained by formula (27) and formula (28), wherein (start.x, start.y) (end.x, end.y) is the two-dimensional coordinate of starting point and impact point；

$\left\{\begin{array}{c}\mathrm{\α}=arctan\left(\frac{end.y-start.y}{end.x-start.x}\right)\\ x_1=start.x+r_1*cos(\mathrm{\α}+{\mathrm{\θ}}_1)\\ y_1=start.y+r_1*sin(\mathrm{\α}+{\mathrm{\θ}}_1)\\ z_1=z_1\end{array}\right.---\left(27\right)$

Based on the W tried to achieve₁(x₁,y₁,z₁), obtain follow-up track points W_iAbsolute value coordinate (x_i,y_i,z_i):

$\left\{\begin{array}{c}sum\mathrm{\θ}={\Σ}_{k=1}^{i+1}{\mathrm{\θ}}_k\\ x_i=x_{i-1}+r_{i+1}*\cos (\mathrm{\α}+sum\mathrm{\θ})\\ y_i=y_{i-1}+r_{i+1}*\sin (\mathrm{\α}+sum\mathrm{\θ})\\ z_i=z_i\end{array}\right.---\left(28\right)$

Wherein, i=2 ..., N-1；

Population at individual flight path node discretization described in step S2, specifically refers to:

On the premise of supposing that aircraft flies at a constant speed with speed v, by flight path node discretization, obtain carrying out constraint checking Discretization point sequence dis{S, P₁,P₂,…,P_N-2,E}；With by two continuous print track points P_i(x_i,y_i,z_i) and P_I+1(x_I+1, y_I+1,z_i+1) flight path section l that constitutes_iAs a example by, i=0 ..., N-1, flight path section l_iDiscrete point sequence disp_l_i{P_i1(x_i1,y_i1, z_i1),P_i2(x_i2,y_i2,z_i2),…,P_im(x_im,y_im,z_im) coordinate figure obtained by following formula:

${\mathrm{\θ}}_i=arctan\left(\frac{y_{i+1}-y_i}{x_{i+1}-x_i}\right)---\left(31\right)$

$m=floor\left(\frac{\sqrt{{(x_{i+1}-x_i)}^2+{(y_{i+1}-y_i)}^2+{(z_{i+1}-z_i)}^2}}{\mathrm{\Δ}t*v}\right)---\left(32\right)$

WhereinFor flight path section l_iThe deviation angle of in the vertical direction, θ_iFor flight path section l_iRelative on two dimensional surface In the deviation angle of x-axis forward, Δ t is step sample time, and floor () is downward bracket function；

Population adaptive differential genetic cross and mutation operation in described step S4, form sub-population, particularly as follows:

S401: be evaluated initial population, i.e. calculates the target function value of each individuality in initial population；

S402: carry out mutation operation, to each target individualThe variation being produced its correspondence by formula (33) is individual

$v_i^{G+1}=x_{r1}^G+F*(x_{r2}^G-x_{r3}^G)---\left(33\right)$

Wherein:WithSubscript r₁, r₂And r₃It is the different number randomly choosed, and and object vector's Subscript i is the most different；F is zoom factor, for the magnification level of control deviation variable；

S403: by formula (34), by target individualWith its variation individualityCarry out intersecting and operate, then obtain correspondence Intersect individualAfter the intersection of all target individual has operated, its intersection individuality forms interim population；

$u_{ij}^{G+1}=\left\{\begin{array}{cc}{v}_{ij}^{G+1},& if(rand\≤CR)orj=q\\ x_{ij}^G,& if(rand>CR)andj\≠q\end{array}\right.---\left(34\right)$

Wherein: CR for intersect the factor, be individuality carry out intersect operation control probability, q be one between [1, D] with The integer that machine produces, to guarantee to intersect individualityAt least can be from variation individualityObtain one-dimensional parameter；Letter D refers to one Constant；

The quick non-dominated ranking of the described improvement in step S6, specific as follows:

If a Yu b is two individualities in population, the quicksort operator improved under constraints is used to complete individual a and b The comparison of dominance relation, flow process is as follows:

S601: obtain individual a and individual b each constraint with self related constraint binding occurrence sum a.err_ind with b.err_ind；

S602: judge whether individual a with b breaks a contract, the most respective err_ind value and the magnitude relationship of 0, if a.err_ind Be both greater than 0 with b.err_ind, i.e. a and b violates constraints, is infeasible solution, then go to S603；If only one of which The err_ind of body is more than 0, then go to S604；If two individual err_ind values are equal to 0, then go to S605；

S603: if the err_ind of two individualities is different, go to S603.1；If two individual err_ind are identical, go to S603.2；

S603.1: by the individual principle that the individual domination binding occurrence sum that binding occurrence sum is big is little, if a.err_ind > b.err_ind, the then individual individual b of a domination, otherwise, then the individual individual a of b domination；

S603.2: individual a with b dominance relation not may compare, and does not arranges；

S604: violate the individual principle of constraints by the individual domination not violating constraints, if b.err_ind > 0 And a.ind_err=0, the then individual individual b of a domination；If a.err_ind > 0 and b.ind_err=0, the then individual individual a of b domination；

S605: the quick non-dominated ranking improved；

The quick non-dominated ranking method operator step improved is as follows: in population each individual p have two parameters np and Sp, the individual amount of the individual p of domination during wherein np is population, Sp is the group of individuals arranged by p in population,

SS1: find the individuality of np=0 in population, and be saved in current collection F1；

SS2: interim set H composes sky；For each individual i in current collection F1, investigate the individual collections that it is arranged Si, subtracts 1 by the nj of individual j each in Si, if nj-1=0, then individuality j puts into set H；Now set H is that current layer is individual The individual collections of next layer of set；

SS3:F1 is the 1st domination layer individual collections, individual identical non-dominated ranking irank in giving this layer；

SS4: with H as current collection, returns step SS2, until the sequence of whole population is complete.

If on the object function of jth optimization aim, had | f_j(a)-f_j(b) |≤δ, then it is assumed that individual a and b is at f_jOn Value equal, if individual a and b has such relation on all of target, then individual a with b does not arranges；Wherein, f_j(a) and f_j(b) represent respectively individual a and individual b in jth object function f_jFunctional value；

Crowding distance for the improvement described in step S7 calculates, with gathering individual on population a certain Pareto layer F Distance calculate as a example by, be given improve crowding distance computational methods:

S701: calculate individual number NF of population a certain Pareto layer F；

S702: the crowding distance of each individuality in F is composed null value, it may be assumed that F [i] .dis=0, i=1 ..., popsize； Popsize is number individual in Population Size, i.e. population；

S703: to each object function f_j, j=1 ..n, obtain functional value F [i] .f of each individuality_j, by F [i] .f_j's Individuality in F is ranked up by size；

S704: according to the ordering scenario that F layer is individual, by infinitely great for crowding distance assignment individual for boundary point, by remaining The crowding distance of body is pressedAssignment, i=2 ..., NF-1, i.e. with body one by one The diagonal sum of besieged rectangle represents the crowding distance of this individuality；Wherein, f_kmaxAnd f_kminRepresent kth function respectively Maximum and minima.

Based on described single rack aircraft paths planning method, the present invention also provides for a kind of multi rack aircraft path planning side Method, comprises the steps:

Step A: each frame aircraft is carried out heuristic initialization population:

A1: determine planning space, starting point S and impact point E, utilizes heuristic initialization algorithm to initialize population, and complete Become the individual population pretreatment converted to absolute value coordinate individuality of contrast pole coordinate；

A2: after all aircraft populations all complete heuristic initialization, enters step B；

Step B: to each frame aircraft population:

B1: according to headway v and sampling time step delta t, population at individual is completed track points discretization；

B2: according to aircraft path planning model, utilizes the individual discretization track points of initial population to carry out flight path about Bundle checking and desired value calculate, and calculate institute's Constrained sum err_ind of each individuality；

B3: carry out parent population carrying out genetic manipulation according to adaptive differential operator, form sub-population；

B4: complete offspring individual track points discretization, bring the track points of discretization into problem model, complete with individuality about Bundle, target function value and the calculating of institute's Constrained sum err_ind；

B5: parent is formed with filial generation and mixes group, mixing group implements the quick non-dominated ranking improved, determines mixing group The Pareto layer of individual distribution；

B6: represent individual selection mechanism according to population, the representative selecting this aircraft population from mixing group is individual；

B7: after all aircraft populations all complete the representative individual selection of mixing group, enters step C；

Step C: to each frame aircraft population:

C1: receive the representative individuality that remaining each aircraft mixing group selects, carries out this aircraft mixing all individualities of group Collaborative collision avoidance binding occurrence calculates, and obtains mixing institute's Constrained sum err_all of each individuality in group；

C2: the crowding distance improving the individuality mixing each layer of group calculates and Crowing Mechanism, according to dominance relation And individuality crowding is from mixing the parent population that the individuality composition choosing Population Size number group is new；

C3: judge whether evolutionary process terminates, if meeting end condition, proceeds to C4；Otherwise proceed to step B；

C4: select its flight path represented individual from population in non-dominant layer and be required flight path.

It is an advantage of the current invention that:

1, this method is a kind of conflict Resolution method that multi-machine collaborative is evolved, and operation efficiency is high, it is possible to meet pre-planning road The demand in footpath；

2, the unit under the conditions of this method disclosure satisfy that low altitude airspace free flight and multimachine path planning demand.

3, this method has certain feedback setting, it is possible to also obtain higher while ensureing the safety of path planning Flight efficiency.

Accompanying drawing explanation

Fig. 1 is aircraft path planning model；

Fig. 2 is for violating body of a map or chart flight path schematic diagram；

Fig. 3 A～Fig. 3 D is four kinds of situations that may be present between Duan Yuyi no-fly zone of a flight path；

Fig. 4 is the threatening area i threat coefficient calculations schematic diagram to route segment j；

Fig. 5 is the flow chart of adaptive differential multi-objective Evolutionary Algorithm；

Fig. 6 A and Fig. 6 B is respectively the same flight path represented with contrast pole coordinate and cartesian coordinate；

Fig. 7 is the flight path node discretization schematic diagram of the i-th flight path section；

Fig. 8 is multi-aircraft trajectory planning flow chart based on collaborative difference multi-objective Evolutionary Algorithm；

Fig. 9 is the collision detection schematic diagram between aircraft's flight track.

Detailed description of the invention

Below in conjunction with drawings and Examples, the present invention is described in further detail.

First the several definition to relating to carry out as described below:

1, flight collision and flight collision: the distance between two frame aircraft is less than conflict threshold (during collision threshold), Then think that this two framves aircraft exists collision risk (risk of collision).

2, individual: each individuality all represents a potential solution of problem, and each individuality is all evaluated quality, and with certain Probability likely carries out genetic manipulation, evolves.Individual also referred to as chromosome.

3, chromosome coding: it is treatable that the feasible solution of a problem is transformed into algorithm institute from its solution space by coding exactly The conversion method of search volume.Encode the performance for intelligent optimization algorithm based on genetic algorithm and population diversity impact The biggest.

4, population: maintain the colony's set being made up of a group individuality to cry population in genetic algorithm, individuality in population Number is called Population Size.

5, parent and filial generation: population was parent before genetic manipulation, the new population obtained after genetic manipulation For its filial generation.

6, collaborative multi-target evolution: a kind of Multipurpose Optimal Method based on coevolution.

7, the dominance relation that multi-target evolution is individual: give and minimize multi-objective optimization question, if p and q is evolution colony The individuality that in pop, any two is different, if meet simultaneously: 1. to all sub-goals, p is poor unlike q, i.e. f_k(p)≤f_k(q), k= 1,2,…,m；The most at least there is a sub-goal, make p better than q, { 1,2 .., m} make f i.e. to there is k ∈_k(p)<f_k(q).Then claim p Domination q.P is (non-dominated) of non-dominant, and q is (dominated) arranged.And if only if exist i ∈ 1, 2 .., m}, make f_i(p)<f_i(q), and { 1,2 .., m} make f to there is j ∈_j(p)<f_jQ (), then claim p with q not arrange.

8, Pareto optimal solution: multi-objective problem exists multiple target, owing to there is conflict between target and cannot comparing, Therefore be difficult to search out an optimal solution and make all object functions all reach optimum, it being usually present a disaggregation, this solves concentration Solution cannot compare good and bad for all object functions, i.e. pareto solves.It is defined as (puts up a question and entitled minimizes multiple target Problem): for feasible solution x, and if only if exists feasible solution x* and makes: 1. F (x)≤F (x*), F (x) feeling the pulse with the finger-tip scalar functions vector set； The most at least there is a f_j(x)<f_j(x*), wherein f_jX () is one of object function.

The present invention provides a kind of and is applicable to the paths planning method of aircraft under complicated low latitude, it is assumed that

(1) all aircraft all fly with friction speed at sustained height layer, and flight speed in each comfortable flight course Do not change.(2) every frame aircraft all determines respective beginning and end.(3) low altitude airspace flight environment of vehicle is it is known that include Prohibited flight area position and scope, threat point position and radiation radius.Wherein prohibited flight area is rectangular area, threatens district Territory be with threaten point as the center of circle, the radiation radius border circular areas as radius.(4) general purpose vehicle physical property limits and includes maximum Deflection angle and maximum flying distance.

Described paths planning method, specifically includes following steps:

The first step, sets up aircraft path planning model；

In conjunction with accompanying drawing 1, for the flight environment of vehicle of simulated flight device, the present invention is from complicated low altitude airspace environment, aircraft machine Dynamic performance and the several aspect of flight path performance indications set up aircraft path planning model.

Wherein, the standard can not run counter in the influence factor that flight path generates, such as self mobility etc., it is referred to as to be about Bundle；And some probabilistic influence factors, such as potential threat Regional Risk, if flight track superscale is little, aircraft has been Entirely can fly under undertaking certain risk premise, then it is called target.Every frame aircraft rule of flight in low altitude airspace The path drawn must is fulfilled for environmental constraints and the physical property of himself constraint, and realizes safe and efficient target.

Described environmental constraints is specific as follows:

(1) constraint diagram: aircraft must fly to guarantee flight safety in known body of a map or chart.Article one, flight road The map binding occurrence in footpath can be obtained by following formula:

$c_{map}={\mathrm{\&Sigma;}}_{i=0}^N{c}_i---\left(1\right)$

$c_i=\left\{\begin{array}{cc}1,& \mathrm{if}(x_{\mathrm{map}}^l\&le;x_i\&le;x_{\mathrm{map}}^u)^(y_{\mathrm{map}}^l\&le;y_i\&le;y_{\mathrm{map}}^u)\\ 0,& \mathrm{else}\end{array}\right.---\left(2\right)$

Wherein, C_mapRefer to a flight path in the constraint violation value of ground constraint diagram, C_iRefer to the node on flight path I in map restricted index, C_iWhen=0, node i is positioned at map, does not breaks a contract, c_iWhen=1 not in map, promise breaking.(x_i, y_i) it is the coordinate of path point i, N is path point number,WithSeat for graph region definitely Parameter bound.As in figure 2 it is shown, a certain flight path starting point is S, terminal is E, includes path point W therebetween₁、W₂、W₃、W₄If, There is path point W₁And W₂Outside body of a map or chart, then regard as violating the path point of body of a map or chart constraint.

(2) prohibited flight area: can there is extremely hazardous region known to some in flight environment of vehicle, aircraft once enters Enter, then face considerable risk.Prohibited flight area (PFZ) represents the region that flight environment of vehicle risk is very big, information is unknown, is to fly The off-limits region of row device.It is a rectangular area that the present invention sets prohibited flight area abstract, and its major parameter is four Apex coordinate position, if the number of prohibited flight area is N in flight space territory_P, note i-th no-fly zone is PFZ_i(i=1, 2,...,N_P), its vertex position is respectively

The flight path that one flight path (also referred to as flight path) enters in the present invention all prohibited flight areas (no-fly zone) is long Degree sum is as this flight path binding occurrence in the constraint of no-fly zone, and when binding occurrence is zero, this flight path is no-fly without prejudice to forbidding entering District retrains.Owing to no-fly zone information is given with two dimensional form, the present invention only calculates flight path projection on the two dimensional surface of no-fly zone Length enters no-fly zone length as it.Article one, the prohibited flight area binding occurrence c of flight track_PFzCan be obtained by formula (3):

$c_{PFZ}={\mathrm{\&Sigma;}}_{i=1}^{Np}{\mathrm{\&Sigma;}}_{j=1}^{N-1}{l}_{ji}---\left(3\right)$

$l_{ji}=\left\{\begin{array}{ccc}0,& if& seg(i,j)=0\\ \sqrt{{(x_{ji}^1-x_j)}^2+{(y_{ji}^1-y_j)}^2},& if& seg(i,j)=1\\ \sqrt{{(x_{ji}^1-x_{ji}^2)}^2+{(y_{ji}^1-y_{ji}^2)}^2},& if& seg(i,j)=2\\ \sqrt{{(x_{j-1}-x_j)}^2+{(y_{j-1}-y_j)}^2},& if& seg(i,j)=-1\end{array}\right.---\left(4\right)$

L in its Chinese style (4)_jiFor the jth flight path section of the flight path two-dimentional flight path length in i-th no-fly zone, Np is The number of prohibited flight area.(i j) returns there may be alternately of a Duan Yuyi rectangle prohibited flight area of flight path to seg One of four kinds of situations, as shown in Fig. 3 A～3D:

(a) flight path section j completely outside i-th no-fly zone, now l_ji=0；As shown in Figure 3A, described flight path section refers to Track points W_j‐i(x_j-1,y_j-1) and track points W_j(x_j,y_jBetween) one section.

(b) such as Fig. 3 B, an end points (x of flight path section j_j,y_j) in i-th no-fly zone, i.e. flight path section j is forbidden with i-th There is an intersection point in a certain border of flight rangeNow

(c) such as Fig. 3 C, two end points (x of flight path section j_j-1,y_j-1) and (x_j,y_j) in i-th no-fly zone, i.e. flight path section j It is fully located in i-th no-fly zone, now

D () such as Fig. 3 D, there is intersection point with certain two border of i-th no-fly zone in flight path section jWithI.e. Flight path section j passes through i-th no-fly zone completely, now

(3) terrain obstruction constraint；

The terrain obstruction of flight range is primarily referred to as the height relief of flight range landform, such as mountain peak, the situation in mountain valley.Fly Row device can run into the barriers such as the high mountain that can not leap in flight course, once hits mountain, can produce serious consequences such as falling, institute Should forbid that when carrying out trajectory planning aircraft direct collision is to terrain obstacle.For a certain bar flight path, its landform Obstacles Constraints can be obtained by formula (6):

$c_{terrain}={\&Sigma;}_{i=0}^N{c}_{ter\_i}---\left(6\right)$

$c_{ter\_i}=\left\{\begin{array}{cc}1,& \begin{array}{cc}if& z_i\&le;terrain(x_i,y_i)\end{array}\\ 0,& else\end{array}\right.---\left(7\right)$

Wherein, terrain (x_i,y_i) it is reentry point (x_i,y_i) the geopotentia function at place.

Described own physical constraint is specific as follows:

(1) deflection angle constraint:

Owing to the physical property of common aero vehicle self limits, when it need to change course, deflection angle should be made not surpass Cross its maximum deflection angle θ allowed.Article one, the deflection angle binding occurrence c of flight path_turnCan be obtained by formula (8):

$c_{turn}={\&Sigma;}_{j=1}^{N-1}{c}_j---\left(8\right)$

$c_j=1-\frac{\frac{a_j^T{a}_{j+1}}{\left|a_j\right|\left|a_{j+1}\right|}}{\cos \mathrm{\&theta;}}---\left(9\right)$

Wherein, a_jRepresent the vector of flight path section j determined by jth+1 node and jth node: a_j=(x_j+1-x_j,y_j+1- y_j)。

(2) maximum angle of pitch restricted model；

Similar with maximum deflection angle, affected by the mobility of aircraft self, aircraft is in the vertical plane of flight Can not climb with arbitrarily angled and dive.Article one, the angle of pitch binding occurrence c of flight path_pitchCan be obtained by formula (10):

$c_{pitch}={\mathrm{\&Sigma;}}_{j=1}^{N-1}{c}_{p\_j}---\left(10\right)$

$c_{p\_j}=\left\{\begin{array}{cc}0,& if(\mathrm{\&beta;}-{\tan }^{-1}\&lsqb;\frac{|z_{j+1}-z_j|}{\left|\sqrt{{(x_{j+1}-x_j)}^2{(y_{j+1}-y_j)}^2}\right|}\&rsqb;)<0\\ \mathrm{\&beta;}-{\tan }^{-1}\&lsqb;\frac{|z_{j+1}-z_j|}{\left|\sqrt{{(x_{j+1}-x_j)}^2+{(y_{j+1}-y_j)}^2}\right|}\&rsqb;,& else\end{array}\right.---\left(11\right)$

Wherein, β represents the maximum angle of pitch of aircraft physical property regulation；(x_j,y_j,z_j) and (x_J+1,y_J+1,z_j+1) point Biao Shi jth flight path node and+1 flight path node of jth.

(3) minimum flight path section restricted model；

Can be safer and fly more accurately in order to ensure aircraft, carry out, at aircraft, flight attitudes such as turning round, climb Before and after adjustment, what aircraft should keep a segment distance flies nonstop to state to complete its pose adjustment process.This flies nonstop to state The step distance that at least should reach is known as the minimum flight path section flying distance constraint of aircraft, represents with Minseg herein Minimum flight path section.Article one, the flight path section constraint of flight path can be obtained by formula (12):

$c_{seglen}={\mathrm{\&Sigma;}}_{j=1}^{N-1}{c}_{sl\_j}---\left(12\right)$

$c_{sl\_j}=\left\{\begin{array}{cc}0,& if{l}_j-Minseg\&GreaterEqual;0\\ 1,& else\end{array}\right.---\left(13\right)$

$l_j=\sqrt{{(x_{j+1}-x_j)}^2+{(y_{j+1}-y_j)}^2+{(z_{j+1}-z_j)}^2}---\left(14\right)$

Wherein, l_jRepresent the physical length of jth flight path section, (x_j,y_j,z_j) and (x_j+1,y_j+1,z_j+1) represent jth respectively Flight path node and+1 flight path node of jth.

(4) maximum flying distance (being also ultimate run) constraint: aircraft is in whole flight course, and its voyage can be subject to , therefore there is ultimate run constraint, be designated as maxL in fuel oil and the restriction of time dispensing.The ultimate run binding occurrence of one paths c_TdisCan be obtained by formula (15):

$c_{Tdis}=\left\{\begin{array}{cc}0,& ifMaxL-{\&Sigma;}_{j=1}^{N-1}lj\&GreaterEqual;0\\ 1,& else\end{array}\right.---\left(15\right)$

$l_j=\sqrt{(x_{j+1}-x_j)+(y_{j+1}-y_j)+(z_{j+1}-z_j)}---\left(16\right)$

During for aircraft path planning, the oil consumption of aircraft, flight time and safety coefficient all should be ensured that.The most Target is made up of two conflicting targets, is minimum risk target and shortest path target respectively.

In the present invention, potential threat region is modeled as border circular areas, and based on threatening operating radius, calculating aircraft flight path is subject to Threat index size to all unknown threatening areas.When aircraft is in the range of threatening area, threaten district center to flying The threat index of row device is strong, and aircraft is further away from the center of threat, and threat index is the most weak.Can be with the threat index of each flight path section Sum represents the threat index of whole flight path, minimum risk desired value f of a flight path_riskCan be calculated by formula (17):

$f_{risk}={\mathrm{\&Sigma;}}_{j=1}^{N-1}risk\left({\mathrm{seg}}_j\right)---\left(17\right)$

With by track points w_j-1With track points w_jAs a example by the threat cost of the jth flight path section determined calculates, choose this flight path section In six points calculate i-th potential threat region (such as the threatening area i) threat index to jth flight path section, such as formula (18) Shown in, Fig. 4 is corresponding schematic diagram:

${\mathrm{seg}}_{rishi}\left({\mathrm{seg}}_j\right)=\frac{l_j}{6}*{\&Sigma;}_{i=1}^{Nsite}\&lsqb;P_i\left(w_{j-1}\right)+P_i\left(P_j\right(\frac{1}{6}\left)\right)+P_i\left(P_j\right(\frac{2}{6}\left)\right)+P_i\left(P_j\right(\frac{3}{6}\left)\right)+P_i\left(P_j\right(\frac{4}{6}\left)\right)+P_i\left(P_j\right(\frac{5}{6}\left)\right)\&rsqb;---\left(18\right)$

Wherein, Seg_riski() is function name, and parameter is segj, represent jth flight path section in i-th potential risk district pair Value-at-risk；l_jRepresenting the length of jth flight path section, Nsite is potential threat district number, P_iP () is for threatening the some i prestige to certain point p Side of body probability, can be calculated by formula (19):

$P_i\left(p\right)=\left\{\begin{array}{cc}\frac{K_i}{\left(d_p\right(i\left)\right)},& \begin{array}{cc}if& d_p^i<R_{safei}\end{array}\\ 0,& else\end{array}\right.---\left(19\right)$

$d_p\left(i\right)=\sqrt{{\left(P_j\right(p).x-Riski.x)}^2+{\left(P_j\right(p).y-Riki.y)}^2},p=0,\frac{1}{6},\mathrm{...},\frac{5}{6}---\left(20\right)$

Wherein, R_safeiFor the threat radius of threatening area i, (Riski.x, Riski.y) is the two dimension that i-th threatens target Position, K_iFor the threat intensity of threatening area i, d_pI () is a p distance away from the threatening area center of circle, p represents in jth flight path section TheAlong ent.

2) shortest path target；

On the premise of ensureing flight safety, flight path is short means that the required time is short, oil consumption is few, so meeting about On the premise of the maximum flight path length of bundle, short path is better than long path, it is desirable to flight path length should be the least.Use true herein The ratio of flying distance and the possible length of minimum defines the shortest path length desired value calculating a paths, i.e. utilizes minimum to fly Row distance l_minActual path length is normalized as object function f_lenr, see formula (21):

$f_{lenr}=\frac{{\&Sigma;}_{i=1}^{N-1}\sqrt{{(x_{i+1}-x_i)}^2+{(y_{i+1}-y_i)}^2+{(z_{i+1}-z_i)}^2}}{l_{\min }}---\left(21\right)$

$l_{min}=\sqrt{{(x_{end}-x_{start})}^2+{(y_{end}-y_{start})}^2+{(z_{end}-z_{start})}^2}---\left(22\right)$

Wherein, l_minFor the air line distance of origin-to-destination, (x_start,y_start,z_start) and (x_end,y_end,z_end) be respectively Beginning and end coordinate；(x_i,y_i,z_i) and (x_i+1,y_i+1,z_i+1) it is respectively the coordinate of arbitrary neighborhood 2 between beginning and end.

For described collaborative collision avoidance constraint, when multimachine (M frame aircraft) performs aerial mission, each machine needs meeting Get around no-fly zone on the basis of its flying quality, avoid threatening area etc. and find a quick path of economy, arrive respective Destination implements task, the flight path of each frame aircraft will meet along flight course along the flight path planned for it not with The arbitrary aircraft flown along remaining aircraft's flight track bumps against.For the sake of security, once any two frame aircraft of synchronization Between distance less than a certain secure threshold, then it is assumed that the flight path that this two framves aircraft is corresponding can produce conflict, is unsafe.

The collaborative collision avoidance binding occurrence C of the flight path a of aircraft A_acoFor: the flight path a of aircraft A (sets with remaining all aircraft One has M frame aircraft, then with other M-1 frame aircraft) respective population represent the sum of conflicting of flight path, such as formula (23) institute Show:

$C_{aco}={\mathrm{\&Sigma;}}_{\underset{\&ForAll;A\&NotEqual;B}{B=1}}^M{\mathrm{\&Sigma;}}_{i=0}^{N_{tab}+1}{c}_i---\left(23\right)$

$c_i=\left\{\begin{array}{c}1-\frac{d_i^{ab}}{d_{\min }},if{d}_i^{ab}<d_{\min }\\ 0\end{array}\right.---\left(24\right)$

$d_i^{ab}=\sqrt{{(x_{ai}-x_{bi})}^2+{(y_{ai}-y_{bi})}^2+{(z_{ai}-z_{bi})}^2}---\left(25\right)$

Wherein, b is that the representative flight path of the population of aircraft B is individual, N_tabFor taking the total degree of time step,For aircraft A along each flight path a and aircraft B along flight path b at synchronization each position point { x_ai,y_ai,z_aiAnd { x_bi,y_bi, z_biSpace length.

Second step, based on described aircraft path planning model, uses adaptive differential target evolution algorithm to single rack Aircraft carries out path planning, in conjunction with flow chart 5, specifically comprises the following steps that

S1: determine planning space, starting point S and impact point E, utilizes heuristic initialization algorithm to initialize population, and complete Become the individual population pretreatment converted to absolute value coordinate individuality of contrast pole coordinate.

S2: according to headway v and sampling time step delta t, population at individual is completed track points discretization.

S3: according to the trajectory planning model set up in the first step, utilize the individual discretization track points of initial population Carry out flight path constraint checking and desired value calculates, and calculate constrained binding occurrence sum err_ind of each individuality.

S4: carry out parent population carrying out genetic cross and mutation operation according to adaptive differential operator, form sub-population.

S5: complete sub-population generation of neutrons individuality track points discretization, bring the track points of discretization into problem model, complete With individuality constraint, target function value and the calculating of constrained binding occurrence sum err_ind.

S6: parent is formed with filial generation and mixes group, and mixing group implements the quick non-dominated ranking improved, and determines that mixing group is each The individual Pareto layer being distributed.

S7: the crowding distance that improves of individuality and the Crowing Mechanism to each layer of group of mixing, according to dominance relation and individual Body crowding is chosen the individuality of Population Size number from mixing group and is formed new parent population.

S8: judge whether evolutionary process terminates.If meeting end condition (reaching the algebraically of greatest iteration), proceed to S9； Otherwise proceed to S4.

S9: select individuality in non-dominant layer (ground floor) from population, its flight path represented is required flight path.

Initialization of population described in step S1, specific as follows:

Carry out the premise of trajectory planning with multi-objective Evolutionary Algorithm to determine that and carry out evolutional operation chromosome structure, its structure Directly determine efficiency and the quality of path planning.The present invention use heuristic population algorithm carry out the side of initialization of population Method, uses two kinds of initial population individuality expression waies, i.e. contrast pole codes co-ordinates mode and absolute value codes co-ordinates mode, comprises Contrast pole Coordinate generation initial chromosome and initial population individuality Coordinate Conversion two step:

(S101) by formula (26) acquisition contrast pole coordinate, (r, θ, value z) represents the three-dimensional of its correspondence to flight path node initializing Flight path node, wherein rand () represents and is uniformly distributed value, the lower limit Minseg of r and the upper limitRepresent between adjacent track points The shortest permission distance, the bound of turn_ θ is relevant with aircraft mobility, the θ obtained_iPositive and negative represent next flight path section The yawing moment of more previous flight path section, is just counterclockwise, is negative clockwise.This strategy benefit is the most abundant in initialization procedure Make use of Given information, limit such as flight path section minimum length and maximum deflection angle limits, thus reduce search volume.

$\left\{\begin{array}{c}{r}_i=rand(Minseg,\frac{MaxL}{N_w})\\ {\mathrm{\&theta;}}_i=rand(-turn\_\mathrm{\&theta;},turn\_\mathrm{\&theta;})\\ z_i=rand\left(\min \right(start.z,end.z),\max (start.z,end.z))\end{array}\right.---\left(26\right)$

Wherein, (start.z, end.z) represents the height coordinate value of starting point and impact point.

(S102) initial population individuality coordinate system conversion:

After utilizing contrast pole coordinate method to generate initial population, utilize contrast pole coordinate and the corresponding relation of absolute value coordinate, By formula (27) and formula (28), the individual two-dimentional flight path of the contrast pole coordinate representation to generating carries out coordinate system conversion, will generate Contrast pole coordinate initial population change into the initial population of absolute value coordinate representation.

Fig. 6 illustrates the same path using (a) contrast pole coordinate and (b) cartesian coordinate respectively.First flight path node W₁Absolute value coordinate (x₁,y₁,z₁) can be by its contrast pole coordinate (r₁,θ₁,z₁) obtained, wherein by formula (27) and formula (28) (start.x, start.y) and the two-dimensional coordinate that (end.x, end.y) is starting point and impact point.

$\left\{\begin{array}{c}\mathrm{\&alpha;}=arctan\left(\frac{end.y-start.y}{end.x-start.x}\right)\\ x_1=start.x+r_1*cos(\mathrm{\&alpha;}+{\mathrm{\&theta;}}_1)\\ y_1=start.y+r_1*sin(\mathrm{\&alpha;}+{\mathrm{\&theta;}}_1)\\ z_1=z_1\end{array}\right.---\left(27\right)$

Based on the W tried to achieve₁(x₁,y₁,z₁), follow-up track points W can be obtained_i(i=2 ..., N-1) absolute value coordinate (x_i, y_i,z_i):

$\left\{\begin{array}{c}sum\mathrm{\&theta;}={\&Sigma;}_{k=1}^{i+1}{\mathrm{\&theta;}}_k\\ x_i=x_{i-1}+r_{i+1}*\cos (\mathrm{\&alpha;}+sum\mathrm{\&theta;})\\ y_i=y_{i-1}+r_{i+1}*\sin (\mathrm{\&alpha;}+sum\mathrm{\&theta;})\\ z_i=z_i\end{array}\right.,i=2toN-1---\left(28\right)$

Population at individual flight path node discretization described in step S2, specifically refers to:

Two the most adjacent flight path nodes form a flight path section, and aircraft flies along flight path section successively, so flight path Point in section is also to need to carry out constraint checking.The present invention, will boat on the premise of supposing that aircraft flies at a constant speed with speed v Mark node discretization, obtains carrying out the discretization point sequence dis{S, P of constraint checking₁,P₂,...,P_N-2,E}.With by two continuous print Track points P_i(x_i,y_i,z_i) and P_i+1(x_i+1,y_i+1,z_i+1) flight path section l that constitutes_i(i=0 ..., N-1) as a example by, as it is shown in fig. 7, Flight path section l_iDiscrete point sequence disp_l_i{P_i1(x_i1,y_i1,z_i1),P_i2(x_i2,y_i2,z_i2),...,P_im(x_im,y_im,z_im) Coordinate figure can be obtained by following formula:

${\mathrm{\&theta;}}_i=arctan\left(\frac{y_{i+1}-y_i}{x_{i+1}-x_i}\right)---\left(31\right)$

$m=floor\left(\frac{\sqrt{{(x_{i+1}-x_i)}^2+{(y_{i+1}-y_i)}^2+{(z_{i+1}-z_i)}^2}}{\mathrm{\&Delta;}t*v}\right)---\left(32\right)$

WhereinFor flight path section l_iThe deviation angle of in the vertical direction, θ_iFor flight path section l_iRelative on two dimensional surface In the deviation angle of x-axis forward, Δ t is step sample time, and floor () is downward bracket function.

Population adaptive differential genetic cross and mutation operation in described step S4, form sub-population, particularly as follows:

S401: be evaluated initial population, i.e. calculates the target function value of each individuality in initial population.

S402: carry out mutation operation, to each target individualThe variation being produced its correspondence by formula (33) is individual

$v_i^{G+1}=x_{r1}^G+F*(x_{r2}^G-x_{r3}^G)---\left(33\right)$

Wherein:WithSubscript r₁, r₂And r₃It is the different number randomly choosed, and and object vector's Subscript i is the most different；F is zoom factor, for the magnification level of control deviation variable.

S403: by formula (34), by target individualWith its variation individualityCarry out intersecting and operate, then it is right to obtain The intersection answered is individualAfter the intersection of all target individual has operated, its intersection individuality forms interim population.

$u_{ij}^{G+1}=\left\{\begin{array}{cc}{v}_{ij}^{G+1},& if(rand\&le;CR)orj=q\\ x_{ij}^G,& if(rand>CR)andj\&NotEqual;q\end{array}\right.---\left(34\right)$

Wherein: CR for intersect the factor, be individuality carry out intersect operation control probability, q be one between [1, D] with The integer that machine produces, to guarantee to intersect individualityAt least can be from variation individualityObtain one-dimensional parameter.Letter D refers to one Constant.

The present invention take self adaptation dynamically adjust zoom factor F and the size strategy of intersection factor CR, in population Per generation individuality all uses different zoom factor F and factor CR of intersecting to carry out differential evolution, F and CR is during evolution along with kind Shown in the rule such as formula (35) and (36) that group evolves and changes:

$F=F_{min}+(F_{max}-F_{min})*\exp \left(\frac{-G}{Maxgen-G+0.1}\right)---\left(35\right)$

$CR=\frac{G}{Maxgen}*{\mathrm{CR}}_{max}---\left(36\right)$

Wherein, Maxgen represents maximum evolutionary generation, and G represents current algebraically.

This step obtains after terminating by initial population carries out the progeny population (sub-population) that difference operation obtains.

The quick non-dominated ranking of the described improvement in step S6, specific as follows:

Owing under the complicated low latitude in the present invention, aircraft's flight track planning is that the multiple target under meeting multi-constraint condition is excellent Change problem, in the case of the multiple-objection optimization under this Prescribed Properties, if a Yu b is two individualities in population, the present invention uses The quicksort operator improved under constraints completes the comparison of the dominance relation of individual a Yu b, and flow process is as follows:

S601: obtain individual a and individual b each constraint with self related constraint binding occurrence sum a.err_ind with b.err_ind。

S602: judge whether individual a with b breaks a contract, the most respective err_ind value and the magnitude relationship of 0.If a.err_ind Be both greater than 0 with b.err_ind, i.e. a and b violates constraints, is infeasible solution, then go to S603；If only one of which The err_ind of body is more than 0, then go to S604；If two individual err_ind values are equal to 0, then go to S605.

S603: if the err_ind of two individualities is different, go to S603.1；If two individual err_ind are identical, go to S603.2.

S603.1: by the individual principle that the individual domination binding occurrence sum that binding occurrence sum is big is little, if a.err_ind > b.err_ind, the then individual individual b of a domination, otherwise, then the individual individual a of b domination.

S603.2: individual a with b dominance relation not may compare, and does not arranges.

S604: violate the individual principle of constraints by the individual domination not violating constraints, if b.err_ind > 0 And a.ind_err=0, the then individual individual b of a domination；If a.err_ind > 0 and b.ind_err=0, the then individual individual a of b domination.

S605: the quick non-dominated ranking improved.

Wherein non-dominated ranking method operator step is as follows: (in population, each individual p has two parameters np and Sp, its Middle np is the individual amount of the individual p of domination in population, and Sp is the group of individuals arranged by p in population)

SS1: find the individuality of np=0 in population, and be saved in current collection F1.

SS2: interim set H composes sky.For each individual i in current collection F1, investigate the individual collections that it is arranged Si, subtracts 1 by the nj of individual j each in Si, if (i.e. individual j is only arranged nj-1=0 by individual i, illustrates that individual j is under individual i One layer), then individuality j is put into set H.Now gather the individual collections of next layer that H is current layer individual collections.

SS3:F1 is the 1st domination layer individual collections, individual identical non-dominated ranking irank in giving this layer.

SS4: with H as current collection, returns step SS2, until the sequence of whole population is complete.

The present invention designs a constant parameter for comparing precision δ for improving quick non-dominated ranking operator: if jth (one has a n optimization aim, j=1 ..., n) on object function, have | f_j(a)-f_j(b) |≤δ, then it is assumed that individual a and b is at f_j On value equal, if individual a and b has such relation on all of target, then individual a with b does not arranges.Its In, f_j(a) and f_j(b) represent respectively individual a and individual b in jth object function f_jFunctional value.So algorithm can be separately The individuality being in a disadvantageous position on one function by another individual domination, thus will gradually remove those at certain in search procedure Or the individuality the least but poor in other targets with remaining individual difference in some target, make final result occurs unreasonable Individual probability reduces.

Crowding distance for the improvement described in step S7 calculates, particularly as follows:

As a example by crowding distance individual on population a certain Pareto layer F below calculates, provide the meter improving crowding distance Calculation method:

S701: calculate individual number NF of population a certain Pareto layer F.

S702: the crowding distance of each individuality in F is composed null value, it may be assumed that F [i] .dis=0, i=1 ..., popsize. Popsize is number individual in Population Size, i.e. population.

S703: to each object function f_j(j=1 ..n), obtains functional value F [i] .f of each individuality_j, by F [i] .f_j Size individuality in F is ranked up.

S704: according to the ordering scenario that F layer is individual, by infinitely great for crowding distance assignment individual for boundary point, by remaining The crowding distance of body is pressed(i=2 ..., NF-1) assignment, i.e. with body one by one The diagonal sum of besieged rectangle represents the crowding distance of this individuality；Wherein, f_kmaxAnd f_kminRepresent kth function respectively Maximum and minima.

Based on population sequence and the Crowing Mechanism of each layer of individual crowding in step S7, by quick non-dominated ranking with And after crowding calculates, each individual i has two attributes: non-dominant sequence i_rankWith crowding distance I [i] .distance.Utilize The two attribute, can choose Population Size number purpose individuality from mixing group (2 times of former Population Sizes) and form a new generation and enter The parent population of travelingization, selection standard is: even two individual non-dominated ranking are different, then take the individuality that sequence number is less；As Really two individualities are in same one-level, then take the individuality of crowding distance big (crowding is little).

Based on described single rack aircraft paths planning method, the invention allows for based on the many mesh of cooperative self-adapted difference The multi-aircraft path planning method of mark evolution algorithm, flow chart as shown in Figure 8, specifically comprises the following steps that

Step A: each frame aircraft is carried out heuristic initialization population:

A1: determine planning space, starting point S and impact point E, utilizes heuristic initialization algorithm to initialize population, and complete Become the individual population pretreatment converted to absolute value coordinate individuality of contrast pole coordinate.

A2: after all aircraft populations all complete heuristic initialization, enters step B.

Step B: to each frame aircraft population:

B1: according to headway v and sampling time step delta t, population at individual is completed track points discretization.

B2: describe the trajectory planning model set up in a joint according to unit trajectory planning problem, utilize initial population Individual discretization track points carries out flight path constraint checking and desired value calculates, and calculates institute's Constrained sum of each individuality err_ind。

B3: carry out parent population carrying out genetic manipulation according to adaptive differential operator, form sub-population.

B4: complete offspring individual track points discretization, bring the track points of discretization into problem model, complete with individuality about Bundle, target function value and the calculating of institute's Constrained sum err_ind.

B5: parent is formed with filial generation and mixes group, mixing group implements the quick non-dominated ranking improved, determines mixing group The Pareto layer of individual distribution.

B6: represent individual selection mechanism according to population, the representative selecting this aircraft population from mixing group is individual.

B7: after all aircraft populations all complete the representative individual selection of mixing group, enters step C.

Step C: to each frame aircraft population:

C1: receive the representative individuality that remaining each aircraft mixing group selects, carries out this aircraft mixing all individualities of group Collaborative collision avoidance binding occurrence calculates, and obtains mixing institute's Constrained sum err_all of each individuality in group.

C2: the crowding distance improving the individuality mixing each layer of group calculates and Crowing Mechanism, according to dominance relation And individuality crowding is from mixing the parent population that the individuality composition choosing Population Size number group is new.

C3: judge whether evolutionary process terminates.If meeting end condition, proceed to C4；Otherwise proceed to step B.

C4: select its flight path represented individual from population in non-dominant layer (ground floor) and be required flight path.

The present invention uses multi-species cooperative, will the kinestate of all unmanned planes, Performance Constraints, threat feelings Condition, target etc. put together and carry out large-scale optimizatoin, obtain meet constraint and the unmanned aerial vehicle flight path of spatio-temporal synergy will become higher-dimension, Close coupling, optimizes the problem that difficulty is big, and amount of calculation is huge, and complexity is high.By means of the thought of coevolution, centralized optimization is asked Topic resolves into low-dimensional problem, and each independent population represents different aircraft, parallel evolution respectively, receives during evolution The representative individuality that remaining population selects carries out collaborative constraint and evaluates its population solution.

In above-mentioned steps C1, all individualities work in coordination with the calculating of collision avoidance binding occurrence, calculate according to collision detection mechanism.Count That calculates a flight path exempts from crash restraint, just need to set up the collision detection mechanism between flight path.The present invention proposes more real-time Collision detection mechanism between " time point is to time point " flight path, tells about below by the flight path crash restraint analysis of certain two frame aircraft The collision detection mechanism of the present invention.

As it is shown in figure 9, aircraft A and aircraft B is respectively from respective starting point S_AAnd S_BFly along its flight path a and flight path b To terminal E_AAnd E_B, W_ak,W_bk(k=1 ... .N_w；N_wFor needing the flight path nodes of planning) it is respectively flight path a and the boat of flight path b Mark node, p_ai,p_biIt is respectively certain synchronization aircraft A and aircraft B location point, v on flight path a and flight path b_aWith v_bBeing respectively the flight speed of aircraft, Δ t is a little time step, then aircraft A plans flight path a's and aircraft B The flow process of the collision detection that planning flight path b is carried out is:

S1: according to position calculation operator (Position-Obtain Strategy, SOS), obtains from taking off Moment rises, and every time step Δ t, two airplanes are respective by location point sequence on flight path a and flight path b respectivelyWithP_aAnd P_bFirst prime number N of point sequence_tDuring by flight Between flight time of short that flight path determine.

S2: from the beginning of taking off, every Δ t, calculates at synchronization, along the aircraft A of flight path a flight with along path Distance between the aircraft B of b flight, i.e. starts to i=N from i=0_t, calculate p_aiAnd p_biDistance | p_aip_bi|, it is designated as

S3: judge the most whether two machines along respective track flight produce conflict: ifLess than the safety preset Distance, i.e.The path b then claiming the path a and aircraft B of aircraft A occurs once to conflict, i.e. the boat of aircraft A The flight path b of mark a and aircraft B exempts from crash restraint promise breaking.

Needing in collision detection mechanism between aircraft's flight track to use position calculation operator SOS, its flow process is as follows:

The respective distance of flight path b of the flight path a and aircraft B of SS1: calculating aircraft A, is designated as L respectively_aAnd L_b, and ByObtain each machine flight time T respectively on this two flight path_aAnd T_b。

SS2: calculating aircraft A at flight path a and aircraft B in upper all flight path sections of flight path b Flight time sequence, the most all is arrived the sequence that the flight time of next flight path node constitutes by a flight path node:WithCalculate again and respectively fly Row device arrives the moment sequence of respective track pointsWithWith flight path As a example by a, wherein t_aCan be calculated by formula (37) and try to achieve:

$t_{ai}={\mathrm{\&Sigma;}}_{k=1}^i{T}_{segak},i=1,2,...,N_w---\left(37\right)$

$T_{segai}=\left\{\begin{array}{cc}\frac{\left|W_{ai}{S}_A\right|}{v},& i=1\\ \frac{\left|W_{ai}{W}_{ai-1}\right|}{v},& i=2:N_w\\ \frac{\left|E_A{W}_{ai}\right|}{v},& i=N_w+1\end{array}\right.---\left(38\right)$

SS3: calculate the testing time N of collision detection_t, whereinCelling () is upwards Bracket function.

SS4: determine that flight path a and flight path b carries out the sampling location point sequence of collision detection With

In position calculation operator, need to solve p_ai{x_ai,y_ai,z_ai, its algorithm flow is:

SSS1:begin_t record is from the time of the starting point of each flight path section starting sample point to this flight path section, begin_ It is 0 that t composes initial value；K is flight path section subscript, k=1:N-1.

SSS2: as k≤Nw+1:

1) calculate accumulative sampled point time det, the det=i* Δ t+being_t in kth flight path section, calculate kth flight path The deviation angle of section in the vertical directionWith deviation angle θ relative to x-axis forward on two dimensional surface_k。

2) when det is less than the flight time T of this flight path section (kth section flight path section)_segakTime, i.e. as det≤t_akTime, by formula (39) p is calculated_ai{x_ai,y_ai,z_ai, and i=i+1:

Otherwise (i.e. as det ＞ T_segakTime):

3) judge whether to arrive last sampled point, i.e. work as i=N_tWhen+1, calculate p by formula 40_ai{x_ai,y_ai,z_ai}:

$\left\{\begin{array}{c}{x}_{ai}=W_{ak}.x\\ y_{ai}=W_{ak}.y\\ z_{ai}=W_{ak}.z\end{array}\right.---\left(40\right)$

Otherwise k=k+1, and being_t=det-t_ak。

On the basis of above-mentioned collision detection mechanism, it is possible to the collaborative collision avoidance binding occurrence (formula in computational problem model 23)。

Propose population in the present invention representing selection mechanism is cooperation individuality selection mechanism.

Below as a example by certain aircraft population, introduce population and represent the flow process of individual selection:

S1: select the individuality of all non-dominant layers in population, form set to be selected.

S2: in set to be selected, finds safe flight target f_riskThe individuality of=0, if f_riskThe individuality of=0 has multiple, Then go to S3；Without f_riskThe individuality of=0, the most all individual f_risk＞ 0, then enter S4；If f_riskThe individuality of=0 is only Have one, then this individuality is as representing individuality, and terminates to represent individual selection operation；

S3: to all f_riskThe individuality of=0 carries out flight path length f_lenrComparison, find, set to be selected will have minimum Flight path length min (f_lenr) individuality as representing individuality, and terminate represent individual selection operation.

S4: all f_riskThe individuality of ＞ 0 carries out f_riskSize compare, find the minimum threat index min in set to be selected (f_risk)。

S5: and with min (f_risk) to all f_riskThe individual f of ＞ 0_riskIt is normalized, such as individual i correspondence phase To threat index it is

S6: discard Ref_riskiThe solution of >=V, wherein V is the comparatively safe threshold value set.

S7: carry out flight path length f in the individuality retained in set to be selected_lenrComparison, will set to be selected have Little flight path length min (f_lenr) individuality as representing individuality, and terminate to represent individual selection operation, and terminate to represent individual choosing Select operation.

So, the most just have selected that potential threat index is relatively low and the solution of fair-sized path.Additionally, plant Group represents individual selection operator in addition to evaluating for the collaborative constraint of aircraft's flight track.Non-has been obtained after running through maximum algebraically Join Pareto disaggregation, it is possible to use this operator solves concentration from non-dominant Pareto and picks out last solution, as this aircraft Finally select flight path.

Claims (5)

1. one kind is applicable to the paths planning method of aircraft under complicated low latitude, it is characterised in that: assume (1) all aircraft all Fly with friction speed at sustained height layer, and in each comfortable flight course, flight speed does not changes；(2) every frame aircraft All determine respective beginning and end；(3) low altitude airspace flight environment of vehicle it is known that include prohibited flight area position and scope, Threaten some position and radiation radius；(4) general purpose vehicle physical property limits and includes maximum deflection angle and maximum flying distance；Institute State paths planning method and specifically include following steps:

The first step, sets up aircraft path planning model；

Second step, uses adaptive differential target evolution algorithm that single rack aircraft is carried out path planning, specifically comprises the following steps that

S1: determine planning space, starting point S and impact point E, utilizes heuristic initialization algorithm to initialize population, and completes phase The population pretreatment that to absolute value coordinate individuality convert individual to polar coordinate；

S2: according to headway v and sampling time step delta t, population at individual is completed track points discretization；

S3: according to the aircraft path planning model set up in the first step, utilize the individual discretization flight path of initial population Point carries out flight path constraint checking and desired value calculates, and calculates constrained binding occurrence sum err_ind of each individuality；

S4: carry out parent population carrying out genetic cross and mutation operation according to adaptive differential operator, form sub-population；

S5: complete sub-population generation of neutrons individuality track points discretization, bring the track points of discretization into problem model, completes with individual Body constraint, target function value and the calculating of constrained binding occurrence sum err_ind；

S6: parent is formed with filial generation and mixes group, and mixing group implements the quick non-dominated ranking improved, and determines the mixing each individuality of group The Pareto layer being distributed；

S7: the crowding distance improving the individuality mixing each layer of group and Crowing Mechanism, gather around according to dominance relation and individuality Squeeze degree from mixing the parent population that the individuality composition choosing Population Size number group is new；

S8: judge whether evolutionary process terminates, if reaching the algebraically of greatest iteration, proceeds to S9；Otherwise proceed to S4；

S9: select individuality in non-dominant layer from population, its flight path represented is required flight path；

Initialization population described in step S1, uses two kinds of initial population individuality expression waies, i.e. contrast pole codes co-ordinates sides Formula and absolute value codes co-ordinates mode, comprise contrast pole Coordinate generation initial chromosome and initial population individuality Coordinate Conversion two Step:

(S101) by formula (26) acquisition contrast pole coordinate, (r, θ, value z) represents the Three-dimensional Track of its correspondence to flight path node initializing Node, wherein rand () represents and is uniformly distributed value, the lower limit Minseg of r and the upper limitRepresent between adjacent track points is the shortest Allowing distance, the bound of turn_ θ is relevant with aircraft mobility, the θ obtained_iPositive and negative represent next flight path section earlier above The yawing moment of one flight path section, is just counterclockwise, is negative clockwise；

$\left\{\begin{array}{c}{r}_i=rand(Minseg,\frac{MaxL}{N_w})\\ {\mathrm{\&theta;}}_i=rand(-turn\_\mathrm{\&theta;},turn\_\mathrm{\&theta;})\\ z_i=rand\left(\min \right(start.z,end.z),\max (start.z,end.z\left)\right)\end{array}\right.---\left(26\right)$

Wherein, (start.z, end.z) represents the height coordinate value of starting point and impact point；

(S102) initial population individuality coordinate system conversion:

After utilizing contrast pole coordinate method to generate initial population, utilize contrast pole coordinate and the corresponding relation of absolute value coordinate, pass through The individual two-dimentional flight path of formula (27) and the formula (28) the contrast pole coordinate representation to generating carries out coordinate system conversion, the phase that will generate Polar coordinate initial population is changed into the initial population of absolute value coordinate representation；First flight path node W₁Absolute value coordinate (x₁,y₁,z₁) by its contrast pole coordinate (r₁,θ₁,z₁) obtained by formula (27) and formula (28), wherein (start.x, start.y) and (end.x, end.y) is the two-dimensional coordinate of starting point and impact point；

$\left\{\begin{array}{c}\mathrm{\&alpha;}=\arctan \left(\frac{end.y-start.y}{end.x-start.x}\right)\\ x_1=start.x+r_1*cos(\mathrm{\&alpha;}+{\mathrm{\&theta;}}_1)\\ y_1=start.y+r_1*sin(\mathrm{\&alpha;}+{\mathrm{\&theta;}}_1)\\ z_1=z_1\end{array}\right.---\left(27\right)$

Based on the W tried to achieve₁(x₁,y₁,z₁), obtain follow-up track points W_iAbsolute value coordinate (x_i,y_i,z_i):

$\left\{\begin{array}{c}sum\mathrm{\&theta;}={\&Sigma;}_{k=1}^{i+1}{\mathrm{\&theta;}}_k\\ x_i=x_{i-1}+r_{i+1}*\cos (\mathrm{\&alpha;}+sum\mathrm{\&theta;})\\ y_i=y_{i-1}+r_{i+1}*\sin (\mathrm{\&alpha;}+sum\mathrm{\&theta;})\\ z_i=z_i\end{array}\right.---\left(28\right)$

Wherein, i=2 ..., N-1；

Population at individual flight path node discretization described in step S2, specifically refers to:

On the premise of supposing that aircraft flies at a constant speed with speed v, by flight path node discretization, obtain carrying out constraint checking from Dispersion point sequence dis{S, P₁,P₂,…,P_N-2,E}；With by two continuous print track points P_i(x_i,y_i,z_i) and P_i+1(x_i+1,y_i+1, z_i+1) flight path section l that constitutes_iAs a example by, i=0 ..., N-1, flight path section l_iDiscrete point sequence disp_l_i{P_i1(x_i1,y_i1,z_i1), P_i2(x_i2,y_i2,z_i2),…,P_im(x_im,y_im,z_im) coordinate figure obtained by following formula:

${\mathrm{\&theta;}}_i=\arctan \left(\frac{y_{i+1}-y_i}{x_{i+1}-x_i}\right)---\left(31\right)$

$m=floor\left(\frac{\sqrt{{(x_{i+1}-x_i)}^2+{(y_{i+1}-y_i)}^2+{(z_{i+1}-z_i)}^2}}{\mathrm{\&Delta;}t*v}\right)---\left(32\right)$

WhereinFor flight path section l_iThe deviation angle of in the vertical direction, θ_iFor flight path section l_iOn two dimensional surface relative to x-axis The deviation angle of forward, Δ t is step sample time, and floor () is downward bracket function；

Population adaptive differential genetic cross and mutation operation in described step S4, form sub-population, particularly as follows:

S401: be evaluated initial population, i.e. calculates the target function value of each individuality in initial population；

S402: carry out mutation operation, to each target individualThe variation being produced its correspondence by formula (33) is individual

$v_i^{G+1}=x_{r1}^G+F*(x_{r2}^G-x_{r3}^G)---\left(33\right)$

Wherein:WithSubscript r₁, r₂And r₃It is the different number randomly choosed, and and object vectorSubscript i Also different；F is zoom factor, for the magnification level of control deviation variable；

S403: by formula (34), by target individualWith its variation individualityCarry out intersecting and operate, then obtain the intersection of correspondence IndividualAfter the intersection of all target individual has operated, its intersection individuality forms interim population；

$u_{ij}^{G+1}=\left\{\begin{array}{ccc}{v}_{ij}^{G+1},& if& (rand\&le;CR)orj=q\\ x_{ij}^G,& if& (rand>CR)andj\&NotEqual;q\end{array}\right.---\left(34\right)$

Wherein: CR, for intersecting the factor, is that individuality carries out intersecting the control probability of operation, and q is a random product between [1, D] Raw integer, to guarantee to intersect individualityAt least can be from variation individualityObtain one-dimensional parameter；Letter D often refers to one Number；

The quick non-dominated ranking of the described improvement in step S6, specific as follows:

If a Yu b is two individualities in population, the quicksort operator improved under constraints is used to complete propping up of individual a Yu b Joining the comparison of relation, flow process is as follows:

S601: obtain individual a and individual b each constraint with self related constraint binding occurrence sum a.err_ind with b.err_ind；

S602: judge whether individual a with b breaks a contract, the most respective err_ind value and the magnitude relationship of 0, if a.err_ind and B.err_ind is both greater than 0, i.e. a and b violates constraints, is infeasible solution, then go to S603；If only one of which is individual Err_ind more than 0, then go to S604；If two individual err_ind values are equal to 0, then go to S605；

S603: if the err_ind of two individualities is different, go to S603.1；If two individual err_ind are identical, go to S603.2；

S603.1: by the individual principle that the individual domination binding occurrence sum that binding occurrence sum is big is little, if a.err_ind > B.err_ind, the then individual individual b of a domination, otherwise, then the individual individual a of b domination；

S603.2: individual a with b dominance relation not may compare, and does not arranges；

S604: violate the individual principle of constraints by the individual domination not violating constraints, if b.err_ind > 0 and A.ind_err=0, the then individual individual b of a domination；If a.err_ind > 0 and b.ind_err=0, the then individual individual a of b domination；

S605: the quick non-dominated ranking improved；

The quick non-dominated ranking method operator step improved is as follows: in population, each individual p has two parameters np and Sp, its Middle np is the individual amount of the individual p of domination in population, and Sp is the group of individuals arranged by p in population,

SS1: find the individuality of np=0 in population, and be saved in current collection F1；

SS2: interim set H composes sky；For each individual i in current collection F1, investigate individual collections Si that it is arranged, will In Si, the nj of each individual j subtracts 1, if nj-1=0, then individuality j puts into set H；Now set H is current layer individual collections The individual collections of next layer；

SS3:F1 is the 1st domination layer individual collections, individual identical non-dominated ranking irank in giving this layer；

SS4: with H as current collection, returns step SS2, until the sequence of whole population is complete；

Designing a constant parameter is to compare precision δ for improving quick non-dominated ranking operator: if the mesh of jth optimization aim On scalar functions, have | f_j(a)-f_j(b) |≤δ, then it is assumed that individual a and b is at f_jOn value equal, if individual a and b is all of Have such relation on target, then individual a with b does not arranges；Wherein, f_j(a) and f_jB () represents individual a and individuality respectively B in jth object function f_jFunctional value；

Crowding distance for the improvement described in step S7 calculates, with crowding distance individual on population a certain Pareto layer F As a example by calculating, be given improve crowding distance computational methods:

S701: calculate individual number NF of population a certain Pareto layer F；

S702: the crowding distance of each individuality in F is composed null value, it may be assumed that F [i] .dis=0, i=1 ..., popsize； Popsize is number individual in Population Size, i.e. population；

S703: to each object function f_j, j=1 ..n, obtain functional value F [i] .f of each individuality_j, by F [i] .f_jSize Individuality in F is ranked up；

S704: according to the ordering scenario that F layer is individual, by infinitely great, by remaining individuality for crowding distance assignment individual for boundary point Crowding distance is pressedAssignment, i=2 ..., NF-1, i.e. besieged with body one by one The diagonal sum of rectangle represent the crowding distance of this individuality；Wherein, f_kmaxAnd f_kminRepresent the maximum of kth function respectively Value and minima.

The most according to claim 1 a kind of it is applicable to the paths planning method of aircraft under complicated low latitude, it is characterised in that: Based on population sequence and the Crowing Mechanism of each layer of individual crowding in step S7, by quick non-dominated ranking and crowding After calculating, each individual i has two attributes: non-dominant sequence i_rankWith crowding distance I [i] .distance；The two is utilized to belong to Property, from mixing group, to choose Population Size number purpose individuality form the parent population that a new generation evolves, selection standard is: If two individual non-dominated ranking are different, then take the individuality that sequence number is less；If two individualities are in same one-level, then take crowding distance Big individuality.

3. a multi-aircraft path planning method based on cooperative self-adapted difference multi-objective Evolutionary Algorithm, it is characterised in that bag Include following steps:

Step A: each frame aircraft is carried out heuristic initialization population:

A1: determine planning space, starting point S and impact point E, utilizes heuristic initialization algorithm to initialize population, and completes phase The population pretreatment that to absolute value coordinate individuality convert individual to polar coordinate；

A2: after all aircraft populations all complete heuristic initialization, enters step B；

Step B: to each frame aircraft population:

B1: according to headway v and sampling time step delta t, population at individual is completed track points discretization；

B2: according to aircraft path planning model, utilizes the individual discretization track points of initial population to carry out flight path constraint and tests Card and desired value calculate, and calculate institute's Constrained sum err_ind of each individuality；

B3: carry out parent population carrying out genetic manipulation according to adaptive differential operator, form sub-population；

B4: complete offspring individual track points discretization, bring the track points of discretization into problem model, complete with individual constraint, Target function value and the calculating of institute's Constrained sum err_ind；

B5: parent is formed with filial generation and mixes group, and mixing group implements the quick non-dominated ranking improved, and determines that mixing group is individual The Pareto layer of distribution；

B6: represent individual selection mechanism according to population, the representative selecting this aircraft population from mixing group is individual；

B7: after all aircraft populations all complete the representative individual selection of mixing group, enters step C；

Step C: to each frame aircraft population:

C1: receive the representative individuality that remaining each aircraft mixing group selects, carries out this aircraft mixing all individualities of group and works in coordination with Collision avoidance binding occurrence calculates, and obtains mixing institute's Constrained sum err_all of each individuality in group；

C2: calculate the crowding distance that improves of individuality of each layer of group of mixing and Crowing Mechanism, according to dominance relation and individual Body crowding is chosen the individuality of Population Size number from mixing group and is formed new parent population；

C3: judge whether evolutionary process terminates, if meeting end condition, proceeds to C4；Otherwise proceed to step B；

C4: select its flight path represented individual from population in non-dominant layer and be required flight path.

A kind of multi-aircraft flight path based on cooperative self-adapted difference multi-objective Evolutionary Algorithm the most according to claim 3 is advised The method of drawing, it is characterised in that: in step C1, all individualities work in coordination with the calculating of collision avoidance binding occurrence, count according to collision detection mechanism Calculate；Described collision detection mechanism, it is assumed that aircraft A and aircraft B is respectively from respective starting point S_AAnd S_BAlong its flight path a and Flight path b flies to terminal E_AAnd E_B, W_ak,W_bkThe respectively flight path node of flight path a and flight path b, k=1 ... .N_w；N_wAdvise for needs The flight path nodes drawn, p_ai,p_biIt is respectively certain synchronization aircraft A and aircraft B position residing on flight path a and flight path b Put a little, v_aAnd v_bBeing respectively the flight speed of aircraft, Δ t is a little time step, then to the planning flight path a of aircraft A with fly The flow process of the collision detection that the planning flight path b of row device B is carried out is:

S1: according to position calculation operator, obtains from departure time, every time step Δ t, two airplanes respectively at flight path a and Each by location point sequence on flight path bWith P_aAnd P_bFirst prime number N of point sequence_tDetermined by the flight time of flight time short that flight path；

S2: from the beginning of taking off, every Δ t, calculates at synchronization, along the aircraft A of flight path a flight and flies along path b Distance between the aircraft B of row, i.e. starts to i=N from i=0_t, calculate p_aiAnd p_biDistance | p_aip_bi|, it is designated as

S3: judge the most whether two machines along respective track flight produce conflict: ifLess than preset safety away from From, i.e.The path b then claiming the path a and aircraft B of aircraft A occurs once to conflict, i.e. the flight path a of aircraft A Break a contract with the crash restraint of exempting from of the flight path b of aircraft B.

A kind of multi-aircraft flight path based on cooperative self-adapted difference multi-objective Evolutionary Algorithm the most according to claim 4 is advised The method of drawing, it is characterised in that: the flow process of described position calculation operator is as follows:

The respective distance of flight path b of the flight path a and aircraft B of SS1: calculating aircraft A, is designated as L respectively_aAnd L_b, and byObtain each machine flight time T respectively on this two flight path_aAnd T_b；

SS2: calculating aircraft A at flight path a and aircraft B in upper all flight path sections of flight path b Flight time sequence, the most all is arrived the sequence that the flight time of next flight path node constitutes by a flight path node:WithCalculate again and respectively fly Row device arrives the moment sequence of respective track pointsWithWith flight path As a example by a, wherein t_aCalculated by formula (37) and try to achieve:

$t_{ai}={\mathrm{\&Sigma;}}_{k=1}^i{T}_{segak},i=1,2,...,N_w---\left(37\right)$

$T_{segai}=\left\{\begin{array}{c}\frac{\left|W_{ai}{S}_A\right|}{v},i=1\\ \frac{\left|W_{ai}{W}_{ai-1}\right|}{v},i=2:N_w\\ \frac{\left|E_A{W}_{ai}\right|}{v},i=N_w+1\end{array}\right.---\left(38\right)$

SS3: calculate the testing time N of collision detection_t, whereinCelling () is for rounding up Function；

SS4: determine that flight path a and flight path b carries out the sampling location point sequence of collision detectionWith

Wherein solve p_ai{x_ai,y_ai,z_aiAlgorithm flow be, i=1 ..., N_t:

SSS1:begin_t record is composed from the time of the starting point of each flight path section starting sample point to this flight path section, begin_t Initial value is 0；K is flight path section subscript, k=1:N-1；

SSS2: as k≤Nw+1:

1) calculate accumulative sampled point time det, the det=i* Δ t+being_t in kth flight path section, calculate kth flight path section and exist Deviation angle on vertical directionWith deviation angle θ relative to x-axis forward on two dimensional surface_k；

2) when det is less than the flight time T of this flight path section_segakTime, i.e. as det≤t_akTime, calculate p by formula (39)_ai{x_ai,y_ai, z_ai, and i=i+1:

Otherwise:

3) judge whether to arrive last sampled point, i.e. work as i=N_tWhen+1, calculate p by formula (40)_ai{x_ai,y_ai,z_ai}:

$\left\{\begin{array}{c}{x}_{ai}=W_{ak}.x\\ y_{ai}=W_{ak}.y\\ z_{ai}=W_{ak}.z\end{array}\right.---\left(40\right)$

Otherwise k=k+1, and being_t=det-t_ak。