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
value
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

    • FMECHANICAL ENGINEERING; LIGHTING; HEATING; WEAPONS; BLASTING
    • F41WEAPONS
    • F41GWEAPON SIGHTS; AIMING
    • F41G3/00Aiming or laying means
    • F41G3/22Aiming 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

本发明一种基于平稳滑翔弹道解析解的高超声速再入制导方法,它包括以下几个步骤:1、建立飞行器动力学模型;2、建立约束模型;3、初始下降段制导方法;4、平稳滑翔段制导方法;5、高度调整段制导方法;优点在于:基于三维弹道解析解进行纵程和横程的预测矫正,极大减少了计算量,提高了计算速度,便于实时应用到机载系统上。横向运动由两次倾侧反转点确定,修正第一次反转点位置来控制终端横程误差,修正第二次反转点位置来满足终端速度约束。相比于传统的利用航向误差门限控制横向运动方法,减少了倾侧反转次数,提高了能量利用效率,减低了对控制系统的需求。高度调整段利用比例导引的思想,使得终端倾侧角收敛到0,增强了目标打击效果。

Figure 202010170138

The present invention is a hypersonic reentry guidance method based on the analytical solution of the smooth gliding trajectory, which comprises the following steps: 1. Establishing a dynamic model of the aircraft; 2. Establishing a constraint model; 3. Guidance method for the initial descent stage; 4. Smoothing Guidance method for gliding segment; 5. Guidance method for height adjustment segment; the advantages are: prediction and correction of longitudinal and transverse distance based on 3D ballistic analytical solution, which greatly reduces the amount of calculation, improves the calculation speed, and is convenient for real-time application to airborne systems superior. The lateral movement is determined by two tilt 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 satisfy the terminal speed constraint. Compared with the traditional method of using the heading error threshold to control the lateral motion, the number of inversions is reduced, the energy utilization efficiency is improved, and the demand for the control system is reduced. The height adjustment section uses the idea of proportional guidance to make the terminal tilt angle converge to 0, which enhances the target strike effect.

Figure 202010170138

Description

一种基于平稳滑翔弹道解析解的高超声速再入制导方法A Hypersonic Reentry Guidance Method Based on the Analytical Solution of Smooth Glide Trajectory

技术领域technical field

本发明提供了一种基于平稳滑翔弹道解析解的高超声速再入制导方法,属于航天技术、武器技术领域The invention provides a hypersonic reentry guidance method based on an analytical solution of a smooth gliding trajectory, belonging to the fields of aerospace technology and weapon technology

背景技术Background technique

再入制导律的设计是高超声速再入滑翔飞行器能否完成精确打击任务的核心技术。再入制导的最终目的是在满足路径约束和终端约束的情况下,实现将高超声速再入飞行器成功导引到指定的目标。The design of the reentry guidance law is the core technology for whether the hypersonic reentry glide vehicle can complete the precision strike mission. The ultimate goal of reentry guidance is to successfully guide the hypersonic reentry vehicle to the designated target under the condition of satisfying the path constraints and terminal constraints.

再入飞行器在再入过程中受到严格的路径约束,然而传统算法对路径约束的处理一直比较困难,给轨迹的在线生成和制导律的在线计算带来了困难和挑战。同时,由于极快的滑翔飞行速度对计算时间提出了更高的要求,使得在线制导律指令生成的时间约束很强,对于目标重新定位的任务,要求算法具有较强的自适应性和快速性,以适应大机动飞行任务的要求。The re-entry vehicle is subject to strict path constraints during the re-entry process. However, it has been difficult for traditional algorithms to deal with the path constraints, which brings difficulties and challenges to the online generation of the trajectory and the online calculation of the guidance law. At the same time, because the extremely fast gliding flight speed puts forward higher requirements for calculation time, the time constraint of online guidance law command generation is very strong. For the task of target repositioning, the algorithm requires strong adaptability and rapidity. , to meet the requirements of large maneuvering missions.

发明内容SUMMARY OF THE INVENTION

针对高超声速飞行器再入制导上述问题,提出了一种基于三维弹道解析解的高精度高超滑翔再入制导方法。首先构建辅助地心旋转坐标系,并基于此建立了一个新的纵程、横程和航向角关于能量的动力学模型,该模型有利于对高度非线性、强耦合的动力学模型进行解耦和线化,得到三维弹道解析解,并基于此解析解设计了制导方法。制导方法利用纵程解析解实时规划纵向参考剖面,利用横程解析解修正第一次反转点以消除横程误差,利用几次弹道仿真来修正第二次反转点来保证终端速度约束、修正攻角最小值保证终端高度约束。Aiming at the above problems of hypersonic vehicle reentry guidance, a high-precision hyperglide reentry guidance method based on 3D ballistic analytical solution was proposed. Firstly, an auxiliary geocentric rotation coordinate system is constructed, and based on this, a new dynamic model of longitudinal, transverse and heading angles with respect to energy is established, which is beneficial to decoupling highly nonlinear and strongly coupled dynamic models. And linearization, the three-dimensional ballistic analytical solution is obtained, and the guidance method is designed based on this analytical solution. The guidance method uses the longitudinal analytical solution to plan the longitudinal reference profile in real time, uses the lateral analytical solution to correct the first reversal point to eliminate the lateral error, and uses several ballistic simulations to correct the second inversion point to ensure the terminal speed constraint, Corrected minimum angle of attack to ensure terminal height constraints.

本发明一种基于平稳滑翔弹道解析解的高超声速再入制导方法,其特征在于:它包括以下几个步骤:The present invention is a hypersonic reentry guidance method based on the analytical solution of the smooth gliding trajectory, characterized in that it comprises the following steps:

步骤一:建立飞行器动力学模型;Step 1: Build the aircraft dynamics model;

高超声速飞行器一般是指飞行马赫数大于5,在大气层和跨大气层中实现高超声速飞行的飞行器。旋转地球背景下,采用标准大气模型和平方反比引力场模型,建立高超声速飞行器的六自由度动力学方程如下:A hypersonic vehicle generally refers to an aircraft with a flight Mach number greater than 5 that achieves hypersonic flight in the atmosphere and across the atmosphere. Under the background of the rotating earth, using the standard atmospheric model and the inverse square gravitational field model, the six-degree-of-freedom dynamic equation of the hypersonic vehicle is established as follows:

Figure BDA0002408897380000021
Figure BDA0002408897380000021

Figure BDA0002408897380000022
Figure BDA0002408897380000022

Figure BDA0002408897380000023
Figure BDA0002408897380000023

Figure BDA0002408897380000024
Figure BDA0002408897380000024

Figure BDA0002408897380000025
Figure BDA0002408897380000025

Figure BDA0002408897380000026
Figure BDA0002408897380000026

Figure BDA0002408897380000027
Figure BDA0002408897380000027

其中,d/dt表示对时间求导,t是时间,h是海拔高度,R是射程,λ是经度,φ是纬度,V是飞行器相对于地球的速率,γ是弹道倾角,ψ是飞行器航向角,以当地北向为基准,σ是倾侧角,Re是地球半径,取值6378.137km,ωe是地球自转角速度,g表示重力加速度;L=ρV2SCL/2m为升力加速度,D=ρV2SCD/2m为阻力加速度,CL、CD分别为升力系数与阻力系数。where d/dt is the derivation of 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, and ψ is the heading of the aircraft Angle, based on the local north direction, σ is the tilt angle, Re is the radius of the earth, taking the value 6378.137km, ω e is the angular velocity of the earth's rotation, g is the acceleration of gravity; L=ρV 2 SC L /2m is the acceleration of lift, D= ρV 2 SC D /2m is the drag acceleration, and CL and CD are the lift coefficient and drag coefficient, respectively .

步骤二:建立约束模型;Step 2: Establish a constraint model;

再入飞行器在飞行过程中需要满足热流密度、动压和过载等过程约束,表示如下:The re-entry vehicle needs to meet the process constraints such as heat flux density, dynamic pressure and overload during flight, which are expressed as follows:

Figure BDA0002408897380000028
Figure BDA0002408897380000028

Figure BDA0002408897380000031
Figure BDA0002408897380000031

Figure BDA0002408897380000032
Figure BDA0002408897380000032

参数ρ0表示标称大气密度、e表示自然常数、h表示海拔高度、hS表示标称高度,n表示过载,m表示质量。The parameter ρ 0 represents the nominal atmospheric density, e the natural constant, h the altitude above sea level, h S the nominal altitude, n the overload, and m the mass.

其中,KQ是热流密度系数;

Figure BDA0002408897380000033
表示最大允许热流密度值;n表示允许过载值;nmax表示最大允许过载值;q表示允许动压值;qmax表示最大允许动压值。where K Q is the heat flux density coefficient;
Figure BDA0002408897380000033
Indicates the maximum allowable heat flux value; n indicates the allowable overload value; n max indicates the maximum allowable overload value; q indicates the allowable dynamic pressure value; q max indicates the maximum allowable dynamic pressure value.

此外,由于飞行控制系统的能力制约,控制量(攻角α和倾侧角σ)的变化率也满足约束:

Figure BDA0002408897380000034
In addition, due to the capability constraints of the flight control system, the rate of change of the control variables (angle of attack α and angle of roll σ) also satisfies the constraints:
Figure BDA0002408897380000034

除了过程约束,飞行器还应满足终端约束,有:VTAEM=2000m/s,HTAEM=25km,Rtm=100km,σTAEM≈0°,γTAEM≈0°andΔψTAEM≈0°.其中下标TAEM表示终端量。VTAEM、HTAEM、σTAEM、γTAEM、ΔψTAEM分别表示终端速度、终端高度、终端倾侧角、终端弹道倾角、终端航向角误差。Rtm表示弹目距离。In addition to the process constraints, the aircraft should also satisfy the terminal constraints, as follows: V TAEM = 2000m/s, H TAEM = 25km, R tm = 100km, σ TAEM ≈0°, γ TAEM ≈0° and Δψ TAEM ≈0°. The subscripts TAEM stands for terminal volume. V TAEM , H TAEM , σ TAEM , γ TAEM , and Δψ TAEM represent terminal velocity, terminal altitude, terminal roll angle, terminal ballistic inclination angle, and terminal heading angle error, respectively. R tm represents the projectile distance.

步骤三:初始下降段制导方法;Step 3: Guidance method for the initial descent segment;

根据再入弹道的特性,将其划分为三部分:初始下降段,平稳滑翔段和高度调整段。设计初始下降段是为了让弹道更平缓地过渡到的平稳滑翔段。全程飞行弹道绝大部分均处于平稳滑翔段,它不但决定了飞行终端射程、速度和高度,还决定了飞行器是否能满足多种过程约束。高度调整段是弹道的最后一个阶段,关系到此弹道能否严格满足终端约束。飞行器在此段飞行速度和飞行高度都较低,对设计高度调整段的制导方法提出了较为严苛的要求。According to the characteristics of the re-entry trajectory, it is divided into three parts: the initial descent section, the smooth glide section and the height adjustment section. The initial descent segment is designed to allow for a smoother glide segment into which the trajectory transitions more smoothly. Most of the whole flight trajectory is in the smooth gliding segment, which not only determines the range, speed and altitude of the flight terminal, but also determines whether the aircraft can meet various process constraints. The height adjustment segment is the last stage of the trajectory, which is related to whether the trajectory can strictly meet the terminal constraints. The flight speed and flight altitude of the aircraft in this section are relatively low, which puts forward more stringent requirements for the guidance method of the design height adjustment section.

在初始下降段,控制量的变化类对飞行影响不大,所以这里设计固定大小的攻角和倾侧角值以避免高度变化。为了降低热流密度,增加射程,弹道的起始点应该设置的高一些。因此这里选择最大攻角αmax=20deg和零倾侧角进行飞行。当弹道倾角变化率

Figure BDA0002408897380000035
时,意味着升力可以满足平稳滑翔需求,所以此时攻角可以平缓地从最大攻角过渡到参考攻角剖面αbsl。In the initial descent stage, the variation of the control quantity has little effect on the flight, so the angle of attack and the angle of inclination are designed with a fixed size to avoid the altitude change. In order to reduce the heat flux density and increase the range, the starting point of the trajectory should be set higher. Therefore, the maximum angle of attack α max = 20deg and zero roll angle are selected for flight here. When the rate of change of ballistic inclination
Figure BDA0002408897380000035
, which means that the lift can meet the requirement of smooth gliding, so the angle of attack can smoothly transition from the maximum angle of attack to the reference angle of attack profile α bsl .

在这个阶段,指令攻角和倾侧角设计如下:At this stage, the commanded attack and roll angles are designed as follows:

Figure BDA0002408897380000041
Figure BDA0002408897380000041

σcmd=0 (12)σ cmd = 0 (12)

其中,in,

Δγ=γ-γSG (13)Δγ=γ-γ SG (13)

Figure BDA0002408897380000042
Figure BDA0002408897380000042

其中,γSG表示平稳滑翔弹道倾角,kγ表示反馈系数且这里取值为5。当Δγ=0时,飞行器进入平稳滑翔阶段。参数αcmd为指令攻角;σcmd为指令倾侧角;Δγ为当前弹道倾角与平稳滑翔弹道倾角的偏差;Δγ1为当弹道倾角变化率

Figure BDA0002408897380000043
为0的时候的弹道倾角与平稳滑翔弹道倾角的偏差值。Among them, γ SG represents the inclination angle of the smooth gliding trajectory, and k γ represents the feedback coefficient and the value here is 5. When Δγ=0, the aircraft enters the smooth gliding phase. The parameter α cmd is the commanded angle of attack; σ cmd is the commanded inclination angle; Δγ is the deviation between the current ballistic inclination and the stable gliding ballistic inclination; Δγ 1 is the change rate of the current ballistic inclination
Figure BDA0002408897380000043
When it is 0, the deviation between the ballistic inclination and the smooth gliding ballistic inclination.

步骤四:平稳滑翔段制导方法;Step 4: Guidance method for smooth gliding segment;

平稳滑翔阶段的弹道规划是所提出的制导算法的基础,需要在保证过程约束的同时使飞行器接近目标。该算法基于一个新的三维解析解来确定纵向射程和横向射程剖面。这里首先简要介绍新的弹道解析解。The ballistic planning of the smooth glide phase is the basis of the proposed guidance algorithm, which needs to make the aircraft approach the target while guaranteeing the process constraints. The algorithm is based on a new three-dimensional analytical solution to determine the longitudinal and lateral range profiles. Here we first briefly introduce the new analytical solution for ballistics.

为了迎合再入弹道纵程远大于横程的特性,首先建立了如下两个坐标系,可以直观得表示纵程和横程的意义。In order to cater for the characteristic that the re-entry trajectory is much longer than the horizontal, the following two coordinate systems are first established, which can intuitively represent the meaning of the vertical and horizontal.

定义一个随地球转动的辅助地心旋转参照系(AGR坐标系):原点在地心E,

Figure BDA0002408897380000044
轴指向飞行器初始位置,
Figure BDA0002408897380000045
轴位于过飞行器与目标点的大圆平面内,且垂直于
Figure BDA0002408897380000046
轴,
Figure BDA0002408897380000047
轴可根据右手定则确定。Define an auxiliary geocentric rotation reference system (AGR coordinate system) that rotates with the earth: the origin is at the center of the earth E,
Figure BDA0002408897380000044
The axis points to the initial position of the aircraft,
Figure BDA0002408897380000045
The axis lies in the great circle plane passing through the aircraft and the target point, and is perpendicular to the
Figure BDA0002408897380000046
axis,
Figure BDA0002408897380000047
The axis can be determined according to the right-hand rule.

同时,定义当地广义北东下坐标系(GNED坐标系):原点o在飞行器质心M向地面的垂直投影点;

Figure BDA0002408897380000048
轴铅垂向下指向地心,
Figure BDA0002408897380000049
轴指向AGR坐标系中
Figure BDA00024088973800000410
指向的“北”方,
Figure BDA00024088973800000411
轴由右手定则确定。At the same time, define the local generalized north-east lower coordinate system (GNED coordinate system): the vertical projection point of the origin o on the ground from the center of mass M of the aircraft;
Figure BDA0002408897380000048
The axis points vertically downward to the center of the earth,
Figure BDA0002408897380000049
The axis points in the AGR coordinate system
Figure BDA00024088973800000410
pointing to the "north" side,
Figure BDA00024088973800000411
The axis is determined by the right-hand rule.

定义纵程为

Figure BDA00024088973800000412
定义横程为
Figure BDA00024088973800000413
其中
Figure BDA00024088973800000414
Figure BDA00024088973800000415
分别表示AGR坐标系中的广义经度和广义纬度。Define the vertical as
Figure BDA00024088973800000412
Define the traverse as
Figure BDA00024088973800000413
in
Figure BDA00024088973800000414
and
Figure BDA00024088973800000415
represent the generalized longitude and generalized latitude in the AGR coordinate system, respectively.

定义能量为Define energy as

Figure BDA0002408897380000051
Figure BDA0002408897380000051

定义纵向升阻比Lcosσ/D=L1/D,横向升阻比Lsinσ/D=L2/D。在AGR坐标轴中的这个广义经纬度坐标系统有利于线化和解耦复杂的动力学模型,获得如下三维弹道解析解。Define the vertical lift-drag ratio Lcosσ/D=L 1 /D, and the lateral lift-drag ratio Lsinσ/D=L 2 /D. This generalized latitude-longitude coordinate system in the AGR coordinate axis facilitates linearization and decoupling of complex dynamic models to obtain the following three-dimensional ballistic analytical solutions.

Figure BDA0002408897380000052
Figure BDA0002408897380000052

Figure BDA0002408897380000053
Figure BDA0002408897380000053

Figure BDA0002408897380000054
Figure BDA0002408897380000054

其中,

Figure BDA0002408897380000055
R*=Re+H;in,
Figure BDA0002408897380000055
R * = Re +H;

Figure BDA0002408897380000056
hij,i=1,2;j=0,1,α1和pij,i=1,2;j=0,1,2,3,4均为常值系数。E,E0分别表示当前能量与初始能量值;xE表示能量积分变量;L1表示纵向平面内的升力值;xC0表示初始横程值;
Figure BDA0002408897380000057
表示GNED坐标系下的初始航向角误差;
Figure BDA0002408897380000058
表示GNED坐标系下的航向角误差。
Figure BDA0002408897380000056
h ij ,i=1,2; j=0,1,α 1 and p ij ,i=1,2;j=0,1,2,3,4 are constant coefficients. E, E 0 represent the current energy and initial energy value respectively; x E represents the energy integral variable; L 1 represents the lift value in the longitudinal plane; x C0 represents the initial transverse value;
Figure BDA0002408897380000057
Represents the initial heading angle error in the GNED coordinate system;
Figure BDA0002408897380000058
Indicates the heading angle error in the GNED coordinate system.

由于该解析解具有很高的精度,可将其应用于制导方法推导中,以减少计算量提高普适性。Due to the high accuracy of the analytical solution, it can be applied to the derivation of the guidance method to reduce the computational complexity and improve the universality.

1.基准攻角设计1. Baseline angle of attack design

因为解析解是以能量为自变量设计的,为更方便地利用解析解,该制导方法中所有的参考剖面均同样采用能量作为自变量。参考攻角剖面αbsl设计为:Because the analytical solution is designed with energy as the independent variable, in order to use the analytical solution more conveniently, all reference profiles in this guidance method also use energy as the independent variable. The reference angle of attack profile α bsl is designed as:

Figure BDA0002408897380000061
Figure BDA0002408897380000061

其中,Eα=-5.55×107J/kg位于平稳滑翔段与最后的高度调整段的交界点附近。为发挥出飞行器的最大能力,设计攻角为最大升阻比对应攻角,即α1=10deg。设计攻角为能量的二次函数是为了使得攻角可以从α1平缓地过渡到α2。在这个制导方法中,首次预设α2=6deg,后面会再应用弹道仿真去对其进行修正,以严格满足终端高度约束。当攻角剖面设计完成时,相应的参考升阻比剖面(L/D)bsl也随之确定。Among them, E α =-5.55×10 7 J/kg is located near the junction of the smooth glide segment and the final height adjustment segment. In order to exert the maximum capability of the aircraft, the design angle of attack is the angle of attack corresponding to the maximum lift-drag ratio, that is, α 1 =10deg. The angle of attack is designed to be a quadratic function of energy so that the angle of attack can transition smoothly from α 1 to α 2 . In this guidance method, α 2 =6deg is preset for the first time, and then ballistic simulation is applied to correct it to strictly satisfy the terminal height constraint. When the design of the angle of attack profile is completed, the corresponding reference lift-to-drag ratio profile (L/D) bsl is also determined.

2.基准升阻比剖面设计2. Benchmark lift-drag ratio profile design

为满足射程要求,设计纵向升阻比剖面如下In order to meet the range requirements, the designed longitudinal lift-drag ratio profile is as follows

Figure BDA0002408897380000062
Figure BDA0002408897380000062

其中,令终端纵程解析解xDf等于剩余射程,即可求得L1/D1。设计L1/D2=L/D(ETAEM),其中ETAEM为终端能量,这样当xE=ETAEM时,(L1/D)bsl|E=E_TAEM=(LcosσTAEM/D)bsl=L/D(ETAEM),可以使得终端的倾侧角σTAEM为0。通过设计式(20)的形式,可在飞行过程中控制倾侧角近似为常值。Wherein, let the terminal longitudinal range analytical solution x Df be equal to the remaining range, then L 1 /D 1 can be obtained. Design L 1 /D 2 =L/D(E TAEM ), where E TAEM is the terminal energy, such that when x E =E TAEM , (L 1 /D) bsl | E = E_TAEM =(Lcosσ TAEM /D) bsl =L/D(E TAEM ), the inclination angle σ TAEM of the terminal can be made zero. Through the form of design formula (20), the roll angle can be controlled to be approximately constant during the flight.

因此,为求解L1/D1,首先需要计算剩余射程xDf。在GER坐标系下,定义η为从地球中心指向飞行器的矢量

Figure BDA0002408897380000063
与从地球中心指向飞行器的矢量
Figure BDA0002408897380000064
之间的夹角。则剩余射程可表示为:Therefore, in order to solve L 1 /D 1 , the remaining range x Df needs to be calculated first. In the GER coordinate system, η is defined as the vector pointing from the center of the earth to the aircraft
Figure BDA0002408897380000063
Vector with aircraft pointing from the center of the earth
Figure BDA0002408897380000064
the angle between. Then the remaining range can be expressed as:

xDf=Reη-STAEM (21)x Df = R e η-S TAEM (21)

其中,STAEM表示再入段结束时到目标的距离。Among them, S TAEM represents the distance to the target at the end of the re-entry segment.

求解L1/D1的过程可描述如下:The process of solving L 1 /D 1 can be described as follows:

当E≥Eα时,有When E≥E α , there is

xD(Eα,E)+xD(ETAEM,Eα)=xDf (22)x D (E α ,E)+x D (E TAEM ,E α )=x Df (22)

其中,in,

Figure BDA0002408897380000071
Figure BDA0002408897380000071

其中,in,

Figure BDA0002408897380000072
Figure BDA0002408897380000072

Figure BDA0002408897380000073
Figure BDA0002408897380000073

并且有,and there is,

Figure BDA0002408897380000074
Figure BDA0002408897380000074

其中,in,

Figure BDA0002408897380000075
Figure BDA0002408897380000075

c2=(a1-c1d0)/d1 (28)c 2 =(a 1 -c 1 d 0 )/d 1 (28)

c3=a0-c2d0 (29)c 3 =a 0 -c 2 d 0 (29)

Figure BDA0002408897380000076
Figure BDA0002408897380000076

Figure BDA0002408897380000077
Figure BDA0002408897380000077

Figure BDA0002408897380000078
Figure BDA0002408897380000078

求解式(22)可得Solving equation (22), we can get

Figure BDA0002408897380000079
Figure BDA0002408897380000079

其中,in,

Figure BDA0002408897380000081
Figure BDA0002408897380000081

Figure BDA0002408897380000082
Figure BDA0002408897380000082

当E<Eα,时,飞行器主要位于高度调整阶段,不再更新L1/D1的值。When E<E α , the aircraft is mainly in the altitude adjustment stage, and the value of L 1 /D 1 is no longer updated.

3.基准横向升阻比和基准倾侧角设计3. Design of reference lateral lift-drag ratio and reference roll angle

横向升阻比的模值可由之前规划好的(L/D)bsl和纵向升阻比(L1/D)bsl决定,为

Figure BDA0002408897380000083
为满足终端横程约束,设计横向升阻比剖面如下:The modulus value of the lateral lift-drag ratio can be determined by the previously planned (L/D) bsl and the longitudinal lift-drag ratio (L 1 /D) bsl , which is
Figure BDA0002408897380000083
In order to meet the terminal transverse constraints, the designed lateral lift-drag ratio profile is as follows:

Figure BDA0002408897380000084
Figure BDA0002408897380000084

其中,EBR1和EBR2分别表示第一次和第二次倾侧反转点。sgn为表示初始倾侧反转方向的符号变量。得到横向升阻比剖面后,参考倾侧角剖面σbsl也随之确定,表示如下:Among them, E BR1 and E BR2 represent the first and second tilt reversal points, respectively. sgn is a symbolic variable representing the direction of the initial roll reversal. After the lateral lift-drag ratio profile is obtained, the reference dip angle profile σ bsl is also determined, which is expressed as follows:

Figure BDA0002408897380000085
Figure BDA0002408897380000085

其中,(L/D)real表示实际的升阻比剖面。Among them, (L/D) real represents the actual lift-drag ratio profile.

4.倾侧反转点修正算法4. Tilt reversal point correction algorithm

在这个算法中,基于横程解析解设计倾侧反转点以消除终端横程误差。这里安排两次倾侧反转点。首先设计并修正第一次反转点EBR1。先令第二次反转点EBR2=Eα,则终端横程xC(ETAEM,E)仅仅是EBR1的函数。从式(17)中可以看出,终端横程xCf可表示为EBR1的分段函数如下:In this algorithm, the roll reversal point is designed based on the analytical solution of the lateral run to eliminate the terminal run error. There are two tilt reversal points arranged here. First design and correct the first reversal point E BR1 . Shilling the second reversal point E BR2 =E α , the terminal stroke x C (E TAEM ,E) is only a function of E BR1 . It can be seen from equation (17) that the terminal traverse x Cf can be expressed as a piecewise function of E BR1 as follows:

Figure BDA0002408897380000086
Figure BDA0002408897380000086

其中,F、G为自定义函数,表达式如下Among them, F and G are self-defined functions, and the expressions are as follows

Figure BDA0002408897380000091
Figure BDA0002408897380000091

Figure BDA0002408897380000092
Figure BDA0002408897380000092

其中,参数xE2,xE1表达能量上下界;k为1-5的常数;P1表示多项式拟合系数。Among them, the parameters x E2 , x E1 express the upper and lower bounds of the energy; k is a constant from 1 to 5; P 1 represents the polynomial fitting coefficient.

对式(38)求EBR1的导数,可得Taking the derivative of EBR1 with respect to equation (38), we can get

Figure BDA0002408897380000093
Figure BDA0002408897380000093

其中,E=EBR1处的右极限用上标+表示,E=EBR1处的左极限用上标-表示。因为L2/Dbsl在E=EBR1处不连续,所以有

Figure BDA0002408897380000094
i=0,1,2,3,4。为了使终端横程误差为0,即需要满足xCf(EBR1)=0。因此这个算法的本质在于求解非线性方程xCf(EBR1)=0的根,用牛顿法求解方法如下:Among them, the right limit at E=E BR1 is represented by a superscript + , and the left limit at E=E BR1 is represented by a superscript - . Since L 2 /D bsl is discontinuous at E=E BR1 , we have
Figure BDA0002408897380000094
i=0,1,2,3,4. In order to make the terminal traverse error 0, it is necessary to satisfy x Cf (E BR1 )=0. Therefore, the essence of this algorithm is to solve the root of the nonlinear equation x Cf (E BR1 )=0. The Newton method is used to solve it as follows:

Figure BDA0002408897380000095
Figure BDA0002408897380000095

参数

Figure BDA0002408897380000096
Figure BDA0002408897380000097
并不一样,
Figure BDA0002408897380000098
表示EBR1的k次幂,
Figure BDA0002408897380000099
表示第k次迭代的EBR1值。求解停止设置为
Figure BDA00024088973800000910
因此xCf(EBR1)随着EBR1单调变化,此非线性方程的解通过简单的几次迭代即可得到。parameter
Figure BDA0002408897380000096
and
Figure BDA0002408897380000097
not the same,
Figure BDA0002408897380000098
represents E BR1 to the power of k,
Figure BDA0002408897380000099
Represents the value of E BR1 for the k-th iteration. Solve Stop is set to
Figure BDA00024088973800000910
Therefore x Cf (E BR1 ) varies monotonically with E BR1 , and the solution to this nonlinear equation can be obtained in a few simple iterations.

5.指令攻角与倾侧角设计5. Design of command angle of attack and tilt angle

为保证参考剖面的跟踪精度,这里采用了弹道阻尼控制技术来抑制再入弹道存在的长周期、弱阻尼的振荡。设计指令攻角αcmd与指令倾侧角σcmd如下:In order to ensure the tracking accuracy of the reference profile, the ballistic damping control technology is used to suppress the long-period and weakly damped oscillation of the re-entry trajectory. The design command angle of attack α cmd and command tilt angle σ cmd are as follows:

αcmd=αbsl+cosσbslKγSG-γ) (43)α cmdbsl +cosσ bsl K γSG -γ) (43)

Figure BDA00024088973800000911
Figure BDA00024088973800000911

其中,Kγ是反馈系数,取值为5;γSG是平稳滑翔弹道倾角的近似解,可以表示如下Among them, K γ is the feedback coefficient, the value is 5; γ SG is the approximate solution of the inclination angle of the smooth gliding trajectory, which can be expressed as follows

Figure BDA0002408897380000101
Figure BDA0002408897380000101

其中,S为飞行器参考面积;CL为升力系数;ρ为空气密度。Among them, S is the reference area of the aircraft; C L is the lift coefficient; ρ is the air density.

平稳滑翔段需要满足过程约束,故将式(8)-(10)的过程约束转变为倾侧角约束,表示如下The smooth gliding segment needs to satisfy the process constraints, so the process constraints of equations (8)-(10) are transformed into the inclination angle constraints, which are expressed as follows

Figure BDA0002408897380000102
Figure BDA0002408897380000102

其中,Lmax为最大升力,表达式如下:Among them, L max is the maximum lift, and the expression is as follows:

Figure BDA0002408897380000103
Figure BDA0002408897380000103

Figure BDA0002408897380000104
Figure BDA0002408897380000104

Figure BDA0002408897380000105
Figure BDA0002408897380000105

其中,Hmin对应滑翔高度的最低边界,可以通过过程约束求得;常值系数kσ=-50。因此,为满足过程约束,倾侧角需要满足下列条件:Among them, H min corresponds to the lowest boundary of the gliding height, which can be obtained through process constraints; the constant coefficient k σ =-50. Therefore, in order to satisfy the process constraints, the tilt angle needs to satisfy the following conditions:

cmd|≤σmax (50)cmd |≤σ max (50)

步骤五:高度调整段制导方法;Step 5: Guidance method for height adjustment segment;

1.指令攻角与第二次反转点修正。1. Command angle of attack and second reversal point correction.

(1)指令攻角α2修正算法( 1 ) Correction algorithm for command angle of attack α2

为了满足终端高度约束,这里采用一次弹道仿真来修正指令攻角α2。根据弹道仿真所得到的终端高度Hf和终端约束高度HTAEM之间的偏差大小,来进一步修正得到最后的指令攻角α2。之后的变量中,下标f表示由弹道仿真计算得到的终端值,下标TAEM表示终端约束值。In order to satisfy the terminal height constraint, a ballistic simulation is used here to correct the command angle of attack α 2 . According to the deviation between the terminal height H f obtained by the ballistic simulation and the terminal constraint height H TAEM , the final command angle of attack α 2 is obtained by further correction. In the following variables, the subscript f represents the terminal value calculated by the ballistic simulation, and the subscript TAEM represents the terminal constraint value.

由于γff

Figure BDA0002408897380000106
可以近似为0,因此可分别假设cosγf=1,cosσf=1和
Figure BDA0002408897380000107
所以根据式(6)可以得到:Since γ ff and
Figure BDA0002408897380000106
can be approximated to 0, so it can be assumed that cosγ f =1, cosσ f =1 and
Figure BDA0002408897380000107
So according to formula (6), we can get:

Figure BDA0002408897380000111
Figure BDA0002408897380000111

其中,CL(est)bsl(Ef))表示在弹道仿真中的终端攻角对应的升力系数,qf表示终端动压。Among them, C L(est)bsl (E f )) represents the lift coefficient corresponding to the terminal angle of attack in the ballistic simulation, and q f represents the terminal dynamic pressure.

根据终端约束,有According to the terminal constraints, there are

Figure BDA0002408897380000112
Figure BDA0002408897380000112

其中,

Figure BDA0002408897380000113
表示修正后的指令攻角,φf、ψf、Vf、Hf分别表示弹道仿真计算得到的终端纬度、航向角、速度、高度;qTAEM、φTAEM、ψTAEM、VTAEM、HTAEM分别表示终端约束的动压值、纬度值、航向角值、速度值、高度值。in,
Figure BDA0002408897380000113
Represents the corrected commanded angle of attack, φ f , ψ f , V f , and H f respectively represent the terminal latitude, heading angle, speed, and altitude calculated by ballistic simulation; q TAEM , φ TAEM , ψ TAEM , V TAEM , H TAEM Respectively represent the dynamic pressure value, latitude value, heading angle value, speed value, and altitude value of the terminal constraint.

在φf=φTAEM和ψf=ψTAEM的假设下,线性化

Figure BDA0002408897380000114
通过令式(51)等于式(52)可以得到:Under the assumptions of φ f = φ TAEM and ψ f = ψ TAEM , the linearization
Figure BDA0002408897380000114
By making equation (51) equal to equation (52), we can get:

Figure BDA0002408897380000115
Figure BDA0002408897380000115

其中,Ef表示弹道仿真计算得到的终端能量值;

Figure BDA0002408897380000116
表示升力系数对攻角的导数;因此,式(53)中求解出的α2可带入式(19)中作为基准攻角进行制导控制。Among them, E f represents the terminal energy value calculated by ballistic simulation;
Figure BDA0002408897380000116
represents the derivative of the lift coefficient with respect to the angle of attack; therefore, α 2 solved in Equation (53) can be brought into Equation (19) as the reference angle of attack for guidance control.

(2)第二次倾侧反转点EBR2修正算法(2) Correction algorithm for the second inclination reversal point E BR2

高度调整阶段设计采用比例导引律进行制导。通过对该制导算法的分析可知,如果减小EBR2,制导律指令将会产生更大的倾侧角以消除最后一次倾侧反转时产生的更大的航向误差,从而使纵向升阻比减小,降低终端速度Vf。因此,Vf可视为关于EBR2的单调函数。为保证终端约束,这个问题可视为求解非线性方程Vf(EBR2)=VTAEM的解的问题。这里采用割线法进行求解,通过多个弹道仿真对Vf进行预测,并根据预测值与终端约束的偏差调整EBR2的值。The height adjustment stage is designed with proportional guidance law for guidance. Through the analysis of the guidance algorithm, it can be seen that if the E BR2 is reduced, the guidance law command will generate a larger pitch angle to eliminate the larger heading error generated during the last roll reversal, thereby reducing the longitudinal lift-to-drag ratio. , reducing the terminal velocity V f . Therefore, V f can be regarded as a monotonic function with respect to E BR2 . To ensure terminal constraints, this problem can be viewed as a problem of solving the solution of the nonlinear equation V f (E BR2 )=V TAEM . Here, the secant method is used to solve the problem, and V f is predicted through multiple ballistic simulations, and the value of E BR2 is adjusted according to the deviation of the predicted value from the terminal constraint.

Figure BDA0002408897380000117
Figure BDA0002408897380000117

EBR2结束迭代修正的条件为求出的偏差满足

Figure BDA0002408897380000121
E BR2 The condition for ending iterative correction is that the obtained deviation satisfies
Figure BDA0002408897380000121

2.基准倾侧角2. Baseline tilt angle

在第二个倾侧反转以后,飞行器进入高度调整阶段。在此阶段,基准攻角仍由式(19)确定,基准倾侧角则由比例导引制导律控制。视线角速率

Figure BDA0002408897380000122
可以表示为After the second roll reversal, the aircraft enters the altitude adjustment phase. At this stage, the reference angle of attack is still determined by equation (19), and the reference tilt angle is controlled by the proportional guidance guidance law. line-of-sight rate
Figure BDA0002408897380000122
It can be expressed as

Figure BDA0002408897380000123
Figure BDA0002408897380000123

然后,通过比例导引获得的横向机动加速度指令aL2如下Then, the lateral maneuvering acceleration command a L2 obtained by proportional guidance is as follows

Figure BDA0002408897380000124
Figure BDA0002408897380000124

其中,kPN为横向加速度机动系数,取值为

Figure BDA0002408897380000125
Figure BDA0002408897380000126
表示第二次倾侧反转时的纵程。假设飞行器近似平稳滑翔,则纵向平面内的合力加速度为0。因此,纵向机动加速度指令aL1为Among them, k PN is the lateral acceleration maneuvering coefficient, and its value is
Figure BDA0002408897380000125
Figure BDA0002408897380000126
Indicates the longitudinal distance at the second roll reversal. Assuming that the aircraft glides approximately smoothly, the resultant acceleration in the longitudinal plane is zero. Therefore, the longitudinal maneuvering acceleration command a L1 is

Figure BDA0002408897380000127
Figure BDA0002408897380000127

所以可得基准倾侧角为Therefore, the reference tilt angle can be obtained as

Figure BDA0002408897380000128
Figure BDA0002408897380000128

3.指令攻角与倾侧角3. Command angle of attack and tilt angle

为了确保终端速度约束,这里通过调整攻角来跟踪以能量为自变量的剩余滑翔距离参考剖面

Figure BDA0002408897380000129
指令攻角与倾侧角可以表示如下:In order to ensure the terminal speed constraint, the remaining glide distance reference profile with energy as the independent variable is tracked by adjusting the angle of attack.
Figure BDA0002408897380000129
The commanded attack angle and pitch angle can be expressed as follows:

Figure BDA00024088973800001210
Figure BDA00024088973800001210

Figure BDA00024088973800001211
Figure BDA00024088973800001211

其中,kα的值可根据任务及时调整,一般可选为2π/(1.8×107);

Figure BDA00024088973800001212
表示剩余射程,值可以通过之前的最后一次弹道仿真获得。Among them, the value of k α can be adjusted in time according to the task, and generally can be selected as 2π/(1.8×10 7 );
Figure BDA00024088973800001212
Indicates the remaining range, the value can be obtained from the last ballistic simulation before.

为了满足过程约束,像前文一样,倾侧角仍然需要满足下列要求:To satisfy the process constraints, as before, the roll angle still needs to satisfy the following requirements:

cmd|≤σmax (61)cmd |≤σ max (61)

通过上述五个步骤,可得到了一种基于平稳滑翔弹道解析解的高超声速再入制导方法。Through the above five steps, a hypersonic reentry guidance method based on the analytical solution of the smooth gliding trajectory can be obtained.

本发明一种基于平稳滑翔弹道解析解的高超声速再入制导方法,优点在于:The present invention is a hypersonic reentry guidance method based on the analytical solution of the smooth gliding trajectory, and has the advantages of:

(1)基于三维弹道解析解进行纵程和横程的预测矫正,极大减少了计算量,提高了计算速度,便于实时应用到机载系统上。(1) The prediction and correction of the longitudinal and transverse distances based on the three-dimensional ballistic analytical solution greatly reduces the amount of calculation, improves the calculation speed, and is convenient for real-time application to airborne systems.

(2)横向运动由两次倾侧反转点确定,修正第一次反转点位置来控制终端横程误差,修正第二次反转点位置来满足终端速度约束。相比于传统的利用航向误差门限控制横向运动方法,该方法减少了倾侧反转次数,提高了能量利用效率,减低了对控制系统的需求。(2) The lateral movement is determined by two tilt 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 satisfy the terminal speed constraint. Compared with the traditional method of using the heading error threshold to control the lateral motion, the method reduces the number of roll reversals, improves the energy utilization efficiency, and reduces the demand for the control system.

(3)高度调整段利用比例导引的思想,使得终端倾侧角收敛到0,增强了目标打击效果。(3) The height adjustment section uses the idea of proportional guidance to make the terminal tilt angle converge to 0, which enhances the target strike effect.

附图说明Description of drawings

图1是地心赤道旋转坐标系(GER)和当地北东下坐标系(NED)示意图。Figure 1 is a schematic diagram of the geocentric equatorial rotating coordinate system (GER) and the local northeast lower coordinate system (NED).

图2是辅助地心旋转参照系(AGR)与辅助当地北东下坐标系(GNED)示意图。FIG. 2 is a schematic diagram of an auxiliary geocentric rotation frame of reference (AGR) and an auxiliary local northeast lower coordinate system (GNED).

图3是五种目标状态的经纬度曲线。Figure 3 is the latitude and longitude curves of five target states.

图4是五种目标状态的高度-能量曲线。Figure 4 is a height-energy curve for five target states.

图5是五种目标状态的速度-能量曲线。Figure 5 is a velocity-energy curve for five target states.

图6是五种目标状态的攻角-能量曲线。Figure 6 is an angle of attack-energy curve for five target states.

图7是五种目标状态的倾侧角-能量曲线。FIG. 7 is a tilt angle-energy curve for five target states.

图8是五种目标状态的弹道倾角-能量曲线。Fig. 8 is the ballistic inclination-energy curve of five target states.

图9是五种目标状态的航向误差-能量曲线。Figure 9 is a heading error-energy curve for five target states.

上述图中,涉及到的符号、代号说明如下:In the above figure, the symbols and codes involved are explained as follows:

图1中,

Figure BDA0002408897380000131
表示地心赤道旋转坐标系(GER),oxyz表示当地北东下坐标系(NED),North表示北向,East表示东向。h是海拔高度,R是射程,λ是经度,φ是纬度,V是飞行器相对于地球的速率,γ是弹道倾角,ψ是飞行器航向角,以当地北向为基准,σ是倾侧角。图2中,
Figure BDA0002408897380000132
表示辅助地心旋转参照系(AGR),
Figure BDA0002408897380000133
表示当地北东下坐标系(NED),Axis表示地球极轴。M是飞行器当前位置点,o是飞行器当前位置点在地球表面的投影点。T是目标点,oT是目标点在地球表面的投影点。ωex表示地球自转角速度在
Figure BDA0002408897380000141
轴的分量,ωey表示地球自转角速度在
Figure BDA0002408897380000142
轴的分量,ωzx表示地球自转角速度在
Figure BDA0002408897380000143
轴的分量。图3中,Longitude表示经度,Latitude表示纬度,T1-T5代表五个不同目标的算例。图4中,Altitude表示高度,Energy表示能量,H表示弹道仿真得到的高度曲线,Hmin表示满足过程约束的最小高度边界。图5中,Velocity表示高度,Energy表示能量。图6中,Angle of Attack表示攻角,Energy表示能量。图7中,Bank Angle表示倾侧角,Energy表示能量。图8中,Flight-pathAngle表示倾侧角,Energy表示能量。图9中,Heading Angle表示航向角,Energy表示能量。In Figure 1,
Figure BDA0002408897380000131
Represents the geocentric equatorial rotation coordinate system (GER), o xyz represents the local north-east lower coordinate system (NED), North represents the north direction, and East represents the east direction. h is the altitude, R is the range, λ is the longitude, φ is the latitude, V is the velocity of the aircraft relative to the earth, γ is the ballistic inclination, ψ is the heading angle of the aircraft, based on the local north direction, and σ is the roll angle. In Figure 2,
Figure BDA0002408897380000132
represents the auxiliary geocentric rotation frame of reference (AGR),
Figure BDA0002408897380000133
Represents the local Northeast Lower Coordinate System (NED), and Axis represents the Earth's polar axis. M is the current position of the aircraft, and o is the projection of the current position of the aircraft on the surface of the earth. T is the target point, o T is the projected point of the target point on the earth's surface. ω ex represents the angular velocity of the earth's rotation in
Figure BDA0002408897380000141
The component of the axis, ω ey represents the angular velocity of the earth's rotation in
Figure BDA0002408897380000142
The component of the axis, ω zx represents the angular velocity of the earth's rotation in
Figure BDA0002408897380000143
Axis components. In Figure 3, Longitude represents longitude, Latitude represents latitude, and T1-T5 represent the calculation examples of five different targets. In Figure 4, Altitude represents altitude, Energy represents energy, H represents the altitude curve obtained by ballistic simulation, and H min represents the minimum altitude boundary that satisfies the process constraints. In Figure 5, Velocity represents height, and Energy represents energy. In Figure 6, Angle of Attack represents the angle of attack, and Energy represents energy. In Fig. 7, Bank Angle represents the bank angle, and Energy represents the energy. In Fig. 8, Flight-pathAngle represents the tilt angle, and Energy represents energy. In Figure 9, Heading Angle represents the heading angle, and Energy represents the energy.

具体实施方式Detailed ways

下面将结合附图和实施案例对本发明作进一步的详细说明。The present invention will be further described in detail below with reference to the accompanying drawings and implementation cases.

针对高超声速飞行器再入制导上述问题,提出了一种基于三维弹道解析解的高精度高超滑翔再入制导方法。首先构建辅助地心旋转坐标系,并基于此建立了一个新的纵程、横程和航向角关于能量的动力学模型,该模型有利于对高度非线性、强耦合的动力学模型进行解耦和线化,得到三维弹道解析解,并基于此解析解设计了制导方法。制导方法利用纵程解析解实时规划纵向参考剖面,利用横程解析解修正第一次反转点以消除横程误差,利用几次弹道仿真来修正第二次反转点来保证终端速度约束、修正攻角最小值保证终端高度约束。Aiming at the above problems of hypersonic vehicle reentry guidance, a high-precision hyperglide reentry guidance method based on 3D ballistic analytical solution was proposed. Firstly, an auxiliary geocentric rotation coordinate system is constructed, and based on this, a new dynamic model of longitudinal, transverse and heading angles with respect to energy is established, which is beneficial to decoupling highly nonlinear and strongly coupled dynamic models. And linearization, the three-dimensional ballistic analytical solution is obtained, and the guidance method is designed based on this analytical solution. The guidance method uses the longitudinal analytical solution to plan the longitudinal reference profile in real time, uses the lateral analytical solution to correct the first reversal point to eliminate the lateral error, and uses several ballistic simulations to correct the second inversion point to ensure the terminal speed constraint, Corrected minimum angle of attack to ensure terminal height constraints.

本发明一种基于平稳滑翔弹道解析解的高超声速再入制导方法,其特征在于:它包括以下几个步骤:The present invention is a hypersonic reentry guidance method based on the analytical solution of the smooth gliding trajectory, characterized in that it comprises the following steps:

步骤一:建立飞行器动力学模型;Step 1: Build the aircraft dynamics model;

如图1所示,建立地心赤道旋转坐标系(GER坐标系):原点在地心E,ze轴垂直于地球赤道表面,指向北极;xe轴和ye轴在赤道平面内,且相互垂直。该坐标系随地球一起转动,其绕ze轴的旋转角速度为地球自转角速度ωeAs shown in Figure 1, establish the geocentric equatorial rotating coordinate system (GER coordinate system): the origin is at the earth's center E, the z e axis is perpendicular to the earth's equatorial surface and points to the north pole; the x e axis and the y e axis are in the equatorial plane, and perpendicular to each other. The coordinate system rotates with the earth, and its rotational angular velocity around the z e axis is the earth's rotation angular velocity ω e .

当地北东下坐标系(NED坐标系):定义原点o在飞行器质心M向地面的垂直投影点;,x轴指向当地北方,y轴指向当地东方,z轴铅垂向下指向地心。Local North-East Lower Coordinate System (NED Coordinate System): Define the vertical projection point of the origin o at the center of mass M of the aircraft to the ground; the x-axis points to the local north, the y-axis points to the local east, and the z-axis points vertically downward to the center of the earth.

旋转地球背景下,采用标准大气模型和平方反比引力场模型,建立高超声速飞行器的六自由度动力学方程如下:Under the background of the rotating earth, using the standard atmospheric model and the inverse square gravitational field model, the six-degree-of-freedom dynamic equation of the hypersonic vehicle is established as follows:

Figure BDA0002408897380000151
Figure BDA0002408897380000151

Figure BDA0002408897380000152
Figure BDA0002408897380000152

Figure BDA0002408897380000153
Figure BDA0002408897380000153

Figure BDA0002408897380000154
Figure BDA0002408897380000154

Figure BDA0002408897380000155
Figure BDA0002408897380000155

Figure BDA0002408897380000156
Figure BDA0002408897380000156

Figure BDA0002408897380000157
Figure BDA0002408897380000157

其中,d/dt表示对时间求导,t是时间,h是海拔高度,R是射程,λ是经度,φ是纬度,V是飞行器相对于地球的速率,γ是弹道倾角,ψ是飞行器航向角,以当地北向为基准,σ是倾侧角,Re是地球半径,取值6378.137km,ωe是地球自转角速度,L=ρV2SCL/2m为升力加速度,D=ρV2SCD/2m为阻力加速度;S为飞行器参考面积;CL为升力系数;ρ为空气密度。where d/dt is the derivation of 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, and ψ is the heading of the aircraft Angle, based on the local north direction, σ is the tilt angle, Re is the radius of the earth, taking the value 6378.137km, ω e is the angular velocity of the earth's rotation, L=ρV 2 SC L /2m is the lift acceleration, D=ρV 2 SC D / 2m is the drag acceleration; S is the reference area of the aircraft; C L is the lift coefficient; ρ is the air density.

步骤二:建立约束模型;Step 2: Establish a constraint model;

再入飞行器在飞行过程中需要满足热流密度、动压和过载等过程约束,表示如下:The re-entry vehicle needs to meet the process constraints such as heat flux density, dynamic pressure and overload during flight, which are expressed as follows:

Figure BDA0002408897380000161
Figure BDA0002408897380000161

Figure BDA0002408897380000162
Figure BDA0002408897380000162

Figure BDA0002408897380000163
Figure BDA0002408897380000163

其中,KQ是热流密度系数;

Figure BDA0002408897380000164
表示最大允许热流密度值;nmax表示最大允许过载值;qmax表示最大允许动压值.此外,由于飞行控制系统的能力制约,控制量(攻角α和倾侧角σ)的变化率也满足约束:
Figure BDA0002408897380000165
where K Q is the heat flux density coefficient;
Figure BDA0002408897380000164
Represents the maximum allowable heat flux density value; n max represents the maximum allowable overload value; q max represents the maximum allowable dynamic pressure value. In addition, due to the ability of the flight control system, the rate of change of the control variables (angle of attack α and angle of tilt σ) also satisfies constraint:
Figure BDA0002408897380000165

除了过程约束,飞行器还应满足终端约束,有:VTAEM=2000m/s,HTAEM=25km,Rtm=100km,σTAEM≈0°,γTAEM≈0°andΔψTAEM≈0°.In addition to the process constraints, the aircraft should also satisfy the terminal constraints, which are: V TAEM = 2000m/s, H TAEM = 25km, R tm = 100km, σ TAEM ≈0°, γ TAEM ≈0° and Δψ TAEM ≈0°.

步骤三:初始下降段制导方法;Step 3: Guidance method for the initial descent segment;

根据再入弹道的特性,将其划分为三部分:初始下降段,平稳滑翔段和高度调整段。设计初始下降段是为了让弹道更平缓地过渡到的平稳滑翔段。全程飞行弹道绝大部分均处于平稳滑翔段,它不但决定了飞行终端射程、速度和高度,还决定了飞行器是否能满足多种过程约束。高度调整段是弹道的最后一个阶段,关系到此弹道能否严格满足终端约束。飞行器在此段飞行速度和飞行高度都较低,对设计高度调整段的制导方法提出了较为严苛的要求。According to the characteristics of the re-entry trajectory, it is divided into three parts: the initial descent section, the smooth glide section and the height adjustment section. The initial descent segment is designed to allow for a smoother glide segment into which the trajectory transitions more smoothly. Most of the whole flight trajectory is in the smooth gliding segment, which not only determines the range, speed and altitude of the flight terminal, but also determines whether the aircraft can meet various process constraints. The height adjustment segment is the last stage of the trajectory, which is related to whether the trajectory can strictly meet the terminal constraints. The flight speed and flight altitude of the aircraft in this section are relatively low, which puts forward more stringent requirements for the guidance method of the design height adjustment section.

在初始下降段,控制量的变化类对飞行影响不大,所以这里设计固定大小的攻角和倾侧角值以避免高度变化。为了降低热流密度,增加射程,弹道的起始点应该设置的高一些。因此这里选择最大攻角αmax=20deg和零倾侧角进行飞行。当弹道倾角变化率

Figure BDA0002408897380000166
时,意味着升力可以满足平稳滑翔需求,所以此时攻角可以平缓地从最大攻角过渡到参考攻角αbsl。In the initial descent stage, the variation of the control quantity has little effect on the flight, so the angle of attack and the angle of inclination are designed with a fixed size to avoid the altitude change. In order to reduce the heat flux density and increase the range, the starting point of the trajectory should be set higher. Therefore, the maximum angle of attack α max = 20deg and zero roll angle are selected for flight here. When the rate of change of ballistic inclination
Figure BDA0002408897380000166
, which means that the lift can meet the needs of smooth gliding, so the angle of attack can smoothly transition from the maximum angle of attack to the reference angle of attack α bsl .

在这个阶段,指令攻角和倾侧角设计如下:At this stage, the commanded attack and roll angles are designed as follows:

Figure BDA0002408897380000167
Figure BDA0002408897380000167

σcmd=0 (12)σ cmd = 0 (12)

其中,in,

Δγ=γ-γSG (13)Δγ=γ-γ SG (13)

Figure BDA0002408897380000171
Figure BDA0002408897380000171

其中,γSG表示平稳滑翔弹道倾角,kγ表示反馈系数且这里取值为5。当Δγ=0时,飞行器进入平稳滑翔阶段。Among them, γ SG represents the inclination angle of the smooth gliding trajectory, and k γ represents the feedback coefficient and the value here is 5. When Δγ=0, the aircraft enters the smooth gliding phase.

步骤四:平稳滑翔段制导方法;Step 4: Guidance method for smooth gliding segment;

平稳滑翔阶段的弹道规划是所提出的制导算法的基础,需要在保证过程约束的同时使飞行器接近目标。该算法基于一个新的三维解析解来确定纵向射程和横向射程剖面。这里首先简要介绍新的弹道解析解。The ballistic planning of the smooth glide phase is the basis of the proposed guidance algorithm, which needs to make the aircraft approach the target while guaranteeing the process constraints. The algorithm is based on a new three-dimensional analytical solution to determine the longitudinal and lateral range profiles. Here we first briefly introduce the new analytical solution for ballistics.

如图2所示,为了迎合再入弹道纵程远大于横程的特性,首先建立了如下两个坐标系,可以直观得表示纵程和横程的意义。As shown in Figure 2, in order to cater for the characteristic that the re-entry trajectory is much longer than the horizontal, the following two coordinate systems are first established, which can intuitively represent the meaning of the vertical and horizontal.

定义一个随地球转动的辅助地心旋转参照系(AGR坐标系):原点在地心E,

Figure BDA0002408897380000172
轴指向飞行器初始位置,
Figure BDA0002408897380000173
轴位于过飞行器与目标点的大圆平面内,且垂直于
Figure BDA0002408897380000174
轴,
Figure BDA0002408897380000175
轴可根据右手定则确定。Define an auxiliary geocentric rotation reference system (AGR coordinate system) that rotates with the earth: the origin is at the center of the earth E,
Figure BDA0002408897380000172
The axis points to the initial position of the aircraft,
Figure BDA0002408897380000173
The axis lies in the great circle plane passing through the aircraft and the target point, and is perpendicular to the
Figure BDA0002408897380000174
axis,
Figure BDA0002408897380000175
The axis can be determined according to the right-hand rule.

同时,定义当地广义北东下坐标系(GNED坐标系):原点o在飞行器质心M向地面的垂直投影点;

Figure BDA0002408897380000176
轴铅垂向下指向地心,
Figure BDA0002408897380000177
轴指向AGR坐标系中
Figure BDA0002408897380000178
指向的“北”方,
Figure BDA0002408897380000179
轴由右手定则确定。At the same time, define the local generalized north-east lower coordinate system (GNED coordinate system): the vertical projection point of the origin o on the ground from the center of mass M of the aircraft;
Figure BDA0002408897380000176
The axis points vertically downward to the center of the earth,
Figure BDA0002408897380000177
The axis points in the AGR coordinate system
Figure BDA0002408897380000178
pointing to the "north" side,
Figure BDA0002408897380000179
The axis is determined by the right-hand rule.

定义纵程为

Figure BDA00024088973800001710
定义横程为
Figure BDA00024088973800001711
其中
Figure BDA00024088973800001712
Figure BDA00024088973800001713
分别表示AGR坐标系中的广义经度和广义纬度。Define the vertical as
Figure BDA00024088973800001710
Define the traverse as
Figure BDA00024088973800001711
in
Figure BDA00024088973800001712
and
Figure BDA00024088973800001713
represent the generalized longitude and generalized latitude in the AGR coordinate system, respectively.

定义能量为Define energy as

Figure BDA00024088973800001714
Figure BDA00024088973800001714

在AGR坐标轴中的这个广义经纬度坐标系统有利于线化和解耦复杂的动力学模型,获得如下三维弹道解析解。This generalized latitude-longitude coordinate system in the AGR coordinate axis facilitates linearization and decoupling of complex dynamic models to obtain the following three-dimensional ballistic analytical solutions.

Figure BDA0002408897380000181
Figure BDA0002408897380000181

Figure BDA0002408897380000182
Figure BDA0002408897380000182

Figure BDA0002408897380000183
Figure BDA0002408897380000183

其中,

Figure BDA0002408897380000184
R*=Re+H;in,
Figure BDA0002408897380000184
R * = Re +H;

Figure BDA0002408897380000185
hij,i=1,2;j=0,1,α1和pij,i=1,2;j=0,1,2,3,4均为常值系数。由于该解析解具有很高的精度,可将其应用于制导方法推导中,以减少计算量提高普适性。
Figure BDA0002408897380000185
h ij ,i=1,2; j=0,1,α 1 and p ij ,i=1,2;j=0,1,2,3,4 are constant coefficients. Due to the high accuracy of the analytical solution, it can be applied to the derivation of the guidance method to reduce the computational complexity and improve the universality.

6.基准攻角设计6. Baseline angle of attack design

因为解析解是以能量为自变量设计的,为更方便地利用解析解,该制导方法中所有的参考剖面均同样采用能量作为自变量。参考攻角剖面设计为:Because the analytical solution is designed with energy as the independent variable, in order to use the analytical solution more conveniently, all reference profiles in this guidance method also use energy as the independent variable. The reference angle of attack profile is designed as:

Figure BDA0002408897380000186
Figure BDA0002408897380000186

其中,Eα=-5.55×107J/kg位于平稳滑翔段与最后的高度调整段的交界点附近。为发挥出飞行器的最大能力,设计攻角为最大升阻比对应攻角,即α1=10deg。设计攻角为能量的二次函数是为了使得攻角可以从α1平缓地过渡到α2。在这个制导方法中,首次预设α2=6deg,后面会再应用弹道仿真去对其进行修正,以严格满足终端高度约束。当攻角剖面设计完成时,相应的升阻比剖面L/Dbsl也随之确定。Among them, E α =-5.55×10 7 J/kg is located near the junction of the smooth glide segment and the final height adjustment segment. In order to exert the maximum capability of the aircraft, the design angle of attack is the angle of attack corresponding to the maximum lift-drag ratio, that is, α 1 =10deg. The angle of attack is designed to be a quadratic function of energy so that the angle of attack can transition smoothly from α 1 to α 2 . In this guidance method, α 2 =6deg is preset for the first time, and then ballistic simulation is applied to correct it to strictly satisfy the terminal height constraint. When the design of the angle of attack profile is completed, the corresponding lift-drag ratio profile L/D bsl is also determined.

7.基准升阻比剖面设计7. Benchmark lift-drag ratio profile design

为满足射程要求,设计纵向升阻比剖面如下In order to meet the range requirements, the designed longitudinal lift-drag ratio profile is as follows

Figure BDA0002408897380000191
Figure BDA0002408897380000191

其中,令纵程解析解等于剩余射程,即可求得L1/D1。设计L1/D2=L/D(ETAEM)是为了使得终端的倾侧角为0。通过设计式(20)的形式,可在飞行过程中控制倾侧角近似为常值。Wherein, let the longitudinal analytical solution equal to the remaining range, then L 1 /D 1 can be obtained. L 1 /D 2 =L/D(E TAEM ) is designed so that the tilt angle of the terminal is zero. Through the form of design formula (20), the roll angle can be controlled to be approximately constant during the flight.

因此,为求解L1/D1,首先需要计算剩余射程xDf。在GER坐标系下,定义η为从地球中心指向飞行器的矢量

Figure BDA0002408897380000192
与从地球中心指向飞行器的矢量
Figure BDA0002408897380000193
之间的夹角。则剩余射程可表示为:Therefore, in order to solve L 1 /D 1 , the remaining range x Df needs to be calculated first. In the GER coordinate system, η is defined as the vector pointing from the center of the earth to the aircraft
Figure BDA0002408897380000192
Vector with aircraft pointing from the center of the earth
Figure BDA0002408897380000193
the angle between. Then the remaining range can be expressed as:

xDf=Reη-STAEM (21)x Df = R e η-S TAEM (21)

其中,STAEM表示再入段结束时到目标的距离。Among them, S TAEM represents the distance to the target at the end of the re-entry segment.

求解L1/D1的过程可描述如下:The process of solving L 1 /D 1 can be described as follows:

当E≥Eα时,有When E≥E α , there is

xD(Eα,E)+xD(ETAEM,Eα)=xDf (22)x D (E α ,E)+x D (E TAEM ,E α )=x Df (22)

其中,in,

Figure BDA0002408897380000194
Figure BDA0002408897380000194

其中,in,

Figure BDA0002408897380000195
Figure BDA0002408897380000195

Figure BDA0002408897380000196
Figure BDA0002408897380000196

并且有,and there is,

Figure BDA0002408897380000197
Figure BDA0002408897380000197

其中,in,

Figure BDA0002408897380000198
Figure BDA0002408897380000198

c2=(a1-c1d0)/d1 (28)c 2 =(a 1 -c 1 d 0 )/d 1 (28)

c3=a0-c2d0 (29)c 3 =a 0 -c 2 d 0 (29)

Figure BDA0002408897380000201
Figure BDA0002408897380000201

Figure BDA0002408897380000202
Figure BDA0002408897380000202

Figure BDA0002408897380000203
Figure BDA0002408897380000203

求解式(22)可得Solving equation (22), we can get

Figure BDA0002408897380000204
Figure BDA0002408897380000204

其中,in,

Figure BDA0002408897380000205
Figure BDA0002408897380000205

Figure BDA0002408897380000206
Figure BDA0002408897380000206

当E<Eα,时,飞行器主要位于高度调整阶段,不再更新L1/D1的值。When E<E α , the aircraft is mainly in the altitude adjustment stage, and the value of L 1 /D 1 is no longer updated.

8.基准横向升阻比和基准倾侧角设计8. Design of reference lateral lift-drag ratio and reference roll angle

横向升阻比的模值可由之前规划好的L/Dbsl和纵向升阻比L1/Dbs决定,为

Figure BDA0002408897380000207
为满足终端横程约束,设计横向升阻比剖面如下:The modulus value of the lateral lift-drag ratio can be determined by the previously planned L/D bsl and the longitudinal lift-drag ratio L 1 /D bs , which is
Figure BDA0002408897380000207
In order to meet the terminal transverse constraints, the designed lateral lift-drag ratio profile is as follows:

Figure BDA0002408897380000208
Figure BDA0002408897380000208

其中,EBR1和EBR2分别表示第一次和第二次倾侧反转点。sgn为表示初始倾侧反转方向的符号变量。得到横向升阻比剖面后,倾侧角剖面也随之确定,表示如下:Among them, E BR1 and E BR2 represent the first and second tilt reversal points, respectively. sgn is a symbolic variable representing the direction of the initial roll reversal. After the lateral lift-drag ratio profile is obtained, the inclination angle profile is also determined, which is expressed as follows:

Figure BDA0002408897380000211
Figure BDA0002408897380000211

其中,L/Dreal表示实际的升阻比剖面。Among them, L/D real represents the actual lift-drag ratio profile.

9.倾侧反转点修正算法9. Tilt reversal point correction algorithm

在这个算法中,基于横程解析解设计倾侧反转点以消除终端横程误差。这里安排两次倾侧反转点。首先设计并修正第一次反转点EBR1。先令第二次反转点EBR2=Eα,则终端横程xC(ETAEM,E)仅仅是EBR1的函数。从式(17)中可以看出,终端横程可表示为EBR1的分段函数如下:In this algorithm, the roll reversal point is designed based on the analytical solution of the lateral run to eliminate the terminal run error. There are two tilt reversal points arranged here. First design and correct the first reversal point E BR1 . Shilling the second reversal point E BR2 =E α , the terminal stroke x C (E TAEM ,E) is only a function of E BR1 . It can be seen from equation (17) that the terminal traverse can be expressed as a piecewise function of E BR1 as follows:

Figure BDA0002408897380000212
Figure BDA0002408897380000212

其中,in,

Figure BDA0002408897380000213
Figure BDA0002408897380000213

Figure BDA0002408897380000214
Figure BDA0002408897380000214

对式(38)求EBR1的导数,可得Taking the derivative of EBR1 with respect to equation (38), we can get

Figure BDA0002408897380000215
Figure BDA0002408897380000215

其中,E=EBR1处的右极限用上标+表示,E=EBR1处的左极限用上标-表示。因为L2/Dbsl在E=EBR1处不连续,所以有

Figure BDA0002408897380000216
i=0,1,2,3,4。为了使终端横程误差为0,即需要满足xCf(EBR1)=0。因此这个算法的本质在于求解非线性方程xCf(EBR1)=0的根,用牛顿法求解方法如下:Among them, the right limit at E=E BR1 is represented by a superscript + , and the left limit at E=E BR1 is represented by a superscript - . Since L 2 /D bsl is discontinuous at E=E BR1 , we have
Figure BDA0002408897380000216
i=0,1,2,3,4. In order to make the terminal traverse error 0, it is necessary to satisfy x Cf (E BR1 )=0. Therefore, the essence of this algorithm is to solve the root of the nonlinear equation x Cf (E BR1 )=0. The Newton method is used to solve it as follows:

Figure BDA0002408897380000217
Figure BDA0002408897380000217

求解停止设置为

Figure BDA0002408897380000218
因此xCf(EBR1)随着EBR1单调变化,此非线性方程的解通过简单的几次迭代即可得到。Solve Stop is set to
Figure BDA0002408897380000218
Therefore x Cf (E BR1 ) varies monotonically with E BR1 , and the solution to this nonlinear equation can be obtained in a few simple iterations.

10.指令攻角与倾侧角设计10. Design of command angle of attack and tilt angle

为保证参考剖面的跟踪精度,这里采用了弹道阻尼控制技术来抑制再入弹道存在的长周期、弱阻尼的振荡。设计指令攻角与倾侧角如下:In order to ensure the tracking accuracy of the reference profile, the ballistic damping control technology is used to suppress the long-period and weakly damped oscillation of the re-entry trajectory. The design command attack and roll angles are as follows:

αcmd=αbsl+cosσbslKγSG-γ) (43)α cmdbsl +cosσ bsl K γSG -γ) (43)

Figure BDA0002408897380000221
Figure BDA0002408897380000221

其中,Kγ是反馈系数,取值为5;γSG是平稳滑翔弹道倾角的近似解,可以表示如下Among them, K γ is the feedback coefficient, the value is 5; γ SG is the approximate solution of the inclination angle of the smooth gliding trajectory, which can be expressed as follows

Figure BDA0002408897380000222
Figure BDA0002408897380000222

平稳滑翔段需要满足过程约束,故将式(8)-(10)的过程约束转变为倾侧角约束,表示如下The smooth gliding segment needs to satisfy the process constraints, so the process constraints of equations (8)-(10) are transformed into the inclination angle constraints, which are expressed as follows

Figure BDA0002408897380000223
Figure BDA0002408897380000223

其中,in,

Figure BDA0002408897380000224
Figure BDA0002408897380000224

Figure BDA0002408897380000225
Figure BDA0002408897380000225

Figure BDA0002408897380000226
Figure BDA0002408897380000226

其中,Hmin对应滑翔高度的最低边界,可以通过过程约束求得;常值系数kσ=-50。因此,为满足过程约束,倾侧角需要满足下列条件:Among them, H min corresponds to the lowest boundary of the gliding height, which can be obtained through process constraints; the constant coefficient k σ =-50. Therefore, in order to satisfy the process constraints, the tilt angle needs to satisfy the following conditions:

cmd|≤σmax (50)cmd |≤σ max (50)

步骤五:高度调整段制导方法;Step 5: Guidance method for height adjustment segment;

4.指令攻角与第二次反转点修正。4. Command angle of attack and second reversal point correction.

(3)指令攻角α2修正算法( 3 ) Correction algorithm of command angle of attack α2

为了满足终端高度约束,这里采用一次弹道仿真来修正指令攻角α2。根据弹道仿真所得到的终端高度Hf和终端约束高度HTAEM之间的偏差大小,来进一步修正得到最后的指令攻角α2。之后的变量中,下标f表示由弹道仿真计算得到的终端值,下标TAEM表示终端约束值。In order to satisfy the terminal height constraint, a ballistic simulation is used here to correct the command angle of attack α 2 . According to the deviation between the terminal height H f obtained by the ballistic simulation and the terminal constraint height H TAEM , the final command angle of attack α 2 is obtained by further correction. In the following variables, the subscript f represents the terminal value calculated by the ballistic simulation, and the subscript TAEM represents the terminal constraint value.

由于γff

Figure BDA0002408897380000231
可以近似为0,因此可分别假设cosγf=1,cosσf=1和
Figure BDA0002408897380000232
所以根据式(6)可以得到:Since γ ff and
Figure BDA0002408897380000231
can be approximated to 0, so it can be assumed that cosγ f =1, cosσ f =1 and
Figure BDA0002408897380000232
So according to formula (6), we can get:

Figure BDA0002408897380000233
Figure BDA0002408897380000233

其中,CL(est)bsl(Ef))表示在弹道仿真中的终端攻角对应的升力系数,qf表示终端动压。Among them, C L(est)bsl (E f )) represents the lift coefficient corresponding to the terminal angle of attack in the ballistic simulation, and q f represents the terminal dynamic pressure.

根据终端约束,有According to the terminal constraints, there are

Figure BDA0002408897380000234
Figure BDA0002408897380000234

其中,

Figure BDA0002408897380000235
表示修正后的指令攻角。in,
Figure BDA0002408897380000235
Indicates the modified command angle of attack.

在φf=φTAEM和ψf=ψTAEM的假设下,线性化

Figure BDA0002408897380000236
通过令式(51)等于式(52)可以得到:Under the assumptions of φ f = φ TAEM and ψ f = ψ TAEM , the linearization
Figure BDA0002408897380000236
By making equation (51) equal to equation (52), we can get:

Figure BDA0002408897380000237
Figure BDA0002408897380000237

因此,式(53)中求解出的α2可带入式(19)中作为基准攻角进行制导控制。Therefore, α 2 solved in Equation (53) can be brought into Equation (19) as a reference angle of attack for guidance control.

(4)第二次倾侧反转点EBR2修正算法(4) Second inclination reversal point E BR2 correction algorithm

高度调整阶段设计采用比例导引律进行制导。通过对该制导算法的分析可知,如果减小EBR2,制导律指令将会产生更大的倾侧角以消除最后一次倾侧反转时产生的更大的航向误差,从而使纵向升阻比减小,降低终端速度Vf。因此,Vf可视为关于EBR2的单调函数。为保证终端约束,这个问题可视为求解非线性方程Vf(EBR2)=VTAEM的解的问题。这里采用割线法进行求解,通过多个弹道仿真对Vf进行预测,并根据预测值与终端约束的偏差调整EBR2的值。The height adjustment stage is designed with proportional guidance law for guidance. Through the analysis of the guidance algorithm, it can be seen that if the E BR2 is reduced, the guidance law command will generate a larger pitch angle to eliminate the larger heading error generated during the last roll reversal, thereby reducing the longitudinal lift-to-drag ratio. , reducing the terminal velocity V f . Therefore, V f can be regarded as a monotonic function with respect to E BR2 . To ensure terminal constraints, this problem can be viewed as a problem of solving the solution of the nonlinear equation V f (E BR2 )=V TAEM . Here, the secant method is used to solve the problem, and V f is predicted through multiple ballistic simulations, and the value of E BR2 is adjusted according to the deviation of the predicted value from the terminal constraint.

Figure BDA0002408897380000241
Figure BDA0002408897380000241

EBR2结束迭代修正的条件为求出的偏差满足

Figure BDA0002408897380000242
E BR2 The condition for ending iterative correction is that the obtained deviation satisfies
Figure BDA0002408897380000242

5.基准倾侧角5. Reference tilt angle

在第二个倾侧反转以后,飞行器进入高度调整阶段。在此阶段,基准攻角仍由式(19)确定,基准倾侧角则由比例导引制导律控制。视线角速率可以表示为After the second roll reversal, the aircraft enters the altitude adjustment phase. At this stage, the reference angle of attack is still determined by equation (19), and the reference tilt angle is controlled by the proportional guidance guidance law. The line-of-sight angular rate can be expressed as

Figure BDA0002408897380000243
Figure BDA0002408897380000243

然后,通过比例导引获得的横向机动加速度指令如下Then, the lateral maneuvering acceleration command obtained by proportional guidance is as follows

Figure BDA0002408897380000244
Figure BDA0002408897380000244

其中,

Figure BDA0002408897380000245
Figure BDA0002408897380000246
表示第二次倾侧反转时的纵程。假设飞行器近似平稳滑翔,则纵向平面内的合力加速度为0。因此,有in,
Figure BDA0002408897380000245
Figure BDA0002408897380000246
Indicates the longitudinal distance at the second roll reversal. Assuming that the aircraft glides approximately smoothly, the resultant acceleration in the longitudinal plane is zero. Therefore, there is

Figure BDA0002408897380000247
Figure BDA0002408897380000247

所以可得基准倾侧角为Therefore, the reference tilt angle can be obtained as

Figure BDA0002408897380000248
Figure BDA0002408897380000248

6.指令攻角与倾侧角6. Command angle of attack and roll angle

为了确保终端速度约束,这里通过调整攻角来跟踪以能量为自变量的剩余滑翔距离参考剖面

Figure BDA0002408897380000249
指令攻角与倾侧角可以表示如下:In order to ensure the terminal speed constraint, the remaining glide distance reference profile with energy as the independent variable is tracked by adjusting the angle of attack.
Figure BDA0002408897380000249
The commanded attack angle and pitch angle can be expressed as follows:

Figure BDA00024088973800002410
Figure BDA00024088973800002410

Figure BDA00024088973800002411
Figure BDA00024088973800002411

其中,kα的值可根据任务及时调整,一般可选为2π/(1.8×107);

Figure BDA00024088973800002412
的值可以通过之前的最后一次弹道仿真获得。Among them, the value of k α can be adjusted in time according to the task, and generally can be selected as 2π/(1.8×10 7 );
Figure BDA00024088973800002412
The value of can be obtained from the last ballistic simulation before.

为了满足过程约束,像前文一样,倾侧角仍然需要满足下列要求:To satisfy the process constraints, as before, the roll angle still needs to satisfy the following requirements:

cmd|≤σmax (61)cmd |≤σ max (61)

通过上述五个步骤,可得到了一种基于平稳滑翔弹道解析解的高超声速再入制导方法。Through the above five steps, a hypersonic reentry guidance method based on the analytical solution of the smooth gliding trajectory can be obtained.

为了证明该制导方法的普适性,这里考虑了不同目标位置的5个算例。用T1到T5分别表示五个不同的目标点,它们分别设置为(-77°,0°)(T1),(-34°,-32.5°)(T2),(5°,-50°)(T3),(65°,-33°)(T4)以及(118°,0°)(T5)。初始条件均为h0=80km,V0=7200m/s,λ0=0°,θ0=45°以及γ0=0°。ψ0指向目标,由目标位置而决定。终端约束条件均为VTAEM=2000m/s,HTAEM=25km以及Rtm=100km。仿真结果见图3至图9。In order to prove the universality of the guidance method, five examples with different target positions are considered here. Denote five different target points with T1 to T5, which are set as (-77°, 0°) (T1), (-34°, -32.5°) (T2), (5°, -50°) (T3), (65°, -33°) (T4) and (118°, 0°) (T5). The initial conditions were all h 0 =80 km, V 0 =7200 m/s, λ 0 =0°, θ 0 =45° and γ 0 =0°. ψ 0 points to the target and is determined by the target position. The terminal constraints are all V TAEM =2000m/s, H TAEM =25km and R tm =100km. The simulation results are shown in Figure 3 to Figure 9.

在图3中,红色的圆圈表示距目标100km的范围。由图易见,所有的弹道都落在了圆圈上,证明该制导方法可以满足终端距离约束。表1展示了其他终端射程、速度、高度、倾侧角、弹道倾角和航向误差等状态量。In Figure 3, the red circle represents the range of 100km from the target. It is easy to see from the figure that all the ballistic trajectories fall on the circle, which proves that the guidance method can satisfy the terminal distance constraint. Table 1 shows other state quantities such as terminal range, speed, altitude, bank angle, ballistic inclination angle, and heading error.

表1仿真结果Table 1 Simulation results

Figure BDA0002408897380000251
Figure BDA0002408897380000251

从图4可以看出,除了两个向西飞的算例以外,所有的算例均满足了终端过程约束。这是因为当飞行器向西飞行时,所受科氏力方向向下,使得滑翔高度降低,热流密度增大。此外,根据图4-图9和表1,我们可以很容易地看出,所有示例均满足Vf≈VTAEM=2000m/s,Hf≈HTAEM=25km,Rtm=100km,σf≈0°,γf≈0°以及Δψf≈0°的终端约束。It can be seen from Figure 4 that all the examples satisfy the terminal process constraints except for the two examples of flying westward. This is because when the aircraft flies westward, the Coriolis force is directed downward, which reduces the gliding height and increases the heat flux density. Furthermore, according to Fig. 4-Fig. 9 and Table 1, we can easily see that all the examples satisfy V f ≈ V TAEM = 2000m/s, H f ≈ H TAEM = 25 km, R tm = 100 km, σ f ≈ Terminal constraints for 0°, γ f ≈ 0° and Δψ f ≈ 0°.

综上所述,通过上述步骤,推导出了本发明方法,即一种基于平稳滑翔弹道解析解的高超声速再入制导方法,案例仿真结果表明本发明方法能够满足再入制导所需的多种约束,且由于本方法基于弹道解析解设计,具有精度高、计算量小的特点,可实时应用到飞行器机载系统中。To sum up, through the above steps, the method of the present invention is deduced, that is, a hypersonic reentry guidance method based on the analytical solution of the smooth gliding trajectory. The simulation results of the case show that the method of the present invention can meet various requirements for reentry guidance. Constraints, and because the method is based on the ballistic analytical solution design, it has the characteristics of high precision and small amount of calculation, and can be applied to the airborne system of the 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 gammaffAnd
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 (16)

* 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 北京航空航天大学 A Guidance Method Applicable to Online No-Fly Zone Avoidance for Aircraft
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
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
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
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
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 北京航空航天大学 A design method for high Mach aircraft maneuver penetration strategy
CN114675545B (en) * 2022-05-26 2022-08-23 中国人民解放军火箭军工程大学 A Re-Entry Collaborative Guidance Method for Hypersonic Vehicles Based on Reinforcement Learning
CN114936531B (en) * 2022-05-31 2024-09-20 西北工业大学 Aircraft hitting area prediction system and method based on fuzzy neural network
CN115454133B (en) * 2022-09-30 2025-05-09 西北工业大学 A missile cascade correction guidance method and system under satellite-missile coordination conditions
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
CN107941087B (en) A reentry guidance method for high lift-drag ratio and ultra-smooth glide based on drag profile
CN108036676B (en) A kind of autonomous reentry guidance method of full directive based on three-dimensional resolution Value of Reentry Vehicle
US20210164783A1 (en) Method for directly planning reentry trajectory in height-velocity profile
US11286065B2 (en) Method for designing reentry trajectory based on flight path angle planning
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
CN111591470B (en) Aircraft precise soft landing closed-loop guidance method adapting to thrust adjustable mode
CN103245257B (en) Guidance law of multi-constraint aircraft based on Bezier curve
CN106054604B (en) Reentry vehicle robust optimal method of guidance based on Model Predictive Control Theory
CN110750850B (en) Three-dimensional profile optimization design method, system and medium under strong constraint complex task condition
CN107121929B (en) Robust Reentry Guidance Method Based on Linear Covariance Model Predictive Control
CN106643341B (en) Power thermal control coupling design method based on quasi-equilibrium gliding principle
Wang et al. Terminal guidance for a hypersonic vehicle with impact time control
CN113093790B (en) Analytical model-based aircraft reentry glide trajectory planning method
CN112256061A (en) Reentry guidance method for hypersonic aircraft under complex environment and task constraint
CN106444822A (en) Space vector field guidance based stratospheric airship&#39;s trajectory tracking control method
Yang et al. Robust entry guidance using multi-segment linear pseudospectral model predictive control
CN115454145B (en) High-speed aircraft cooperative guidance method and system under uncontrollable flight speed conditions
Zhang et al. Small fixed-wing unmanned aerial vehicle path following under low altitude wind shear disturbance
CN113741509B (en) Hypersonic gliding aircraft hold-down section energy management method
CN113093789B (en) Planning method for avoiding trajectory of aircraft no-fly zone based on path point optimization
CN110232215A (en) Method, system and medium for three-dimensional profile layered iterative planning considering maneuvering task requirements
CN115268501B (en) Multi-aircraft collaborative reentry trajectory planning method, system, electronic equipment and medium
CN116560403A (en) Intelligent time collaborative guidance method, system and equipment for hypersonic aircraft

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