CN104392047B - Quick trajectory programming method based on smooth glide trajectory analytic solution - Google Patents

Abstract

The invention discloses a quick trajectory programming method based on a smooth glide trajectory analytic solution. The quick trajectory programming method based on the smooth glide trajectory analytic solution includes that step 1, modeling glide trajectory programming problems; step 2, designing glide trajectory programming variables; step 3, calculating a glide trajectory analytic solution; step 4, designing a glide trajectory terminal speed control scheme; step 5, designing a glide trajectory re-entry corridor regulating proposal; step 6, generating initial values of glide trajectory programming; step 7, designing a glide trajectory programming flow. The quick trajectory programming method based on the smooth glide trajectory analytic solution uses longitudinal maneuvering acceleration proportion coefficients and transverse maneuvering acceleration proportion coefficients as the glide trajectory programming variables so that differential equations of the trajectory inclination angle, trajectory deflection angle, height, longitude and latitude in motion equations do not comprise a speed item. The quick trajectory programming method based on the smooth glide trajectory analytic solution obtains the glide trajectory analytic solution corresponding to a fixed longitudinal maneuvering acceleration proportion coefficient and a fixed transverse maneuvering acceleration proportion coefficient.

Description

A kind of quick trajectory planning method based on steady glide trajectories analytic solutions

Technical field

The present invention relates to a kind of based on the quick trajectory planning method for steadily sliding trajectory analytic solutions, belong to space technology, Weapon technologies field.

Background technology

With the fast development of hypersonic technology, reentry trajectory is formulated for for study hotspot.Gliding section is to reenter bullet The main composition part in road, determines the voyage and maneuvering range of ablated configuration.Quickly planning not only can divide gliding section trajectory The performance of analysis hypersonic aircraft, it may also be used for online trajectory planning and Predictor-corrector guidance, with higher researching value.

It is extremely sensitive to planning variable because the section trajectory that glides has the flight time long, and reenter process nonlinear restriction By force so that the feasible zone of glide trajectories is narrower, the calculating time grown is generally required very much when solving using traditional optimized algorithm.Although The optimized algorithms such as pseudo- spectrometry, SQP are greatly improved in computational efficiency, and it calculates and takes still more than second level, and works as When process constraints touch border, its computational efficiency can be substantially reduced, or even have influence on the convergence of result.In order to improve gliding section bullet Road planning speed, some special constraints be introduced into gliding section trajectory planning, such as resistance curve, heat flow curve, etc. dynamic pressure Curve, equilibrium glide condition etc..The common feature of these constraints is that gliding section trajectory is limited to into certain special steady gliding State, so as to reduce the sensitivity of trajectory planning problem.But due to also needing to carry out numerical integration in trajectory planning, limit The raising of planning speed.In order to break away from the dependence of the integration to trajectory, dynamic inverse is introduced into gliding section trajectory planning.The method is straight Planning geometry ballistic-shaped is connect, and by adjusting ballistic-shaped come meet the constraint requirement.But because the environment before and after section of gliding becomes Change violent so that gliding section overall trajectory is difficult to use low order curve matching.

The content of the invention

The purpose of the present invention is that, by solving steady glide trajectories high accuracy analytic solutions, acquisition is independent of the cunning of trajectory integration The quick planing method of Xiang section trajectory, to analyze hypersonic aircraft performance and Predictor-corrector guidance technical support is provided.

It is a kind of based on the quick trajectory planning method for steadily sliding trajectory analytic solutions, including following step：

Step 1：The problem modeling of gliding section trajectory planning；

Step 2：Gliding section trajectory planning Variational Design；

Step 3：Gliding section trajectory analytic solutions are solved；

Step 4：Gliding section endgame speed-control scheme；

Step 5：Gliding section trajectory reentry corridor Adjusted Option；

Step 6：Gliding section trajectory planning forming initial fields；

Step 7：Gliding section trajectory planning flow scheme design.

It is an advantage of the current invention that：

(1) propose with longitudinal maneuver acceleration proportionality coefficient and crossrange maneuvering acceleration proportionality coefficient as the section orbitron gage that glides The amount of wiping so that do not contain speed term in the equation of motion in the differential equation of trajectory tilt angle, trajectory deflection angle, height, longitude and latitude；

(2) fixed longitudinal maneuver acceleration proportionality coefficient and the corresponding gliding of crossrange maneuvering acceleration proportionality coefficient are obtained Section trajectory analytic solutions, including height analytic solutions, flying distance analytic solutions, trajectory deflection angle analytic solutions, longitude analytic solutions and latitude solution Analysis solution；

(3) the forming initial fields side of longitudinal maneuver acceleration proportionality coefficient and crossrange maneuvering acceleration proportionality coefficient is given Method；

(4) propose to incline with longitudinal maneuver acceleration proportionality coefficient, crossrange maneuvering acceleration proportionality coefficient, horizontal reversion trajectory Angle and initial trajectory inclination angle increment are gliding section trajectory planning variable, wherein：Longitudinal maneuver acceleration proportionality coefficient correspondence is glided Section flying distance；Crossrange maneuvering acceleration proportionality coefficient counterpart terminal speed；Laterally invert trajectory tilt angle counterpart terminal longitude and latitude； The increment corresponding process constraint of initial trajectory inclination angle, it is separate between above-mentioned four.

(5) present invention determines trajectory using glide trajectories analytic solutions rule, and integrates and be used only for terminal velocity correction so that Trajectory planning has very fast speed.

Description of the drawings

Fig. 1 is steady glide trajectories planning modeling flow chart；

Fig. 2 is great circle coordinate system；

Fig. 3 isWith V_fRelation；

Fig. 4 is impacts of the Δ γ to process constraints；

Fig. 5 is spheric flying distance estimations；

Fig. 6 is k_εWith s_fRelation；

Fig. 7 is the quick planing method flow process of gliding section trajectory；

Fig. 8 is Longitudinal Trajectory analytic solutions precision test；

Fig. 9 is horizontal trajectory analytic solutions precision test；

Figure 10 is Bell analytic solutions；

Figure 11 is range coverage border trajectory angle of attack family of curves；

Figure 12 is range coverage border trajectory angle of heel family of curves；

Figure 13 is the horizontal ballistic curve race in range coverage border；

Figure 14 is range coverage border Longitudinal Trajectory family of curves；

Figure 15 is range coverage border ballistic velocity family of curves；

Figure 16 is range coverage border trajectory heat flow density family of curves；

Figure 17 is range coverage border trajectory dynamic pressure family of curves；

Figure 18 is range coverage border trajectory overload curves race；

Figure 19 is gliding section trajectory planning range coverage；

Figure 20 is Floorplanning trajectory；

Figure 21 is planning trajectory angle of attack curve；

Figure 22 is planning trajectory tilt angular curve.

Specific embodiment

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

The present invention is on the basis of the coupling condition of the analysis equation of motion, it is proposed that a kind of fixing longitudinal maneuver acceleration Proportionality coefficient and the trajectory pattern that crossrange maneuvering acceleration proportionality coefficient is control variable.Under the trajectory pattern, speed and other Decouple between state variable, so as to obtain the ballistic-shaped differential equation of decoupling, and obtained using the differential equation High accuracy three-dimensional ballistic-shaped analytic solutions.On the basis of this, the quick planing method of gliding section trajectory is constructed：First with bullet Road shape analytic solutions cook up the trajectory for meeting position constraint；Then terminal velocity is met about by adjusting crossrange maneuvering size Beam；Meet process constraints requirement finally by adjustment initial trajectory inclination angle.

The present invention is a kind of based on the quick trajectory planning method for steadily sliding trajectory analytic solutions, planing method modeling procedure As shown in figure 1, planning process is as shown in fig. 7, comprises following step：

Step 1：The problem modeling of gliding section trajectory planning

For the ease of theory analysis, it is assumed that the earth is spherosome, the impact of earth rotation is not considered, then half speed coordinate system Under Three Degree Of Freedom equation of particle motion it is as follows,

In formula, h is the height of aircraft, and r is the distance from the earth's core to aircraft barycenter, and relation between the two is h= r-R₀, wherein R₀For earth radius；θ, φ are respectively longitude and latitude；S is flying distance；V is the speed of the relative earth；γ is Local trajectory tilt angle；ψ is trajectory deflection angle.WithRespectively height, flying distance, longitude, latitude The derivative of degree, trajectory tilt angle, trajectory deflection angle and velocity versus time.σ is angle of heel；G=μ/r²For local gravitational acceleration, μ is The related constant of terrestrial gravitation acceleration；L and D are respectively lift acceleration and drag acceleration, and its expression formula is,

L=ρ V²S_refC_L/(2m) (2)

D=ρ V²S_refC_D/(2m) (3)

In formula, m is vehicle mass；S_refFor pneumatic area of reference；C_LAnd C_DRespectively lift coefficient and resistance coefficient, with Mach number Ma is related to angle of attack；ρ is atmospheric density, generally adopts exponential atmosphere model, as follows,

ρ=ρ_seae^-βh (4)

In formula, ρ_seaFor sea-level atmosphere density；β is exponential atmosphere model constants, generally takes β=1/7200m^-1。

Also need to consider the process constraints such as heat flow density, dynamic pressure, overload in the trajectory planning of gliding section, it is as follows

In formula,For heat flow density,For maximum heat flow density, k is constant coefficient；Q is dynamic pressure, q_maxFor max-Q；n For total overload, n_maxFor maximum total overload, g₀=9.81m/s²。

Step 2：Gliding section trajectory planning Variational Design

Introduce longitudinal maneuver acceleration a_εWith crossrange maneuvering acceleration a_β, expression formula is as follows,

a_ε=Lcos σ+(V²/r-g)cosγ (6)

a_β=Lsin σ/cos γ+(V²/r)cosγsinψtanφ (7)

In formula, a_εAnd a_βRespectively longitudinal maneuver acceleration and crossrange maneuvering acceleration.Formula (6) and formula (7) are brought into respectively Formula (1), then the trajectory tilt angle differential equation and the trajectory deflection angle differential equation can turn to,

Assume a_ε＜ 0, then the γ monotone decreasings of section of gliding, therefore independent variable in the equation of motion can be replaced with γ.By formula (1) divided by the trajectory tilt angle differential equation in formula (8), obtain

In formula,WithRespectively height, flying distance, trajectory deflection angle, longitude, The differential of latitude and speed to γ.Wherein, the first behavior gliding section lengthwise movement differential equation, the second behavior gliding section is laterally transported The dynamic differential equation, the third line is velocity differentials equation.

By formula (9) if can be seen that a_εAnd a_βAnd V²Be directly proportional, then can be in the lengthwise movement differential equation and transverse movement Do not show in the differential equation and contain speed term；Simultaneously because the differential equation of θ and φ is related to sin ψ and cos ψ respectively, to carry out Analytical Integration, then the expression formula of ψ is more simple better, therefore a_εAnd a_βBetween be also directly proportional.In addition, by formula (7) as can be seen that being Make the angle of heel amplitude of variation of planning trajectory less, a_βAlso need to be directly proportional to ρ.Generally speaking, a of gliding section can be set_εWith a_βShape Formula is as follows,

In formula, k_εAnd k_βRespectively longitudinal maneuver acceleration proportionality coefficient and crossrange maneuvering acceleration proportionality coefficient.By formula (10) bring formula (9) into, and take sin γ=γ, cos γ=1 (in gliding section γ in a small amount), then have

From formula (11) to formula (12), as long as given k_εAnd k_β, so that it may Analytical Solution is carried out to the above-mentioned differential equation.Cause This is with k_εAnd k_βAs planning variable, then gliding section trajectory planning problem can be converted into the problem for solving nonlinear equation.

Step 3：Gliding section trajectory analytic solutions are solved；

If the initial trajectory inclination angle of gliding section is γ₀, elemental height is h₀, initial trajectory drift angle is ψ₀, initial longitude and latitude Spend for respectively θ₀And φ₀；The terminal of gliding section is highly h_f, for given k_εAnd k_β, then formula (11) to (15) formula can be entered Row Analytical Solution.

Firstth, height analytic solutions

Bringing formula (4) into formula (11) can obtain,

Formula (16) integration can be obtained,

In formula, γ₀For initial trajectory inclination angle；ρ₀For section initial atmosphere density of gliding.Formula (17) is further arranged can be obtained,

In formula, C₁For a constant, value is

Secondth, flying distance analytic solutions

Bringing formula (18) into formula (12) can obtain,

If initial flight is apart from s₀=0, then formula (19) integration can obtain,

S=f_s(γ)-f_s(γ₀) (20)

In formula, f_s(γ) it is by the function of formula (19) indefinite integral acquisition, its expression formula and C₁Positive negative correlation, following institute Show,

In formula, x is independent variable.

3rd, road drift angle analytic solutions

Further, can be obtained by formula (13) integration,

ψ=ψ₀+(k_β/k_ε)(γ-γ₀) (21)

In formula, ψ₀For initial trajectory drift angle.When the section trajectory that glides is located at equator, the value of ψ is near pi/2；If Trajectory not under the line near, then can set up the great circle coordinate system that gliding section beginning and end is determined, then under great circle coordinate system Trajectory deflection angle ψ is also near pi/2.(as shown in Fig. 2 zero is located at the earth's core, x_e1Axle points to gliding starting point, y by the earth's core_e1Position In the great circle determined by the earth's core, gliding starting point, gliding terminal, and 90deg, z are less than with the angle in initial velocity direction_e1By the right side Hand rule is determined.)

4th, latitude analytic solutions

By cos ψ near pi/2 5 rank Taylor expansion, can obtain,

Take y=pi/2-ψ, wherein y is integration intermediate variable, and bring formula (18), formula (21) and formula (22) into formula (15) can ,

In formula, C₂For a constant, C is met₂=γ₀+k_ε(π/2-ψ₀)/k_β.Take,

In formula, C_a0、C_a1And C_a2It is constant.Bringing formula (24) into formula (23) can obtain,

Further take,

In formula, C_b1、C_b2、C_b3And C_b4It is constant.Formula (25) can be turned to using formula (26),

Formula (27) integration can be obtained,

φ=φ₀+f_φ(y)-f_φ(y₀) (28)

In formula, φ₀For initial latitude；y₀Related to initial trajectory drift angle, value is y₀=pi/2-ψ₀；f_φY () is and y phases The function of pass, it is as follows

In formula, g_φY () is the function related to y, as follows

5th, longitude analytic solutions

Formula (18), formula (21) and formula (29) sets forth height, trajectory deflection angle and the latitude of gliding section trajectory with trajectory The analytic solutions of change of pitch angle rule, they are substituted into respectively formula (14) can obtain,

As can be seen that the right of formula (30) is only related to γ, but complexity is higher, it is difficult to Analytical Solution, therefore available Guass-Legendre quadrature formulas directly obtain result, and expression is as follows,

In formula,γ_iFor [γ, γ₀] in i-th Gaussian node；A_iFor quadrature coefficient；N is Gaussian node Sum.

6th, analytic solutions precision analysis

In superincumbent analytic solutions, trajectory deflection angle analytic solutions (formula 21) are without any it is assumed that belonging to Exact Solutions；Height analytic solutions The source of error of (formula 18) launches in sin γ first order Taylors, usual γ ＜ 2deg, therefore the relative error of height analytic solutions is about For 10^-4；Flying distance analytic solutions (formula 20) are simultaneously not introduced into new error term, thus relative error with it is highly identical；Latitude is parsed The five rank Taylor expansions of cos ψ are introduced in solution (formula 29), if the excursion of ψ is [0, π], is then launched error and is to the maximum 0.0047, therefore the relative error of latitude analytic solutions is between 0.1% to 1%；And used latitude solution when solving longitude analytic solutions Analysis solution, and also there is quadrature error, therefore its relative error is maximum, 1% or so.

Step 4：Gliding section endgame speed-control scheme；

In step 3, given k is given_εAnd k_βUnder the conditions of gliding section trajectory analytic solutions, but do not provide the solution of speed Analysis solution.In order to obtain the size of terminal velocity, need to carry out numerical integration to speed using above-mentioned analytic solutions.

If negative specific energy e=μ/r-V²/ 2, then can be obtained by formula (1),

In formula,To bear the derivative of specific energy.Can be obtained by formula (8), formula (10) and formula (32),

By formula (33) as can be seen that C_DIt is the key of speed solution.By in the case of, C_DCan be write as C_LWith the function of e, such as Shown in lower,

C_D=f_CD(C_L,e) (34)

In formula, f_CDFor lift coefficient and the relation function of resistance coefficient.The key of solution formula (33) integration is solution C_L, by Formula (6), (7) formula and formula (10) can be obtained,

In addition, can be obtained by formula (4) and formula (18),

C can be tried to achieve by (35), formula (36) and formula (37)_L, by C_LBringing formula (34) into can obtain C_DSize.

Finally, by r and C_DExpression formula bring formula (33) into, and be integrated using Runge-Kutta method, terminal can be obtained and born Specific energy, so as to solve the velocity magnitude of terminal.It is pointed out that being apparent from C by formula (33)_DMore big then terminal velocity is less, And C_LMore big then C_DIt is bigger, therefore adjust C_LCan be with control terminal speed.Due to k_εUsually one in a small amount, also determines range Size, therefore cannot be used for speed regulation；Thus k can only be adopted_βTo adjust terminal velocity.Fig. 3 givesWith V_fRelation, It can be seen that linearisation degree between the two is very high, be conducive to numerical solution during trajectory planning.

When the trajectory of planning is comprising tilt reversal point, then the k before and after reversal point_βNeed to meet following condition of contact,

In formula, k_β-And k_β+Crossrange maneuvering acceleration proportionality coefficient respectively before and after reversal point；ρ_c、ψ_cAnd φ_cRespectively γ_cThe atmospheric density at place, trajectory deflection angle and latitude.

Step 5：Gliding section trajectory reentry corridor Adjusted Option；

Using position of the initial trajectory Inclination maneuver reentry trajectory in corridor.As shown in Figure 4, initial trajectory tilt angle gamma₀ Bigger, then the height of reentry trajectory is higher, and distance is constrained by maximum heat flow density, max-Q constraint and maximum overload constraint are true Fixed reentry corridor lower boundary is more remote.But in view of the flatness of the gliding section angle of attack, initial trajectory inclination angle value is as follows,

γ₀=γ^*+Δγ (39)

In formula, Δ γ is the initial trajectory inclination angle component for adjusting trajectory reentry corridor position；γ * are to meet angle of attack flatness The initial trajectory inclination angle component of requirement, meets,

γ^*=-2g/ (K_NβV₀ ²) (40)

K in formula_N=Lcos σ/D, are longitudinal lift-drag ratio.It is pointed out that due to the gliding section trajectory of present invention planning Meet steady gliding condition, thus heat flow density, dynamic pressure and overload change are more gentle, peak value is also less, therefore Δ γ is usual It is zero.

Step 6：Gliding section trajectory planning forming initial fields；

Before gliding section trajectory planning is carried out, need to s_f、γ^*、k_εAnd k_βInitial value estimated.

Firstth, s_fEstimation

Under spheric coordinate systems, (θ₀,φ₀) and (θ_f,φ_f) distance be (as shown in Figure 5),

In formula, s_sFor (θ₀,φ₀) to (θ_f,φ_f) spherical distance；WithMeet respectively,

Deviate the situation in face of penetrating in view of initial trajectory drift angle, terminal flying distance is,

In formula,WhereinFor section starting point meridian and the angle for penetrating face place great circle of gliding.

Secondth, γ^*Estimation

In gliding section, there is following relation in longitudinal lift-drag ratio with terminal flying distance,

Bring formula (42) into formula (40), can be in the hope of γ^*Valuation be,

3rd, k_εEstimation

If Δ γ=0, then γ₀=γ^*.From formula (17) and formula (20), γ_f、k_εAnd s_fMeet following relation,

In formula, γ_fAnd ρ_fRespectively terminal trajectory tilt angle and terminal atmospheric density.From formula (44), s is given_f, then deposit In unique k_εCorrespond to (as shown in Figure 6) therewith.The s that formula (41) is obtained_fBring into by formula (44) is solved and obtain k_εValuation.

4th, k_βEstimation

As shown in Figure 2, in k_βWhen less,With V_fLinear relationship it is preferable.Therefore in estimation k_βPrinciple during initial value is peaceful It is little not big, it is iterated amendment after being so conducive to.Known by formula (21), in known (γ-γ₀)/k_εUnder conditions of, k_βOnly with ψ_f-ψ₀It is related.Assume ψ_f-ψ₀It is only used for correcting ψ₁₀Impact, then according to plane circular arc trajectory assume have,

ψ_f-ψ₀=2 ψ₁₀

So as to k can be obtained_βInitial valuation be,

k_β=2 ψ₁₀C_N1/(γ-γ₀) (45)

Step 7：Gliding section trajectory planning flow scheme design；

Gliding section trajectory planning task is the suitable k of searching_ε、k_β, Δ γ and angle of heel reversal point γ_cSo that ablated configuration Device from gliding starting point (θ₀,φ₀,r₀,V₀,ψ₀) extremely gliding terminal (θ_f,φ_f,r_f,V_f).Aforementioned four planning variable correspond to respectively Different constraint requirements, hardly couple each other, wherein：k_εFor the flying distance of adjustment gliding section；k_βFor adjusting cunning The speed of Xiang segment endpoint；γ_cFor the final position of adjustment gliding section；Δ γ will be used to adjust reentry trajectory height, to meet Process constraints are required.Using gliding section ballistic-shaped analytic solutions, can will gliding section trajectory planning be split into ballistic-shaped planning and Terminal velocity corrects two parts.Wherein, ballistic-shaped planning need not carry out trajectory integration, solve terminal location restricted problem； The velocity correction of terminal end then needs to be integrated velocity differentials equation, solves the problems, such as that terminal velocity is constrained and process constraints.It is complete Whole gliding section trajectory flow process is as shown in fig. 7, be described as follows：

Firstth, according to step 6, s is obtained_f、γ^*、k_εAnd k_βPlanning initial value, and assume Δ γ=0 and γ_c=γ₀。

Secondth, according to given Δ γ, γ^*、γ_c、k_εAnd k_β, using step 3 θ (γ are calculated_f) and φ (γ_f).Wherein, γ_cFor reversal point trajectory tilt angle；θ(γ_f) and φ (γ_f) it is planning endgame longitude and latitude.

3rd, γ is corrected_cSo that terminal latitude meets and requires, the method for amendment is Newton iteration method, the termination bar of iteration Part is | φ (γ_f)-φ_f| ＜ φ_limit.φ in formula_limitPrecision is planned for latitude, it is considered to the parsing ballistic computation essence in step 3 Degree, can use φ_limit=0.01deg.If | φ (γ_f)-φ_f|≥φ_limitThen turn to the second step in step 7.

4th, s is corrected_f(or k_ε) terminal longitude is met require, modification method is Newton iteration method, iteration ends Condition is | θ (γ_f)-θ_f| ＜ θ_limit.θ in formula_limitPrecision is planned for longitude, it is contemplated that the parsing ballistic computation essence in step 3 Degree, can use θ_limit=0.1deg.If | θ (γ_f)-θ_f|≥θ_limitThen turn to the second step in step 7.

5th, according to Δ γ, γ^*、γ_c、k_εAnd k_β, using step 4 integration the section endgame speed V (γ that glides is obtained_f), And using Newton iteration method amendment k_βSo that | V (γ_f)-V_f| ＜ V_limit.V in formula_limitFor speed planning precision, it is contemplated that step Rapid 3 parsing ballistic computation precision, can use V_limit=10m/s.If | V (γ_f)-V_f|≥V_limitThen turn to the second step in step 7.

6th, judged whether to meet process constraints according to the integration trajectory of previous step, and Δ γ is corrected as follows,

In formula, h is practical flight height；h_limitIt is by the constraint of maximum heat flow density, max-Q constraint and maximum overload The reentry corridor height lower boundary that constraint determines；min(h-h_limit) ＜ 0 represents that gliding section trajectory exceeds reentry corridor, need to increase Big Δ γ, and turn to the second step in step 7；min(h-h_limit) >=0 represents that gliding section trajectory meets reentry corridor requirement； γ_limFor min (h-h_limit) trajectory tilt angle when taking extreme value；Δγ^(old)With Δ γ^(new)It is before respectively correcting and revised first Beginning trajectory tilt angle correction；h(γ_lim)、s(γ_lim) and h_limit(γ_lim) be respectively trajectory tilt angle and take γ_limWhen height, flight Distance and height lower boundary.

Case study on implementation：

In order to check the precision of Analytical Solution algorithm, from CAV as computation model, numerical simulation effect is carried out.Emulation Using without slewing circle spherical model and exponential atmosphere model, parameter setting such as form 1.Simulation computer CPU is Core i3-2120, Internal memory 4GB, simulated environment is Matlab.

The simulation parameter of form 1 is arranged

First, analytic solutions precision effect

Before trajectory planning is carried out, ballistic-shaped analytic solutions precision is verified first, and carried out with Bell analytic solutions Contrast.Fig. 8 and Fig. 9 sets forth the contrast of Longitudinal Trajectory analytic solutions and horizontal trajectory analytic solutions and trajectory integral result, can To find out, ballistic-shaped analytic solutions are almost completely superposed with trajectory integral result, with high precision.Specifically, adopt Analytic solutions directly predict gliding segment endpoint position, and its height error is less than 1m；Trajectory deflection angle error is less than 0.01deg；Great circle is sat Latitude error under mark system is less than 0.1deg, and longitude error is less than 0.3deg；Demonstrate the error analyses in text.Figure 10 gives The horizontal trajectory analytic solutions that Bell is given, it can be seen that Bell analytic solutions only have higher precision in less scope, when Flying distance farther out when, horizontal journey deviation is larger.This is due to earth curvature is have ignored in Bell analytic solutions to horizontal trajectory Affect, and solved using plane coordinate system, bring larger solution error.

Second, the analysis of trajectory planning range coverage

After the precision for demonstrating ballistic-shaped analytic solutions, the range coverage for analyzing trajectory planning is also needed.Can solving During up to region, it is assumed that trajectory deflection angle is fixed value, thus terminal location is mainly adjusted by tilting reversion, therefore without reversion Trajectory is the border of trajectory planning range coverage.Figure 11 to Figure 18 gives border ballistic curve race.Can by Figure 11 and Figure 12 Know, plan the angle of attack curve and tilt angular curve unusual light of trajectory and amplitude of variation is less；Glide as shown in Figure 13 section trajectory The vertical journey coverage of the range coverage of planning is about 11000 kms, and horizontal journey coverage is about 6700 kms；By Figure 14 and Figure 15 Understand that the terminal height and terminal velocity of planning trajectory meet to require；Understand that planning trajectory meets process constraints by Figure 16 to Figure 18 Require.When target is in range coverage, terminal longitude and latitude is adjusted due to only increasing angle of heel reversal point, because This its angle of attack curve, tilt curve, Longitudinal Trajectory curve and process constraints curve trend by with Figure 11 to Figure 18 in be given As a result it is similar to.

On the basis of Figure 13, Figure 19 further obtains the range coverage that gliding section plans trajectory.With steady gliding most Big range coverage is compared, smaller, the maximum range loss about 81.7km of the range coverage that the present invention is obtained, minimum range Loss is about 302.8km, and maximum horizontal journey loss is about 923.9km.This is because the present invention is from fixed k_βAs horizontal machine Dynamic planning variable so that angle of heel has larger change (as shown in figure 11) in the case of high maneuver, and this can increase horizontal journey machine Dynamic energy loss, so as to reduce crossrange maneuvering ability.When range is less or horizontal journey is larger, larger horizontal machine is needed It is dynamic, thus there is the loss of the range coverage shown in Figure 18.

3rd, gliding section trajectory is quickly planned

The initial value and final value of trajectory planning as previously shown, to Figure 19 given by the gliding section trajectory planning flow process be given using Fig. 5 The impact point in range coverage for going out carries out trajectory planning, and the result for obtaining is as shown in form 1 and Figure 20 to Figure 22.By form 2 As can be seen that can be planned within the time of 0.03s or so based on the gliding section trajectory planning method of feedback linearization analytic solutions Go out the reentry trajectory that a meet the constraint is required, and endgame has higher positional precision.Deficiency is planning trajectory Terminal velocity error it is larger, this be due to terminal velocity be according to geometry trajectory, using dynamic inversion integration obtain, make The little deviation obtained in ballistic-shaped is amplified by the integration of solving speed.Figure 20 gives the shape of planning trajectory；Figure 21 angle of attack curves for having provided out planning trajectory, it can be seen that angle of attack variation is relatively smooth, amplitude of variation is less than 4deg；Figure 22 gives Go out tilt angular curve, it is equally also relatively smooth.

The gliding section trajectory planning result of form 2

Claims (1)

1. a kind of based on the quick trajectory planning method for steadily sliding trajectory analytic solutions, including following step：

Step 1：The problem modeling of gliding section trajectory planning

If the earth is spherosome, earth rotation is not considered, then the Three Degree Of Freedom equation of particle motion under half speed coordinate system is as follows It is shown,

$\begin{array}{c}\begin{array}{cc}\stackrel{\·}{h}=V\sin \mathrm{\γ}& \stackrel{\·}{\mathrm{\θ}}=V\cos \mathrm{\γ}\sin \mathrm{\ψ}/\left(r\cos \mathrm{\φ}\right)\\ \stackrel{\·}{s}=V{\mathrm{cos\γR}}_0/r& \stackrel{\·}{\mathrm{\φ}}=V\cos \mathrm{\γ}\cos \mathrm{\ψ}/r\\ \stackrel{\·}{\mathrm{\γ}}=\[L\cos \mathrm{\σ}+(V^2/r-g)\cos \mathrm{\γ}\]/V& \stackrel{\·}{\mathrm{\ψ}}=\[L\sin \mathrm{\σ}/\cos \mathrm{\γ}+(V^2/r)\cos \mathrm{\γ}\sin \mathrm{\ψ}\tan \mathrm{\φ}\]/V\end{array}\\ \stackrel{\·}{V}=-D-g\sin \mathrm{\γ}\end{array}---\left(1\right)$

In formula, h for aircraft height, r is the distance from the earth's core to aircraft barycenter, h=r-R₀, wherein R₀For earth radius； θ, φ are respectively longitude and latitude；S is flying distance；V is the speed of the relative earth；γ is local trajectory tilt angle；ψ is that trajectory is inclined Angle；WithRespectively height, flying distance, longitude, latitude, trajectory tilt angle, trajectory deflection angle and speed Derivative of the degree to the time；σ is angle of heel；G=μ/r²For local gravitational acceleration, μ is the related constant of terrestrial gravitation acceleration； L and D are respectively lift acceleration and drag acceleration, and its expression formula is,

L=ρ V²S_refC_L/(2m) (2)

D=ρ V²S_refC_D/(2m) (3)

In formula, m is vehicle mass；S_refFor pneumatic area of reference；C_LAnd C_DRespectively lift coefficient and resistance coefficient, with Mach Ma is related to angle of attack for number；ρ is atmospheric density, as follows,

ρ=ρ_seae^-βh (4)

In formula, ρ_seaFor sea-level atmosphere density；β is exponential atmosphere model constants；

Process constraints are as follows

$\begin{array}{ccc}\stackrel{\·}{Q}=k\sqrt{\mathrm{\ρ}}{V}^{3.15}\≤{\stackrel{\·}{Q}}_{max}& q=\frac{1}{2}{\mathrm{\ρV}}^2\≤q_{\max }& n=\sqrt{L^2+D^2}/g_0\≤n_{\max }\end{array}---\left(5\right)$

In formula,For heat flow density,For maximum heat flow density, k is constant coefficient；Q is dynamic pressure, q_maxFor max-Q；N is total Overload, n_maxFor maximum total overload, g₀=9.81m/s²；

Step 2：Gliding section trajectory planning Variational Design

Introduce longitudinal maneuver acceleration a_εWith crossrange maneuvering acceleration a_β, expression formula is as follows,

a_ε=Lcos σ+(V²/r-g)cosγ (6)

a_β=Lsin σ/cos γ+(V²/r)cosγsinψtanφ (7)

In formula, a_εAnd a_βRespectively longitudinal maneuver acceleration and crossrange maneuvering acceleration；Bring formula (6) and formula (7) into formula respectively (1), then the trajectory tilt angle differential equation and the trajectory deflection angle differential equation are,

$\begin{array}{cc}\stackrel{\·}{\mathrm{\γ}}=a_{\mathrm{\ϵ}}/V& \stackrel{\·}{\mathrm{\ψ}}=a_{\mathrm{\β}}/V\end{array}---\left(8\right)$

Assume a_ε＜ 0, the then γ monotone decreasings of section of gliding, by independent variable in the equation of motion γ is replaced with；By formula (1) divided by formula (8) In the trajectory tilt angle differential equation, obtain

$\begin{array}{c}\begin{array}{cc}\frac{dh}{d\mathrm{\γ}}=\frac{V^2\sin \mathrm{\γ}}{a_{\mathrm{\ϵ}}}& \frac{ds}{d\mathrm{\γ}}=\frac{R_0{V}^2\cos \mathrm{\γ}}{{\mathrm{ra}}_{\mathrm{\ϵ}}}\end{array}\\ \begin{array}{ccc}\frac{d\mathrm{\ψ}}{d\mathrm{\γ}}=\frac{a_{\mathrm{\β}}}{a_{\mathrm{\ϵ}}}& \frac{d\mathrm{\θ}}{d\mathrm{\γ}}=\frac{V^2\cos \mathrm{\γ}\sin \mathrm{\ψ}}{{\mathrm{ra}}_{\mathrm{\ϵ}}\cos \mathrm{\φ}}& \frac{d\mathrm{\φ}}{d\mathrm{\γ}}=\frac{V^2\cos \mathrm{\γ}\cos \mathrm{\ψ}}{{\mathrm{ra}}_{\mathrm{\ϵ}}}\end{array}\\ \frac{dV}{d\mathrm{\γ}}=-\frac{(D+g\sin \mathrm{\γ})V}{a_{\mathrm{\ϵ}}}\end{array}---\left(9\right)$

In formula,WithRespectively height, flying distance, trajectory deflection angle, longitude, latitude With differential of the speed to γ；Wherein, the first behavior gliding section lengthwise movement differential equation, the second behavior gliding section transverse movement is micro- Divide equation, the third line is velocity differentials equation；

If a of gliding section_εWith a_βForm is as follows,

$\begin{array}{cc}{a}_{\mathrm{\ϵ}}=\frac{{\mathrm{\ρV}}^2{S}_{ref}{k}_{\mathrm{\ϵ}}}{2m}& a_{\mathrm{\β}}=\frac{{\mathrm{\ρV}}^2{S}_{ref}{k}_{\mathrm{\β}}}{2m}\end{array}---\left(10\right)$

In formula, k_εAnd k_βRespectively longitudinal maneuver acceleration proportionality coefficient and crossrange maneuvering acceleration proportionality coefficient；By formula (10) band Enter formula (9), and take sin γ=γ, cos γ=1, then have

$\frac{dh}{d\mathrm{\γ}}=\frac{2m}{{\mathrm{\ρS}}_{ref}{k}_{\mathrm{\ϵ}}}\mathrm{\γ}---\left(11\right)$

$\frac{ds}{d\mathrm{\γ}}=\frac{2m}{{\mathrm{\ρS}}_{ref}{k}_{\mathrm{\ϵ}}}\frac{R_0}{r}---\left(12\right)$

$\frac{d\mathrm{\ψ}}{d\mathrm{\γ}}=\frac{k_{\mathrm{\β}}}{k_{\mathrm{\ϵ}}}---\left(13\right)$

$\frac{d\mathrm{\θ}}{d\mathrm{\γ}}=\frac{2m}{{\mathrm{\ρS}}_{ref}{k}_{\mathrm{\ϵ}}}\frac{sin\mathrm{\ψ}}{rcos\mathrm{\φ}}---\left(14\right)$

$\frac{d\mathrm{\φ}}{d\mathrm{\γ}}=\frac{2m}{{\mathrm{\ρS}}_{ref}{k}_{\mathrm{\ϵ}}}\frac{cos\mathrm{\ψ}}{r}---\left(15\right)$

With k_εAnd k_βAs planning variable, gliding section trajectory planning problem is converted into the problem for solving nonlinear equation；

Step 3：Gliding section trajectory analytic solutions are solved；

If the initial trajectory inclination angle of gliding section is γ₀, elemental height is h₀, initial trajectory drift angle is ψ₀, initial longitude and latitude be Respectively θ₀And φ₀；The terminal of gliding section is highly h_f, for given k_εAnd k_β, parsing is carried out to formula (11) to (15) formula and is asked Solution；

Firstth, height analytic solutions

Bringing formula (4) into formula (11) can obtain,

$e^{-\mathrm{\β}h}dh=\frac{2m}{{\mathrm{\ρ}}_{sea}{S}_{ref}{k}_{\mathrm{\ϵ}}}\mathrm{\γ}d\mathrm{\γ}---\left(16\right)$

Formula (16) integration can be obtained,

${\mathrm{\γ}}^2-{\mathrm{\γ}}_0^2=\frac{S_{ref}{k}_{\mathrm{\ϵ}}}{m\mathrm{\β}}({\mathrm{\ρ}}_0-\mathrm{\ρ})---\left(17\right)$

In formula, γ₀For initial trajectory inclination angle；ρ₀For section initial atmosphere density of gliding；Formula (17) is further arranged can be obtained,

$\frac{m}{{\mathrm{\ρS}}_{ref}{k}_{\mathrm{\ϵ}}}=\frac{1}{\mathrm{\β}(C_1-{\mathrm{\γ}}^2)}---\left(18\right)$

In formula, C₁For a constant, value is

Secondth, flying distance analytic solutions

Bringing formula (18) into formula (12) can obtain,

$\frac{ds}{d\mathrm{\γ}}=\frac{-2R_0}{\mathrm{\β}r({\mathrm{\γ}}^2-C_1)}---\left(19\right)$

If initial flight is apart from s₀=0, then formula (19) integration can obtain,

S=f_s(γ)-f_s(γ₀) (20)

In formula, f_s(γ) it is by the function of formula (19) indefinite integral acquisition, its expression formula and C₁Positive negative correlation, it is as follows,

$f_s\left(x\right)=\left\{\begin{array}{cc}-\frac{2R_0}{\mathrm{\&beta;}r\sqrt{-C_1}}arctan\left(\frac{x}{\sqrt{-C_1}}\right)& C_1<0\\ 2R_0/\left(\mathrm{\&beta;}rx\right)& C_1=0\\ -\frac{R_0}{\mathrm{\&beta;}r\sqrt{C_1}}ln\left(\frac{x-\sqrt{C_1}}{x+\sqrt{C_1}}\right)& C_1>0\end{array}\right.$

In formula, x is independent variable；

3rd, road drift angle analytic solutions

Further, can be obtained by formula (13) integration,

ψ=ψ₀+(k_β/k_ε)(γ-γ₀) (21)

In formula, ψ₀For initial trajectory drift angle；

4th, latitude analytic solutions

By cos ψ near pi/2 5 rank Taylor expansion, can obtain,

$cos\mathrm{\&psi;}=(\frac{\mathrm{\&pi;}}{2}-\mathrm{\&psi;})-\frac{1}{6}{(\frac{\mathrm{\&pi;}}{2}-\mathrm{\&psi;})}^3+\frac{1}{120}{(\frac{\mathrm{\&pi;}}{2}-\mathrm{\&psi;})}^5---\left(22\right)$

Taking y=pi/2-ψ, wherein y is integration intermediate variable, and bring formula (18), formula (21) and formula (22) into formula (15) to obtain,

$\frac{d\mathrm{\&phi;}}{dy}=\left(\frac{k_{\mathrm{\&epsiv;}}}{{\mathrm{r\&beta;k}}_{\mathrm{\&beta;}}}\right)\frac{2y-y^3/3+y^5/60}{{\&lsqb;C_2-(k_{\mathrm{\&epsiv;}}/k_{\mathrm{\&beta;}})y\&rsqb;}^2-C_1}---\left(23\right)$

In formula, C₂For a constant, C is met₂=γ₀+k_ε(π/2-ψ₀)/k_β；Take,

$\begin{array}{c}\begin{array}{cc}{C}_{a0}=k_{\mathrm{\&beta;}}/\left(60{\mathrm{r\&beta;k}}_{\mathrm{\&epsiv;}}\right)& C_{a1}=2C_2{k}_{\mathrm{\&beta;}}/k_{\mathrm{\&epsiv;}}\end{array}\\ C_{a2}=(C_2^2-C_1){(k_{\mathrm{\&beta;}}/k_{\mathrm{\&epsiv;}})}^2\end{array}---\left(24\right)$

In formula, C_a0、C_a1And C_a2It is constant；Bringing formula (24) into formula (23) can obtain,

$\frac{d\mathrm{\&phi;}}{dy}=\frac{C_{a0}(y^5-20y^3+120y)}{y^2-C_{a1}y+C_{a2}}---\left(25\right)$

Further take,

$\begin{array}{c}{C}_{b1}=C_{a1}^2-C_{a2}-20\\ C_{b2}=C_{a1}^3-2C_{a1}{C}_{a2}-2C_{a1}\\ C_{b3}=C_{a1}^4-3C_{a1}^2{C}_{a2}-20C_{a1}^2+C_{a2}^2+20C_{a2}+120\\ C_{b4}=-C_{a1}^3{C}_{a2}+2C_{a1}{C}_{a2}^2+20C_{a1}{C}_{a2}\end{array}---\left(26\right)$

In formula, C_b1、C_b2、C_b3And C_b4It is constant；Formula (25) can be turned to using formula (26),

$\begin{array}{c}\frac{d\mathrm{\&phi;}}{dy}=C_{a0}(y^3+C_{a1}{y}^{2}y+C_{b1}y+C_{b2})\\ +\frac{(C_{a0}{C}_{b3}/2)(2y-C_{a1})}{y^2-C_{a1}y+C_{a2}}+\frac{C_{a0}(C_{b3}{C}_{a1}/2+C_{b4})}{y^2-C_{a1}y+C_{a2}}\end{array}---\left(27\right)$

Formula (27) integration can be obtained,

φ=φ₀+f_φ(y)-f_φ(y₀) (28)

In formula, φ₀For initial latitude；y₀Related to initial trajectory drift angle, value is y₀=pi/2-ψ₀；f_φY () is the letter related to y Number, it is as follows

$\begin{array}{c}{f}_{\mathrm{\&phi;}}\left(y\right)=C_{a0}(\frac{y^4}{4}+\frac{C_{a1}{y}^3}{3}+\frac{C_{b1}{y}^2}{2}+C_{b2}y)\\ +\frac{C_{a0}{C}_{b3}}{2}\ln (y^2-C_{a1}y+C_{a2})\\ +C_{a0}(C_{b3}{C}_{a1}/2+C_{b4})g_{\mathrm{\&phi;}}\left(y\right)\end{array}---\left(29\right)$

In formula, g_φY () is the function related to y, as follows

$g_{\mathrm{\&phi;}}\left(y\right)=\left\{\begin{array}{cc}\frac{2}{\sqrt{\left|C_{a1}^2-4C_{a2}\right|}}arctan\left(\frac{2y-C_{a1}}{\sqrt{\left|C_{a1}^2-4C_{a2}\right|}}\right)& C_{a1}^2<4C_{a2}\\ -\frac{2}{2y-C_{a1}}& C_{a1}^2=4C_{a2}\\ \frac{1}{\sqrt{\left|C_{a1}^2-4C_{a2}\right|}}ln\left(\frac{2y-C_{a1}-\sqrt{\left|C_{a1}^2-4C_{a2}\right|}}{2y-C_{a1}+\sqrt{\left|C_{a1}^2-4C_{a2}\right|}}\right)& C_{a1}^2>4C_{a2}\end{array}\right.$

5th, longitude analytic solutions

Formula (18), formula (21) and formula (29) sets forth height, trajectory deflection angle and the latitude of gliding section trajectory with trajectory tilt angle The analytic solutions of Changing Pattern, they are substituted into respectively formula (14) can obtain,

$\frac{d\mathrm{\&theta;}}{d\mathrm{\&gamma;}}=\frac{2sin\&lsqb;{\mathrm{\&psi;}}_0+(k_{\mathrm{\&beta;}}/k_{\mathrm{\&epsiv;}})(\mathrm{\&gamma;}-{\mathrm{\&gamma;}}_0)\&rsqb;}{\mathrm{\&beta;}r(C_1-{\mathrm{\&gamma;}}^2)cos\&lsqb;{\mathrm{\&phi;}}_0+f_{\mathrm{\&phi;}}\left(x\right)-f_{\mathrm{\&phi;}}\left(x_0\right)\&rsqb;}---\left(30\right)$

Result is directly obtained using Guass-Legendre quadrature formulas, expression is as follows,

$\mathrm{\&theta;}={\mathrm{\&theta;}}_0+\frac{\mathrm{\&gamma;}-{\mathrm{\&gamma;}}_0}{2}\underset{i=1}{\overset{n}{\&Sigma;}}{A}_i{\stackrel{\&CenterDot;}{\mathrm{\&theta;}}}_1\left({\mathrm{\&gamma;}}_i\right)---\left(31\right)$

In formula,γ_iFor i-th Gaussian node in [γ, γ 0]；A_iFor quadrature coefficient；N is that Gaussian node is total Number；

Step 4：Gliding section endgame speed-control scheme；

If negative specific energy e=μ/r-V²/ 2, then can be obtained by formula (1),

$\stackrel{\&CenterDot;}{e}=VD---\left(32\right)$

In formula,To bear the derivative of specific energy；Can be obtained by formula (8), formula (10) and formula (32),

$\frac{de}{d\mathrm{\&gamma;}}=\frac{2(\mathrm{\&mu;}/r-e)C_D}{k_{\mathrm{\&epsiv;}}}---\left(33\right)$

Found out by formula (33), C_DIt is the key of speed solution；By in the case of, C_DCan be write as C_LWith the function of e, it is as follows,

C_D=f_CD(C_L,e) (34)

In formula, f_CDFor lift coefficient and the relation function of resistance coefficient；The key of solution formula (33) integration is solution C_L, by formula (6), (7) formula and formula (10) can be obtained,

$C_{L}cos\mathrm{\&sigma;}=\frac{m\&lsqb;g/(\mathrm{\&mu;}/r-e)-2/r\&rsqb;cos\mathrm{\&gamma;}}{{\mathrm{\&rho;S}}_{ref}}+k_{\mathrm{\&epsiv;}}---\left(35\right)$

$C_{L}sin\mathrm{\&sigma;}=k_{\mathrm{\&beta;}}-\frac{2m{\left(cos\mathrm{\&gamma;}\right)}^{2}sin\mathrm{\&psi;}tan\mathrm{\&phi;}}{{\mathrm{\&rho;rS}}_{ref}}---\left(36\right)$

In addition, can be obtained by formula (4) and formula (18),

$r=R_0-\frac{1}{\mathrm{\&beta;}}ln\&lsqb;\frac{m\mathrm{\&beta;}(C_1-{\mathrm{\&gamma;}}^2)}{{\mathrm{\&rho;}}_0{S}_{ref}{k}_{\mathrm{\&epsiv;}}}\&rsqb;---\left(37\right)$

C can be tried to achieve by (35), formula (36) and formula (37)_L, by C_LBringing formula (34) into can obtain C_DSize；

Finally, by r and C_DExpression formula bring formula (33) into, and be integrated using Runge-Kutta method, obtain the negative specific energy of terminal, So as to solve the velocity magnitude of terminal；

When the trajectory of planning is comprising tilt reversal point, then the k before and after reversal point_βNeed to meet following condition of contact,

$k_{\mathrm{\&beta;}-}+k_{\mathrm{\&beta;}+}=\frac{4m{\left({\mathrm{cos\&gamma;}}_c\right)}^2{\mathrm{sin\&psi;}}_c{\mathrm{tan\&phi;}}_c}{{\mathrm{\&rho;}}_c{\mathrm{rS}}_{ref}}---\left(38\right)$

In formula, k_β-And k_β+Crossrange maneuvering acceleration proportionality coefficient respectively before and after reversal point；ρ_c、ψ_cAnd φ_cRespectively γ_cPlace Atmospheric density, trajectory deflection angle and latitude；

Step 5：Gliding section trajectory reentry corridor Adjusted Option；

Initial trajectory inclination angle value is as follows,

γ₀=γ^*+Δγ (39)

In formula, Δ γ is the initial trajectory inclination angle component for adjusting trajectory reentry corridor position；γ^*Require to meet angle of attack flatness Initial trajectory inclination angle component, meet,

${\mathrm{\&gamma;}}^{*}=-2g/\left(K_N{\mathrm{\&beta;V}}_0^2\right)---\left(40\right)$

K in formula_N=L cos σ/D, are longitudinal lift-drag ratio；

Step 6：Gliding section trajectory planning forming initial fields；

Before gliding section trajectory planning is carried out, need to s_f、γ^*、k_εAnd k_βInitial value estimated；

Firstth, s_fEstimation

Under spheric coordinate systems, (θ₀,φ₀) and (θ_f,φ_f) distance be,

$s_s=R_{0}arccos\&lsqb;{\left({\stackrel{\&RightArrow;}{r}}_0\right)}^T{\stackrel{\&RightArrow;}{r}}_f\&rsqb;$

In formula, s_sFor (θ₀,φ₀) to (θ_f,φ_f) spherical distance；WithMeet respectively,

$\begin{array}{cc}{\stackrel{\&RightArrow;}{r}}_0=\left[\begin{array}{c}{\mathrm{cos\&phi;}}_0{\mathrm{cos\&theta;}}_0\\ {\mathrm{cos\&phi;}}_0{\mathrm{sin\&theta;}}_0\\ {\mathrm{sin\&phi;}}_0\end{array}\right]& {\stackrel{\&RightArrow;}{r}}_f=\left[\begin{array}{c}{\mathrm{cos\&phi;}}_f{\mathrm{cos\&theta;}}_f\\ {\mathrm{cos\&phi;}}_f{\mathrm{sin\&theta;}}_f\\ {\mathrm{sin\&phi;}}_f\end{array}\right]\end{array}$

Then terminal flying distance is,

$s_f=\frac{{\mathrm{\&psi;}}_{10}{s}_s}{{\mathrm{sin\&psi;}}_{10}}---\left(41\right)$

In formula,WhereinFor section starting point meridian and the angle for penetrating face place great circle of gliding；

Secondth, γ^*Estimation

In gliding section, there is following relation in longitudinal lift-drag ratio with terminal flying distance,

$K_N=\frac{2s_f}{rln\&lsqb;(gr-V_f^2)/(gr-V_0^2)\&rsqb;}---\left(42\right)$

Bring formula (42) into formula (40), try to achieve γ^*Valuation be,

${\mathrm{\&gamma;}}^{*}=-\frac{grln\&lsqb;(gr-V_f^2)/(gr-V_0^2)\&rsqb;}{s_f{\mathrm{\&beta;V}}_0^2}---\left(43\right)$

3rd, k_εEstimation

If Δ γ=0, then γ₀=γ^*；From formula (17) and formula (20), γ_f、k_εAnd s_fMeet following relation,

$\begin{array}{c}{\mathrm{\&gamma;}}_f^2-{\mathrm{\&gamma;}}_0^2=\frac{S_{ref}{k}_{\mathrm{\&epsiv;}}}{m\mathrm{\&beta;}}({\mathrm{\&rho;}}_0-{\mathrm{\&rho;}}_f)\\ s_f=f_s\left({\mathrm{\&gamma;}}_f\right)-f_s\left({\mathrm{\&gamma;}}_0\right)\end{array}---\left(44\right)$

In formula, γ_fAnd ρ_fRespectively terminal trajectory tilt angle and terminal atmospheric density；From formula (44), s is given_f, then exist only One k_εCorrespond to therewith；The s that formula (41) is obtained_fBring into by formula (44) is solved and obtain k_εValuation；

4th, k_βEstimation

Assume ψ_f-ψ₀It is only used for correcting ψ₁₀Impact, then according to plane circular arc trajectory assume have,

ψ_f-ψ₀=2 ψ₁₀

k_βInitial valuation be,

k_β=2 ψ₁₀C_N1/(γ-γ₀) (45)

Step 7：Gliding section trajectory planning flow scheme design；

Specially：

Firstth, according to step 6, s is obtained_f、γ^*、k_εAnd k_βPlanning initial value, and assume Δ γ=0 and γ_c=γ₀；

Secondth, according to given Δ γ, γ^*、γ_c、k_εAnd k_β, using step 3 θ (γ are calculated_f) and φ (γ_f)；Wherein, γ_cFor anti- Turning point trajectory tilt angle；θ(γ_f) and φ (γ_f) it is planning endgame longitude and latitude；

3rd, γ is corrected_cSo that terminal latitude meets and requires, the method for amendment is Newton iteration method, and the end condition of iteration is | φ(γ_f)-φ_f| ＜ φ_limit；φ in formula_limitPrecision is planned for latitude；If | φ (γ_f)-φ_f|≥φ_limitThen turn to step 7 In second step；

4th, s is corrected_fOr k_εSo that terminal longitude meets requiring, modification method is Newton iteration method, and stopping criterion for iteration is | θ(γ_f)-θ_f| ＜ θ_limit；θ in formula_limitPrecision is planned for longitude；If | θ (γ_f)-θ_f|≥θ_limitThen turn to the in step 7 Two steps；

5th, according to Δ γ, γ^*、γ_c、k_εAnd k_β, using step 4 integration the section endgame speed V (γ that glides is obtained_f), and adopt With Newton iteration method amendment k_βSo thatIn formulaFor speed planning precision；If | V (γ_f)-V_f|≥ V_limitThen turn to the second step in step 7；

6th, judged whether to meet process constraints according to the integration trajectory of previous step, and Δ γ is corrected as follows,

${\mathrm{\&Delta;\&gamma;}}^{\left(new\right)}=\left\{\begin{array}{cc}{\mathrm{\&Delta;\&gamma;}}^{\left(old\right)}-\frac{h\left({\mathrm{\&gamma;}}_{\lim }\right)-h_{\lim it}\left({\mathrm{\&gamma;}}_{\lim }\right)}{s\left({\mathrm{\&gamma;}}_{\lim }\right)}& \min (h-h_{\lim it})<0\\ {\mathrm{\&Delta;\&gamma;}}^{\left(old\right)}& \min (h-h_{\lim it})\&GreaterEqual;0\end{array}\right.$

In formula, h is practical flight height；h_limitIt is by the constraint of maximum heat flow density, max-Q constraint and maximum overload constraint It is determined that reentry corridor height lower boundary；min(h-h_limit) ＜ 0 represents that gliding section trajectory exceeds reentry corridor, need to increase Δ γ, and turn to the second step in step 7；min(h-h_limit) >=0 represents that gliding section trajectory meets reentry corridor requirement；γ_limFor min(h-h_limit) trajectory tilt angle when taking extreme value；Δγ^(old)With Δ γ^(new)With revised initial trajectory before respectively correcting Inclination correction amount；h(γ_lim)、s(γ_lim) and h_limit(γ_lim) be respectively trajectory tilt angle and take γ_limWhen height, flying distance and Height lower boundary.