气动参数在线辨识在火箭助推段制导中的应用

康继芝， 张士峰， 胡铖

(国防科技大学 空天科学学院，湖南 长沙 410073)

助推段主要在大气层中，气动力是主要作用力之一。精确地计算火箭受到的气动力，对提高火箭制导精度至关重要。由于气动环境的复杂多变，很难直接建立气动力精确计算模型。目前使用的解决办法是针对特定型号火箭外型，结合空气动力学理论计算和风洞实验校正，给出计算气动力所需的图表、曲线等[1]。工程上，通常提供气动力系数的曲线或插值表，以此来计算火箭飞行过程中的气动力。但是实际飞行过程气动环境复杂，标准气动参数与实际参数存在偏差。这种偏差在制导过程中如果没有处理，一定程度上会影响制导精度，最终增大终端误差。因此，利用辨识算法对气动参数进行实时在线辨识，并将辨识结果用于修正制导指令，对于提高火箭的制导精度有一定研究意义。

摄动制导方法要求在射前计算参考弹道并装载，在飞行过程中将实际弹道参数与参考弹道的偏离量视为小扰动量，再基于摄动理论计算制导指令，使真实弹道始终保持在标准弹道附近[2]。由于运动参数具有可观测性，摄动制导一般选择运动参数扰动(如位置坐标、速度分量、角度值等)作为偏差量来计算制导指令。但是火箭真实飞行过程中的参数偏差量远不止这些运动参数[1]。本文拟将辨识得到的气动系数偏差量视为小扰动量，对传统摄动制导方法进行改进，提高制导精度。

现有的参数估计方法包括递推最小二乘法、递推极大似然法、Wiener滤波和Kalman滤波等，其中Kalman滤波是处理随机估计应用最广泛的方法[3]。为扩展适用范围，Kalman滤波出现了多种改进形式，包括扩展Kalman滤波、无迹Kalman滤波[4-5]、UD分解滤波[6]等。在飞行器气动参数辨识上，文献[7-8]分别应用了具有指数遗忘的最小二乘法和径向基神经网络方法分别对滑翔飞行器和旋翼飞行器进行气动参数辨识，达到较好效果。文献[9-12]研究了扩展Kalman滤波器在飞行器气动参数辨识上的应用，使得扩展Kalman滤波方法得到了广泛应用。本文将基于文献[9]中的基本思想，构建气动参数在线辨识模型，用于修正火箭助推段摄动制导指令，提高末端精度。

1 动力学模型与摄动制导

1.1 火箭助推段动力学模型

为了便于后续辨识系统方程的推导，本文拟引入速度坐标系下的动力学方程[13]。由于火箭助推段飞行时间较短，可假设地球为不旋转均质圆球；并假设火箭为对称圆柱体形状，不考虑倾侧角控制。故在速度坐标系下的动力学方程简化为：

(1)

式中：v表示速度大小；Pe表示轴向推力大小；α为攻角；β为侧滑角；ρ为大气密度；Sm为火箭参考面积；CD、CL为阻力系数和升力系数；θ为速度倾角；Re为地球平均半径；h为飞行高度；m为火箭质量。CD、CL均由标准值与扰动量组成：

(2)

式中：CDd、CLd分别为阻力系数和升力系数的扰动量，该扰动量正是本文中要辨识的参量。

1.2 摄动制导模型

摄动制导的目的是将火箭的实际飞行轨迹保持在标准弹道附近。在实际飞行过程中，导航系统计算得到火箭实际位置，制导系统将其与存储在箭载计算机中的标准弹道进行对比，经制导方程计算出制导指令。为便于执行机构设计，摄动制导通常分为纵向制导和横向导引，本文仅考虑纵平面内的制导。不考虑姿态变化的动态过程，制导指令俯仰角可简化为[14]：

φ(t)=φpr(t)+uφ

(3)

式中：φpr(t)为t时刻对应的程序俯仰角；uφ为导引控制量，由真实状态和参考状态的偏差计算得到。导引控制量由导引目标量和导引系数组成，即uφ=KδXδXf。本文选择终端速度倾角偏差δθf作为导引目标量。摄动制导通常选择可测量的状态量偏差来计算导引目标量：

(4)

但是根据摄动理论，“扰动”定义为实际弹道飞行条件与标准弹道飞行条件之间的偏差[1]。因此，扰动量除了起始点位置速度外，还包括气温、压强、大气密度等环境参数和质量、秒耗量、气动力系数等弹道参数。本文拟将辨识得到的气动参数偏差作为新的小扰动量。在不考虑气动参数偏差的情况下，终端速度倾角θf是速度大小v、速度倾角θ和高度h的函数，即：

θf=θf(v,θ,h)

(5)

故其摄动方程如式(4)所示。当考虑气动参数扰动时，式(5)可写为：

θf=θf(v,θ,h,CD,CL)

(6)

对式(6)求微分：

(7)

式中:δCD、δCL是气动参数的等时偏差;在标准模型中，CD是马赫数Ma和高度h的函数，CL是马赫数Ma和攻角α的函数，即：

(8)

而Ma=v/a，声速a又是高度的函数。攻角α=φpr-θ，是速度倾角θ的函数。所以在标准情况下，气动参数的等时偏差δCD、δCL可化为：

(9)

式中：

(10)

由式(4)可知，标准状况下气动参数的等时偏差可表示为速度大小、高度和速度倾角的等时偏差，将式(9)代入式(4)，得到的结果和式(4)等价。但是考虑气动参数扰动时，气动参数等时偏差除了来源于速度大小、高度和速度倾角等时偏差外，还来源于其他未知因素引起的气动参数扰动。结合式(9)可扩展为：

(11)

将式(11)代入式(7)中，相同的等时偏差项可以合并，得：

(12)

而标准弹道的气动参数扰动量均为0，所以式(12)中δCDd=CDd，δCLd=CLd。式化(12)为：

(13)

式(13)就是考虑气动参数偏差的摄动方程。

2 气动参数在线辨识预测模型

2.1 扩展Kalman滤波辨识模型

CDd、CLd为气动参数扰动量，具有一定的随机性，难以用精确模型进行描述。而一阶高斯马尔可夫过程符合大多数物理过程，且数学描述较为简单。因此采用一阶高斯马尔可夫过程对气动参数扰动进行描述[9]：

(14)

Y=h(X)+V

(15)

式中：Y为观测量；V为观测噪声向量。火箭惯性敏感器件可以测量除地球引力外的其余力产生的加速度。因此观测量选择为轴向和法向的视加速度，即Y=[Nx,Ny]T。其中：

(16)

只考虑加速度计的随机漂移，故ωNx、ωNy为高斯白噪声。

扩展Kalman滤波的基本原理如下：

1)预测方程：

(17)

式中：Qk+1为系统噪声矩阵；Φk,k+1为tk到tk+1时刻的状态转移矩阵，其表达式为：

(18)

式中:F(X)为状态方程f(X)的雅克比矩阵。

2)更新方程：

(19)

大多数的人，第一次接触武侠，都是在很年轻的时候。这也许是因为，年轻的人，甚至有点中二的人，总是血气方刚的，于是就更容易受到武侠小说里那种自由、热血、张扬和少年意气的影响。

2.2 限定记忆递推最小二乘法

扩展Kalman滤波方法能辨识出当前时刻气动参数扰动量，但是在计算下一时间段内的制导指令时，需要下一时间段内的气动参数扰动量。因此，需要通过已经辨识出来的结果对下一时间段内的气动参数扰动量进行预测。由于气动参数扰动量具有一定的随机性，其与其他状态量的函数关系不明确。因此，建立扰动量的时变多项式模型，再利用参数估计方法对多项式参数进行估计：

(20)

参数估计方法有递推最小二乘法、遗忘因子递推最小二乘法、限定记忆递推最小二乘法等[15]。递推最小二乘法虽然算法简单，但是随着观测数据的增多，会出现“数据饱和”现象，使新数据对估计结果的修正作用减弱。遗忘因子递推最小二乘法通过对旧数据增加“遗忘因子”以减弱数据的作用，而限定记忆递推最小二乘法通过不断地剔除老数据，增加新数据，始终选择固定长度的最新数据(即“数据矩形窗”)进行参数估计。限定记忆递推最小二乘法在一定程度上完全避免了数据饱和现象，更充分地体现了新数据的修正作用。但是为了保证预测结果的合理性，每次预测计算采用的数据量不能过少。本文选择限定记忆递推最小二乘法作为参数估计方法[16]。该方法的计算过程：

(21)

2.3 气动参数在线辨识模型

参数在线辨识的关键是实时性，即辨识结果及时应用于下一阶段的制导计算。但在实际飞行过程中，由于存在计算时延，辨识结果不能保持时刻更新。因此，本文考虑采用周期性迭代的方法，来实现在线气动参数实时更新：首先基于实时观测数据，采用扩展Kalman滤波方法对气动参数扰动进行辨识；然后再基于已有的辨识结果，采用限定记忆递推最小二乘法，对下一制导周期的气动参数扰动量进行预测，将预测结果用于下一制导周期的计算。依次循环，便构建了“辨识—预测—计算”迭代更新的在线辨识模型，如图1所示。

图1 气动参数在线辨识在摄动制导中的应用流程Fig.1 Flow chart of the aerodynamic parameter online identification in perturbation guidance

3 仿真结果与分析

本文将通过数值仿真对以上方法进行验证。仿真条件如下：对象是某型号的运载火箭，飞行阶段为起飞至一级分离，起飞方式为垂直起飞，状态初值为：[v0,θ0,h0]=[0,π/2,0]，飞行标准程序如图2所示。

图2 数值仿真标准弹道程序俯仰角Fig.2 Program pitch angle in numerical simulation

3.1 扩展Kalman滤波辨识效果

为了验证本文设计扩展Kalman滤波的有效性，基于上述对象，设计算例对Kalman滤波器的辨识效果进行验证。

在滤波参数设计中，假设测量噪声仅存在随机漂移，且保持常值：

Rk=diag(10-6,10-6)

(22)

Qk=diag(0,0,0,10-4,10-4,10-10,10-10)

(23)

系统状态初值：

[v0,θ0,h0,CDd0,CLd0,ξD0,ξD0]=

[0,π/2,0,0.001,0.001,0,0]

状态误差协方差矩阵的初值为：

P0=diag(0,0,0,10-2,10-2,10-6,10-6)

(24)

分别设置扰动量真实值为0.1、-0.1，利用本文中的扩展Kalman滤波器进行辨识，辨识结果如图3所示。从图3可以发现，滤波器的收敛时间约为10 s，收敛速度较快；辨识误差保持在10%以内，精度较高。

图3 扩展Kalman滤波器辨识曲线Fig.3 The performance curves of extended Kalman filter

3.2 气动参数辨识在摄动制导中的应用

为了验证基于气动参数在线辨识的摄动制导方案的有效性，设计仿真对比实验。分别设置多组真实扰动值，计算传统摄动制导方案和基于气动参数在线辨识的摄动制导方案的终端状态偏差。如表1所示。

表1 气动参数扰动下的终端状态偏差Table 1 Deviation of terminal state under the aerodynamic parameter disturbance

对表1中的结果进行分析，新的摄动制导方案相对于传统摄动制导方案，终端速度倾角偏差和高度偏差有所降低，终端速度偏差改变不大。从能量的角度分析，气动参数扰动量直接影响气动力的大小，进而改变气动力消耗的能量大小。而终端速度直接与终端能量相关，而仅靠制导方法是无法完全弥补能量损耗的，因此新的摄动制导方案对于终端速度偏差没有明显的改善效果。综合以上分析，基于气动参数在线辨识的摄动制导方案相对于传统摄动制导，可以一定程度上减小终端偏差，提高制导精度。

为了进一步验证新的摄动制导方案的可靠性，设计Monte Carlo仿真实验。将气动参数扰动量CDd、CLd设置为随机扰动，均服从零值正态分布，3δ偏差CDd为0.15，CLd为0.015，进行500次Monte Carlo打靶，结果如图4、5所示。图4为终端高度-终端速度倾角剖面图，由图可知：基于气动参数在线辨识的摄动制导方法终端高度-速度倾角散布明显小于传统摄动制导：终端高度散布从约2 000 m减小到约1 000 m，终端速度倾角散布从0.016°减小到0.006°。图5为高度-时间曲线散布对比，由图可知基于气动参数辨识的摄动制导的高度-时间曲线散布明显小于传统摄动制导。综上所述，基于气动参数在线辨识的摄动制导方案不仅能提高终端精度，还能让真实弹道更好地保持在标准弹道附近。

图4 Monte Carlo打靶终端状态高度-速度倾角剖面图Fig.4 Terminal altitude-velocity profile in Monte Carlo simulation

图5 Monte Carlo打靶高度-时间曲线Fig.5 Altitude-velocity curves in Monte Carlo simulation

4 结论

1) 采用扩展Kalman滤波方法进行气动参数辨识，收敛较快，辨识精度较高；

2) 针对基于气动参数在线辨识的摄动制导方法，设计对比仿真，发现该方法相对于传统摄动制导方法，能有效降低终端状态偏差；并针对随机气动扰动，进行Monte Carlo打靶仿真，发现改进的摄动制导方案能有效地提高弹道收敛性，使弹道更好地保持在标准弹道附近。