CN111306989B - Hypersonic velocity reentry guidance method based on steady glide trajectory analytic solution - Google Patents

Hypersonic velocity reentry guidance method based on steady glide trajectory analytic solution Download PDF

Info

Publication number
CN111306989B
CN111306989B CN202010170138.1A CN202010170138A CN111306989B CN 111306989 B CN111306989 B CN 111306989B CN 202010170138 A CN202010170138 A CN 202010170138A CN 111306989 B CN111306989 B CN 111306989B
Authority
CN
China
Prior art keywords
angle
terminal
taem
attack
representing
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202010170138.1A
Other languages
Chinese (zh)
Other versions
CN111306989A (en
Inventor
陈万春
张晚晴
余文斌
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Beihang University
Original Assignee
Beihang University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Beihang University filed Critical Beihang University
Priority to CN202010170138.1A priority Critical patent/CN111306989B/en
Publication of CN111306989A publication Critical patent/CN111306989A/en
Application granted granted Critical
Publication of CN111306989B publication Critical patent/CN111306989B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • F—MECHANICAL ENGINEERING; LIGHTING; HEATING; WEAPONS; BLASTING
    • F41—WEAPONS
    • F41G—WEAPON SIGHTS; AIMING
    • F41G3/00—Aiming or laying means
    • F41G3/22—Aiming or laying means for vehicle-borne armament, e.g. on aircraft

Landscapes

  • Engineering & Computer Science (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • General Engineering & Computer Science (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Feedback Control In General (AREA)

Abstract

The invention relates to a hypersonic reentry guidance method based on a steady gliding trajectory analytic solution, which comprises the following steps: 1. establishing an aircraft dynamics model; 2. establishing a constraint model; 3. an initial descent stage guidance method; 4. a steady glide phase guidance method; 5. a height adjustment section guidance method; has the advantages that: the longitudinal and transverse courses are predicted and corrected based on the three-dimensional trajectory analytic solution, so that the calculation amount is greatly reduced, the calculation speed is improved, and the method is convenient to apply to an airborne system in real time. The transverse movement is determined by twice tilting reversal points, the position of the first reversal point is corrected to control the terminal traverse error, and the position of the second reversal point is corrected to meet the terminal speed constraint. Compared with the traditional method for controlling the transverse motion by using the course error threshold, the method has the advantages of reducing the number of times of tilting and reversing, improving the energy utilization efficiency and reducing the requirements on a control system. The height adjusting section utilizes the idea of proportional guidance, so that the terminal roll angle is converged to 0, and the target hitting effect is enhanced.

Description

Hypersonic velocity reentry guidance method based on steady glide trajectory analytic solution
Technical Field
The invention provides a hypersonic velocity reentry guidance method based on steady glide trajectory analytic solution, belonging to the technical field of space technology and weapons
Background
The design of the reentry guidance law is a core technology of whether the hypersonic reentry gliding aircraft can finish the accurate hitting task. The ultimate goal of reentry guidance is to achieve successful guidance of the hypersonic reentry vehicle to a specified target while satisfying path constraints and terminal constraints.
The reentry aircraft is subjected to strict path constraint in the reentry process, however, the traditional algorithm is always difficult to process the path constraint, and difficulties and challenges are brought to the online generation of the trajectory and the online calculation of the guidance law. Meanwhile, the extremely fast glide flight speed puts forward higher requirements on the calculation time, so that the time constraint generated by the online guidance law instruction is very strong, and for the task of target relocation, the algorithm is required to have stronger adaptivity and rapidity so as to adapt to the requirements of the large-maneuvering flight task.
Disclosure of Invention
Aiming at the problems of reentry guidance of the hypersonic aircraft, a high-precision hypersonic reentry guidance method based on a three-dimensional trajectory analytic solution is provided. Firstly, an auxiliary geocentric rotation coordinate system is constructed, a new kinetic model of the longitudinal course, the transverse course and the course angle with respect to the energy is established on the basis of the auxiliary geocentric rotation coordinate system, the model is favorable for decoupling and linearizing the highly nonlinear and strongly coupled kinetic model to obtain a three-dimensional ballistic analytic solution, and a guidance method is designed on the basis of the analytic solution. The guidance method utilizes a longitudinal analytic solution to plan a longitudinal reference profile in real time, utilizes a transverse analytic solution to correct a first inversion point to eliminate a transverse error, utilizes trajectory simulation for several times to correct a second inversion point to ensure terminal speed constraint and corrects the minimum value of an attack angle to ensure terminal height constraint.
The invention discloses a hypersonic reentry guidance method based on a steady glide trajectory analytic solution, which is characterized by comprising the following steps: the method comprises the following steps:
the method comprises the following steps: establishing an aircraft dynamics model;
the hypersonic flight vehicle generally refers to an aircraft with flight Mach number larger than 5 and capable of realizing hypersonic flight in an atmospheric layer and a trans-atmospheric layer. Under the background of the rotating earth, a standard atmosphere model and an inverse square gravitational field model are adopted to establish a six-degree-of-freedom kinetic equation of the hypersonic aerocraft as follows:
Figure BDA0002408897380000021
Figure BDA0002408897380000022
Figure BDA0002408897380000023
Figure BDA0002408897380000024
Figure BDA0002408897380000025
Figure BDA0002408897380000026
Figure BDA0002408897380000027
where d/dt represents the derivative over time, t is time, h is altitude, R is range, λ is longitude, φ is latitude, V is the velocity of the aircraft relative to the earth, γ is the ballistic inclination, ψ is the aircraft heading angle, based on local north, σ is the roll angle, R is the altitude, and R is the altitude, whereeIs the radius of the earth, and takes the value of 6378.137km, omegaeIs the angular velocity of rotation of the earth, g represents the acceleration of gravity; l ═ ρ V2SCLThe lift acceleration is given by/2 m, and D is equal to rho V2SCD(2 m) is the acceleration of resistance, CL、CDRespectively lift coefficient and drag coefficient.
Step two: establishing a constraint model;
the reentry aircraft needs to meet the process constraints of heat flux density, dynamic pressure, overload and the like in the flight process, and the process constraints are expressed as follows:
Figure BDA0002408897380000028
Figure BDA0002408897380000031
Figure BDA0002408897380000032
parameter p0Denotes nominal atmospheric density, e denotes natural constant, h denotes altitude, hSIndicating nominal height, n overload, and m mass.
Wherein, KQIs the heat flux density coefficient;
Figure BDA0002408897380000033
representing a maximum allowable heat flow density value; n represents an allowable overload value; n ismaxRepresents a maximum allowable overload value; q represents an allowable dynamic pressure value; q. q.smaxIndicating the maximum allowable dynamic pressure value.
Furthermore, due to the capability constraints of the flight control system, the rate of change of the control quantities (angle of attack α and roll angle σ) also satisfies the constraints:
Figure BDA0002408897380000034
in addition to process constraints, the aircraft should also satisfy terminal constraints, with: vTAEM＝2000m/s,HTAEM＝25km,Rtm＝100km,σTAEM≈0°,γTAEM≈0°andΔψTAEM0 deg. where the subscript TAEM denotes the terminal amount. VTAEM、HTAEM、σTAEM、γTAEM、ΔψTAEMRespectively representing the terminal speed, the terminal height, the terminal roll angle, the terminal trajectory inclination angle and the terminal course angle error. RtmIndicating the shot-eye distance.
Step three: an initial descent stage guidance method;
according to the characteristics of reentry trajectory, the reentry trajectory is divided into three parts: an initial descent section, a steady glide section and a height adjustment section. The initial descent segment is designed to allow a more gradual transition of the trajectory to a smooth glide segment. The whole flight trajectory is mostly in a stable gliding section, which not only determines the range, speed and height of the flight terminal, but also determines whether the aircraft can meet various process constraints. The height adjustment segment is the last phase of the trajectory, as to whether the trajectory can strictly meet the terminal constraints. The flying speed and flying height of the aircraft are low in this section, and severe requirements are provided for a guidance method for designing a height adjusting section.
In the initial descent, the variation of the controlled variable has little influence on the flight, so that the fixed values of the angle of attack and the roll angle are designed to avoid altitude variations. To decrease the heat flux density and increase the range, the starting point of the trajectory should be set higher. The maximum angle of attack α is therefore chosen heremaxFly 20deg and zero roll angle. Rate of change of ballistic inclination
Figure BDA0002408897380000035
At this time, it means that the lift force can satisfy the requirement of smooth gliding, so that the attack angle can be smoothly transited from the maximum attack angle to the reference attack angle section alphabsl。
At this stage, the command angle of attack and roll are designed as follows:
Figure BDA0002408897380000041
σcmd＝0 (12)
wherein,
Δγ＝γ-γSG (13)
Figure BDA0002408897380000042
wherein, γSGRepresenting the angle of inclination, k, of the smooth gliding trajectoryγRepresenting the feedback coefficient and here taking the value 5. When Δ γ is 0, the aircraft enters a smooth glide phase. Parameter alphacmdIs an instruction angle of attack; sigmacmdIs a commanded roll angle; delta gamma is the deviation of the current trajectory inclination angle and the steady glide trajectory inclination angle; delta gamma1As the rate of change of ballistic inclination
Figure BDA0002408897380000043
Deviation value of trajectory inclination angle at 0 from the trajectory inclination angle of smooth glide.
Step four: a steady glide phase guidance method;
ballistic planning in the steady glide phase is the basis of the proposed guidance algorithm, which requires approaching the aircraft to the target while guaranteeing process constraints. The algorithm determines the longitudinal range and lateral range profiles based on a new three-dimensional analytical solution. Here, a new ballistic analytical solution is first briefly introduced.
In order to meet the characteristic that the longitudinal stroke of the reentry trajectory is far greater than the transverse stroke, two coordinate systems are firstly established, and the meanings of the longitudinal stroke and the transverse stroke can be intuitively represented.
Defining an auxiliary earth-centered rotational reference system (AGR coordinate system) rotating with the earth: the origin is at the geocentric E,
Figure BDA0002408897380000044
the shaft points to the initial position of the aircraft,
Figure BDA0002408897380000045
the axis is in the plane of a great circle passing through the aircraft and the target point and perpendicular to the aircraft
Figure BDA0002408897380000046
The shaft is provided with a plurality of axial holes,
Figure BDA0002408897380000047
the axes may be determined according to the right hand rule.
Meanwhile, a local generalized north east down coordinate system (GNED coordinate system) is defined: the vertical projection point of the origin o from the mass center M of the aircraft to the ground;
Figure BDA0002408897380000048
the axis is vertically downwards directed to the center of the earth,
Figure BDA0002408897380000049
axis-oriented AGR coordinate system
Figure BDA00024088973800000410
The direction of the "north" direction is,
Figure BDA00024088973800000411
the axes are determined by the right hand rule.
Define the longitudinal extent as
Figure BDA00024088973800000412
Define the course as
Figure BDA00024088973800000413
Wherein
Figure BDA00024088973800000414
And
Figure BDA00024088973800000415
respectively representing a generalized longitude and a generalized latitude in an AGR coordinate system.
Defining the energy as
Figure BDA0002408897380000051
Defining the longitudinal lift-to-drag ratio Lcos sigma/D ═ L1D, transverse lift-drag ratio Lsin sigma/D ═ L2and/D. The generalized longitude and latitude coordinate system in the AGR coordinate axis is beneficial to linearizing and decoupling a complex dynamic model, and the following three-dimensional ballistic analysis solution is obtained.
Figure BDA0002408897380000052
Figure BDA0002408897380000053
Figure BDA0002408897380000054
Wherein,
Figure BDA0002408897380000055
R*＝Re+H；
Figure BDA0002408897380000056
hij,i＝1,2；j＝0,1,α1and pijI is 1, 2; j is 0,1,2,3,4 are constant coefficients. E, E0Respectively representing the current energy and the initial energyAn initial energy value; x is the number ofERepresenting an energy integral variable; l is1Representing lift values in a longitudinal plane; x is the number ofC0Representing an initial traverse value;
Figure BDA0002408897380000057
representing an initial course angle error under a GNED coordinate system;
Figure BDA0002408897380000058
indicating the heading angle error in the GNED coordinate system.
Since the analytic solution has high precision, the analytic solution can be applied to guidance method derivation to reduce calculation amount and improve universality.
1. Design of reference angle of attack
Because the analytic solution is designed by taking energy as an independent variable, all reference profiles in the guidance method also adopt energy as the independent variable so as to more conveniently utilize the analytic solution. Reference angle of attack profile alphabslThe design is as follows:
Figure BDA0002408897380000061
wherein E isα＝-5.55×107J/kg is located near the intersection of the steady glide section and the last height adjustment section. To maximize the capability of the aircraft, the angle of attack is designed to be the angle of attack corresponding to the maximum lift-to-drag ratio, i.e., alpha 110 deg. The quadratic function of the energy of the angle of attack is designed so that the angle of attack can be varied from a1Smoothly transits to alpha2. In this guidance method, α is preset for the first time2The ballistic simulation is later applied to correct it to strictly meet the terminal height constraint, 6 deg. When the design of the attack angle section is finished, the corresponding reference lift-drag ratio section (L/D)bslAs determined accordingly.
2. Baseline lift-to-drag ratio profile design
In order to meet the requirement of range, the longitudinal lift-drag ratio profile is designed as follows
Figure BDA0002408897380000062
Wherein, let the terminal longitudinal analysis solve xDfEqual to the remaining range, then obtain L1/D1. Design L1/D2＝L/D(ETAEM) In which ETAEMIs terminal energy, so when xE＝ETAEMWhen (L)1/D)bsl|E＝E_TAEM＝(LcosσTAEM/D)bsl＝L/D(ETAEM) Can make the terminal tilt angle sigmaTAEMIs 0. By designing the form of equation (20), the roll angle can be controlled to be approximately constant during flight.
Therefore, to solve for L1/D1First, the remaining range x needs to be calculatedDf. In the GER coordinate system, η is defined as a vector pointing from the center of the earth to the aircraft
Figure BDA0002408897380000063
Vector directed to aircraft from earth center
Figure BDA0002408897380000064
The included angle therebetween. The remaining range can be expressed as:
xDf＝Reη-STAEM (21)
wherein S isTAEMIndicating the distance to the target at the end of the reentry segment.
Solving for L1/D1The process of (a) can be described as follows:
when E is more than or equal to EαWhen there is
xD(Eα,E)+xD(ETAEM,Eα)＝xDf (22)
Wherein,
Figure BDA0002408897380000071
wherein,
Figure BDA0002408897380000072
Figure BDA0002408897380000073
and is provided with a plurality of groups of the following components,
Figure BDA0002408897380000074
wherein,
Figure BDA0002408897380000075
c2＝(a1-c1d0)/d1 (28)
c3＝a0-c2d0 (29)
Figure BDA0002408897380000076
Figure BDA0002408897380000077
Figure BDA0002408897380000078
the formula (22) is solved to obtain
Figure BDA0002408897380000079
Wherein,
Figure BDA0002408897380000081
Figure BDA0002408897380000082
when E < EαIn the meantime, the aircraft is mainly in the altitude adjustment phase and does not update L any more1/D1The value of (c).
3. Reference lateral lift-drag ratio and reference roll angle design
The modulus of the transverse lift-drag ratio can be well planned from the previous (L/D)bslAnd longitudinal lift-to-drag ratio (L)1/D)bslIs determined as
Figure BDA0002408897380000083
In order to satisfy the terminal traverse constraint, the transverse lift-drag ratio profile is designed as follows:
Figure BDA0002408897380000084
wherein E isBR1And EBR2Respectively, the first and second roll reversal points. sgn is a sign variable representing the initial roll reversal direction. After obtaining the transverse lift-drag ratio section, referring to the roll angle section sigmabslIt was also determined that this is expressed as follows:
Figure BDA0002408897380000085
wherein, (L/D)realShowing the actual lift-to-drag ratio profile.
4. Roll reversal point correction algorithm
In this algorithm, the roll reversal point is designed based on a course analytic solution to eliminate dead-end course errors. Here two tilting reversal points are arranged. The first reversal point E is first designed and correctedBR1. First order second reversal point EBR2＝EαThen the terminal traverse xC(ETAEME) is only EBR1As a function of (c). As can be seen from equation (17), the terminal traverse xCfWatch capable of showingShown as EBR1The piecewise function of (c) is as follows:
Figure BDA0002408897380000086
f, G is a custom function, and the expression is as follows
Figure BDA0002408897380000091
Figure BDA0002408897380000092
Wherein, the parameter xE2,xE1Expressing upper and lower energy bounds; k is a constant of 1 to 5; p1Representing polynomial fitting coefficients.
Obtaining E from equation (38)BR1Derivative of (2) can be obtained
Figure BDA0002408897380000093
Wherein E ═ EBR1Upper mark for right limit of+Is represented by E ═ EBR1Upper mark for left limit of department-And (4) showing. Because L is2/DbslWhen E is equal to EBR1Is not continuous in place, so there is
Figure BDA0002408897380000094
i is 0,1,2,3, 4. In order to make the terminal traverse error 0, x needs to be satisfiedCf(EBR1) 0. The essence of this algorithm is therefore to solve the non-linear equation xCf(EBR1) Root 0, the newton method is as follows:
Figure BDA0002408897380000095
parameter(s)
Figure BDA0002408897380000096
And
Figure BDA0002408897380000097
and not the same as the above-mentioned case,
Figure BDA0002408897380000098
represents EBR1To the power of k in the first order,
Figure BDA0002408897380000099
e representing the kth iterationBR1The value is obtained. Solution stop set to
Figure BDA00024088973800000910
Thus xCf(EBR1) With EBR1Monotonously changing, the solution of this non-linear equation can be obtained by a few simple iterations.
5. Design of command attack angle and inclination angle
To ensure the tracking accuracy of the reference profile, ballistic damping control techniques are employed herein to suppress long-period, weakly damped oscillations present in the reentry trajectory. Design command angle of attack alphacmdAnd the command roll angle sigmacmdThe following were used:
αcmd＝αbsl+cosσbslKγ(γSG-γ) (43)
Figure BDA00024088973800000911
wherein, KγIs a feedback coefficient, and takes the value of 5; gamma raySGIs an approximate solution to the angle of inclination of the smooth glide trajectory, which can be expressed as follows
Figure BDA0002408897380000101
Wherein S is the aircraft reference area; cLIs the coefficient of lift; ρ is the air density.
The steady glide phase needs to satisfy the process constraint, so the process constraint of equations (8) - (10) is converted into the roll angle constraint, which is expressed as follows
Figure BDA0002408897380000102
Wherein L ismaxFor maximum lift, the expression is as follows:
Figure BDA0002408897380000103
Figure BDA0002408897380000104
Figure BDA0002408897380000105
wherein HminThe lowest boundary of the corresponding glide height can be obtained through process constraint; constant coefficient kσ-50. Thus, to satisfy process constraints, the roll angle needs to satisfy the following conditions:
|σcmd|≤σmax (50)
step five: a height adjustment section guidance method;
1. and correcting the command attack angle and the second inversion point.
(1) Angle of attack alpha of instruction2Correction algorithm
In order to satisfy the terminal height constraint, a ballistic simulation is adopted to correct the command attack angle alpha2. Terminal height H obtained from trajectory simulationfAnd terminal constraint height HTAEMThe deviation magnitude therebetween, to further correct the final command attack angle alpha2. In the following variables, the subscript f denotes the terminal value calculated from the ballistic simulation, and the subscript TAEM denotes the terminal constraint value.
Due to gammaf,σfAnd
Figure BDA0002408897380000106
can be approximated to 0, and thus can assume cos γ, respectivelyf＝1,cosσ f1 and
Figure BDA0002408897380000107
therefore, according to equation (6):
Figure BDA0002408897380000111
wherein, CL(est)(αbsl(Ef) Q) represents the lift coefficient, q, corresponding to the terminal angle of attack in a ballistic simulationfIndicating the terminal dynamic pressure.
According to terminal constraints, have
Figure BDA0002408897380000112
Wherein,
Figure BDA0002408897380000113
indicating the corrected command angle of attack, phif、ψf、Vf、HfRespectively representing terminal latitude, course angle, speed and height obtained by trajectory simulation calculation; q. q.sTAEM、φTAEM、ψTAEM、VTAEM、HTAEMRespectively representing the dynamic pressure value, the latitude value, the course angle value, the speed value and the height value of the terminal constraint.
At phif＝φTAEMAnd psif＝ψTAEMUnder the assumption of (1), linearizing
Figure BDA0002408897380000114
By making equation (51) equal to equation (52), it can be obtained:
Figure BDA0002408897380000115
wherein,Efrepresenting a terminal energy value obtained by ballistic simulation calculation;
Figure BDA0002408897380000116
representing the derivative of the lift coefficient with respect to the angle of attack; therefore, α solved in the formula (53)2The control device can be taken into a formula (19) to be used as a reference attack angle for guidance control.
(2) Second roll reversal point EBR2Correction algorithm
And in the height adjusting stage, guidance is performed by adopting a proportional guidance law. Through the analysis of the guidance algorithm, if E is reducedBR2The guidance law command will generate a larger roll angle to eliminate a larger course error generated during the last roll reversal, so that the longitudinal lift-drag ratio is reduced, and the terminal speed V is reducedf. Thus, VfCan be regarded as relating to EBR2Is a monotonic function of (a). To guarantee the terminal constraints, this problem can be considered as solving a non-linear equation Vf(EBR2)＝VTAEMThe problem of (2). The solution is carried out by adopting a secant method, and V is simulated by a plurality of trajectoriesfMaking a prediction and adjusting E according to the deviation of the predicted value from the terminal constraintBR2The value of (c).
Figure BDA0002408897380000117
EBR2The condition for ending the iterative correction is that the calculated deviation satisfies
Figure BDA0002408897380000121
2. Reference roll angle
After the second roll reversal, the aircraft enters a height adjustment phase. At this stage, the reference angle of attack is still determined by equation (19), and the reference roll angle is controlled by the proportional guidance law. Angular rate of view
Figure BDA0002408897380000122
Can be expressed as
Figure BDA0002408897380000123
Then, the lateral maneuvering acceleration command a obtained by the proportional guidanceL2As follows
Figure BDA0002408897380000124
Wherein k isPNIs a transverse acceleration maneuvering coefficient and takes the value as
Figure BDA0002408897380000125
Figure BDA0002408897380000126
The longitudinal stroke at the time of the second tilt reversal is shown. Assuming that the aircraft is approximately gliding smoothly, the resultant acceleration in the longitudinal plane is 0. Thus, the longitudinal maneuvering acceleration command aL1Is composed of
Figure BDA0002408897380000127
So that the reference roll angle can be obtained as
Figure BDA0002408897380000128
3. Command angle of attack and roll angle
To ensure terminal velocity constraints, the energy-independent residual glide distance reference profile is tracked by adjusting the angle of attack
Figure BDA0002408897380000129
The commanded angle of attack and roll angle may be expressed as follows:
Figure BDA00024088973800001210
Figure BDA00024088973800001211
wherein k isαThe value of (A) can be adjusted in time according to the task, and can be selected as 2 pi/(1.8 × 10)7)；
Figure BDA00024088973800001212
Representing the remaining range, the value can be obtained by the last previous ballistic simulation.
To meet the process constraints, the roll angle still needs to satisfy the following requirements, as before:
|σcmd|≤σmax (61)
through the five steps, the hypersonic reentry guidance method based on the steady gliding trajectory analytic solution can be obtained.
The invention relates to a hypersonic reentry guidance method based on steady glide trajectory analytic solution, which has the advantages that:
(1) the longitudinal and transverse courses are predicted and corrected based on the three-dimensional trajectory analytic solution, so that the calculation amount is greatly reduced, the calculation speed is improved, and the method is convenient to apply to an airborne system in real time.
(2) The transverse movement is determined by twice tilting reversal points, the position of the first reversal point is corrected to control the terminal traverse error, and the position of the second reversal point is corrected to meet the terminal speed constraint. Compared with the traditional method for controlling the transverse motion by using the course error threshold, the method reduces the tilting and reversing times, improves the energy utilization efficiency and reduces the requirement on a control system.
(3) The height adjusting section utilizes the idea of proportional guidance, so that the terminal roll angle is converged to 0, and the target hitting effect is enhanced.
Drawings
Fig. 1 is a schematic view of the equatorial rotation coordinate system (GER) of the earth's center and the eastern coordinate system (NED) of the north and east of the local area.
Fig. 2 is a schematic diagram of an auxiliary earth rotation reference system (AGR) and an auxiliary north east down coordinate system (GNED).
FIG. 3 is a latitude and longitude plot for five target states.
Fig. 4 is a height-energy curve for five target states.
Fig. 5 is a velocity-energy curve for five target states.
Fig. 6 is an angle of attack versus energy curve for five target states.
Fig. 7 is a roll angle versus energy curve for five target states.
Figure 8 is a ballistic dip-energy curve for five target states.
FIG. 9 is a heading error-energy curve for five target states.
In the above figures, the symbols and symbols are as follows:
in the context of figure 1 of the drawings,
Figure BDA0002408897380000131
representing the equatorial rotation coordinate system (GER), oxyzRepresents the local North-East down coordinate system (NED), North represents the North direction, East represents the East direction. h is altitude, R is range, λ is longitude, φ is latitude, V is velocity of the vehicle relative to the earth, γ is ballistic inclination, ψ is vehicle heading angle, σ is roll angle, based on local north. In the context of figure 2, it is shown,
Figure BDA0002408897380000132
represents an assisted earth-rotation reference frame (AGR),
Figure BDA0002408897380000133
representing the north, east and down-the-earth coordinate system (NED), Axis represents the earth polar Axis. M is the current position point of the aircraft, and o is the projection point of the current position point of the aircraft on the earth surface. T is the target point, oTIs the projected point of the target point on the surface of the earth. OmegaexIndicating the rotational angular velocity of the earth
Figure BDA0002408897380000141
Component of the axis, ωeyIndicating the rotational angular velocity of the earth
Figure BDA0002408897380000142
Component of the axis, ωzxIndicating the rotational angular velocity of the earth
Figure BDA0002408897380000143
The component of the axis. In FIG. 3, Longitude represents Longitude, Latitude represents Latitude, and T1-T5 represent examples of five different objects. In fig. 4, Altitude represents height, Energy represents Energy, H represents a height curve obtained by ballistic simulation, and H representsminRepresenting a minimum height boundary that satisfies the process constraints. In fig. 5, Velocity represents altitude and Energy represents Energy. In FIG. 6, Angle of attach represents Angle of Attack and Energy represents Energy. In FIG. 7, Bank Angle represents roll Angle and Energy represents Energy. In FIG. 8, Flight-pathAngle represents roll angle and Energy represents Energy. In fig. 9, Heading Angle represents a Heading Angle, and Energy represents Energy.
Detailed Description
The invention will be further explained in detail with reference to the drawings and the embodiments.
Aiming at the problems of reentry guidance of the hypersonic aircraft, a high-precision hypersonic reentry guidance method based on a three-dimensional trajectory analytic solution is provided. Firstly, an auxiliary geocentric rotation coordinate system is constructed, a new kinetic model of the longitudinal course, the transverse course and the course angle with respect to the energy is established on the basis of the auxiliary geocentric rotation coordinate system, the model is favorable for decoupling and linearizing the highly nonlinear and strongly coupled kinetic model to obtain a three-dimensional ballistic analytic solution, and a guidance method is designed on the basis of the analytic solution. The guidance method utilizes a longitudinal analytic solution to plan a longitudinal reference profile in real time, utilizes a transverse analytic solution to correct a first inversion point to eliminate a transverse error, utilizes trajectory simulation for several times to correct a second inversion point to ensure terminal speed constraint and corrects the minimum value of an attack angle to ensure terminal height constraint.
The invention discloses a hypersonic reentry guidance method based on a steady glide trajectory analytic solution, which is characterized by comprising the following steps: the method comprises the following steps:
the method comprises the following steps: establishing an aircraft dynamics model;
as shown in fig. 1, the geocentric equatorial rotation coordinate system (GER coordinate system) is established: with origin at the groundHeart E, zeThe axis is perpendicular to the earth's equatorial surface, pointing to the north pole; x is the number ofeAxis and yeThe axes are in the equatorial plane and perpendicular to each other. The coordinate system rotates with the earth, which rotates around zeThe angular velocity of rotation of the shaft is the rotational angular velocity omega of the earthe。
North east local coordinate system (NED coordinate system): defining a vertical projection point of an origin o to the ground at the mass center M of the aircraft; the x axis points to the local north, the y axis points to the local east, and the z axis points to the geocentric vertically and downwardly.
Under the background of the rotating earth, a standard atmosphere model and an inverse square gravitational field model are adopted to establish a six-degree-of-freedom kinetic equation of the hypersonic aerocraft as follows:
Figure BDA0002408897380000151
Figure BDA0002408897380000152
Figure BDA0002408897380000153
Figure BDA0002408897380000154
Figure BDA0002408897380000155
Figure BDA0002408897380000156
Figure BDA0002408897380000157
wherein d isWhere/dt represents the derivative over time, t is time, h is altitude, R is range, λ is longitude, φ is latitude, V is the velocity of the aircraft relative to the earth, γ is the ballistic inclination, ψ is the aircraft heading angle, σ is the roll angle, R is the velocity of the aircraft relative to the local north, andeis the radius of the earth, and takes the value of 6378.137km, omegaeIs the angular velocity of rotation of the earth, L ═ ρ V2SCLThe lift acceleration is given by/2 m, and D is equal to rho V2SCDThe/2 m is resistance acceleration; s is the aircraft reference area; cLIs the coefficient of lift; ρ is the air density.
Step two: establishing a constraint model;
the reentry aircraft needs to meet the process constraints of heat flux density, dynamic pressure, overload and the like in the flight process, and the process constraints are expressed as follows:
Figure BDA0002408897380000161
Figure BDA0002408897380000162
Figure BDA0002408897380000163
wherein, KQIs the heat flux density coefficient;
Figure BDA0002408897380000164
representing a maximum allowable heat flow density value; n ismaxRepresents a maximum allowable overload value; q. q.smaxIn addition, due to the capability constraints of the flight control system, the rate of change of the control quantities (angle of attack α and roll angle σ) also satisfies the constraint:
Figure BDA0002408897380000165
in addition to process constraints, the aircraft should also satisfy terminal constraints, with: vTAEM＝2000m/s,HTAEM＝25km,Rtm＝100km,σTAEM≈0°,γTAEM≈0°andΔψTAEM≈0°.
Step three: an initial descent stage guidance method;
according to the characteristics of reentry trajectory, the reentry trajectory is divided into three parts: an initial descent section, a steady glide section and a height adjustment section. The initial descent segment is designed to allow a more gradual transition of the trajectory to a smooth glide segment. The whole flight trajectory is mostly in a stable gliding section, which not only determines the range, speed and height of the flight terminal, but also determines whether the aircraft can meet various process constraints. The height adjustment segment is the last phase of the trajectory, as to whether the trajectory can strictly meet the terminal constraints. The flying speed and flying height of the aircraft are low in this section, and severe requirements are provided for a guidance method for designing a height adjusting section.
In the initial descent, the variation class of the control quantity has little influence on the flight, so that fixed values of the angle of attack and the roll angle are designed here to avoid altitude variations. To decrease the heat flux density and increase the range, the starting point of the trajectory should be set higher. The maximum angle of attack α is therefore chosen heremaxFly 20deg and zero roll angle. Rate of change of ballistic inclination
Figure BDA0002408897380000166
At this time, it means that the lift force can satisfy the requirement of smooth gliding, so that the attack angle can be smoothly transited from the maximum attack angle to the reference attack angle alphabsl。
At this stage, the command angle of attack and roll are designed as follows:
Figure BDA0002408897380000167
σcmd＝0 (12)
wherein,
Δγ＝γ-γSG (13)
Figure BDA0002408897380000171
wherein, γSGRepresenting the angle of inclination, k, of the smooth gliding trajectoryγRepresenting the feedback coefficient and here taking the value 5. When Δ γ is 0, the aircraft enters a smooth glide phase.
Step four: a steady glide phase guidance method;
ballistic planning in the steady glide phase is the basis of the proposed guidance algorithm, which requires approaching the aircraft to the target while guaranteeing process constraints. The algorithm determines the longitudinal range and lateral range profiles based on a new three-dimensional analytical solution. Here, a new ballistic analytical solution is first briefly introduced.
As shown in fig. 2, in order to meet the characteristic that the reentry trajectory has a greater longitudinal range than transverse range, two coordinate systems are first established, and the meanings of the longitudinal range and the transverse range can be visually expressed.
Defining an auxiliary earth-centered rotational reference system (AGR coordinate system) rotating with the earth: the origin is at the geocentric E,
Figure BDA0002408897380000172
the shaft points to the initial position of the aircraft,
Figure BDA0002408897380000173
the axis is in the plane of a great circle passing through the aircraft and the target point and perpendicular to the aircraft
Figure BDA0002408897380000174
The shaft is provided with a plurality of axial holes,
Figure BDA0002408897380000175
the axes may be determined according to the right hand rule.
Meanwhile, a local generalized north east down coordinate system (GNED coordinate system) is defined: the vertical projection point of the origin o from the mass center M of the aircraft to the ground;
Figure BDA0002408897380000176
the axis is vertically downwards directed to the center of the earth,
Figure BDA0002408897380000177
axis-oriented AGR coordinate system
Figure BDA0002408897380000178
The direction of the "north" direction is,
Figure BDA0002408897380000179
the axes are determined by the right hand rule.
Define the longitudinal extent as
Figure BDA00024088973800001710
Define the course as
Figure BDA00024088973800001711
Wherein
Figure BDA00024088973800001712
And
Figure BDA00024088973800001713
respectively representing a generalized longitude and a generalized latitude in an AGR coordinate system.
Defining the energy as
Figure BDA00024088973800001714
The generalized longitude and latitude coordinate system in the AGR coordinate axis is beneficial to linearizing and decoupling a complex dynamic model, and the following three-dimensional ballistic analysis solution is obtained.
Figure BDA0002408897380000181
Figure BDA0002408897380000182
Figure BDA0002408897380000183
Wherein，
Figure BDA0002408897380000184
R*＝Re+H；
Figure BDA0002408897380000185
hij,i＝1,2；j＝0,1,α1And pijI is 1, 2; j is 0,1,2,3,4 are constant coefficients. Since the analytic solution has high precision, the analytic solution can be applied to guidance method derivation to reduce calculation amount and improve universality.
6. Design of reference angle of attack
Because the analytic solution is designed by taking energy as an independent variable, all reference profiles in the guidance method also adopt energy as the independent variable so as to more conveniently utilize the analytic solution. The reference angle of attack profile is designed as:
Figure BDA0002408897380000186
wherein E isα＝-5.55×107J/kg is located near the intersection of the steady glide section and the last height adjustment section. To maximize the capability of the aircraft, the angle of attack is designed to be the angle of attack corresponding to the maximum lift-to-drag ratio, i.e., alpha 110 deg. The quadratic function of the energy of the angle of attack is designed so that the angle of attack can be varied from a1Smoothly transits to alpha2. In this guidance method, α is preset for the first time2The ballistic simulation is later applied to correct it to strictly meet the terminal height constraint, 6 deg. When the design of the attack angle section is finished, the corresponding lift-drag ratio section L/DbslAs determined accordingly.
7. Baseline lift-to-drag ratio profile design
In order to meet the requirement of range, the longitudinal lift-drag ratio profile is designed as follows
Figure BDA0002408897380000191
Wherein, the length range resolution is equal to the remaining range, so as to obtain L1/D1. Design L1/D2＝L/D(ETAEM) In order to make the terminal tilt angle 0. By designing the form of equation (20), the roll angle can be controlled to be approximately constant during flight.
Therefore, to solve for L1/D1First, the remaining range x needs to be calculatedDf. In the GER coordinate system, η is defined as a vector pointing from the center of the earth to the aircraft
Figure BDA0002408897380000192
Vector directed to aircraft from earth center
Figure BDA0002408897380000193
The included angle therebetween. The remaining range can be expressed as:
xDf＝Reη-STAEM (21)
wherein S isTAEMIndicating the distance to the target at the end of the reentry segment.
Solving for L1/D1The process of (a) can be described as follows:
when E is more than or equal to EαWhen there is
xD(Eα,E)+xD(ETAEM,Eα)＝xDf (22)
Wherein,
Figure BDA0002408897380000194
wherein,
Figure BDA0002408897380000195
Figure BDA0002408897380000196
and is provided with a plurality of groups of the following components,
Figure BDA0002408897380000197
wherein,
Figure BDA0002408897380000198
c2＝(a1-c1d0)/d1 (28)
c3＝a0-c2d0 (29)
Figure BDA0002408897380000201
Figure BDA0002408897380000202
Figure BDA0002408897380000203
the formula (22) is solved to obtain
Figure BDA0002408897380000204
Wherein,
Figure BDA0002408897380000205
Figure BDA0002408897380000206
when E < EαIn the meantime, the aircraft is mainly in the altitude adjustment phase and does not update L any more1/D1The value of (c).
8. Reference lateral lift-drag ratio and reference roll angle design
The module value of the transverse lift-drag ratio can be well planned by the previously planned L/DbslAnd a longitudinal lift-to-drag ratio L1/DbsIs determined as
Figure BDA0002408897380000207
In order to satisfy the terminal traverse constraint, the transverse lift-drag ratio profile is designed as follows:
Figure BDA0002408897380000208
wherein E isBR1And EBR2Respectively, the first and second roll reversal points. sgn is a sign variable representing the initial roll reversal direction. After obtaining the transverse lift-drag ratio profile, the roll angle profile is also determined as follows:
Figure BDA0002408897380000211
wherein, L/DrealShowing the actual lift-to-drag ratio profile.
9. Roll reversal point correction algorithm
In this algorithm, the roll reversal point is designed based on a course analytic solution to eliminate dead-end course errors. Here two tilting reversal points are arranged. The first reversal point E is first designed and correctedBR1. First order second reversal point EBR2＝EαThen the terminal traverse xC(ETAEME) is only EBR1As a function of (c). As can be seen from equation (17), the terminal traverse can be represented as EBR1The piecewise function of (c) is as follows:
Figure BDA0002408897380000212
wherein,
Figure BDA0002408897380000213
Figure BDA0002408897380000214
obtaining E from equation (38)BR1Derivative of (2) can be obtained
Figure BDA0002408897380000215
Wherein E ═ EBR1Upper mark for right limit of+Is represented by E ═ EBR1Upper mark for left limit of department-And (4) showing. Because L is2/DbslWhen E is equal to EBR1Is not continuous in place, so there is
Figure BDA0002408897380000216
i is 0,1,2,3, 4. In order to make the terminal traverse error 0, x needs to be satisfiedCf(EBR1) 0. The essence of this algorithm is therefore to solve the non-linear equation xCf(EBR1) Root 0, the newton method is as follows:
Figure BDA0002408897380000217
solution stop set to
Figure BDA0002408897380000218
Thus xCf(EBR1) With EBR1Monotonously changing, the solution of this non-linear equation can be obtained by a few simple iterations.
10. Design of command attack angle and inclination angle
To ensure the tracking accuracy of the reference profile, ballistic damping control techniques are employed herein to suppress long-period, weakly damped oscillations present in the reentry trajectory. The design command attack angle and roll angle are as follows:
αcmd＝αbsl+cosσbslKγ(γSG-γ) (43)
Figure BDA0002408897380000221
wherein, KγIs a feedback coefficient, and takes the value of 5; gamma raySGIs an approximate solution to the angle of inclination of the smooth glide trajectory, which can be expressed as follows
Figure BDA0002408897380000222
The steady glide phase needs to satisfy the process constraint, so the process constraint of equations (8) - (10) is converted into the roll angle constraint, which is expressed as follows
Figure BDA0002408897380000223
Wherein,
Figure BDA0002408897380000224
Figure BDA0002408897380000225
Figure BDA0002408897380000226
wherein HminThe lowest boundary of the corresponding glide height can be obtained through process constraint; constant coefficient kσ-50. Thus, to satisfy process constraints, the roll angle needs to satisfy the following conditions:
|σcmd|≤σmax (50)
step five: a height adjustment section guidance method;
4. and correcting the command attack angle and the second inversion point.
(3) Angle of attack alpha of instruction2Correction algorithm
In order to satisfy the terminal height constraint, a ballistic simulation is adopted to correct the command attack angle alpha2. Terminal height H obtained from trajectory simulationfAnd terminal constraint height HTAEMThe deviation magnitude therebetween, to further correct the final command attack angle alpha2. In the following variables, the subscript f denotes the terminal value calculated from the ballistic simulation, and the subscript TAEM denotes the terminal constraint value.
Due to gammaf,σfAnd
Figure BDA0002408897380000231
can be approximated to 0, and thus can assume cos γ, respectivelyf＝1,cosσ f1 and
Figure BDA0002408897380000232
therefore, according to equation (6):
Figure BDA0002408897380000233
wherein, CL(est)(αbsl(Ef) Q) represents the lift coefficient, q, corresponding to the terminal angle of attack in a ballistic simulationfIndicating the terminal dynamic pressure.
According to terminal constraints, have
Figure BDA0002408897380000234
Wherein,
Figure BDA0002408897380000235
representing the corrected commanded angle of attack.
At phif＝φTAEMAnd psif＝ψTAEMUnder the assumption of (1), linearizing
Figure BDA0002408897380000236
By makingEquation (51) equals equation (52) to yield:
Figure BDA0002408897380000237
therefore, α solved in the formula (53)2The control device can be taken into a formula (19) to be used as a reference attack angle for guidance control.
(4) Second roll reversal point EBR2Correction algorithm
And in the height adjusting stage, guidance is performed by adopting a proportional guidance law. Through the analysis of the guidance algorithm, if E is reducedBR2The guidance law command will generate a larger roll angle to eliminate a larger course error generated during the last roll reversal, so that the longitudinal lift-drag ratio is reduced, and the terminal speed V is reducedf. Thus, VfCan be regarded as relating to EBR2Is a monotonic function of (a). To guarantee the terminal constraints, this problem can be considered as solving a non-linear equation Vf(EBR2)＝VTAEMThe problem of (2). The solution is carried out by adopting a secant method, and V is simulated by a plurality of trajectoriesfMaking a prediction and adjusting E according to the deviation of the predicted value from the terminal constraintBR2The value of (c).
Figure BDA0002408897380000241
EBR2The condition for ending the iterative correction is that the calculated deviation satisfies
Figure BDA0002408897380000242
5. Reference roll angle
After the second roll reversal, the aircraft enters a height adjustment phase. At this stage, the reference angle of attack is still determined by equation (19), and the reference roll angle is controlled by the proportional guidance law. The line-of-sight angular rate may be expressed as
Figure BDA0002408897380000243
Then, the lateral maneuver acceleration command obtained by the proportional steering is as follows
Figure BDA0002408897380000244
Wherein,
Figure BDA0002408897380000245
Figure BDA0002408897380000246
the longitudinal stroke at the time of the second tilt reversal is shown. Assuming that the aircraft is approximately gliding smoothly, the resultant acceleration in the longitudinal plane is 0. Thus, there are
Figure BDA0002408897380000247
So that the reference roll angle can be obtained as
Figure BDA0002408897380000248
6. Command angle of attack and roll angle
To ensure terminal velocity constraints, the energy-independent residual glide distance reference profile is tracked by adjusting the angle of attack
Figure BDA0002408897380000249
The commanded angle of attack and roll angle may be expressed as follows:
Figure BDA00024088973800002410
Figure BDA00024088973800002411
wherein k isαThe value of (A) can be adjusted in time according to the task, and can be selected as 2 pi/(1.8 × 10)7)；
Figure BDA00024088973800002412
The value of (c) can be obtained by the last previous ballistic simulation.
To meet the process constraints, the roll angle still needs to satisfy the following requirements, as before:
|σcmd|≤σmax (61)
through the five steps, the hypersonic reentry guidance method based on the steady gliding trajectory analytic solution can be obtained.
To demonstrate the universality of the guidance method, 5 examples of different target positions are considered here. Five different target points are indicated by T1 to T5, respectively, which are set at (-77 °, 0 °) (T1), (-34 °, -32.5 °) (T2), (5 °, -50 °) (T3), (65 °, -33 °) (T4), and (118 °, 0 °) (T5), respectively. The initial conditions are all h0＝80km,V0＝7200m/s,λ0＝0°,θ045 ° and γ0＝0°。ψ0The pointing direction is determined by the target position. The terminal constraint conditions are all VTAEM＝2000m/s,HTAEM25km and R tm100 km. The simulation results are shown in fig. 3 to 9.
In fig. 3, a red circle represents a range of 100km from the target. As is apparent from the figure, all the trajectories fall on the circles, and the guidance method is proved to meet the terminal distance constraint. Table 1 shows the state quantities of range, speed, altitude, roll angle, trajectory inclination and course error of other terminals.
TABLE 1 simulation results
Figure BDA0002408897380000251
It can be seen from fig. 4 that all but two westward-fly examples satisfy the terminal process constraints. This is because when the aircraft flies westward, the aircraft is subjected to a downward Coriolis force, so that the glide height is reducedThe heat flux density increases. Furthermore, from fig. 4-9 and table 1, we can readily see that all examples satisfy Vf≈VTAEM＝2000m/s,Hf≈HTAEM＝25km,Rtm＝100km,σf≈0°,γf0DEG and [ Delta ] phifTerminal constraint of 0 deg..
In conclusion, the method is deduced through the steps, namely the hypersonic velocity reentry guidance method based on the steady glide trajectory analytic solution, case simulation results show that the method can meet various constraints required by reentry guidance, and the method is designed based on the trajectory analytic solution, has the characteristics of high precision and small calculated amount, and can be applied to an airborne system of an aircraft in real time.

Claims (2)

1. A hypersonic velocity reentry guidance method based on steady glide trajectory analytic solution is characterized in that: the method comprises the following steps:
the method comprises the following steps: establishing an aircraft dynamics model;
under the background of the rotating earth, a standard atmosphere model and an inverse square gravitational field model are adopted to establish a six-degree-of-freedom kinetic equation of the hypersonic aerocraft as follows:
Figure FDA0003267701800000011
Figure FDA0003267701800000012
Figure FDA0003267701800000013
Figure FDA0003267701800000014
Figure FDA0003267701800000015
Figure FDA0003267701800000016
Figure FDA0003267701800000017
where d/dt represents the derivative over time, t is time, h is altitude, R is range, λ is longitude, φ is latitude, V is the velocity of the aircraft relative to the earth, γ is the ballistic inclination, ψ is the aircraft heading angle, based on local north, σ is the roll angle, R is the altitude, and R is the altitude, whereeIs the radius of the earth, and takes the value of 6378.137km, omegaeIs the angular velocity of rotation of the earth, g represents the acceleration of gravity; l ═ ρ V2SCLThe lift acceleration is given by/2 m, and D is equal to rho V2SCD(2 m) is the acceleration of resistance, CL、CDRespectively a lift coefficient and a drag coefficient;
step two: establishing a constraint model;
the reentry vehicle needs to meet the process constraints of heat flux density, dynamic pressure and overload during flight, as follows:
Figure FDA0003267701800000021
Figure FDA0003267701800000022
Figure FDA0003267701800000023
parameter p0Denotes nominal atmospheric density, e denotes natural constant, h denotes altitude, hSRepresenting nominal height, n representing overload, m representing mass;
wherein, KQIs the heat flux density coefficient;
Figure FDA0003267701800000024
representing a maximum allowable heat flow density value; n ismaxRepresents a maximum allowable overload value; q represents an allowable dynamic pressure value; q. q.smaxRepresents the maximum allowable dynamic pressure value;
furthermore, due to the capability constraints of the flight control system, the rate of change of the control quantity also satisfies the constraint:
Figure FDA0003267701800000025
Figure FDA0003267701800000026
in addition to process constraints, the aircraft should also satisfy terminal constraints, with: vTAEM＝2000m/s,HTAEM＝25km,Rtm＝100km,σTAEM≈0°,γTAEM≈0°andΔψTAEM0 °, where subscript TAEM denotes the terminal amount; vTAEM、HTAEM、σTAEM、γTAEM、ΔψTAEMRespectively representing the terminal speed, the terminal height, the terminal roll angle, the terminal trajectory inclination angle and the terminal course angle error; rtmRepresenting the shot-eye distance;
step three: an initial descent stage guidance method;
in the initial descent phase, the maximum angle of attack α is selectedmaxFlying at 20deg and zero roll angle; rate of change of ballistic inclination
Figure FDA0003267701800000028
By time, it is meant that lift meets the smooth gliding requirement, so that the angle of attack transitions smoothly from the maximum angle of attack to the reference angle of attack αbsl；
At this stage, the command angle of attack and the command roll angle are designed as follows:
Figure FDA0003267701800000027
σcmd＝0 (12)
wherein,
Δγ＝γ-γSG (13)
Figure FDA00032677018000000315
wherein, γSGRepresenting the angle of inclination, k, of the smooth gliding trajectoryγRepresents the feedback coefficient and takes the value of 5; when the delta gamma is 0, the aircraft enters a stable gliding stage; parameter alphacmdIs an instruction angle of attack; sigmacmdIs a commanded roll angle; delta gamma is the deviation of the current trajectory inclination angle and the steady glide trajectory inclination angle; delta gamma1As the rate of change of ballistic inclination
Figure FDA0003267701800000031
A deviation value of the trajectory inclination angle when it is 0 from the steady glide trajectory inclination angle;
step four: a steady glide phase guidance method;
ballistic planning in the steady glide phase is the basis of the proposed guidance algorithm, and the aircraft needs to approach the target while process constraints are ensured;
firstly, establishing the following two coordinate systems, and intuitively obtaining the meanings of representing the longitudinal course and the transverse course;
defining an auxiliary earth center rotation reference system (AGR coordinate system) rotating along with the earth: the origin is at the geocentric E,
Figure FDA0003267701800000032
the shaft points to the initial position of the aircraft,
Figure FDA0003267701800000033
the axis is located in a large circular plane passing through the aircraft and the target point and is vertical to the aircraftIn a direction perpendicular to
Figure FDA0003267701800000034
The shaft is provided with a plurality of axial holes,
Figure FDA0003267701800000035
the axis is determined according to the right-hand rule;
meanwhile, local generalized north, east and down coordinate systems, i.e. GNED coordinate systems, are defined: the vertical projection point of the origin o from the mass center M of the aircraft to the ground;
Figure FDA0003267701800000036
the axis is vertically downwards directed to the center of the earth,
Figure FDA0003267701800000037
axis-oriented AGR coordinate system
Figure FDA0003267701800000038
The direction of the "north" direction is,
Figure FDA0003267701800000039
the axis is determined by the right hand rule;
define the longitudinal extent as
Figure FDA00032677018000000310
Define the course as
Figure FDA00032677018000000311
Wherein
Figure FDA00032677018000000312
And
Figure FDA00032677018000000313
respectively representing a generalized longitude and a generalized latitude in an AGR coordinate system;
defining the energy as
Figure FDA00032677018000000314
Defining the longitudinal lift-to-drag ratio Lcos sigma/D ═ L1D, transverse lift-drag ratio Lsin sigma/D ═ L2D; the generalized longitude and latitude coordinate system in the AGR coordinate axis is beneficial to linearizing and decoupling a complex dynamic model, and the following three-dimensional ballistic analysis solution is obtained;
Figure FDA0003267701800000041
Figure FDA0003267701800000042
Figure FDA0003267701800000043
wherein,
Figure FDA0003267701800000044
Figure FDA0003267701800000045
and pijI is 1, 2; j is 0,1,2,3 and 4 are constant coefficients; e, E0Respectively representing the current energy and the initial energy value; x is the number ofERepresenting an energy integral variable; l is1Representing lift values in a longitudinal plane; x is the number ofC0Representing an initial traverse value;
Figure FDA0003267701800000046
representing an initial course angle error under a GNED coordinate system;
Figure FDA0003267701800000047
representing the course angle error under the GNED coordinate system;
4.1 datum Angle of attack design
All reference profiles in the guidance method also adopt energy as independent variables; reference angle of attack alphabslThe design is as follows:
Figure FDA0003267701800000048
wherein E isα＝-5.55×107J/kg is positioned near the junction point of the stable gliding section and the last height adjusting section; alpha is alpha1And alpha2Two reference angle of attack parameters; to maximize the capacity of the aircraft, design alpha1For maximum lift-to-drag ratio corresponding to angle of attack, i.e. alpha110 deg; the quadratic function of the energy of the angle of attack is designed so that the angle of attack is from alpha1Smoothly transits to alpha2(ii) a In this guidance method, α is preset for the first time2If the terminal height is 6deg, trajectory simulation is applied to correct the terminal height to strictly meet the terminal height constraint; when the design of the attack angle section is finished, the corresponding reference lift-drag ratio section (L/D)bslIs also determined accordingly;
4.2 benchmark lift-drag ratio profile design
In order to meet the requirement of the range, a standard lift-drag ratio profile is designed as follows:
Figure FDA0003267701800000051
wherein, let the terminal longitudinal analysis solve xDfEqual to the remaining range, i.e. obtaining L1/D1(ii) a Design L1/D2＝L/D(ETAEM) In which ETAEMIs terminal energy, so when xE＝ETAEMWhen the temperature of the water is higher than the set temperature,
(L1/D)bsl|E＝E_TAEM＝(LcosσTAEM/D)bsl＝L/D(ETAEM) So that the terminal is inclined at an angle σTAEMIs 0; by means of the form of the design formula (20), the terminal roll angle is controlled to be approximately constant in the flying process;
therefore, to solve for L1/D1First, the remaining range x needs to be calculatedDf(ii) a In the GER coordinate system, η is defined as a vector pointing from the center of the earth to the aircraft
Figure FDA0003267701800000052
Vector directed to aircraft from earth center
Figure FDA0003267701800000053
The included angle between them; the remaining range is then expressed as:
xDf＝Reη-STAEM (21)
wherein S isTAEMIndicating the distance to the target at the end of the reentry segment;
solving for L1/D1The process of (a) is described as follows:
when E is more than or equal to EαWhen there is
xD(Eα,E)+xD(ETAEM,Eα)＝xDf (22)
Wherein,
Figure FDA0003267701800000054
wherein,
Figure FDA0003267701800000055
Figure FDA0003267701800000056
and is provided with a plurality of groups of the following components,
Figure FDA0003267701800000057
wherein,
Figure FDA0003267701800000061
c2＝(a1-c1d0)/d1 (28)
c3＝a0-c2d0 (29)
Figure FDA0003267701800000062
Figure FDA0003267701800000063
Figure FDA0003267701800000064
the formula (22) is solved to obtain
Figure FDA0003267701800000065
Wherein,
Figure FDA0003267701800000066
Figure FDA0003267701800000067
when E < EαWhen the aircraft is in the altitude adjustment stage, L is not updated any more1/D1A value of (d);
4.3 reference transverse lift-drag ratio and reference roll angle design
The module value of the reference transverse lift-drag ratio is well planned from the previous (L/D)bslAnd longitudinal lift-to-drag ratio (L)1/D)bslIs determined as
Figure FDA0003267701800000068
In order to meet the terminal traverse constraint, the design reference transverse lift-drag ratio is as follows:
Figure FDA0003267701800000071
wherein E isBR1And EBR2Respectively representing a first and a second roll reversal point; sgn is a sign variable representing the initial roll reversal direction; after obtaining the reference transverse lift-drag ratio, the reference roll angle sigmabslIt was also determined that this is expressed as follows:
Figure FDA0003267701800000072
wherein, (L/D)realRepresenting the actual lift-to-drag ratio profile;
4.4 roll reversal Point correction Algorithm
Designing a roll reversal point based on a traverse analytic solution to eliminate a terminal traverse error; arranging two tilting reversal points; the first reversal point E is first designed and correctedBR1(ii) a First order second reversal point EBR2＝EαThen the terminal traverse xC(ETAEME) is only EBR1A function of (a); as seen from the formula (17), the terminal traverse xCfIs represented by EBR1The piecewise function of (c) is as follows:
Figure FDA0003267701800000073
f, G is a custom function, and the expression is as follows
Figure FDA0003267701800000074
Figure FDA0003267701800000075
Wherein, the parameter xE2,xE1Expressing upper and lower energy bounds; k is a constant of 1 to 5; p1Representing a polynomial fitting coefficient;
obtaining E from equation (38)BR1Derivative of (2) can be obtained
Figure FDA0003267701800000076
Wherein E ═ EBR1Upper mark for right limit of+Is represented by E ═ EBR1Upper mark for left limit of department-Represents; because L is2/DbslWhen E is equal to EBR1Is not continuous in place, so there is
Figure FDA0003267701800000077
Figure FDA0003267701800000078
In order to make the terminal traverse error 0, x needs to be satisfiedCf(EBR1) 0; the essence of this algorithm is therefore to solve the non-linear equation xCf(EBR1) Root 0, the newton method is as follows:
Figure FDA0003267701800000081
parameter(s)
Figure FDA0003267701800000082
And
Figure FDA0003267701800000083
and not the same as the above-mentioned case,
Figure FDA0003267701800000084
represents EBR1To the power of k in the first order,
Figure FDA0003267701800000085
e representing the kth iterationBR1A value; solution stop set to
Figure FDA0003267701800000086
Thus xCf(EBR1) With EBR1Monotonously changing, and the solution of the nonlinear equation is obtained by simple iterations;
4.5 Command attack Angle and Command roll Angle design
In order to ensure the tracking precision of the reference profile, a trajectory damping control technology is adopted to inhibit the oscillation of long-period and weak damping existing in a reentry trajectory; design command angle of attack alphacmdAnd the command roll angle sigmacmdThe following were used:
αcmd＝αbsl+cosσbslKγ(γSG-γ) (43)
Figure FDA0003267701800000087
wherein, KγIs a feedback coefficient, and takes the value of 5; gamma raySGIs an approximate solution to the angle of inclination of the smooth glide trajectory, as shown below
Figure FDA0003267701800000088
Wherein S is the aircraft reference area; cLIs the coefficient of lift; ρ is the air density;
the smooth glide phase needs to satisfy the process constraint, so the process constraint of equation (8) -10 is converted into the command roll angle constraint, which is expressed as follows
Figure FDA0003267701800000089
Wherein L ismaxFor maximum lift, the expression is as follows:
Figure FDA00032677018000000810
Figure FDA00032677018000000811
Figure FDA0003267701800000091
wherein HminThe lowest boundary corresponding to the glide height is obtained through process constraint; constant coefficient kσ-50; thus, to satisfy the process constraints, the following conditions need to be satisfied for the commanded roll angle:
|σcmd|≤σmax (50)
step five: a height adjustment section guidance method;
5.1 correcting the command attack angle and the second reversal point;
(1) reference angle of attack parameter alpha2Correction algorithm
In order to meet the height constraint of the terminal, a primary trajectory simulation is adopted to correct a reference attack angle parameter alpha2(ii) a Terminal height H obtained from trajectory simulationfAnd terminal constraint height HTAEMThe deviation between the two parameters is further corrected to obtain the final reference attack angle parameter alpha2(ii) a In the following variables, subscript f represents a terminal value obtained by ballistic simulation calculation, and subscript TAEM represents a terminal constraint value;
due to gammaf,σfAnd
Figure FDA0003267701800000092
is approximately 0, so cos γ is assumed separatelyf＝1,cosσf1 and
Figure FDA0003267701800000093
therefore, it is obtained from equation (6):
Figure FDA0003267701800000094
wherein, CL(est)(αbsl(Ef) Q) represents the lift coefficient, q, corresponding to the terminal angle of attack in a ballistic simulationfRepresenting terminal dynamic pressure; according to terminal constraints, have
Figure FDA0003267701800000095
Wherein,
Figure FDA0003267701800000096
indicating the corrected command angle of attack, phif、ψf、Vf、HfRespectively representing terminal latitude, course angle, speed and height obtained by trajectory simulation calculation; q. q.sTAEM、φTAEM、ψTAEM、VTAEM、HTAEMRespectively representing a dynamic pressure value, a latitude value, a course angle value, a speed value and a height value of terminal constraint; at phif＝φTAEMAnd psif＝ψTAEMUnder the assumption of (1), linearizing
Figure FDA0003267701800000097
By making equation (51) equal to equation (52), we obtain:
Figure FDA0003267701800000101
wherein E isfRepresenting a terminal energy value obtained by ballistic simulation calculation;
Figure FDA0003267701800000102
representing the derivative of the lift coefficient with respect to the angle of attack; therefore, α solved in the formula (53)2Taking the belt-in type (19) as a reference attack angle to conduct guidance control;
(2) second roll reversal point EBR2Correction algorithm
In the height adjustment stage, guidance is designed by adopting a proportional guidance law; through the analysis of the guidance algorithm, if E is reducedBR2The guidance law command will generate a larger roll angle to eliminate a larger course error generated during the last roll reversal, so that the longitudinal lift-drag ratio is reduced, and the terminal speed V is reducedf(ii) a Thus, VfIs regarded as relating to EBR2A monotonic function of (a); to guarantee the terminal constraints, this problem is treated as solving the nonlinear equation Vf(EBR2)＝VTAEMThe problem of the solution of (a); solving by adopting a secant method, and simulating V by a plurality of trajectoriesfMaking a prediction and adjusting E according to the deviation of the predicted value from the terminal constraintBR2A value of (d);
Figure FDA0003267701800000103
EBR2the condition for ending the iterative correction is that the calculated deviation satisfies
Figure FDA0003267701800000104
5.2 reference roll angle
After the second roll reversal, the aircraft enters a height adjustment phase; at this stage, the reference angle of attack is still determined by equation (19), and the reference roll angle is controlled by the proportional guidance law; angular rate of view
Figure FDA0003267701800000105
Is shown as
Figure FDA0003267701800000106
Then, the lateral maneuvering acceleration command a obtained by the proportional guidanceL2As follows
Figure FDA0003267701800000107
Wherein k isPNIs a transverse acceleration maneuvering coefficient and takes the value as
Figure FDA0003267701800000108
Figure FDA0003267701800000109
Represents the longitudinal stroke at the time of the second tilt reversal; assuming that the aircraft glides approximately smoothly, the resultant acceleration in the longitudinal plane is 0; thus, the longitudinal maneuvering acceleration command aL1Is composed of
Figure FDA0003267701800000111
So that the reference roll angle is
Figure FDA0003267701800000112
5.3 Command Angle of attack and Command roll
To ensure terminal velocity constraints, the energy-independent residual glide distance reference profile is tracked by adjusting the angle of attack
Figure FDA0003267701800000113
The commanded angle of attack and commanded roll angle are expressed as follows:
Figure FDA0003267701800000114
Figure FDA0003267701800000115
wherein k isαThe value of (2) is selected as 2 pi/(1.8 × 10) according to task timely adjustment7)；
Figure FDA0003267701800000116
Representing the remaining range, obtained by the last trajectory simulation before;
to meet the process constraints, as before, the commanded roll angle still needs to meet the following requirements:
|σcmd|≤σmax (61)
through the five steps, the hypersonic reentry guidance method based on the steady glide trajectory analytic solution is obtained.
2. The hypersonic re-entry guidance method based on the smooth gliding trajectory analytic solution as claimed in claim 1, characterized in that: the hypersonic flight vehicle is an aircraft with a flight Mach number larger than 5 and capable of realizing hypersonic flight in an atmosphere layer and a trans-atmosphere layer.
CN202010170138.1A 2020-03-12 2020-03-12 Hypersonic velocity reentry guidance method based on steady glide trajectory analytic solution Active CN111306989B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010170138.1A CN111306989B (en) 2020-03-12 2020-03-12 Hypersonic velocity reentry guidance method based on steady glide trajectory analytic solution

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010170138.1A CN111306989B (en) 2020-03-12 2020-03-12 Hypersonic velocity reentry guidance method based on steady glide trajectory analytic solution

Publications (2)

Publication Number Publication Date
CN111306989A CN111306989A (en) 2020-06-19
CN111306989B true CN111306989B (en) 2021-12-28

Family

ID=71149646

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010170138.1A Active CN111306989B (en) 2020-03-12 2020-03-12 Hypersonic velocity reentry guidance method based on steady glide trajectory analytic solution

Country Status (1)

Country Link
CN (1) CN111306989B (en)

Families Citing this family (15)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111674573B (en) * 2020-08-11 2020-11-10 北京控制与电子技术研究所 Non-parallel gravitational field deep space impact control method and system based on proportional guidance
CN112009698B (en) * 2020-08-28 2021-11-23 北京空间技术研制试验中心 Return flight energy management method for hypersonic cruise aircraft
CN112162567B (en) * 2020-09-09 2022-05-10 北京航空航天大学 Avoidance guidance method suitable for online no-fly zone of aircraft
CN112507466B (en) * 2020-12-21 2023-05-09 中国人民解放军96901部队23分队 Method for calculating glide trajectory energy-dependent degradation solution considering earth rotation rate
CN112861249B (en) * 2020-12-21 2023-03-31 中国人民解放军96901部队23分队 Method for calculating degradation solution of gliding trajectory along with speed change based on attack angle and resistance
CN112507467B (en) * 2020-12-21 2023-05-12 中国人民解放军96901部队23分队 Method for calculating descending order solution of glide trajectory along with speed change based on resistance and lift-drag ratio
CN112861250B (en) * 2020-12-21 2023-03-28 中国人民解放军96901部队23分队 Method for calculating degradation solution of gliding trajectory along with energy change based on attack angle and resistance
CN112507465B (en) * 2020-12-21 2023-04-14 中国人民解放军96901部队23分队 Gliding trajectory energy change degradation solution calculation method based on resistance and lift-drag ratio
CN113375634B (en) * 2021-04-30 2022-10-14 北京临近空间飞行器系统工程研究所 Altitude measurement method based on atmospheric model and aircraft normal overload combination
CN113741509B (en) * 2021-07-30 2024-02-27 北京航空航天大学 Hypersonic gliding aircraft hold-down section energy management method
CN114265433B (en) * 2021-12-24 2024-02-02 北京航空航天大学 Transverse guidance method and system for meeting large-journey maneuver
CN114610057B (en) * 2022-02-25 2024-06-28 北京航空航天大学 Design method for maneuver burst prevention strategy of high Mach aircraft
CN114675545B (en) * 2022-05-26 2022-08-23 中国人民解放军火箭军工程大学 Hypersonic aircraft reentry cooperative guidance method based on reinforcement learning
CN114936531B (en) * 2022-05-31 2024-09-20 西北工业大学 Aircraft hitting area prediction system and method based on fuzzy neural network
CN117724335B (en) * 2023-12-18 2024-07-19 苏州星幕航天科技有限公司 Analysis algorithm of amplitude change relation of inertial coordinate system ZEM guidance state

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8489258B2 (en) * 2009-03-27 2013-07-16 The Charles Stark Draper Laboratory, Inc. Propulsive guidance for atmospheric skip entry trajectories
CN104035335A (en) * 2014-05-27 2014-09-10 北京航空航天大学 High accuracy longitudinal and cross range analytical prediction method based smooth gliding reentry guidance method
CN104656450A (en) * 2015-01-20 2015-05-27 北京航空航天大学 Method for designing hypersonic aircraft to re-enter into trajectory in smooth gliding mode
CN105550402A (en) * 2015-12-07 2016-05-04 北京航空航天大学 Attack angle or inclination angle frequency conversion based design method for hypersonic steady maneuver gliding trajectory
CN108036676A (en) * 2017-12-04 2018-05-15 北京航空航天大学 A kind of autonomous reentry guidance method of full directive based on three-dimensional resolution Value of Reentry Vehicle

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8489258B2 (en) * 2009-03-27 2013-07-16 The Charles Stark Draper Laboratory, Inc. Propulsive guidance for atmospheric skip entry trajectories
CN104035335A (en) * 2014-05-27 2014-09-10 北京航空航天大学 High accuracy longitudinal and cross range analytical prediction method based smooth gliding reentry guidance method
CN104656450A (en) * 2015-01-20 2015-05-27 北京航空航天大学 Method for designing hypersonic aircraft to re-enter into trajectory in smooth gliding mode
CN105550402A (en) * 2015-12-07 2016-05-04 北京航空航天大学 Attack angle or inclination angle frequency conversion based design method for hypersonic steady maneuver gliding trajectory
CN108036676A (en) * 2017-12-04 2018-05-15 北京航空航天大学 A kind of autonomous reentry guidance method of full directive based on three-dimensional resolution Value of Reentry Vehicle

Also Published As

Publication number Publication date
CN111306989A (en) 2020-06-19

Similar Documents

Publication Publication Date Title
CN111306989B (en) Hypersonic velocity reentry guidance method based on steady glide trajectory analytic solution
CN109508030B (en) Collaborative analysis reentry guidance method considering multi-forbidden flight area constraint
CN108036676B (en) A kind of autonomous reentry guidance method of full directive based on three-dimensional resolution Value of Reentry Vehicle
CN104035335B (en) Steady glide reentry guidance method based on the longitudinal and transverse journey analytical Prediction method of high accuracy
CN107941087B (en) A kind of superb steady gliding reentry guidance method of high lift-drag ratio based on resistance profiles
CN109740198B (en) Analytic prediction-based three-dimensional reentry guidance method for gliding aircraft
Xie et al. Highly constrained entry trajectory generation
Li et al. Improved artificial potential field based lateral entry guidance for waypoints passage and no-fly zones avoidance
Yu et al. Analytical entry guidance for coordinated flight with multiple no-fly-zone constraints
Zhang et al. Entry trajectory planning based on three-dimensional acceleration profile guidance
CN110015446B (en) Semi-analytic Mars entry guidance method
CN109062241B (en) Autonomous full-shot reentry guidance method based on linear pseudo-spectrum model predictive control
Pan et al. 3D guidance for hypersonic reentry gliders based on analytical prediction
CN106444822A (en) Space vector field guidance based stratospheric airship&#39;s trajectory tracking control method
CN113900448B (en) Aircraft prediction correction composite guidance method based on sliding mode interference observer
CN105022403B (en) The defining method of longitudinal TRAJECTORY CONTROL gain of glide vehicle
CN115657711B (en) Hypersonic aircraft avoidance guidance method based on deep reinforcement learning
CN105718660B (en) The a wide range of Maneuver Ballistic Trajectory three-dimensional envelope calculation method of near space
CN105354380B (en) The glide trajectories rapid correction method of compensation is influenced towards perturbation factors
CN113835442B (en) Hypersonic gliding aircraft linear pseudo-spectrum reentry guidance method and hypersonic gliding aircraft linear pseudo-spectrum reentry guidance system
CN113093789B (en) Planning method for avoiding trajectory of aircraft no-fly zone based on path point optimization
CN113093776A (en) Method and device for determining off-orbit parameters of spacecraft
Yomchinda et al. Autonomous control and path planning for autorotation of unmanned helicopters
CN115924127A (en) Ground-hitting spacecraft space-atmosphere inner orbit full-section control method
CN113759966B (en) Terminal guidance method with controllable terminal speed in three-dimensional space

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant