CN109508030A - A kind of collaboration parsing reentry guidance method considering more no-fly zone constraints - Google Patents

Abstract

A kind of collaboration parsing reentry guidance method for considering more no-fly zone constraints of the present invention, comprising steps of 1, reentry guidance problem describes；2, the flight time analytic solutions based on rotation earth model；3, consider the parsing reentry guidance method of more no-fly zones and arrival time constraint；4, multi-aircraft arrival time cooperative approach.The invention has the advantages that: (1) the case where proposing the Exact of the hypersonic glide trajectories flight time based on rotation earth model, predicting that error is maintained within 3%, be suitable for non-constant value longitudinal direction lift resistance ratio section.(2) it can solve the problems, such as to guide multiple hypersonic glide vehicles in more no-fly zone environment while reach the collaboration flight of target based on three-dimensional resolution Value of Reentry Vehicle and flight time analytic solutions.(3) the multiple target iterative numerical programme based on online Ballistic Simulation of Underwater for considering terminal time, speed and requirement for height is devised.Utilization orientation derivative improves quasi-Newton method, reduces Ballistic Simulation of Underwater number.

Description

A kind of collaboration parsing reentry guidance method considering more no-fly zone constraints

Technical field

The present invention relates to a kind of collaborations for considering more no-fly zone constraints to parse reentry guidance method, belongs to space technology, force Device technology, Guidance and control field.

Background technique

Collaboration formation flight technology is of great significance for promoting system fighting efficiency, gathers around and has broad application prospects, Extensive and in-depth research has been obtained in unmanned plane field at present.But it for hypersonic glide vehicle, forms into columns and flies Row is very difficult.Because such aircraft initial velocity is very high but unpowered, need accurately to manage flight by crossrange maneuvering appropriate Energy.But for the aircraft interfered with different sliding states, in various degree, velocity attenuation curve and crossrange maneuvering Amplitude is all different, thus is difficult to form the more stable flight formation in relative position.

In collaboration formation flight technical aspect, common formation form has wedge team, echelon, rank, column and V-arrangement etc., can To realize the complex tasks such as coordinated investigation, defence and attack.In addition, being compiled by generating the ascending air in vortex using leader Team's flight can also effectively save fuel, increase voyage.But restricted by the time-varying characteristics of speed and crossrange maneuvering amplitude, There is presently no the theory and methods that can successfully manage hypersonic glide vehicle collaboration formation flight problem.

Yu W., Chen W., Analytical entry guidance for no-fly-zone avoidance (needle The parsing reentry guidance rule that no-fly zone is evaded) [J] Aerospace Science and Technology, 2018,72:426- 442. propose with high-precision earth rotation compensation model, obtain performance preferably three-dimensional resolution Value of Reentry Vehicle, so The tilt for devising incremental based on analytic solutions afterwards inverts Sequence Planning scheme, realizes with less tilt reversion cost reply The target of a large amount of no-fly zone constraint.Document one is referred to as in following explanation of the invention.

Gaxley C., Atmospheric entry (atmosphere reenters) [M] .The RAND Corporation, 1960. For fixed earth model, the flight time analytic solutions under the conditions of the lift resistance ratio of constant value longitudinal direction are given.In following theory of the invention Document two is referred to as in bright.

Summary of the invention

The purpose of the present invention is to solve there is presently no can successfully manage hypersonic glide vehicle collaboration to compile The problem of team's flight, proposes a kind of collaboration parsing reentry guidance method for considering more no-fly zone constraints, solves in more no-fly zone rings Multiple hypersonic glide vehicles are guided in border while reaching the collaboration flight problem of target.

Main contents of the invention are divided into three aspects: (1) for rotation earth model, when having derived high-precision flight Between analytic solutions；(2) it is directed to single aircraft, devises the parsing reentry guidance side for considering more no-fly zones and arrival time constraint Method；(3) multiple aircraft are directed to, is devised and is guided multiple hypersonic glide vehicles in more no-fly zone environment while reaching mesh Target cooperates with flight scenario.Wherein, when parsing the design of reentry guidance method, controlling terminal speed and arrival time are by longitudinal direction Lift resistance ratio control, therefore be a drive lacking problem in control, it is devised herein using vertical journey analytic solutions and time resolution solution A kind of longitudinal reference section planing method that can take into account energy management requirements and arrival time requirement, can effectively solve this and ask Topic.

A kind of collaboration parsing reentry guidance method for considering more no-fly zone constraints of the present invention, including the following steps:

Step 1: the description of reentry guidance problem

Here particle dynamics model is reentered using the three-dimensional under the spheroidal rotation earth, wherein utilizing longitude λ, latitude φ Aircraft position information is described with height H, and describes velocity information with rate V, trajectory tilt angle γ and course angle ψ.

High flying speed causes power thermal environment extremely severe.Therefore, in order to guarantee the subsystems of aircraft just Often work, ablated configuration track needs to meet heat flow densityIncoming flow dynamic pressure q, overload n constraint.In view of flight control system ability It is limited, it needs to limit the angle of attack and tilt angular region and change rate.In addition to conventional constraint, multiple circles are further considered here The constraint of shape no-fly zone.

Re-entry flight is S in the horizontal distance of aircraft to target_TAEMWhen terminate.Desired terminal height is at this time H_TAEM, terminal velocity V_TAEM, terminal course error | Δ ψ_TAEM| and terminal angle of heel | σ_TAEM| meet constraint condition.In addition, this In also specify desired terminal juncture be t_TAEM。

Step 2: the flight time analytic solutions based on rotation earth model

In order to cope with arrival time constraint, flight time analytic solutions are derived here.Energy definition is as follows

Wherein, V is speed, and H is height, R_eIt is earth mean radius, μ is gravitational constant.Energy is to the derivative of time t

G is acceleration of gravity in formula.In the case where rotating earth background, there is the nonlinear equation of following complexity

Wherein, γ is trajectory tilt angle, and D is resistance, and m is quality, ω_eIt is rotational-angular velocity of the earth, φ is latitude, and ψ is boat To angle.The 3rd, the right in formula (4)4th It is the inertia force component as caused by earth rotation.It will be considered as herein along the inertia force component of directional velocity as additional drag, it is as follows

Defining equivalent drag isFormula (3)~(5) are substituted into formula (2), and take its inverse, can be obtained

Defining longitudinal lift resistance ratio isWherein L₁=Lcos (σ) is component of the lift in fore-and-aft plane, rule It drawsFor

Wherein a₀, a₁, a₂For quadratic polynomial coefficient.

Now withSection estimates that the equivalent drag under steady gliding state (is denoted as).The trajectory tilt angle time leads Number is as follows

Assuming thatIt can estimate steady longitudinal lift that glides, it is as follows

Wherein, Δ L₁It is the component of inertia force vertically, is considered as additional longitudinal lift, formula is as follows

At the same time, using above-mentioned it is assumed that Δ D can also be reduced to

ThenIt can be estimated by following formula

Formula (12) are substituted into formula (6), and can be obtained using E replacement V

Due to H < < R_e, therefore H can be enabled to take the median H of gliding height^*, and then enable R^*=R_e+H^*.Although inertia force is small Amount, but when speed is close to the first universal speed, inertia force is possible to will cause above formula denominator to be zero, to occur unusual. In order to avoid the above-mentioned singularity of generation, inertia force is regarded as in a small amount here, and carry out single order Taylor expansion, it is as follows

Wherein, simplify symbol h_z1、h_mIt is defined as follows

h_m=2E+ μ/R^* (16)

According to Δ L₁It, can be by h with the expression formula of Δ D_z1It is divided into two part h_PAF1And h_PAF2, as follows

h_z1=h_PAF1+h_PAF2 (17)

Wherein,

h_PAF1=-2R^*Vω_e cos(φ)sin(ψ) (18)

Due to h_PAF2Influence to be much smaller than h_PAF1, therefore, can be roughly using linear function to h_PAF2Approximation is carried out, It is as follows

h_PAF2(E)≈K_PAF2(1)V+k_PAF2(0) (20)

Wherein, coefficient k_PAF2(1)And k_PAF2(0)It is determined by two endpoints, it is as follows

Wherein, h_PAF2(E) formula (19) are substituted by current state to be calculated, and h_PAF2(E_TAEM) it is by desired terminal State substitutes into formula (19) and is calculated.

In fact, hypersonic glide vehicle surrounds some the great circle target for being denoted as gen-eralized equators, aircraft Around gen-eralized equators or so crossrange maneuvering, the azimuth of gen-eralized equatorsCertain weighting that can be regarded as aircraft course angle is flat Mean value, therefore, in the case where subsequent course angular curve is even unknown, the present invention utilizes the azimuth of gen-eralized equatorsSubstitution h_PAF2(E) and h_PAF2(E_TAEM) course angle in formula.So far, in formula (14), other than independent variable E, other parameters are equal For constant value, and then formula (14) can be integrated.It largely derived, arranged, flight time parsing solution's expression can be obtained It is as follows

Wherein, coefficient k_t(1)、k_t(2)、k_t(3)、k_t(4)、k_t(5)、k_t(6)、k_t(7)Expression formula it is as follows

Step 3: considering the parsing reentry guidance method of more no-fly zones and arrival time constraint

This step is directed to single aircraft, and design considers the reentry guidance method of more no-fly zones and arrival time constraint.According to Ballistic trajectory character, reentering process, one is divided into three phases: descending branch, steady glide phase and height adjusting stage.

S31: descending branch

Cause heat flow density excessive in order to avoid falling into dense atmosphere, descending branch aircraft with maximum can with the angle of attack, Zero angle of heel glides.When lift is enough that aircraft is supported steadily to glide, the angle of attack is smoothly transitted into the benchmark angle of attack, enters immediately flat Steady glide phase.

S32: steady glide phase

Steady glide phase is longest, most important and most complicated ablated configuration stage, the Celestial Guidance Scheme benefit in this stage With three-dimensional resolution Value of Reentry Vehicle and flight time the analytic solutions ginseng that quickly planning meets no-fly zone online and the arrival time constrains Examine trajectory.The guidance process design in equilibrium glide stage is as follows:

One S321, design benchmark angle of attack section, and determine corresponding benchmark lift resistance ratio section；

S322, consider energy management requirements and arrival time constraint, rationally design longitudinal lift resistance ratio of parametrization Section, while can determine the lateral lift resistance ratio of parametrizationSection；

S323, the influence in order to compensate for earth rotation,WithOn the basis of section, equivalent vertical, horizontal is established Lift resistance ratio (With) section；

S324, it is required according to energy management and arrival time, utilizes the vertical journey analytic solutions and the step 2 in document one The flight time analytic solutions arrived solve longitudinal profile parameter；

S325, every 60s, using the horizontal journey analytic solutions in vertical journey analytic solutions and document one according to no-fly zone and terminal location Constraint updates primary tilt reversion sequence.

S326, the equivalent longitudinal lift resistance ratio section of benchmark tilt angle tracking, and tilt reversion according to plan are adjusted；

S327, in order to inhibit trajectory to vibrate, by trajectory damping feedback be introduced into the benchmark angle of attack and angle of heel, to obtain It guidances command.

S328, it is instructed according to process constraints limitation angle of heel；

S329, it is repeated the above process since step S322, until steady glide phase terminates.

S33: height adjusting stage

(remember that the corresponding aircraft energy of this node is E in last time rollback point_BR(nTR)) after, aircraft entry altitude Adjusting stage, but (be denoted as corresponding aircraft energy is E to some node before last time rollback point_ST), aircraft The preparation of height adjusting stage is begun to.Here the present invention devises a kind of multiple target number based on online Ballistic Simulation of Underwater It is worth iteration programme.In this scheme, utilization orientation derivative improves quasi-Newton method iterative algorithm, so that it is imitative to reduce trajectory True number can greatly improve iterative convergence speed.The process of the multiple target iterative numerical programme are as follows:

S331, terminal time and the SOT state of termination are predicted using Ballistic Simulation of Underwater；

S332, judge whether terminal velocity or terminal time meet the requirements；If then jump procedure S333, step is otherwise jumped Rapid S334；

S333, required according to terminal velocity and terminal time, using Quasi-Newton iterative method adjust longitudinal lift resistance ratio section and Last time rollback point, jump procedure S331；

S334, judge whether terminal height meets the requirements；If so, jump procedure S336, otherwise jump procedure S335；

S335, benchmark angle of attack parameter, jump procedure S331 are corrected according to terminal height error；

S336, iteration planning is completed.

Work as E=E_STWhen, this method of guidance plans that longitudinal lift resistance ratio, residue fly using this multiple target iterative numerical programme The reference section of row distance and flight time relative energy, and subsequent track is really finely tuned, to meet terminal time, speed, height Degree requires.Due to only needing to be emulated for several times to subsequent a bit of trajectory, so not will cause biggish computation burden.When E_st> E > E_BR(nTR)When, this method of guidance is tracked in E=E_STThe longitudinal direction that moment utilizes the programme planning of multiple target iterative numerical to obtain The reference section of lift resistance ratio, remaining flying distance and flight time relative energy, and work as E_BR(nTR)When > E, then proportion of utilization Guidance law determines benchmark angle of heel, to eliminate course error, and tracks flight time and remaining flying distance by the fine tuning angle of attack Reference section, to guarantee to meet terminal time and rate request.

Step 4: multi-aircraft arrival time cooperative approach

In the case of one, the time that aircraft reaches the destination is influenced by two key factors: launch time and being reentered To the control of time in flight course, wherein the former for adjusting the arrival time on a large scale, and the latter is dry for overcoming It disturbs, precise fine-adjustment is carried out to the time.But due to the present invention it is not intended that motors in boost phase penetration and push section flight course, herein Launch time is substituted using the initial time of reentry stage, substitutes the arrival time using the terminal time of reentry stage.In actual conditions In, then need to comprehensively consider this three sections flight course.Although the uncertain factor during boosting can cause reentry stage to originate The deviation at moment and initial state, but utilize method of guidance described above that can effectively overcome these in ablated configuration partially Poor factor.

Common time coordination problem has: same target is synchronously arrived at, it is asynchronous to reach same target, and synchronously arrive at different mesh Mark, asynchronous arrival different target and the combinatorial problem between them.Since the countermeasure of these problems is similar, below with Cooperative approach is introduced for the problem of synchronously arriving at same target.

Assuming that there is n_HGVA aircraft, is denoted as HGV_i, i=1,2 ..., n_HGV.After given launch point and target point, root According to booster rocket performance and selected powered phase guidance scheme, the reentry stage initial state of aircraft can be determined.Later, according only to Energy management requirements determine corresponding longitudinal lift resistance ratio sectional parameter using vertical journey analytic solutions.In turn, offline Ballistic Simulation of Underwater is utilized It can predict the flight time t of all aircraft_EF(HGVi), i=1,2 ..., n_HGV, when therefrom selecting a time longest flight Between t_EF(max), and reserved launch preparation time t_preWith assisting flight time t_boost, can determine desired arrival time t_TAEM= t_pre+t_boost+t_EF(max), and then utilize t_TAEM-t_EF(HGVi)It can determine the reentry stage initial time of each aircraft.Later, root According to selected boostphase guidance method, corresponding launch time can be determined.

The present invention has the advantages that

(1) Exact of the hypersonic glide trajectories flight time based on rotation earth model, prediction have been invented The case where error is maintained within 3%, is applicable to non-constant value longitudinal direction lift resistance ratio section.And traditional prediction technique is not due to examining The case where considering geocyclic influence, and being only applicable to constant value longitudinal direction lift resistance ratio, error is larger；

(2) it can meet no-fly zone based on three-dimensional resolution Value of Reentry Vehicle and flight time analytic solutions, online quickly planning, support Up to the reference trajectory of the multiple constraints such as time, final energy, solves and guide multiple hypersonic glidings in more no-fly zone environment Aircraft reaches the collaboration flight problem of target simultaneously.

(3) it devises and considers that the multiple target numerical value based on online Ballistic Simulation of Underwater of terminal time, speed and requirement for height changes For programme.Meanwhile in order to mitigate computation burden, utilization orientation derivative improves quasi-Newton method, reduces Ballistic Simulation of Underwater time Number.

Detailed description of the invention

Fig. 1 is workflow schematic diagram of the present invention；

Fig. 2 is benchmark lift resistance ratio section curve；

Fig. 3 is the multiple target iterative numerical programme flow chart of height adjusting stage；

Fig. 4 is the floor projection of aerial vehicle trajectory；

Fig. 5 is height-rate curve of aircraft；

Fig. 6 is the angle of attack-time graph of aircraft；

Fig. 7 is angle of heel-time graph of aircraft；

Fig. 8 is remaining flying distance-energy curve of aircraft；

Fig. 9 is the flight time-energy curve of aircraft.

Figure 10 a, b are the aircraft reentry trajectories of Monte Carlo emulation；

Figure 11 a, b are aircraft altitude-rate curves of Monte Carlo emulation；

Figure 12 a, b are the Aircraft Angle of Attack curves of Monte Carlo emulation；

Figure 13 a, b are the aircraft tilt angular curves of Monte Carlo emulation；

Figure 14 a, b are the aircraft terminal speed and terminal course error distribution situation of Monte Carlo emulation；

Figure 15 a, b are the aircraft arrival time distribution situations of Monte Carlo emulation.

Specific embodiment

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

Main contents of the invention are divided into three aspects: (1) for rotation earth model, when having derived high-precision flight Between analytic solutions；(2) it is directed to single aircraft, devises the parsing reentry guidance side for considering more no-fly zones and arrival time constraint Method；(3) multiple aircraft are directed to, is devised and is guided multiple hypersonic glide vehicles in more no-fly zone environment while reaching mesh Target cooperates with flight scenario.Workflow schematic diagram of the present invention is as shown in Figure 1, wherein when parsing the design of reentry guidance method, Controlling terminal speed and arrival time are controlled by longitudinal lift resistance ratio, therefore are a drive lacking problem in control, herein benefit A kind of longitudinal ginseng that can take into account energy management requirements and arrival time requirement is devised with vertical journey analytic solutions and time resolution solution Section planing method is examined, can effectively solve the problems, such as this.

Whole process including the following steps:

Step 1: the description of reentry guidance problem

Here particle dynamics model is reentered using the three-dimensional under the spheroidal rotation earth, wherein utilizing longitude λ, latitude φ Aircraft position information is described with height H, and describes velocity information with rate V, trajectory tilt angle γ and course angle ψ.

High flying speed causes power thermal environment extremely severe.Therefore, in order to guarantee the subsystems of aircraft just Often work, ablated configuration track needs to meet heat flow densityIncoming flow dynamic pressure q, overload n constraint.In view of flight control system ability It is limited, it needs to limit the angle of attack and tilt angular region and change rate.In addition to conventional constraint, multiple circles are further considered here The constraint of shape no-fly zone.

Re-entry flight is S in the horizontal distance of aircraft to target_TAEMWhen terminate.Desired terminal height is at this time H_TAEM, terminal velocity V_TAEM, terminal course error | Δ ψ_TAEM| and terminal angle of heel | σ_TAEM| meet constraint condition.In addition, this In also specify desired terminal juncture be t_TAEM。

Step 2: the flight time analytic solutions based on rotation earth model

In order to cope with arrival time constraint, flight time analytic solutions are derived here.Energy definition is as follows

Wherein, V is speed, and H is height, R_eIt is earth mean radius, μ is gravitational constant.Energy is to the derivative of time t

G is acceleration of gravity in formula.In the case where rotating earth background, there is the nonlinear equation of following complexity

Wherein, γ is trajectory tilt angle, and D is resistance, and m is quality, ω_eIt is rotational-angular velocity of the earth, φ is latitude, and ψ is boat To angle.The 3rd, the right in formula (4)4th It is the inertia force component as caused by earth rotation.It will be considered as herein along the inertia force component of directional velocity as additional drag, it is as follows

Defining equivalent drag isFormula (33)~(35) are substituted into formula (32), and take its inverse, can be obtained

Defining longitudinal lift resistance ratio isWherein L₁=Lcos (σ) is component of the lift in fore-and-aft plane, and Plan that longitudinal lift resistance ratio is as follows,

Wherein a₀, a₁, a₂For quadratic polynomial coefficient.

Now withSection estimates that the equivalent drag under steady gliding state (is denoted as).The trajectory tilt angle time leads Number is as follows

Assuming that 0 He of γ ≈It can estimate steady longitudinal lift L that glides_1(SG), as follows

Wherein, Δ L₁It is the component of inertia force vertically, is considered as additional longitudinal lift, formula is as follows

At the same time, using above-mentioned it is assumed that Δ D can also be reduced to

The then equivalent drag under steady gliding stateIt can be estimated by following formula

Formula (42) are substituted into formula (36), and can be obtained using E replacement V

Due to H < < R_e, therefore H can be enabled to take the median H of gliding height^*, and then enable average the earth's core away from R^*=R_e+H^*.Although Inertia force is a small amount of, but when speed is close to the first universal speed, inertia force is possible to will cause above formula denominator to be zero, thus Occur unusual.In order to avoid the above-mentioned singularity of generation, inertia force is regarded as in a small amount here, and carry out single order Taylor expansion, It is as follows

Wherein, simplify symbol h_z1、h_mIt is defined as follows

h_m=2E+ μ/R^* (46)

According to Δ L₁It, can be by h with the expression formula of Δ D_z1It is divided into two part h_PAF1And h_PAF2, as follows

h_z1=h_PAF1+h_PAF2 (47)

Wherein,

h_PAF1=-2R^*Vω_e cos(φ)sin(ψ) (48)

Due to h_PAF2Influence to be much smaller than h_PAF1, therefore, can be roughly using linear function to h_PAF2Approximation is carried out, It is as follows

h_PAF2(E)≈k_PAF2(1)V+k_PAF2(0) (50)

Wherein, coefficient k_PAF2(1)And k_PAF2(0)It is determined by two endpoints, it is as follows

Wherein, h_PAF2(E) formula (49) are substituted by current state to be calculated, and h_PAF2(E_TAEM) it is by desired terminal State substitutes into formula (49) and is calculated.

In fact, hypersonic glide vehicle surrounds some the great circle target for being denoted as gen-eralized equators, aircraft Around gen-eralized equators or so crossrange maneuvering, the azimuth of gen-eralized equatorsCertain weighting that can be regarded as aircraft course angle is flat Mean value, therefore, in the case where subsequent course angular curve is even unknown, the present invention utilizes the azimuth of gen-eralized equatorsSubstitution h_PAF2(E) and h_PAF2(E_TAEM) course angle in formula.So far, in formula (44), other than independent variable E, other parameters are equal For constant value, and then formula (44) can be integrated.It is as follows that flight time parsing solution's expression can be obtained through derivation

Wherein, coefficient k_t(1)、k_t(2)、k_t(3)、k_t(4)、k_t(5)、k_t(6)、k_t(7)Expression formula it is as follows

The Comparative result analysis such as embodiment of flight time analytic solutions and conventional method and numerical value Ballistic Simulation of Underwater of the invention Shown in one.

Step 3: considering the parsing reentry guidance method of more no-fly zones and arrival time constraint

This step is directed to single aircraft, and design considers the reentry guidance method of more no-fly zones and arrival time constraint.According to Ballistic trajectory character, reentering process, one is divided into three phases: descending branch, steady glide phase and height adjusting stage.Descending branch Elemental height is higher, and atmosphere is thin, and aircraft quickly falls height.In steady glide phase, atmospheric density is moderate, and lift is enough to balance Gravity, gliding height gently decline.Have benefited from high-lift and high speed, aircraft this stage gliding distance up to public affairs up to ten thousand In.When aircraft distance objective is close enough, the entry altitude adjusting stage.Aircraft is obtained by suitably adjusting the angle of attack at this time It is expected that flying height.Compared to method of guidance (document one) before, main here there are two the improvement of aspect: (1) according to energy Buret reason requires and the arrival time requires, and the planning side of reference section is improved using flight time analytic solutions and vertical journey analytic solutions Method；(2) it in order to which the track to the height adjusting stage carries out precise fine-adjustment, according to terminal time, speed, highly constrained, devises Iterative numerical programme and corresponding tracking strategy.The Celestial Guidance Scheme in these three stages is discussed in detail in below step 4~6.

Step 4: the design of descending branch method of guidance

Cause heat flow density excessive in order to avoid falling into dense atmosphere, descending branch aircraft with maximum can with the angle of attack, Zero angle of heel glides.When lift is enough that aircraft is supported steadily to glide, the angle of attack is smoothly transitted into the benchmark angle of attack, enters immediately flat Steady glide phase.

Step 5: steady glide phase method of guidance design

Steady glide phase is longest, most important and most complicated ablated configuration stage, the Celestial Guidance Scheme benefit in this stage With three-dimensional resolution Value of Reentry Vehicle and flight time the analytic solutions ginseng that quickly planning meets no-fly zone online and the arrival time constrains Examine trajectory.The guidance process in equilibrium glide stage is as follows:

S51 designs a benchmark angle of attack section, and determines corresponding benchmark lift resistance ratio section；

S52 considers energy management requirements and arrival time constraint, rationally designs longitudinal lift resistance ratio of parametrizationIt cuts open Face, while can determine the lateral lift resistance ratio of parametrizationSection；

S53 in order to compensate for earth rotation influence,WithOn the basis of section, equivalent vertical, horizontal liter is established Resistance than (With) section；

S54 is required according to energy management and arrival time, using in document one vertical journey analytic solutions and the step 2 obtain Flight time analytic solutions solve longitudinal profile parameter；

S55 is every 60s, using the horizontal journey analytic solutions in vertical journey analytic solutions and document one according to no-fly zone and terminal location Constraint updates primary tilt reversion sequence.

S56 adjusts the equivalent longitudinal lift resistance ratio section of benchmark tilt angle tracking, and tilt reversion according to plan；

Trajectory damping feedback is introduced into the benchmark angle of attack and angle of heel, in order to inhibit trajectory to vibrate to be made by S57 Lead instruction.

S58 is according to process constraints limitation angle of heel instruction；

S59 is repeated the above process since step S52, until steady glide phase terminates.

(1) the benchmark angle of attack and benchmark lift resistance ratio

In order to play aircraft maximum capacity, following benchmark angle of attack is taken_bslSection

Wherein, E_α=-5.55 × 10⁷J/kg is located at the handing over to the next shift near a little of steady gliding and height adjusting stage.E_TAEMIt is again Enter the desired final energy of section.α₁=10 ° be lift resistance ratio maximum when angle of attack value.In the height adjusting stage, in order to it is expected Terminal height, the benchmark angle of attack is gradually decrease to angle of attack₂.Here, α₂6 ° are tentatively set as, and adjusts rank in entry altitude Duan Qianhui carries out precise fine-adjustment to it.In order to maintain gliding, artificial limitation α under big disturbed condition_bsl≥5°.With the benchmark angle of attack α_bslCorresponding benchmark lift resistance ratioSection is as shown in Figure 2.

(2) reference section is improved

Since flight energy management and arrival time control intercouple, mutually restrict, merely with longitudinal lift resistance ratio section Controlling terminal energy and time are not easy to.It is largely explored, present invention discover that longitudinal lift resistance ratio section of following formEffective adjusting to the arrival time can be realized under the premise of guaranteeing to meet final energy requirement.

Wherein, E represents the energy at current time, x_ERepresent the energy of any time.Function f_L1/D(1)(x_E)、f_L1/D(2) (x_E)、f_L1/D(3)(x_E) be defined as follows

WithIt is the corresponding longitudinal lift resistance ratio of three energy nodes and section Three parameters, wherein the first two parameter is for management energy and controls the flight time, can hereinafter be asked using analytic solutions Solution.In order to enable angle of heel finally converges to zero, third parameter is enabledWherein,It is the corresponding ideal lift resistance ratio of the expectation SOT state of termination, coefficientFor mending roughly It repays and pneumatically draws inclined influence,It is the currently practical lift resistance ratio obtained using Pneumatic Identification technology, andIt is The corresponding ideal lift resistance ratio of current state.Hereinafter, when being not specific to the corresponding functional value of a certain energy node, for convenience of rising See the writing for saving independent variable, such asIt is abbreviated as

After determining benchmark lift resistance ratio section and longitudinal lift resistance ratio section, the cross of parametrization can be determined by geometrical relationship To lift resistance ratio sectionIt is as follows

Wherein, n_RIt is the reversion number having occurred and that, sgn is sign function, for determining the initial tilt side of aircraft To,It isModulus value, estimated by following formula

In order to compensate for the influence of earth rotation, inertia force and aerodynamic force group are combined into equivalent aerodynamic force.Pass through reasonable analysis The changing rule of inertia force, above-mentionedWithOn the basis of section, following inverse proportion function combining form can establish Equivalent longitudinal lift resistance ratio sectionLateral lift resistance ratio section

Here,

Wherein, h_mAs shown in formula (46),WithIt is the approximation public affairs of the inertia force group item of some complexity Formula is as follows

The coefficient k of approximate formula_h1(0)、k_h1(1)、k_h2(0)、k_h2(1)、k_h3(0)、k_h3(1)、k_h3(2)、k_h4(0)、k_h4(1)By current shape State and the desired SOT state of termination determine that calculation method is as follows.

1) design factor k_h1(0)And k_h1(1)

Wherein

It is the current longitudinal lift resistance ratio obtained in a upper guidance period, is desired terminal lift resistance ratio, φ_TAEMIt is mesh Target latitude,It is the gen-eralized equators azimuth course angle of aircraft,It is the gen-eralized equators azimuth course angle of aircraft In the value of terminal point.

2) design factor k_h2(0)And k_h2(1)

Wherein,

3) design factor k_h3(0)、k_h3(1)And k_h3(2)

4) design factor k_h4(0)And k_h4(1)

Wherein

(3) longitudinal profile parameter calculation

Here, the present invention is required according to energy management and arrival time, utilizes vertical journey analytic solutions and flight time analytic solutions To determine the parameter in formula (62)With

Remember that A is longitudinal lift resistance ratio multinomial of one form, i.e. the coefficient array of formula (37) is as follows

It then using the vertical journey analytic solutions formula in document one, can obtain under longitudinal lift resistance ratio control herein, aircraft is from energy The vertical journey of node E0 to energy node E is

Wherein, function f_xD(1)(E, E₀, A), f_xD(2)(E, E₀, A) and f_xD(3)(E, E₀, A) and it is defined as follows

Remember A₁、A₂、A₃And A₄Multinomial coefficient array respectively relevant to formula (62)-(63) is as follows

Then for longitudinal direction lift resistance ratio section shown in formula (62), predictable terminal indulges journey and is

Meanwhile for longitudinal direction lift resistance ratio section shown in formula (62), it is as follows that terminal time can be predicted using formula (53)

The then parameter in formula (62)WithIt can be obtained by solving following linear equation in two unknowns group

Wherein, s_goIt is remaining flying distance, t_TAEMThe desired arrival time.

But it is on the one hand limited to the restriction of energy management requirements, on the other hand fails the influence for fully considering crossrange maneuvering, The solution of equation group (90) can cause angle of heel substantially to vibrate, and interfere no-fly zone to evade, be unfavorable for practical application.Therefore, of the invention It proposes to overcome the time deviation as caused by crossrange maneuvering by substantially adjusting longitudinal lift resistance ratio section.Through exploring, to longitudinal liter Resistance is modified than sectional parameter, is enabled

Wherein,WithIt is the ideally corresponding longitudinal profile parameter of nominal trajectory, WithIt is then the parameter for needing to be modified according to the actual situation.Formula (91) are substituted into (88) and (89), and are arranged It can obtain

Wherein, x_D(bsl)(E_TAEM, E) and t (E_TAEM, E) and it is the corresponding part of nominal trajectory, and Δ x_D(E_TAEM, E) and Δ t (E_TAEM, E) and it is then amendment to actual conditions, formula is as follows

Under nominal case,WithThe vertical journey of aircraft should be made to meet end conswtraint, i.e.,

x_D(bsl)(E_TAEM, E) and=s_go-S_TAEM (97)

Then enableFormula (93) substitution formula (97) can be solved

Note that following integral property is utilized in the expression formula in above formula denominator

On the other hand,WithFor eliminating terminal time error delta t_f, but additional boat cannot be caused Journey.It can be obtained using formula (94) and (96)

It solves

Wherein, it is abbreviated symbol D₀Be expressed as follows

In order to predict terminal time error delta t_f, prior segregation reasons nominal trajectory is needed, and obtain two reference sections, That is range section s_go(ref)(E) and time section t_(ref)(E), wherein longitudinal lift resistance ratio section of nominal trajectory is according only to energy Management requires to determine, and terminal time is then controlled by launch time.But the section of actual conditions may be s_go(E) and t (E).Cause This prediction terminal time error delta t_fFirst component part be Δ t_f(1)=t (E)-t_(ref)(E), second component part be then It is by error in the voyage Δ s_go=s_go(E)-s_go(ref)(E) Δ t caused by_f(2).It is derived under benchmark control now by Δ s_goCause Time error Δ t_f(2).Formula (98) substitution formula (95) can be obtained

It can be obtained using difference formula by Δ s_goCaused time error Δ t_f(2)

So far, the work of longitudinal profile parameter calculation is completed.

(4) tilt reversion Sequence Planning

Longitudinal reference section is tracked since hypersonic glide vehicle is mainly adjusted by angle of heel modulus value, so true After fixed longitudinal direction lift resistance ratio section, angle of heel modulus value section determines substantially.Therefore, aircraft is only capable of by changing tilt side in due course To control crossrange maneuvering, to evade no-fly region and finally to arrive at target.The present invention is tilted using analysis iteration programme planning Invert sequence.Due to preferentially adjusting existing rollback point reply no-fly zone constraint, reversion number needed for this strategy is less, can sufficiently send out The crossrange maneuvering ability for waving aircraft copes with a large amount of no-fly zone constraints.Strategic process is described below.

1) in order to anti-interference, artificially expand all no-fly zone radiuses so that aircraft and true no-fly zone boundary it Between retain a safe distance.In view of the cumulative effect of ballistic error, radius increment is set to off is arrived in aircraft here The linear function of no-fly zone distance, and can be gradually reduced with aircraft close to target；

2) it is constrained according only to terminal location, plans two rollback points using horizontal journey analytic solutions；

3) since nearest no-fly zone, using the analytic solutions search track points nearest apart from no-fly zone, to judge that this is no-fly Whether area's constraint is met.If it is, continuing to judge next no-fly zone, otherwise, execute 4)；

4) judging the no-fly zone is before being located at penultimate rollback point, is between most latter two rollback point or last After rollback point.If it is before being located at penultimate rollback point, then execute 5)；If it is positioned at most latter two reversion Between point, then execute 6)；After being located at last time rollback point, then execute 7)；

If 5) this no-fly zone is located at before secondary rollback point second from the bottom, show this no-fly zone apart from terminal farther out, because This can not consider that terminal location constrains.Here rollback point is adjusted first with no-fly zone evasion tactics, copes with current no-fly zone Then constraint adjusts subsequent rollback point using the horizontal process control strategy 1 of terminal, to guarantee to meet terminal location requirement.Jump execution 8).Hereinafter can to no-fly zone evasion tactics, the horizontal process control strategy 1 of terminal and 6) in the horizontal process control strategy 2 of terminal carry out Explanation；

If 6) no-fly zone is located between last rollback point twice, the reversion sequence of planning tilt at this time needs to take into account no-fly zone Constraint and terminal location requirement.Here equally first with no-fly zone evasion tactics reply no-fly zone constraint, but due at this time Distance objective is closer, subsequent rollback point can be planned using the horizontal process control strategy 2 of terminal later, to meet the horizontal range request of terminal.It jumps Turn to execute 8)；

7) if no-fly zone is located at after last time rollback point, since no-fly zone is too close apart from terminal, ignore this here No-fly zone constraint, and retain original reversion point sequence.Jump execution 8)；

8) it determines next no-fly zone, returns to and 3) continue to judge whether no-fly zone constraint is met, until having detected There is no-fly zone constraint.

Illustrate the principle of the horizontal process control strategy 1 of no-fly zone evasion tactics, terminal and the horizontal process control strategy 2 of terminal below.

1) no-fly zone evasion tactics

The reply no-fly zone constraint in two kinds of situation of this strategy: current no-fly zone is located at before secondary rollback point second from the bottom；When Preceding no-fly zone constraint is located between last rollback point twice.

In the first case, farther out due to distance objective, without considering during handling current no-fly zone and constraining Terminal location constraint.The present invention preferentially adjusts existing rollback point using analytic solutions, so that aircraft is from close to the one of no-fly zone Side bypasses no-fly zone.If it fails, then further attempting to get around from the other side of no-fly zone, or increase a rollback point newly Cope with no-fly zone constraint.After successfully evading current no-fly zone, since reversion sequence adjusts the track before changing, so Whether the no-fly zone constraint before needing to check still meets.If conditions are not met, then further judge which threat level is higher, To decide whether to adopt new reversion point sequence.

In second situation, since distance objective is closer, rollback point Sequence Planning needs to take into account no-fly zone and evades at this time It is constrained with terminal location.The present invention preferentially adjusts existing rollback point reply no-fly zone constraint, but the difference is that will make to fly as far as possible Row device bypasses no-fly zone by the side of close-target along no-fly zone.Subsequent process is same as above.

2) the horizontal process control strategy 1 of terminal

The case where this strategy is located at before secondary rollback point second from the bottom primarily directed to current no-fly zone.Utilizing no-fly zone After evasion tactics achievement copes with current no-fly zone constraint, since track adjusts, subsequent reversion point sequence is unsatisfactory for terminal position Set constraint.Therefore, the horizontal process control strategy 1 of terminal will abandon subsequent rollback point, and be constrained according to terminal location, plan again most It inverts twice afterwards.

3) the horizontal process control strategy 2 of terminal

The case where this strategy is located between last rollback point twice primarily directed to the constraint of current no-fly zone.Equally, at After function calls no-fly zone evasion tactics, subsequent rollback point is abandoned.Since distance objective is closer at this time, adjust first here most A rollback point reply terminal location constraint afterwards, then judges whether the rollback point meets steady glide phase and height adjustment rank The hand-over condition of section.If conditions are not met, then increasing 1-2 reversion by the way that analysis crossrange maneuvering rule is appropriate.

(5) the benchmark angle of heel of steady glide phase

Using the geometrical relationship between longitudinal lift resistance ratio and total lift resistance ratio, benchmark angle of heel can be determined.Note n-th_RIt is secondary to incline It is E that side, which inverts corresponding energy,_BR(nR).Work as E_BR(nR)+ΔE≥E≥E_BR(nR+1)When+Δ E, benchmark angle of heel σ_bslFor

Wherein, it is appropriate in advance to make tilt reversion by offset Δ E, so that compensation response as caused by roll rate limit is stagnant Afterwards.

(6) the instruction angle of attack and angle of heel of steady glide phase

In order to inhibit trajectory to vibrate, the three-dimensional trajectory damping control technology of introducing is as follows

α_cmd=α_bsl+cos(σ_bsl)k_γ(γ_SG-γ) (106)

Wherein, α_cmdAnd σ_cmdIt is the instruction angle of attack and angle of heel, k_γIt is feedback gain, γ_SGIt is that steady glide trajectories incline Angle can be determined by current state.

It is constrained to meet heat flow density, incoming flow dynamic pressure and overload etc., process constraints can be converted to angle of heel constraint, into And limit angle of heel instruction size.

Step 6: the design of height adjusting stage method of guidance

(remember that the corresponding aircraft energy of this node is E in last time rollback point_BR(nTR)) after, aircraft entry altitude Adjusting stage, but (be denoted as corresponding aircraft energy is E to some node before last time rollback point_ST), aircraft The preparation of height adjusting stage is begun to.Here the present invention devises a kind of multiple target number based on online Ballistic Simulation of Underwater It is worth iteration programme.Work as E=E_STWhen, using multiple target iterative numerical programme obtain longitudinal lift resistance ratio, remaining flight away from From and flight time relative energy reference section, and the subsequent track of precise fine-adjustment, to meet terminal time, speed, height It is required that.Due to only needing to be emulated for several times to subsequent a bit of trajectory, so not will cause biggish computation burden.Work as E_st > E > E_BR(nTR)When, this method of guidance is tracked in E=E_STWhen by what multiple target iterative numerical programme obtained longitudinal rise resistance Than, remaining flying distance and the reference section of flight time relative energy, and work as E_BR(nTR)When > E, then proportion of utilization guides It restrains and determines benchmark angle of heel, to eliminate course error, and flight time and remaining flying distance reference are tracked by the fine tuning angle of attack Section, to guarantee to meet terminal time and rate request.The specific Celestial Guidance Scheme of height adjusting stage is as follows:

(1) multiple target iterative numerical programme

Multiple target iterative numerical programme considers terminal time, speed and highly constrained, wherein terminal time and speed Degree is by parameterAnd E_BR(nTR)Control, terminal height is by benchmark angle of attack parameter alpha₂Control.As shown in figure 3, more Target value iteration programme process is as follows:

S611, terminal time and the SOT state of termination are predicted using Ballistic Simulation of Underwater；

S612, judge whether terminal velocity or terminal time meet the requirements；If then jump procedure S613, step is otherwise jumped Rapid S614；

S613, required according to terminal velocity and terminal time, using Quasi-Newton iterative method adjust longitudinal lift resistance ratio section and Last time rollback point, jump procedure S611；

S614, judge whether terminal height meets the requirements；If so, jump procedure S616, otherwise jump procedure S615；

S615, benchmark angle of attack parameter, jump procedure S611 are corrected according to terminal height error；

S616, iteration planning is completed.

In step S613, for different situations, countermeasure of the invention is as follows:

If 1) predict that discovery terminal velocity and terminal time are greater than desired value through Ballistic Simulation of Underwater, can suitably increaseReduceAnd E_BR(nTR)；

If 2) terminal velocity is greater than desired value and terminal time is less than desired value, can suitably reduceAnd E_BR(nTR)；

If 3) terminal velocity is less than desired value and terminal time is greater than desired value, can suitably increaseAnd E_BR(nTR)；

If 4) terminal velocity and terminal time are respectively less than desired value, can suitably reduceIncreaseAnd E_BR(nTR)。

Due toWith E_BR(nTR)Synchronization variation can incite somebody to actionAs an auxiliary adjustment means, and design such as Lower relational expression

Wherein,With Δ E_BR(nTR)It respectively indicatesWith E_BR(nTR)Variable quantity, k_f=2 × 10^-7To adjust Save coefficient.If above formula indicates E_BR(nTR)Change 0.1 × 10⁷KJ/kg, thenIt is corresponding to change 0.2.One is reduced in this way Dimension is searched for, helps to reduce algorithm complexity.

Meet terminal time, speed and highly constrained below with Newton iteration method searchAnd E_BR(nTR).It is fixed Adopted vector x and vector function G

Wherein,WithRespectively parameterAnd E_BR(nTR)It is corresponding The terminal velocity and time that online Ballistic Simulation of Underwater obtains.Solve equation G (x)=0 newton iteration formula be

Wherein, subscript " (k) " represents kth time iteration, matrix J_GIt is the Jacobian matrix of vector function G, such as formula (81) It is shown, andIt is J_GInverse matrix.

Due toWithIt needs to obtain by integral, J_GIt can not Analytical Solution. If solved using conventional difference method, 3 trajectory integrals are needed to be implemented in each step, it is negative that this will cause biggish calculating Load.

In order to mitigate computation burden, the present invention proposes a kind of J based on directional derivative_GSolution scheme.Definition vector p^(k)= x^(k)-x^(k-1), definition vector q^(k)It is perpendicular to p^(k), mould a length of 0.01 a small vector.Then there is the approximation of following directional derivative Calculation formula

Wherein, | | p^(k)| | and | | q^(k)| | respectively indicate vector p^(k)And q^(k)Mould it is long.Utilization orientation derivative and partial derivative Between relationship can obtain

Wherein, θ_x(p^(k)) and θ_y(p^(k)) it is vector p^(k)Respectively with the angle of x-axis, y-axis, and θ_x(q^(k)) and θ_y(q^(k)) then It is vector q^(k)Respectively with the angle of x-axis, y-axis.It can be solved using formula (113)

In current iteration, due to only needing to calculate G (x^(k)+q^(k)) and G (x^(k)), only need 2 trajectories of execution imitative here Very, so that calculation amount be greatly reduced.When due to every single-step iteration,Variable quantityThe order of magnitude be about 0.1~1, and Δ E_BR(nTR)The order of magnitude be 10⁶, above-mentioned equation is usually morbid state, but can use exchange entry and solve this and ask Topic.

Benchmark angle of attack sectional parameter α then is finely tuned using following formula (85) in step S615₂.Due to formula accuracy compared with Height only needs to correct here primary.α₂The knots modification Δ α needed₂For

Wherein, C_Lf(est)It is the terminal lift coefficient obtained by Pneumatic Identification technology,Be estimation lifting line it is oblique Rate, S_refIt is aircraft area of reference, E_f、q_f、V_f、H_f、φ_f、ψ_fIt is the final energy obtained after line Ballistic Simulation of Underwater respectively, dynamic Pressure, speed, height, latitude, course angle value, q_TAEM、V_TAEMIt is desired terminal dynamic pressure and velocity amplitude respectively.

(2) the instruction angle of attack and angle of heel of height adjusting stage

Work as E_st> E > E_BR(nTR)When, aircraft is in the preparation stage before the height adjusting stage.At this time Celestial Guidance Scheme with Longitudinal lift resistance ratio section, remaining flying distance section and the flight time that track is obtained by multiple target iterative numerical programme Section is denoted as respectivelys_go(ref)(E)、t_ref(E).Benchmark angle of heel σ_bslCalculation formula it is as follows

Wherein,Here longitudinal lift resistance ratioMainly byComposition, but in order to guarantee that energy management and terminal time require under noisy condition, design and track consideration The modified remaining flying distance section of time error, it is as follows

s_go(ref)(2)(E)=s_go(ref)(E)-k_sgoV(t-t_ref(E)) (117)

Wherein, k_sgo=0.2s_go/s_go(ST)。s_go(ST)It is to execute multiple target iterative numerical programme moment corresponding residue Flying distance.

Hereafter, formula (116) are substituted into formula (107), and combines formula (106), so that it may instruct the angle of attack and angle of heel.

As E < E_BR(nTR)When, the aircraft entry altitude adjusting stage, still such as formula (61) are shown for the benchmark angle of attack, and base Then proportional guidance law determines quasi- angle of heel.Sight azimuth rate is

Wherein, Δ ψ is course error.Proportional guidance law generate need be with crossrange maneuvering acceleration

Wherein, Δ L₂It is inertia force along L₂The component in direction.Section 2-Δ L on the right of above formula₂/ m is to compensate for the earth certainly The influence turned.Initial angle of heel saturation in order to prevent, enables effectively guiding compare k here_PNIt is tapered to remaining flying distance from 2 4.On the other hand, in steady gliding, lift acceleration a in fore-and-aft plane_L1With gravity, centrifugal force and balance of shaking force, I.e.

Then benchmark angle of heel is

In turn, it designs the angle of attack and angle of heel instruction is as follows

α_cmd=α_bsl+k_α[s_go-s_go(ref)(2)(E)] (122)

Step 7: the design of multi-aircraft arrival time cooperative approach

Step 3-6 completes parsing guidance side for single aircraft, considering more no-fly zones and terminal time constraint Method design, studies arrival time cooperative approach further directed to multiple aircraft here.In the case of one, aircraft arrives at purpose The time on ground is influenced by two key factors: the control during launch time and ablated configuration to the time, wherein the former uses In being adjusted on a large scale to the arrival time, and the latter is for overcoming interference, to time progress precise fine-adjustment.But due to this The flight course it is not intended that motors in boost phase penetration and pushing section is invented, when utilizing the initial time substitution transmitting of reentry stage in the present invention Between, the arrival time is substituted using the terminal time of reentry stage.Note: needing to comprehensively consider this three sections of flight in a practical situation Process.Although the uncertain factor during boosting can cause the deviation of reentry stage initial time and initial state, utilize Previously described method of guidance can effectively overcome these bias factors in ablated configuration.

Common time coordination problem has: same target is synchronously arrived at, it is asynchronous to reach same target, and synchronously arrive at different mesh Mark, asynchronous arrival different target and the combinatorial problem between them.Since the countermeasure of these problems is similar, below with Cooperative approach is introduced for the problem of synchronously arriving at same target.

Assuming that there is n_HGVA aircraft, is denoted as HGV respectively_i, i=1,2 ..., n_HGV.Given launch point and target point it Afterwards, according to booster rocket performance and selected powered phase guidance scheme, the reentry stage initial state of aircraft can be determined.Later, Same formula (98) enables longitudinal lift resistance ratio sectional parameterAccording only to energy management requirements, benefit Corresponding longitudinal lift resistance ratio sectional parameter is determined with vertical journey analytic solutions, as shown in formula (94).Here, in subscript " HGV (i) " generation, refers to I-th of aircraft.

In turn, the flight time t of each aircraft can be predicted using offline Ballistic Simulation of Underwater_EF(HGVi), i=1,2 ..., n_HGV, the time longest one is therefrom selected, and reserved transmitting prepares and the assisting flight time, when can determine desired arrival Between t_TAEM, and then utilize t_TAEM-t_EF(HGVi)It can determine the reentry stage initial time of each aircraft.Later, according to selected Boostphase guidance method can determine corresponding launch time.

Embodiment:

Embodiment one

The present embodiment carries out flight time analytic solutions and conventional method and numerical value results of trajectory simulation of the invention and compares, and tests Demonstrate,prove the precision of prediction of flight time analytic solutions of the invention.

Flight time formula in document two under the conditions of the lift resistance ratio of constant value longitudinal direction is

The initial longitude λ of aircraft is set₀=0deg, initial latitude φ₀=50deg, primary power E₀=-3.8602 × 10⁴KJ/kg and final energy E_f=-5.5 × 10⁴kJ/kg.Consider 5 different headings: ψ₀=100deg, 180deg ,- 100deg,20deg,-20deg.The longitudinal lift resistance ratio being arranged in document two isLateral lift resistance ratio is

Simulation result is as shown in table 1.As can be seen that being influenced by earth rotation from simulation result, aircraft is along different The flight time in direction is not identical, but the analytic solutions of document two fail to consider the time difference as caused by earth rotation.It is logical It crosses and is compared with results of trajectory simulation it can be seen that the influence due to analytic solutions of the present invention to earth rotation is compensated, Its precision is higher, and prediction error is maintained within 3%.In computational efficiency, the calculating time-consuming of time resolution solution at least compares trajectory Emulate small 5 orders of magnitude.

Table 1

Embodiment two

The present embodiment is ideal glitch-free situation, and solution guides multiple hypersonic glidings in more no-fly zone environment and flies Row device (V1, V2, V3) emits from different location but reaches the collaboration flight problem of same target simultaneously.It is arranged in flight range 64 radiuses are the round no-fly zones of 200km, and the primary condition of aircraft is as shown in table 1, and aiming spot is longitude λ_T= 130deg, latitude φ_T=-20deg.Using the multi-aircraft arrival time cooperative approach of step 7, t can be chosen_TAEM=2900s, And determine each aircraft reenters initial time, respectively 109.29s, 358.69s and 751.62s.

Table 2

Re-entry flight is S in the horizontal distance of aircraft to target_TAEMIt is terminated when=50km.Desired terminal is high at this time Degree is H_TAEM=25km, terminal velocity V_TAEM=2000m/s, terminal course error meet | Δ ψ_TAEM|≤5deg and terminal are inclined Side angle meets | σ_TAEM|≤30deg。

Fig. 4 illustrates aircraft and has successfully evaded all no-fly zones under method of guidance of the present invention control, and arrives at target. Fig. 5 is height-rate curve, wherein H_minCurve is the height lower bound determined by process constraints.It can be seen from the figure that due to More severe in high speed stage power thermal environment, aircraft is easy to fly close to height lower bound.Furthermore it can also be seen that terminal velocity Highly meet demanding terminal.Fig. 6 illustrates the angle of attack rule of aircraft, it can be seen that three aircraft are with different Reenter initial time.Fig. 7 illustrates tilt angle sections.It is corresponding with reference to remaining flying distance section that Fig. 8 illustrates three kinds of situations. It is corresponding with reference to flight time section that Fig. 9 illustrates three kinds of situations.

Embodiment three

The present embodiment is by Monte Carlo simulating, verifying method of guidance of the present invention in the case where dummy vehicle has and draws inclined situation Robustness.Wherein, Aerodynamic Coefficient uses lower linear such as to draw inclined model

Wherein, δ_CLAnd δ_CDIt is that lift coefficient and resistance coefficient draw inclined percentage respectively, changes with angle of attack and Mach number Ma. δ_CL0、It is relevant zero-mean normal distribution random parameter.Wind-force and atmospheric density are taken the photograph Movable model is as follows

δ_ρ=k_ρδ_ρ(max) (130)

Wherein,It is the wind speed along east-west direction,It is wind speed along the north-south direction, δ_ρIt is that atmospheric density draws inclined percentage Than.V_wind(max)It is maximum possible wind speed, with height change, in high-altitude up to 170m/s.δ_ρ(max)It is atmospheric density maximum possible Inclined percentage is drawn, in high-altitude up to 50%.And k_ρIt is corresponding zero-mean normal distribution random parameter.It is above-mentioned with The statistical property of machine parameter and primary condition error is as shown in table 3.

Table 3

Here consider to emit the scene that two aircraft (V1, V2) carry out collaboration flight, nominal item from two different locations Part is as shown in table 4, and target position is longitude λ_T=120deg and latitude φ_T=0deg.No-fly zone distribution is as shown in Figure 10.Setting It is expected that the arrival time is 2900s, then it can have determined that reentry stage rises under nominal case using multi-aircraft arrival time cooperative approach Time beginning is 426.68s and 404.29s.Note that being remained required for nominal trajectory and Guidance Law due to being prior off-line calculation Remaining flying distance, flight time reference section are obtained based on ideal aerodynamic model and ideal atmosphere model, do not consider, It can not consider the inclined data of drawing that any on-line measurement is arrived.

Table 4

Figure 10 illustrates the reentry trajectory of aircraft, therefrom as it can be seen that there are in uncertain noises environment, this method of guidance The directing aircraft that can succeed evades no-fly zone, and smoothly arrives at target.Figure 11 illustrates height-rate curve, wherein H_minIt is The height lower bound determined by process constraints such as heat flow densities.As seen from the figure, this method of guidance can guarantee flight safety.Figure 12-13 Illustrate angle of attack curve and tilt angular curve.Figure 14 illustrates the distribution situation of terminal velocity and terminal course error, is wanting Within the scope of asking.Figure 15 has counted arrival time distribution situation, accounting of the arrival time error of aircraft V1 within ± 5s It is 79.01%；Accounting of the arrival time error of aircraft V2 within ± 5s is 74.36%.

Claims (3)

1. a kind of collaboration parsing reentry guidance method for considering more no-fly zone constraints, it is characterised in that: this method includes following several A step:

Step 1: the description of reentry guidance problem

Three-dimensional under the earth is rotated using spheroidal and reenters particle dynamics model, wherein utilizing longitude λ, latitude φ and height H Aircraft position information is described, and describes velocity information with rate V, trajectory tilt angle γ and course angle ψ；

In order to guarantee that the subsystems of aircraft work normally, ablated configuration track needs to meet heat flow densityIncoming flow dynamic pressure Q, overload n constraint；In view of flight control system ability is limited, need to limit the angle of attack and tilt angular region and change rate；It removes Conventional constraint, it is also necessary to the case where considering multiple round no-fly zones constraints；

Re-entry flight is S in the horizontal distance of aircraft to target_TAEMWhen terminate；Desired terminal height is H at this time_TAEM, eventually End speed is V_TAEM, terminal course error | Δ ψ_TAEM| and terminal angle of heel | σ_TAEM| meet constraint condition；In addition, also referring to here Fixed desired terminal juncture is t_TAEM；

Step 2: the flight time analytic solutions based on rotation earth model,

In order to cope with arrival time constraint, flight time analytic solutions are derived here；Energy definition is as follows

Wherein, V is speed, and H is height, R_eIt is earth mean radius, μ is gravitational constant；Energy is to the derivative of time t

G is acceleration of gravity in formula；In the case where rotating earth background, there is the nonlinear equation of following complexity

Wherein, γ is trajectory tilt angle, and D is resistance, and m is quality, ω_eIt is rotational-angular velocity of the earth, φ is latitude, and ψ is course angle；

The 3rd, the right in formula (4)4th It is the inertia force component as caused by earth rotation；The present invention will be considered as along the inertia force component of directional velocity as additional drag, such as Under

Defining equivalent drag isFormula (3)~(5) are substituted into formula (2), and take its inverse, can be obtained

Defining longitudinal lift resistance ratio isWherein L₁=Lcos (σ) is component of the lift in fore-and-aft plane, planningFor

Wherein a₀,a₁,a₂For quadratic polynomial coefficient；

Now withSection estimates that the equivalent drag under steady gliding state (is denoted as)；Trajectory tilt angle time-derivative is such as Under

Assuming that 0 He of γ ≈It can estimate steady longitudinal lift that glides, it is as follows

Wherein, Δ L₁It is the component of inertia force vertically, is considered as additional longitudinal lift, formula is as follows

At the same time, using above-mentioned it is assumed that Δ D can also be reduced to

ThenIt can be estimated by following formula

Formula (12) are substituted into formula (6), and can be obtained using E replacement V

Due to H < < R_e, therefore H can be enabled to take the median H of gliding height^*, and then enable R^*=R_e+H^*；Although inertia force be it is a small amount of, It is when speed is close to the first universal speed, inertia force is possible to will cause above formula denominator to be zero, to occur unusual；In order to keep away Exempt from that above-mentioned singularity occurs, here regard inertia force in a small amount as, and carries out single order Taylor expansion, it is as follows

Wherein, simplify symbol h_z1、h_mIt is defined as follows

h_m=2E+ μ/R^* (16)

According to Δ L₁It, can be by h with the expression formula of Δ D_z1It is divided into two part h_PAF1And h_PAF2, as follows

h_z1=h_PAF1+h_PAF2 (17)

Wherein,

h_PAF1=-2R^*Vω_ecos(φ)sin(ψ) (18)

Due to h_PAF2Influence to be much smaller than h_PAF1, therefore, can be roughly using linear function to h_PAF2Approximation is carried out, it is as follows

h_PAF2(E)≈k_PAF2(1)V+k_PAF2(0) (20)

Wherein, coefficient k_PAF2(1)And k_PAF2(0)It is determined by two endpoints, it is as follows

Wherein, h_PAF2(E) formula (19) are substituted by current state to be calculated, and h_PAF2(E_TAEM) it is by the desired SOT state of termination Formula (19) are substituted into be calculated；

In fact, hypersonic glide vehicle surrounds some the great circle target for being denoted as gen-eralized equators, aircraft is surrounded Gen-eralized equators or so crossrange maneuvering, the azimuth of gen-eralized equatorsIt can be regarded as certain weighted average at aircraft course angle Value, therefore, in the case where subsequent course angular curve is even unknown, the present invention utilizes the azimuth of gen-eralized equatorsSubstitute h_PAF2 (E) and h_PAF2(E_TAEM) course angle in formula；So far, in formula (14), other than independent variable E, other parameters are normal Value, and then formula (14) can be integrated；It derived, arranged, it is as follows that flight time parsing solution's expression can be obtained:

Wherein, coefficient k_t(1)、k_t(2)、k_t(3)、k_t(4)、k_t(5)、k_t(6)、k_t(7)Expression formula it is as follows

Step 3: considering the parsing reentry guidance method of more no-fly zones and arrival time constraint

This step is directed to single aircraft, and design considers the reentry guidance method of more no-fly zones and arrival time constraint；According to trajectory Feature, the process of reentering are generally divided into three phases: descending branch, steady glide phase and height adjusting stage；

S31: descending branch

Cause heat flow density excessive in order to avoid falling into dense atmosphere, descending branch aircraft can be inclined with maximum with the angle of attack, zero Side angle glides；When lift is enough that aircraft is supported steadily to glide, the angle of attack is smoothly transitted into the benchmark angle of attack, enters immediately steady sliding The Xiang stage；

S32: steady glide phase

Steady glide phase is longest, most important and most complicated ablated configuration stage, and the Celestial Guidance Scheme in this stage utilizes three Dimension resolution Value of Reentry Vehicle and the flight time analytic solutions reference bullet that quickly planning meets no-fly zone online and the arrival time constrains Road；

S33: height adjusting stage

After last time rollback point, the aircraft entry altitude adjusting stage, but certain before last time rollback point A node, aircraft have begun to the preparation of height adjusting stage；Wherein, the corresponding aircraft of note last time rollback point Energy is E_BR(nTR), it is E that some node before last time rollback point, which is denoted as corresponding aircraft energy,_ST；

Design a kind of multiple target iterative numerical programme based on online Ballistic Simulation of Underwater；In this scheme, utilization orientation derivative Improve quasi-Newton method iterative algorithm；

Work as E=E_STWhen, it is subsequent using the multiple target iterative numerical programme precise fine-adjustment based on online Ballistic Simulation of Underwater Track, to meet terminal time, speed, requirement for height；Work as E_st> E > E_BR(nTR)When, the method for the present invention tracking is planned in iteration The reference section of longitudinal lift resistance ratio of middle acquisition, remaining flying distance and flight time relative energy, and work as E_BR(nTR)> E When, proportion of utilization guidance law determines benchmark angle of heel, to eliminate course error, and tracks the flight time by the fine tuning angle of attack and remains Remaining flying distance reference section, to guarantee to meet terminal time and rate request；

Step 4: multi-aircraft arrival time cooperative approach

Assuming that there is n_HGVA aircraft, is denoted as HGV_i, i=1,2 ..., n_HGV；After given launch point and target point, according to boosting Rocket performance and selected powered phase guidance scheme, can determine the reentry stage initial state of aircraft；Later, according only to energy pipe Reason requires, and determines corresponding longitudinal lift resistance ratio sectional parameter using vertical journey analytic solutions；It in turn, can be pre- using offline Ballistic Simulation of Underwater Survey the flight time t of all aircraft_EF(HGVi), i=1,2 ..., n_HGV, therefrom select the time longest flight time t_EF(max), and reserved launch preparation time t_preWith assisting flight time t_boost, can determine desired arrival time t_TAEM= t_pre+t_boost+t_EF(max), and then utilize t_TAEM-t_EF(HGVi)It can determine the reentry stage initial time of each aircraft；Later, root According to selected boostphase guidance method, corresponding launch time can be determined.

2. a kind of collaboration parsing reentry guidance method for considering more no-fly zone constraints according to claim 1, feature exist In: the step S32: the guidance process of steady glide phase is specific as follows:

One S321, design benchmark angle of attack section, and determine corresponding benchmark lift resistance ratio section；

S322, consider energy management requirements and arrival time constraint, rationally design longitudinal lift resistance ratio of parametrizationIt cuts open Face, while can determine the lateral lift resistance ratio of parametrizationSection；

S323, the influence in order to compensate for earth rotation,WithOn the basis of section, establishes equivalent vertical, horizontal and rise resistance Than (With) section；

S324, required according to energy management and arrival time, using in document one vertical journey analytic solutions and the step 2 obtain Flight time analytic solutions solve longitudinal profile parameter；

S325, every 60s, constrained using the horizontal journey analytic solutions in vertical journey analytic solutions and document one according to no-fly zone and terminal location Update primary tilt reversion sequence；

S326, the equivalent longitudinal lift resistance ratio section of benchmark tilt angle tracking, and tilt reversion according to plan are adjusted；

S327, in order to inhibit trajectory to vibrate, by trajectory damping feedback be introduced into the benchmark angle of attack and angle of heel, to be guided Instruction；

S328, it is instructed according to process constraints limitation angle of heel；

S329, it is repeated the above process since step S322, until steady glide phase terminates.

3. a kind of collaboration parsing reentry guidance method for considering more no-fly zone constraints according to claim 1, feature exist In: the multiple target iterative numerical programme described in step S33 based on online Ballistic Simulation of Underwater, process are specific as follows:

S331, terminal time and the SOT state of termination are predicted using Ballistic Simulation of Underwater；

S332, judge whether terminal velocity or terminal time meet the requirements；If then jump procedure S333, otherwise jump procedure S334；

S333, it is required according to terminal velocity and terminal time, using the longitudinal lift resistance ratio section of Quasi-Newton iterative method adjustment and finally Rollback point, jump procedure S331；

S334, judge whether terminal height meets the requirements；If so, jump procedure S336, otherwise jump procedure S335；

S335, benchmark angle of attack parameter, jump procedure S331 are corrected according to terminal height error；

S336, iteration planning is completed.