基于阻力加速度倒数剖面的再入轨迹规划与制导方法

2019-01-18黄汉斌梁禄扬杨业

航空学报 2018年12期

黄汉斌，梁禄扬，杨业

北京航天自动控制研究所,北京 100854

以航天飞机[1-2]为代表的升力式飞行器具有巨大优势，高升阻比的气动布局，使其可以充分利用气动力实现大航程和高精度飞行任务。再入段[3]通常是约束最复杂、耗时最长、飞行距离最远的飞行段，再入飞行轨迹规划与制导算法是实现再入精确飞行的核心技术[3-4]。因此，再入轨迹规划与制导算法受到广泛关注。

升力式飞行器再入轨迹规划与制导算法[4]，分为标准轨迹跟踪制导和数值预测校正制导两种，其中数值预测校正制导算法需要迭代数值积分，由于计算效率不高，不能在线应用，而传统的标准轨迹制导算法基于离线规划的轨迹，不够灵活，因此需要研究在线轨迹规划与制导算法。由于再入过程中阻力加速度可以准确测量，同时可以克服模型等多种误差带来的影响，因此采用阻力加速度剖面作为标准轨迹[3]具有一定的优势。在生成阻力加速度剖面时，需要满足多种约束，因此需要建立阻力加速度约束走廊，同时由于再入终端时间不确定，采用能量作为阻力加速度自变量[5]，描述阻力加速度走廊，将高度和速度信息包含在能量中。为了获得满足待飞航程要求的阻力加速度-能量剖面，文献[5-11]采用多段线性函数将阻力加速度剖面设定为参数化的形式，进而利用优化的方法对参数进行求解，以得到标准轨迹；但线性函数转折点不光滑，不利于跟踪制导，文献[12] 则重新定义能量高度，采用分段三次函数描述剖面，利用简单实时积分的方法迭代获得光滑标准轨迹。阻力加速度标准轨迹必须满足阻力加速度走廊约束，文献[13-14]直接将阻力加速度剖面设定为走廊上下界的函数来进行优化，以使得阻力加速度剖面始终在走廊内。为了加快优化速度，简化优化过程，Mease等[15]采用降阶的动力学模型，在进行纵向阻力加速度剖面设计时，采用三段线性阻力加速度剖面进行优化，该方法在一定程度上加快了优化速度。

获得标准轨迹之后，Dukeman[16]以能量作为自变量，建立小扰动线性化模型，通过反馈航程、高度和航迹倾角的偏差设计LQR控制器对阻力加速度进行跟踪；Hanson[17]、Leavitt[18]等利用阻力加速度的动态误差特性，根据标准轨迹和实际阻力加速度求出制导律，实现纵向制导。横向制导，通过改变倾侧角符号[16-18]将航向角误差限制在走廊内。

标准轨迹生成时，优化的过程通常利用数值积分计算待飞航程，且初值随机选取，因而计算时间较长，Yang等[19]将阻力加速度-能量剖面分成多段，分段采用三次样条函数拟合，然后利用卡尔丹分解，推导出待飞航程的解析计算公式，大大提高了计算效率，但是卡尔丹分解推导过程复杂，本文对该解析航程预测方法进行了改进，消除卡尔丹分解过程，同时，基于该方法，提出了一种满足终端和过程约束的阻力加速度轨迹在线生成与更新的新方法，并将其与跟踪制导相结合，完成了再入轨迹在线规划和跟踪制导的闭环仿真。

本文首先利用过程约束、平衡滑翔约束建立阻力加速度倒数-能量走廊，将飞行器再入终端约束加入走廊中，选取阻力加速度倒数上下界的多个点，用三次样条函数对其进行拟合。然后，利用简化的航程计算公式，得到待飞航程的上下界，当待飞航程在约束走廊对应的航程上下界以内时，满足待飞航程和约束的阻力加速度倒数剖面一定在走廊上下界之间，因此，阻力加速度剖面可以利用待飞航程以及待飞航程的上下界的近似线性关系生成阻力加速度剖面，并利用二分法对剖面进一步优化，从而得到阻力加速度剖面，之后根据待飞航程偏差，利用当前剖面来更新阻力加速度剖面。最后，利用阻力加速度的动态误差特性建模，对其进行跟踪。横侧向制导根据航向角偏差，改变倾侧角符号，来消除横向误差。

1 运动模型

1.1 有量纲运动模型

假设地球为匀质圆球，忽略地球自转，微分方程的6个变量为高度H，经度λ，纬度φ，速度Vd，航迹倾角γ，航向角ψ，建立式(1)中以时间t为自变量的运动模型。其中，re为地球半径，g为重力加速度，控制量为倾侧角σ和攻角α，升力加速度和阻力加速度采用式(2)和式(3)计算：

(1)

(2)

(3)

式中：m、Sref分别为飞行器质量和参考面积；ρ为大气密度；CL、CD分别为升力系数和阻力系数。

1.2 阻力加速度倒数约束

(4)

(5)

(6)

(7)

由于飞行器初始速度、高度，以及终端速度、高度给定，则初始能量、阻力加速度，以及终端能量、阻力加速度即给定，阻力加速度走廊需要满足初始和终端约束，因此需要对阻力加速度走廊进行处理。

利用同样的方法对倒数第2个能量点进行处理，得到图1中由N个节点(红色点)组成的阻力加速度能量的上下界。

为了简化计算待飞航程，将N个节点构成的无量纲阻力加速度走廊转换成图2中N个无量纲阻力加速度倒数节点三次样条插值后的走廊。得到式(8)中三次样条函数描述的无量纲阻力加速度倒数约束走廊：

(8)

图1 无量纲阻力加速度走廊Fig.1 Dimensionless drag acceleration corridor

图2 无量纲阻力加速度倒数走廊Fig.2 Inverse dimensionless drag acceleration corridor

2 再入制导算法

基于阻力加速度倒数剖面的标准轨迹跟踪制导算法分为2个步骤：① 生成阻力加速度标准轨迹；② 计算跟踪制导律。其中在线生成并更新满足各种约束的阻力加速度剖面是本文算法的核心。本文算法结构框图如图3所示。

图3 本文算法结构框图Fig.3 Structure diagram of proposed algorithm

2.1 在线轨迹生成与更新

2.1.1 待飞航程预测

标准轨迹制导算法的核心之一在于轨迹的生成，轨迹生成方法主要依据一般是待飞航程的预测，待飞航程定义为当前点地心矢径到目标点地心矢径在地球表面所构成的大圆弧的长度，计算公式为

(9)

在平衡滑翔条件下航迹倾角γ较小，cosγ≈1，同时由于航向角偏差Δψ较小，cos Δψ≈1，采用能量作为待飞航程微分公式的自变量，则可以推导出待飞航程计算公式为

(10)

则将式(8)代入式(10)，可以得到待飞航程上下界的计算公式为

Smax=

(11)

Smin=

(12)

2.1.2 阻力加速度剖面在线解析生成算法

利用阻力加速度倒数剖面，通过式(10)计算待飞航程，则满足待飞航程要求的阻力加速度倒数剖面与阻力加速度倒数上下界近似成线性关系，则可以利用这种关系，首先得到一条相近的阻力加速度剖面，然后利用二分法得到精确的阻力加速度剖面。

待飞航程与其上下界对应的比例系数为

(13)

因此，利用式(14)得到近似满足待飞航程的阻力加速度倒数剖面上的N个节点。

然后式(15)利用三次样条函数插值，得到倒数剖面。

i=0,1,…,N-1

(14)

i=1,2,…,N-1

(15)

对与上下边界有比例系数η的节点进行插值后，节点之间的点不一定成比例系数，利用式(11)或者式(12)估计的结果不一定等于待飞航程，因此，需要进一步迭代求取剖面，这里采用二分法进行迭代，迭代过程如图4所示。

经过迭代，得到满足待飞航程要求的阻力加速度倒数剖面，该剖面与阻力加速度倒数上下界仍然成线性关系，该比例系数仍然采用η描述。得到图5中的阻力加速度倒数剖面。

同时，阻力加速度倒数剖面可以无损转化为阻力加速度剖面，式(16)给出转化方式，由此便得到了满足航程要求的阻力加速度剖面。

图4 倒数剖面迭代过程Fig.4 Reciprocal profile iteration process

图5 无量纲阻力加速度倒数剖面Fig.5 Inverse dimensionless drag acceleration profile

i=1,2,…,N-1

(16)

2.1.3 剖面更新算法

阻力加速度剖面依据当前点到目标点所构成的纵向平面进行规划，由于控制量倾侧角以及模型误差和气动误差的存在，飞行器存在横向机动来消耗能量使得飞行器跟踪阻力加速度剖面，所以飞行器最终不能够到达终点(图6)，因此需要在一定的周期下对阻力加速度剖面进行更新。

由式(10)～式(12)可知，待飞航程在阻力加速度倒数走廊内近似成线性分布，因此可以根据待飞航程的偏差更新η，从而得到新的阻力加速度剖面。由当前飞行器位置到目标之间的待飞航程S(e)，同时，根据阻力加速度-能量剖面以及待飞航程计算公式，可以得到当前能量点处，沿规划的阻力加速度剖面的剩余待飞航程Sego(e)，由于线性关系的存在，令

(17)

则满足剩余待飞航程的阻力加速度倒数剖面与上下界之间的比例系数为

η(j+1)=η(j)+Δη

(18)

由式(14)、式(15)和式(18)，可得到满足待飞航程和约束的新的阻力加速度剖面。

图6 阻力加速度剖面Fig.6 Drag acceleration profile

2.2 阻力加速度跟踪制导

2.2.1 纵向制导

纵向制导对规划的阻力加速度剖面进行跟踪，再入过程中通常采用与速度相关的固定的攻角剖面，仅通过倾侧角进行跟踪制导，利用阻力加速度动态误差特性推导制导律[14]。动态误差特性公式为

(19)

式中：ξ为阻尼系数；ωn为自然频率；k1为积分系数，实际阻力加速度一阶导D′和D″二阶导采用文献[18]得到的推导公式，代入式(19)，得到纵向制导律的表达式为

|σ|=

(20)

(21)

2.2.2 横向制导

侧向制导旨在消除侧向误差，不去跟踪某一标准轨迹，因此，仅需要将航向角偏差的限定在以速度为函数的范围内，即当航向角偏差超出边界，改变倾侧角符号来实现横向制导。

航向角偏差的计算公式为

Δψ=ψ-az

(22)

式中：ψ为当前航向角；az为目标点关于当前点的方位角，计算公式如式(23)所示:

(23)

式中：λf为终端时刻经度；φf为终端时刻纬度。

航向角偏差的上下界±Δψup，采用插值函数给出，如图7和表1所示，当航向角偏差超出边界，倾侧角反向，使得航向角偏差减小。应当注意的是，由于倾侧角反向时，有速度限制，并不能立即使得航向角偏差减小，因此在设计航向角偏差边界时，要留有一定的余量。

图7 航向角偏差走廊Fig.7 Azimuth deviation corridor

3 仿真结果

3.1 仿真模型及环境

对纵向制导律即倾侧角的大小，设定上界，|σ|≤90°，倾侧角反向速度小于15 (°)/s，攻角剖面采用图8中与速度相关的固定剖面。

仿真环境：Intel Core i7-4790、Windows XP系统、MATLAB R2013b。

图8 攻角剖面Fig.8 Angle of attack profile

表1 航向角误差上下界插值表Table 1 Azimuth deviation corridor interpolation table

Vd/(km·s-1)Δψup/(°)87775144162140.64

3.2 标准任务仿真

利用表2中给出的标准任务下飞行器再入初始和终端要求，根据表中条件可以算出初始和终端能量与阻力加速度。

仿真过程中采用的阻力加速度剖面更新周期为1 s，再入段终点为终端能量管理段的起点，即当距离目标小于终端能量管理段航程STAEM时，进入终端能量管理段。选取STAEM=100 km，式(24)为待飞航程计算公式：

S=re·arccos[sinφsinφf+

cosφcosφfcos(λf-λ)]-STAEM

(24)

表2 标准任务初始和终端要求Table 2 Standard task initial and end requirements

仿真结果显示，首次轨迹规划周期小于0.023 s，每更新一次轨迹所需时间小于0.005 s，(因仿真环境不同，存在一定差异)，而以文献[19]中未改进的待飞航程预测方法为基础来生成轨迹的时间多次仿真均大于0.09 s，由于本文方法省去了卡尔丹分解过程，减少了待飞航程计算时的矩阵求逆和解方程的次数，轨迹生成时间得到明显减少，满足工程在线应用要求。

图9给出了阻力加速度走廊以及阻力加速度指令和实际阻力加速度的曲线，起始段有一个小的波动，是由于阻力加速度指令变大，而高空大气密度较小，初始航迹倾角不理想，因此跟踪阻力加速度指令，需要调整过程；在能量大约为0.92时，实际阻力加速度有一个小的波动，这是由于在此时速度Vd=3 km/s，攻角指令在此处非线性，对于升力系数和阻力系数产生影响而造成的，调整由非线性造成的影响需要一定的时间。

图10为制导给出的倾侧角指令，初始倾侧角直接达到最大值，为使得实际阻力加速度尽快调节到与跟踪指令一致，倾侧角直接达到了最大限幅值。倾侧角反向，是由于航向角达到了调整的边界，在横向制导律的作用下，倾侧角改变符号。