CN104392047A - 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 solution

Technical field

The present invention relates to a kind of quick trajectory planning method based on steadily sliding trajectory analytic solution, belonging to spationautics, weapon technologies field.

Background technology

Along with the fast development of hypersonic technology, reentry trajectory is formulated in order to study hotspot.Glide section is the main composition part of reentry trajectory, determines voyage and the maneuvering range of ablated configuration.Glide section trajectory plans the performance not only can analyzing hypersonic aircraft fast, also can be used for online trajectory planning and Predictor-corrector guidance, has higher researching value.

Due to glide section trajectory, to have the flight time long, extremely sensitive to planning variable, and reenter solve nonlinear and be about beam intensity, makes the feasible zone of glide trajectories narrower, often need the computing time grown very much when adopting traditional optimized algorithm to solve.Although the optimized algorithms such as pseudo-spectrometry, SQP are greatly improved in counting yield, it calculates consuming time still more than second level, and when process constraints touches border, its counting yield can reduce greatly, even has influence on the convergence of result.In order to improve the speed of glide section trajectory planning, some special constraints are introduced into glide section trajectory planning, as resistance curve, etc. heat flow curve, etc. dynamic pressure curve, equilibrium glide condition etc.The common feature of these constraints glide section trajectory is limited to certain special steady glide state, thus reduce the sensitivity of trajectory planning problem.But owing to also needing when trajectory planning to carry out numerical integration, limit the raising of planning speed.In order to break away from the dependence of the integration to trajectory, dynamic inverse is introduced into glide section trajectory planning.The method directly plans geometry ballistic-shaped, and meets constraint requirements by adjustment ballistic-shaped.But because the environmental change before and after glide section is violent, make glide section overall trajectory be difficult to use low order curve.

Summary of the invention

The object of the invention is by solving steady glide trajectories high precision analytic solution, obtaining the quick planing method of glide section trajectory not relying on trajectory integration, for analysis hypersonic aircraft performance and Predictor-corrector guidance provide technical support.

Based on the quick trajectory planning method steadily sliding trajectory analytic solution, comprise following step:

Step 1: glide section trajectory planning problem modeling;

Step 2: glide section trajectory planning Variational Design;

Step 3: glide section trajectory analytic solution solve;

Step 4: glide section endgame speed-control scheme;

Step 5: glide section trajectory reentry corridor Adjusted Option;

Step 6: glide section trajectory planning forming initial fields;

Step 7: glide section trajectory planning flow scheme design.

The invention has the advantages that:

(1) to propose with longitudinal maneuver acceleration scale-up factor and crossrange maneuvering acceleration scale-up factor as glide section trajectory planning variable, make trajectory tilt angle in the equation of motion, trajectory deflection angle, highly, in the differential equation of longitude and latitude not containing speed term;

(2) obtain fixing longitudinal maneuver acceleration scale-up factor and glide section trajectory analytic solution corresponding to crossrange maneuvering acceleration scale-up factor, comprise height analytic solution, flying distance analytic solution, trajectory deflection angle analytic solution, longitude analytic solution and latitude analytic solution;

(3) the forming initial fields method of longitudinal maneuver acceleration scale-up factor and crossrange maneuvering acceleration scale-up factor is given;

(4) propose with longitudinal maneuver acceleration scale-up factor, crossrange maneuvering acceleration scale-up factor, laterally reversion trajectory tilt angle and initial trajectory inclination angle increment for glide section trajectory planning variable, wherein: longitudinal maneuver acceleration scale-up factor corresponding glide section flying distance; Crossrange maneuvering acceleration scale-up factor counterpart terminal speed; Laterally reversion trajectory tilt angle counterpart terminal longitude and latitude; The increment corresponding process constraint of initial trajectory inclination angle, separate between above-mentioned four.

(5) the present invention adopts glide trajectories analytic solution to advise to determine trajectory, and integration only corrects for terminal velocity, makes trajectory planning have speed quickly.

Accompanying drawing explanation

Fig. 1 is steady glide trajectories planning modeling process flow diagram;

Fig. 2 is great circle coordinate system;

Fig. 3 is with V _frelation;

Fig. 4 is the impact of Δ γ on 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 glide section trajectory;

Fig. 8 is Longitudinal Trajectory analytic solution precision tests;

Fig. 9 is horizontal trajectory analytic solution precision test;

Figure 10 is Bell analytic solution;

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 glide section trajectory planning range coverage;

Figure 20 is Floorplanning trajectory;

Figure 21 is planning trajectory angle of attack curve;

Figure 22 is planning trajectory angle of heel curve.

Embodiment

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

The present invention is on the basis of coupling condition analyzing the equation of motion, and proposing a kind of is the trajectory pattern of control variable with fixing longitudinal maneuver acceleration scale-up factor and crossrange maneuvering acceleration ratio.Under this trajectory pattern, decoupling zero between speed and other state variables, thus the ballistic-shaped differential equation obtaining decoupling zero, and utilize this differential equation to obtain high accuracy three-dimensional ballistic-shaped analytic solution.On this basis, construct the quick planing method of glide section trajectory: first utilize ballistic-shaped analytic solution to cook up the trajectory meeting position constraint; Then terminal velocity constraint is met by adjustment crossrange maneuvering size; Process constraints requirement is met finally by adjustment initial trajectory inclination angle.

The present invention is a kind of quick trajectory planning method based on steadily sliding trajectory analytic solution, and as shown in Figure 1, planning process as shown in Figure 7, comprises following step to planing method modeling procedure:

Step 1: glide section trajectory planning problem modeling

For the ease of theoretical analysis, suppose that the earth is spherosome, do not consider the impact of earth rotation, then the Three Degree Of Freedom equation of particle motion under half speed coordinate system is as follows,

$\begin{array}{cc}\stackrel{\·}{h}=V\sin \mathrm{\γ}& \stackrel{\·}{\mathrm{\θ}}=V\cos \mathrm{\γ}\sin \mathrm{\ψ}/\left(r\cos \mathrm{\φ}\right)\end{array}$

$\begin{array}{cc}\stackrel{\·}{s}=V\cos \mathrm{\γ}{R}_0/r& \stackrel{\·}{\mathrm{\φ}}=V\cos \mathrm{\γ}\cos \mathrm{\ψ}/r\end{array}$

$\begin{array}{cc}\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}----\left(1\right)$

$\stackrel{\·}{V}=-D-g\sin \mathrm{\γ}$

In formula, h is the height of aircraft, and r is the distance from the earth's core to aircraft barycenter, and pass is between the two 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. with be respectively the derivative of height, flying distance, longitude, latitude, trajectory tilt angle, trajectory deflection angle and velocity versus time.σ is angle of heel; G=μ/r ²for local gravitational acceleration, μ is the constant that terrestrial gravitation acceleration is relevant; L and D is 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 _dbe respectively lift coefficient and resistance coefficient, relevant with angle of attack to Mach number Ma; ρ is atmospheric density, usually adopts exponential atmosphere model, as follows,

ρ＝ρ _seae ^-βh(4)

In formula, ρ _seafor sea-level atmosphere density; β is exponential atmosphere model constants, usually gets β=1/7200m ^-1.

Also need to consider the process constraints such as heat flow density, dynamic pressure, overload in the trajectory planning of glide section, 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 current 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: glide section trajectory planning Variational Design

Introduce longitudinal maneuver acceleration a _εwith crossrange maneuvering acceleration a _β, expression formula is as follows,

a _ε＝L cosσ+(V ²/r-g)cosγ (6)

a _β＝L sin σ/cosγ+(V ²/r)cosγsinψtanφ (7)

In formula, a _εand a _βbe respectively longitudinal maneuver acceleration and crossrange maneuvering acceleration.Bring formula (6) and formula (7) into formula (1) respectively, then the trajectory tilt angle differential equation and the trajectory deflection angle differential equation can turn to,

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

Suppose a _ε< 0, then the γ monotone decreasing of the section of glide, therefore can replace with γ by independent variable in the equation of motion.By formula (1) divided by the trajectory tilt angle differential equation in formula (8),

$\begin{array}{cc}\frac{\mathrm{dh}}{\mathrm{d\&gamma;}}=\frac{V^2\sin \mathrm{\&gamma;}}{a_{\mathrm{\&epsiv;}}}& \frac{\mathrm{ds}}{\mathrm{d\&gamma;}}=\frac{R_0{V}^2\cos \mathrm{\&gamma;}}{r{a}_{\mathrm{\&epsiv;}}}\end{array}$

$\begin{array}{ccc}\frac{\mathrm{d\&psi;}}{\mathrm{d\&gamma;}}=\frac{a_{\mathrm{\&beta;}}}{a_{\mathrm{\&epsiv;}}}& \frac{\mathrm{d\&theta;}}{\mathrm{d\&gamma;}}=\frac{V^2\cos \mathrm{\&gamma;}\sin \mathrm{\&psi;}}{r{a}_{\mathrm{\&epsiv;}}\cos \mathrm{\&phi;}}& \frac{\mathrm{d\&phi;}}{\mathrm{d\&gamma;}}=\frac{V^2\cos \mathrm{\&gamma;}\cos \mathrm{\&psi;}}{r{a}_{\mathrm{\&epsiv;}}}\end{array}---\left(9\right)$

$\frac{\mathrm{dV}}{\mathrm{d\&gamma;}}=-\frac{(D+g\sin \mathrm{\&gamma;})V}{a_{\mathrm{\&epsiv;}}}$

In formula, with be respectively the differential equation of height, flying distance, trajectory deflection angle, longitude, latitude and speed.Wherein, the first behavior glide section lengthwise movement differential equation, the second behavior glide section transverse movement differential equation, the third line is velocity differentials equation.

As can be seen from formula (9), if a _εand a _βall and V ²be directly proportional, then can be all not aobvious containing speed term in the lengthwise movement differential equation and the transverse movement differential equation; Simultaneously because the differential equation of θ with φ is relevant with cos ψ to sin ψ 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, as can be seen from formula (7), in order to make the angle of heel amplitude of variation of planning trajectory less, a _βalso need to be directly proportional to ρ.Generally speaking, a of glide section can be established _εwith a _βform is as follows,

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

In formula, k _εand k _βbe respectively longitudinal maneuver acceleration scale-up factor and crossrange maneuvering acceleration scale-up factor.Bring formula (10) into formula (9), and get sin γ=γ, cos γ=1 (at glide section γ in a small amount), then have

$\frac{\mathrm{dh}}{\mathrm{d\&gamma;}}=\frac{2m}{\mathrm{\&rho;}{S}_{\mathrm{ref}}{k}_{\mathrm{\&epsiv;}}}\mathrm{\&gamma;}---\left(11\right)$

$\frac{\mathrm{ds}}{\mathrm{d\&gamma;}}=\frac{2m}{\mathrm{\&rho;}{S}_{\mathrm{ref}}{k}_{\mathrm{\&epsiv;}}}\frac{R_0}{r}---\left(12\right)$

$\frac{\mathrm{d\&psi;}}{\mathrm{d\&gamma;}}=\frac{k_{\mathrm{\&beta;}}}{k_{\mathrm{\&epsiv;}}}---\left(13\right)$

$\frac{\mathrm{d\&theta;}}{\mathrm{d\&gamma;}}=\frac{2m}{\mathrm{\&rho;}{S}_{\mathrm{ref}}{k}_{\mathrm{\&epsiv;}}}\frac{\sin \mathrm{\&psi;}}{r\cos \mathrm{\&phi;}}---\left(14\right)$

$\frac{\mathrm{d\&phi;}}{\mathrm{d\&gamma;}}=\frac{2m}{\mathrm{\&rho;}{S}_{\mathrm{ref}}{k}_{\mathrm{\&epsiv;}}}\frac{\cos \mathrm{\&psi;}}{r}---\left(15\right)$

From formula (11) to formula (12), as long as given k _εand k _β, just can carry out Analytical Solution to the above-mentioned differential equation.Therefore with k _εand k _βas planning variable, then glide section trajectory planning problem can be converted into the problem solving nonlinear equation.

Step 3: glide section trajectory analytic solution solve;

If the initial trajectory inclination angle of glide section is γ ₀, elemental height is h ₀, initial trajectory drift angle is ψ ₀, initial longitude and latitude be for being respectively θ ₀and φ ₀; The terminal height of glide section is h _f, for given k _εand k _β, then Analytical Solution can be carried out to formula (11) to (15) formula.

The first, height analytic solution

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

$e^{-\mathrm{\&beta;h}}\mathrm{dh}=\frac{2m}{{\mathrm{\&rho;}}_{\mathrm{sea}}{S}_{\mathrm{ref}}{k}_{\mathrm{\&epsiv;}}}\mathrm{\&gamma;d\&gamma;}---\left(16\right)$

Can obtain formula (16) integration,

${\mathrm{\&gamma;}}^2-{\mathrm{\&gamma;}}_0^2=\frac{S_{\mathrm{ref}}{k}_{\mathrm{\&epsiv;}}}{m}({\mathrm{\&rho;}}_0-\mathrm{\&rho;})---\left(17\right)$

In formula, γ ₀for initial trajectory inclination angle; ρ ₀for glide section initial atmosphere density.Formula (17) is arranged further and can be obtained,

$\frac{m}{{\mathrm{\&rho;S}}_{\mathrm{ref}}{k}_{\mathrm{\&epsiv;}}}=\frac{1}{\mathrm{\&beta;}(C_1-{\mathrm{\&gamma;}}^2)}---\left(18\right)$

In formula, C ₁be a constant, value is $C_1=S_{\mathrm{ref}}{k}_{\mathrm{\&epsiv;}}{\mathrm{\&rho;}}_0/\left(\mathrm{m\&beta;}\right)+{\mathrm{\&gamma;}}_0^2.$

The second, flying distance analytic solution

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

$\frac{\mathrm{ds}}{\mathrm{d\&gamma;}}=\frac{-2R_0}{\mathrm{\&beta;r}({\mathrm{\&gamma;}}^2-C_1)}---\left(19\right)$

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

s＝f _s(γ)-f _s(γ ₀) (20)

In formula, f _s(γ) be the function obtained by formula (19) indefinite integral, its expression formula and C ₁positive negative correlation, 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.

Three, drift angle, road analytic solution

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

ψ＝ψ ₀+(k _β/k _ε)(γ-γ ₀) (21)

In formula, ψ ₀for initial trajectory drift angle.When glide section trajectory is positioned 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 glide section starting point and terminal determine, then the trajectory deflection angle ψ under great circle coordinate system is also near pi/2.(as shown in Figure 2, true origin is positioned at the earth's core, x _e1axle points to glide starting point by the earth's core, y _e1be positioned at the great circle determined by the earth's core, glide starting point, glide terminal, and be less than 90deg with the angle in initial velocity direction, z _e1determined by the right-hand rule.)

Four, latitude analytic solution

By cos ψ 5 rank Taylor expansions near pi/2, 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)$

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

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

In formula, C ₂be a constant, meet C ₂=γ ₀+ k _ε(pi/2-ψ ₀)/k _β.Get,

$\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 _a2be constant.Bring formula (24) into formula (23) can obtain,

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

Get further,

$\begin{array}{c}{C}_{b1}=C_{a1}^2-C_{a2}-20\\ C_{b2}=C_{a1}^3-2C_{a1}{C}_{a2}-20C_{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 _b4be constant.Utilize formula (26) formula (25) can be turned to,

$\begin{array}{c}\frac{\mathrm{d\&phi;}}{\mathrm{dy}}=C_{a0}(y^3+C_{a1}{y}^2+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)$

Can obtain formula (27) integration,

φ＝φ ₀+f _φ(y)-f _φ(y ₀) (28)

In formula, φ ₀for initial latitude; y ₀relevant to initial trajectory drift angle, value is y ₀=pi/2-ψ ₀; f _φy () is the function relevant to y, 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 relevant to y, as follows

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

Five, longitude analytic solution

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

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

Can find out, the right of formula (30) is only relevant to γ, but complexity is higher, is difficult to Analytical Solution, and Guass-Legendre quadrature formula therefore can be utilized directly to obtain result, and expression is as follows,

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

In formula, γ _ifor [γ, γ ₀] in i-th Gaussian node; A _ifor quadrature coefficient; N is Gaussian node sum.

Six, analytic solution precision analysis

In superincumbent analytic solution, trajectory deflection angle analytic solution (formula 21), without any hypothesis, belong to Exact Solutions; The source of error of height analytic solution (formula 18) launches in sin γ first order Taylor, usual γ < 2deg, and therefore the relative error of height analytic solution is approximately 10 ^-4; Flying distance analytic solution (formula 20) do not introduce new error term, and therefore relative error is with highly identical; Introduce the five rank Taylor expansions of cos ψ in latitude analytic solution (formula 29), if the variation range of ψ is [0, π], then launch error and be 0.0047 to the maximum, therefore the relative error of latitude analytic solution is between 0.1% to 1%; And when solving longitude analytic solution, employ latitude analytic solution, and there is quadrature error, therefore its relative error is maximum, about 1%.

Step 4: glide section endgame speed-control scheme;

In step 3, given k is given _εand k _βglide section trajectory analytic solution under condition, but do not provide the analytic solution of speed.In order to obtain the size of terminal velocity, need to utilize above-mentioned analytic solution to carry out numerical integration to speed.

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

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

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

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

As can be seen from formula (33), C _dit is the key that speed solves.By in situation, C _dc can be write as _lit is with the function of e, as follows,

C _D＝f _CD(C _L,e) (34)

In formula, f _cDfor the relation function of lift coefficient and resistance coefficient.Solve the key of formula (33) integration for solving C _l, can be obtained by formula (6), (7) formula and formula (10),

$C_L\cos \mathrm{\&sigma;}=\frac{m[g/(\mathrm{\&mu;}/r-e)-2/r]\cos \mathrm{\&gamma;}}{\mathrm{\&rho;}{S}_{\mathrm{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;}{\mathrm{rS}}_{\mathrm{ref}}}---\left(36\right)$

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

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

C can be tried to achieve by (35), formula (36) and formula (37) _l, by C _lbring into formula (34) can C _dsize.

Finally, by r and C _dexpression formula bring formula (33) into, and adopt Runge-Kutta method to carry out integration, terminal can be obtained and bear specific energy, thus solve the velocity magnitude of terminal.It is pointed out that and easily know C by formula (33) _dlarger then terminal velocity is less, and C _llarger then C _dlarger, therefore adjust C _lcan control terminal speed.Due to k _εbe generally one in a small amount, also determine the size of range, therefore can not be used for speed and regulate; Thus k can only be adopted _βadjust terminal velocity.Fig. 3 gives with V _frelation, can find out that linearization degree is between the two very high, be conducive to numerical solution during trajectory planning.

When the trajectory planned comprises canting rollback point, then the k before and after rollback point _βthe following condition of contact of demand fulfillment,

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

In formula, k _β-and k _β+be respectively the crossrange maneuvering acceleration scale-up factor before and after rollback point; ρ _c, ψ _cand φ _cbe respectively γ _cthe atmospheric density located, trajectory deflection angle and latitude.

Step 5: glide section trajectory reentry corridor Adjusted Option;

Adopt the position of initial trajectory Inclination maneuver reentry trajectory in corridor.As shown in Figure 4, initial trajectory tilt angle gamma ₀larger, then the height of reentry trajectory is higher, and it is far away that distance retrains by the constraint of maximum heat current density, max-Q constraint and maximum overload the reentry corridor lower boundary determined.But consider the flatness of the glide section angle of attack, initial trajectory inclination angle value is as follows,

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

In formula, Δ γ is the initial trajectory inclination angle component of adjustment trajectory reentry corridor position; γ ^*for meeting the initial trajectory inclination angle component that angle of attack flatness requires, meet,

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

K in formula _n=L cos σ/D is longitudinal lift-drag ratio.It is pointed out that the glide section trajectory due to the present invention's planning meets steady glide condition, thus heat flow density, dynamic pressure and overload change are comparatively mild, and peak value is also less, and therefore Δ γ is generally zero.

Step 6: glide section trajectory planning forming initial fields;

Before carrying out glide section trajectory planning, need s _f, γ ^*, k _εand k _βinitial value estimate.

The first, s _festimation

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

$s_s=R_0\arccos \left[{\left({\stackrel{\&RightArrow;}{r}}_0\right)}^T{\stackrel{\&RightArrow;}{r}}_f\right]$

In formula, s _sfor (θ ₀, φ ₀) to (θ _f, φ _f) spherical distance; with meet respectively,

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

Consider that the situation in the face of penetrating is departed from initial trajectory drift angle, terminal flying distance is,

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

In formula, wherein for glide section starting point meridian and the angle penetrating place, face great circle.

The second, γ ^*estimation

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

$K_N=\frac{2s_f}{r\ln [(\mathrm{gr}-V_f^2)/(\mathrm{gr}-V_0^2)]}---\left(42\right)$

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

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

Three, k _εestimation

If Δ γ=0, then γ ₀=γ ^*.From formula (17) and formula (20), γ _f, k _εand s _fmeet following relation,

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

s _f＝f _s(γ _f)-f _s(γ ₀)

In formula, γ _fand ρ _fbe respectively terminal trajectory tilt angle and terminal atmospheric density.From formula (44), given s _f, then the k of existence anduniquess _εwith it corresponding (as shown in Figure 6).By the s that formula (41) obtains _fbring formula (44) into solve and can obtain k _εvaluation.

Four, k _βestimation

As shown in Figure 2, at k _βtime less, with V _flinear relationship better.Therefore at estimation k _βprinciple during initial value is rather little not large, carries out iterated revision after being conducive to like this.Known by formula (21), at known (γ-γ ₀)/k _εcondition under, k _βonly with ψ _f-ψ ₀relevant.Suppose ψ _f-ψ ₀only for revising ψ ₁₀impact, then according to plane circular arc trajectory hypothesis have,

ψ _f-ψ ₀＝2ψ ₁₀

Thus can k be obtained _βinitial valuation be,

k _β＝2ψ ₁₀C _N1/(γ-γ ₀) (45)

Step 7: glide section trajectory planning flow scheme design;

Glide section trajectory planning task is find suitable k _ε, k _β, Δ γ and angle of heel rollback point γ _c, make reentry vehicle from glide starting point (θ ₀, φ ₀, r ₀, V ₀, ψ ₀) to the terminal (θ that glides _f, φ _f, r _f, V _f).Above-mentioned four planning variablees correspond to different constraint requirements respectively, are coupled hardly each other, wherein: k _εfor adjusting the flying distance of glide section; k _βfor adjusting the speed of glide segment endpoint; γ _cfor adjusting the final position of glide section; Δ γ will be used for adjusting reentry trajectory height, to meet process constraints requirement.Utilize glide section ballistic-shaped analytic solution, glide section trajectory planning can be split into ballistic-shaped planning and terminal velocity correct two parts.Wherein, ballistic-shaped planning does not need to carry out trajectory integration, solves terminal location restricted problem; The velocity correction of terminal end then needs to carry out integration to velocity differentials equation, solves terminal velocity constraint and process constraints problem.Complete glide section trajectory flow process as shown in Figure 7, is described as follows:

The first, according to step 6, s is obtained _f, γ ^*, k _εand k _βplanning initial value, and suppose Δ γ=0 and γ _c=γ ₀.

The second, according to given Δ γ, γ ^*, γ _c, k _εand k _β, utilize step 3 to calculate θ (γ _f) and φ (γ _f).Wherein, γ _cfor rollback point trajectory tilt angle; θ (γ _f) and φ (γ _f) be planning endgame longitude and latitude.

Three, γ is revised _c, terminal latitude is met the demands, and the method for correction is Newton iteration method, and the end condition of iteration is | φ (γ _f)-φ _f| < φ _limit.φ in formula _limitfor latitude planning precision, consider the parsing ballistic computation precision in step 3, desirable φ _limit=0.01deg.If | φ (γ _f)-φ _f|>=φ _limitthen turn to the second step in step 7.

Four, s is revised _f(or k _ε) terminal longitude is met the demands, modification method is Newton iteration method, and stopping criterion for iteration is | θ (γ _f)-θ _f| < θ _limit.θ in formula _limitfor longitude planning precision, consider the parsing ballistic computation precision in step 3, desirable θ _limit=0.1deg.If | θ (γ _f)-θ _f|>=θ _limitthen turn to the second step in step 7.

Five, according to Δ γ, γ ^*, γ _c, k _εand k _β, utilize step 4 integration active terminal speed V (γ _f), and adopt Newton iteration method correction k _β, make | V (γ _f)-V _f| < V _limit.V in formula _limitfor speed planning precision, consider that step 3 resolves ballistic computation precision, desirable V _limit=10m/s.If | V (γ _f)-V _f|>=V _limitthen turn to the second step in step 7.

6th, the integration trajectory according to previous step judges whether to meet process constraints, and revises as follows Δ γ,

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

In formula, h is practical flight height; h _limitfor being retrained by maximum heat current density, max-Q retrains and maximum overload retrains the reentry corridor height lower boundary determined; Min (h-h _limit) < 0 represents that glide section trajectory exceeds reentry corridor, need to increase Δ γ, and turn to the second step in step 7; Min (h-h _limit)>=0 represents that glide section trajectory meets reentry corridor requirement; γ _limfor min (h-h _limit) trajectory tilt angle when getting extreme value; Δ γ ^(old)with Δ γ ^(new)be respectively before revising and revised initial trajectory inclination correction amount; H (γ _lim), s (γ _lim) and h _limit(γ _lim) be respectively trajectory tilt angle and get γ _limtime height, flying distance and height lower boundary.

Case study on implementation:

In order to check the precision of resolving derivation algorithm, selecting CAV as computation model, carrying out numerical simulation effect.Emulation adopts without slewing circle spherical model and exponential atmosphere model, and optimum configurations is as form 1.Simulation computer CPU is Core i3-2120, internal memory 4GB, and simulated environment is Matlab.

Form 1 simulation parameter is arranged

The first, analytic solution precision effect

Before carrying out trajectory planning, first ballistic-shaped analytic solution precision is verified, and contrast with Bell analytic solution.Fig. 8 and Fig. 9 sets forth the contrast of Longitudinal Trajectory analytic solution and horizontal trajectory analytic solution and trajectory integral result, can find out, ballistic-shaped analytic solution almost overlap completely with trajectory integral result, has high precision.Specifically, adopt analytic solution directly to predict glide segment endpoint position, its height error is less than 1m; Trajectory deflection angle error is less than 0.01deg; Latitude error under great circle coordinate system is less than 0.1deg, and longitude error is less than 0.3deg; Demonstrate the error analysis in literary composition.Figure 10 gives the Bell horizontal trajectory analytic solution provided, and can find out that Bell analytic solution only have higher precision in less scope, when flying distance is far away, horizontal journey deviation is larger.This is owing to have ignored the impact of earth curvature on horizontal trajectory in Bell analytic solution, and adopts plane coordinate system to solve, and brings larger to solve error.

The second, trajectory planning range coverage is analyzed

After the precision demonstrating ballistic-shaped analytic solution, also need the range coverage analyzing trajectory planning.When solving range coverage, assuming that trajectory deflection angle being fixed value, thus regulating terminal location mainly through canting reversion, being therefore the border of trajectory planning range coverage without the trajectory reversed.Figure 11 to Figure 18 gives border ballistic curve race.From Figure 11 and Figure 12, the planning angle of attack curve of trajectory and angle of heel curve unusual light and amplitude of variation is less; The range coverage of the trajectory planning of glide section is as shown in Figure 13 indulged journey coverage and is about 11000 kms, and horizontal journey coverage is about 6700 kms; Met the demands by the terminal height of the known planning trajectory of Figure 14 and Figure 15 and terminal velocity; Process constraints requirement is met by the known planning trajectory of Figure 16 to Figure 18.For the situation of target in range coverage, owing to only increasing angle of heel rollback point to regulate terminal longitude and latitude, therefore the trend of its angle of attack curve, canting curve, Longitudinal Trajectory curve and process constraints curve is by similar with the result that provides in Figure 11 to Figure 18.

On the basis of Figure 13, Figure 19 obtains the range coverage of glide section planning trajectory further.Compared with maximum range coverage of steadily gliding, the range coverage that the present invention obtains smaller, maximum range loss is about 81.7km, and minimum range loss is about 302.8km, and maximum horizontal journey loss is about 923.9km.This is because the present invention selects fixing k _βas crossrange maneuvering planning variable, make angle of heel in high maneuver situation there is larger change (as shown in figure 11), this can increase the motor-driven energy loss of horizontal journey, thus reduces crossrange maneuvering ability.When range is less or horizontal journey is larger, need larger crossrange maneuvering, thus occur the range coverage loss shown in Figure 18.

3rd, glide section trajectory is planned fast

As previously shown, adopt the glide section trajectory planning flow process that Fig. 5 provides, the impact point in the range coverage provide Figure 19 carries out trajectory planning, and the result obtained is as shown in form 1 and Figure 20 to Figure 22 for the initial value of trajectory planning and final value.As can be seen from form 2, the glide section trajectory planning method based on feedback linearization analytic solution can cook up the reentry trajectory that meets constraint requirements within the time of about 0.03s, and endgame has higher positional precision.Deficiency is that the terminal velocity error planning trajectory is comparatively large, and this is according to geometry trajectory, adopts dynamic inversion integration to obtain due to terminal velocity, and the little deviation in ballistic-shaped is amplified by the integration of solving speed.Figure 20 gives the shape of planning trajectory; Figure 21 has provided out the angle of attack curve of planning trajectory, and can find out that angle of attack variation is comparatively smooth, amplitude of variation is less than 4deg; Figure 22 gives angle of heel curve, equally also comparatively smooth.

Form 2 section of glide trajectory planning result

Claims (1)

1., based on the quick trajectory planning method steadily sliding trajectory analytic solution, comprise following step:

Step 1: glide section trajectory planning problem modeling

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

$\begin{array}{cc}\stackrel{.}{h}=V\sin \mathrm{\&gamma;}& \stackrel{.}{\mathrm{\&theta;}}=V\cos \mathrm{\&gamma;}\sin \mathrm{\&psi;}/\left(r\cos \mathrm{\&phi;}\right)\end{array}$

$\begin{array}{cc}\stackrel{.}{s}=V\cos \mathrm{\&gamma;}{R}_0/r& \stackrel{.}{\mathrm{\&phi;}}=V\cos \mathrm{\&gamma;}\cos \mathrm{\&psi;}/r\end{array}$

(1)

$\begin{array}{cc}\stackrel{.}{\mathrm{\&gamma;}}=[L\cos \mathrm{\&sigma;}+(V^2/r-g)\cos \mathrm{\&gamma;}]/V& \stackrel{.}{\mathrm{\&psi;}}=[L\sin \mathrm{\&sigma;}/\cos \mathrm{\&gamma;}+(V^2/r)\cos \mathrm{\&gamma;}\sin \mathrm{\&psi;}\tan \mathrm{\&phi;}]/V\end{array}$

$\stackrel{.}{V}=-D-g\sin \mathrm{\&gamma;}$

In formula, h is the height of aircraft, and 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 trajectory deflection angle; with be respectively the derivative of height, flying distance, longitude, latitude, trajectory tilt angle, trajectory deflection angle and velocity versus time; σ is angle of heel; G=μ/r ²for local gravitational acceleration, μ is the constant that terrestrial gravitation acceleration is relevant; L and D is 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 _dbe respectively lift coefficient and resistance coefficient, relevant with angle of attack to Mach number Ma; ρ is atmospheric density, as follows,

ρ＝ρ _seae ^-βh(4)

In formula, ρ _seafor sea-level atmosphere density; β is exponential atmosphere model constants;

Process constraints is as follows

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

In formula, for heat flow density, for maximum heat current 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: glide 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 _βbe respectively longitudinal maneuver acceleration and crossrange maneuvering acceleration; Bring formula (6) and formula (7) into formula (1) respectively, then the trajectory tilt angle differential equation and the trajectory deflection angle differential equation are,

$\begin{array}{cc}\stackrel{.}{\mathrm{\&gamma;}}=a_{\mathrm{\&epsiv;}}/V& \stackrel{.}{\mathrm{\&psi;}}=a_{\mathrm{\&beta;}}/V---\left(8\right)\end{array}$

Suppose a _ε< 0, then the γ monotone decreasing of the section of glide, replaces with γ by independent variable in the equation of motion; By formula (1) divided by the trajectory tilt angle differential equation in formula (8),

$\begin{array}{cc}\frac{\mathrm{dh}}{\mathrm{d\&gamma;}}=\frac{V^2\sin \mathrm{\&gamma;}}{a_{\mathrm{\&epsiv;}}}& \frac{\mathrm{ds}}{\mathrm{d\&gamma;}}=\frac{R_0{V}^2\cos \mathrm{\&gamma;}}{{\mathrm{ra}}_{\mathrm{\&epsiv;}}}\end{array}$

$\begin{array}{ccc}\frac{\mathrm{d\&psi;}}{\mathrm{d\&gamma;}}=\frac{a_{\mathrm{\&beta;}}}{a_{\mathrm{\&epsiv;}}}& \frac{\mathrm{d\&theta;}}{\mathrm{d\&gamma;}}=\frac{V^2\cos \mathrm{\&gamma;}\sin \mathrm{\&psi;}}{{\mathrm{ra}}_{\mathrm{\&epsiv;}}\cos \mathrm{\&phi;}}& \frac{\mathrm{d\&phi;}}{\mathrm{d\&gamma;}}=\frac{V^2\cos \mathrm{\&gamma;}\cos \mathrm{\&psi;}}{{\mathrm{ra}}_{\mathrm{\&epsiv;}}}\end{array}---\left(9\right)$

$\frac{\mathrm{dV}}{\mathrm{d\&gamma;}}=-\frac{(D+g\sin \mathrm{\&gamma;})V}{a_{\mathrm{\&epsiv;}}}$

In formula, with be respectively the differential equation of height, flying distance, trajectory deflection angle, longitude, latitude and speed; Wherein, the first behavior glide section lengthwise movement differential equation, the second behavior glide section transverse movement differential equation, the third line is velocity differentials equation;

If a of glide section _εwith a _βform is as follows,

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

In formula, k _εand k _βbe respectively longitudinal maneuver acceleration scale-up factor and crossrange maneuvering acceleration scale-up factor; Bring formula (10) into formula (9), and get sin γ=γ, cos γ=1, then have

$\frac{\mathrm{dh}}{\mathrm{d\&gamma;}}=\frac{2m}{\mathrm{\&rho;}{S}_{\mathrm{ref}}{k}_{\mathrm{\&epsiv;}}}\mathrm{\&gamma;}---\left(11\right)$

$\frac{\mathrm{ds}}{\mathrm{d\&gamma;}}=\frac{2m}{\mathrm{\&rho;}{S}_{\mathrm{ref}}{k}_{\mathrm{\&epsiv;}}}\frac{R_0}{r}---\left(12\right)$

$\frac{\mathrm{d\&psi;}}{\mathrm{d\&gamma;}}=\frac{k_{\mathrm{\&beta;}}}{k_{\mathrm{\&epsiv;}}}---\left(13\right)$

$\frac{\mathrm{d\&theta;}}{\mathrm{d\&gamma;}}=\frac{2m}{\mathrm{\&rho;}{S}_{\mathrm{ref}}{k}_{\mathrm{\&epsiv;}}}\frac{\sin \mathrm{\&psi;}}{r\cos \mathrm{\&phi;}}---\left(14\right)$

$\frac{\mathrm{d\&phi;}}{\mathrm{d\&gamma;}}=\frac{2m}{\mathrm{\&rho;}{S}_{\mathrm{ref}}{k}_{\mathrm{\&epsiv;}}}\frac{\cos \mathrm{\&psi;}}{r}---\left(15\right)$

With k _εand k _βas planning variable, glide section trajectory planning problem is converted into the problem solving nonlinear equation;

Step 3: glide section trajectory analytic solution solve;

If the initial trajectory inclination angle of glide section is γ ₀, elemental height is h ₀, initial trajectory drift angle is ψ ₀, initial longitude and latitude be for being respectively θ ₀and φ ₀; The terminal height of glide section is h _f, for given k _εand k _β, Analytical Solution is carried out to formula (11) to (15) formula;

The first, height analytic solution

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

$e^{-\mathrm{\&beta;h}}\mathrm{dh}=\frac{2m}{{\mathrm{\&rho;}}_{\mathrm{sea}}{S}_{\mathrm{ref}}{k}_{\mathrm{\&epsiv;}}}\mathrm{\&gamma;d\&gamma;}---\left(16\right)$

Can obtain formula (16) integration,

${\mathrm{\&gamma;}}^2-{\mathrm{\&gamma;}}_0^2=\frac{S_{\mathrm{ref}}{k}_{\mathrm{\&epsiv;}}}{\mathrm{m\&beta;}}({\mathrm{\&rho;}}_0-\mathrm{\&rho;})---\left(17\right)$

In formula, γ ₀for initial trajectory inclination angle; ρ ₀for glide section initial atmosphere density; Formula (17) is arranged further and can be obtained,

$\frac{m}{\mathrm{\&rho;}{S}_{\mathrm{ref}}{k}_{\mathrm{\&epsiv;}}}=\frac{1}{\mathrm{\&beta;}(C_1-{\mathrm{\&gamma;}}^2)}---\left(18\right)$

In formula, C ₁be a constant, value is $C_1=S_{\mathrm{ref}}{k}_{\mathrm{\&epsiv;}}{\mathrm{\&rho;}}_0/\left(\mathrm{m\&beta;}\right)+{\mathrm{\&gamma;}}_0^2;$

The second, flying distance analytic solution

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

$\frac{\mathrm{ds}}{\mathrm{d\&gamma;}}=\frac{-2R_0}{\mathrm{\&beta;r}({\mathrm{\&gamma;}}^2-C_1)}---\left(19\right)$

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

s＝f _s(γ)-f _s(γ ₀) (20)

In formula, f _s(γ) be the function obtained by formula (19) indefinite integral, its expression formula and C ₁positive negative correlation, 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_0=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;

Three, drift angle, road analytic solution

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

ψ＝ψ ₀+(k _β/k _ε)(γ-γ ₀) (21)

In formula, ψ ₀for initial trajectory drift angle;

Four, latitude analytic solution

By cos ψ 5 rank Taylor expansions near pi/2, 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)$

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

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

In formula, C ₂be a constant, meet C ₂=γ ₀+ k _ε(pi/2-ψ ₀)/k _β; Get,

C _a0＝k _β/(60rβk _ε)C _a1＝2C ₂k _β/k _ε

(24)

$C_{a2}=(C_2^2-C_1){(k_{\mathrm{\&beta;}}/k_{\mathrm{\&epsiv;}})}^2$

In formula, C _a0, C _a1and C _a2be constant; Bring formula (24) into formula (23) can obtain,

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

Get further,

$C_{b1}=C_{a1}^2-C_{a2}-20$

$C_{b2}=C_{a1}^3-2C_{a1}{C}_{a2}-20C_{a1}$

(26)

$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}$

In formula, C _b1, C _b2, C _b3and C _b4be constant; Utilize formula (26) formula (25) can be turned to,

$\begin{array}{c}\frac{\mathrm{d\&phi;}}{\mathrm{dy}}=C_{a0}(y^3+C_{a1}{y}^2+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)$

Can obtain formula (27) integration,

φ＝φ ₀+f _φ(y)-f _φ(y ₀) (28)

In formula, φ ₀for initial latitude; y ₀relevant to initial trajectory drift angle, value is y ₀=pi/2-ψ ₀; f _φy () is the function relevant to y, 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 relevant to y, as follows

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

Five, longitude analytic solution

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

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

Utilize Guass-Legendre quadrature formula directly to obtain result, expression is as follows,

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

In formula, γ _ifor [γ, γ ₀] in i-th Gaussian node; A _ifor quadrature coefficient; N is Gaussian node sum;

Step 4: glide section endgame speed-control scheme;

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

$\stackrel{.}{e}=\mathrm{VD}---\left(32\right)$

In formula, for the derivative of negative specific energy; Can be obtained by formula (8), formula (10) and formula (32),

$\frac{\mathrm{de}}{\mathrm{d\&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 that speed solves; By in situation, C _dc can be write as _lit is with the function of e, as follows,

C _D＝f _CD(C _L,e) (34)

In formula, f _cDfor the relation function of lift coefficient and resistance coefficient; Solve the key of formula (33) integration for solving C _l, can be obtained by formula (6), (7) formula and formula (10),

$C_L\cos \mathrm{\&sigma;}=\frac{m[g/(\mathrm{\&mu;}/r-e)-2/r]\cos \mathrm{\&gamma;}}{\mathrm{\&rho;}{S}_{\mathrm{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;r}{S}_{\mathrm{ref}}}---\left(36\right)$

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

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

C can be tried to achieve by (35), formula (36) and formula (37) _l, by C _lbring into formula (34) can C _dsize;

Finally, by r and C _dexpression formula bring formula (33) into, and adopt Runge-Kutta method to carry out integration, obtain terminal and bear specific energy, thus solve the velocity magnitude of terminal;

When the trajectory planned comprises canting rollback point, then the k before and after rollback point _βthe following condition of contact of demand fulfillment,

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

In formula, k _β-and k _β+be respectively the crossrange maneuvering acceleration scale-up factor before and after rollback point; ρ _c, ψ _cand φ _cbe respectively γ _cthe atmospheric density located, trajectory deflection angle and latitude;

Step 5: glide section trajectory reentry corridor Adjusted Option;

Initial trajectory inclination angle value is as follows,

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

In formula, Δ γ is the initial trajectory inclination angle component of adjustment trajectory reentry corridor position; γ ^*for meeting the initial trajectory inclination angle component that angle of attack flatness requires, meet,

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

K in formula _n=Lcos σ/D is longitudinal lift-drag ratio;

Step 6: glide section trajectory planning forming initial fields;

Before carrying out glide section trajectory planning, need s _f, γ ^*, k _εand k _βinitial value estimate;

The first, s _festimation

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

$s_s=R_0\arccos \left[{\left({\stackrel{\&RightArrow;}{r}}_0\right)}^T{\stackrel{\&RightArrow;}{r}}_f\right]$

In formula, s _sfor (θ ₀, φ ₀) to (θ _f, φ _f) spherical distance; with meet respectively,

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

Then terminal flying distance is,

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

In formula, wherein for glide section starting point meridian and the angle penetrating place, face great circle;

The second, γ ^*estimation

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

$K_N=\frac{2s_f}{r\ln [(\mathrm{gr}-V_f^2)/(\mathrm{gr}-V_0^2)]}---\left(42\right)$

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

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

Three, k _εestimation

If Δ γ=0, then γ ₀=γ ^*; From formula (17) and formula (20), γ _f, k _εand s _fmeet following relation,

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

s _f＝f _s(γ _f)-f _s(γ ₀)

In formula, γ _fand ρ _fbe respectively terminal trajectory tilt angle and terminal atmospheric density; From formula (44), given s _f, then the k of existence anduniquess _εcorresponding with it; By the s that formula (41) obtains _fbring formula (44) into solve and can obtain k _εvaluation;

Four, k _βestimation

Suppose ψ _f-ψ ₀only for revising ψ ₁₀impact, then according to plane circular arc trajectory hypothesis have,

ψ _f-ψ ₀＝2ψ ₁₀

K _βinitial valuation be,

k _β＝2ψ ₁₀C _N1/(γ-γ ₀) (45)

Step 7: glide section trajectory planning flow scheme design;

Be specially:

The first, according to step 6, s is obtained _f, γ ^*, k _εand k _βplanning initial value, and suppose Δ γ=0 and γ _c=γ ₀;

The second, according to given Δ γ, γ ^*, γ _c, k _εand k _β, utilize step 3 to calculate θ (γ _f) and φ (γ _f); Wherein, γ _cfor rollback point trajectory tilt angle; θ (γ _f) and φ (γ _f) be planning endgame longitude and latitude;

Three, γ is revised _c, terminal latitude is met the demands, and the method for correction is Newton iteration method, and the end condition of iteration is | φ (γ _f)-φ _f| < φ _limit; φ in formula _limitfor latitude planning precision; If | φ (γ _f)-φ _f|>=φ _limitthen turn to the second step in step 7;

Four, s is revised _for k _εterminal longitude is met the demands, and modification method is Newton iteration method, and stopping criterion for iteration is | θ (γ _f)-θ _f| < θ _limit; θ in formula _limitfor longitude planning precision; If | θ (γ _f)-θ _f|>=θ _limitthen turn to the second step in step 7;

Five, according to Δ γ, γ ^*, γ _c, k _εand k _β, utilize step 4 integration active terminal speed V (γ _f), and adopt Newton iteration method correction k _β, make | V (γ _f)-V _f| < V _limit; V in formula _limitfor speed planning precision; If | V (γ _f)-V _f|>=V _limitthen turn to the second step in step 7;

6th, the integration trajectory according to previous step judges whether to meet process constraints, and revises as follows Δ γ,

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

In formula, h is practical flight height; h _limitfor being retrained by maximum heat current density, max-Q retrains and maximum overload retrains the reentry corridor height lower boundary determined; Min (h-h _limit) < 0 represents that glide section trajectory exceeds reentry corridor, need to increase Δ γ, and turn to the second step in step 7; Min (h-h _limit)>=0 represents that glide section trajectory meets reentry corridor requirement; γ _limfor min (h-h _limit) trajectory tilt angle when getting extreme value; Δ γ ^(old)with Δ γ ^(new)be respectively before revising and revised initial trajectory inclination correction amount; H (γ _lim), s (γ _lim) and h _limit(γ _lim) be respectively trajectory tilt angle and get γ _limtime height, flying distance and height lower boundary.