两种卡尔曼滤波辨识高速旋转弹丸阻力系数的对比
2018-07-31郑宇程易文俊余春华
郑宇程，易文俊，余春华，管 军
(1.南京理工大学 瞬态物理国家重点实验室， 南京 210094； 2.南京理工大学 理学院， 南京 210094)
火炮在现代战争中具有重要的使命，在对敌火力压制过程中能够起到关键性作用。获取传统高速旋转稳定弹丸准确的气动参数，对于提高火炮射表精度、减小落点散布、增强打击精度具有重要的意义。目前，国内外获取弹丸气动参数的方法主要有三种[1]：第一种是理论计算，第二种是风洞吹风，第三种是利用弹丸的自由飞行数据对弹丸的气动参数进行离线辨识。理论计算法虽然简单，但由于模型中的未建模因素和不确定因素导致计算结果存在误差。风洞吹风法作用于弹丸模型，结果较为准确但成本较高，不能精准模拟高速旋转状态。利用弹丸自由飞行数据辨识弹丸的气动参数，不仅符合实际情况，还能根据辨识结果，及时对弹丸运动状态进行调整，从而提高炮弹的打击精度。
用于参数辨识的方法通常有递推最小二乘法、递推极大似然法、卡尔曼(kalman)滤波法等[2-5]。近年来，许多学者对此进行了研究。史继刚[2]基于粒子群初值选取的牛顿迭代优化算法辨识弹丸的零升阻力系数，该方法基于最大似然准则，相比于卡尔曼滤波法实现更加困难；管军等[3]提出一种新的自适应混沌变异粒子群算法求解该准则下的气动参数最优解，得到弹丸的气动参数，但是在工程上难以实现；Rogers等[4]提出了一种基于证据理论的参数估计方法，偏于理论计算；史金光等[6]利用扩展卡尔曼滤波法对弹道修正弹的阻力和升力符合系数进行了辨识，并且由此对后续弹道进行修正，然而该方法技术要求较高，难以在实际应用中实现；杨靖等[7]利用改进的混杂扩展卡尔曼滤波法对弹丸的阻力系数进行辨识，其对EKF的改进不大。
相比于最小二乘法和极大似然法，卡尔曼滤波法应用最为广泛。在处理非线性问题上，扩展性卡尔曼滤波(EKF)和无迹卡尔曼滤波(UKF)已经成功地应用于各种参数辨识问题[8-11]。本文主要利用扩展卡尔曼滤波和无迹卡尔曼滤波对弹丸的阻力系数进行辨识，并进行了对比研究。为了进行气动参数的辨识，本文将未知的气动参数看成系统的状态量，分别用EKF方法和UKF方法对气动参数进行辨识。为了检验实验的正确性，又分为气动参数不变和气动参数变化两种情况进行研究。
1 弹道模型
参数辨识的本质问题就是如何通过输出数据准确辨识相关输入参数或系统参数。如何选取系统的模型是重要的步骤之一。
阻力系数是传统无控旋转弹丸所有气动参数中最重要的气动参数之一，其直接影响了弹丸的射程和落点散布。本文的主要研究内容是弹丸阻力系数辨识，虽然外弹道学中的弹道模型多种多样，但2D弹道模型已经能够描述阻力系数对弹道性能的绝大部分影响。且2D弹道模型简单、便捷，能够直接揭示阻力系数和弹道参数之间的关系。因此，本文选用2D弹道模型作为系统的模型。
2D弹道模型[12]的方程为
(1)
其中，ρ为当前高度的大气密度，S为参考面积，m为弹体质量。
2 卡尔曼滤波原理
2.1 扩展卡尔曼滤波算法
卡尔曼滤波法最初只用于离散线性系统，后来经过改进，产生了许多种不同的方法。扩展性卡尔曼滤波(EKF)和无迹卡尔曼滤波(UKF)是应用得比较广泛的两种方法。EKF方法在线性化处理的过程中对原系统进行泰勒展开，略去了二阶以及二阶以上的高阶项，只保留了线性项。
扩展卡尔曼滤波的基本原理为：
1) 状态一步预测方程
(2)
2) 求一步预测均方误差
(3)
式中：Qk-1是系统噪声的协方差矩阵，Γk-1是系统噪声驱动阵，Φk,k-1为一步转移矩阵，计算式如下：
Φk,k-1=I+A·Δt
(4)
3) 求滤波增益矩阵
(5)
4) 估计均方误差
Pk=[I-KkHk]Pk,k-1
(6)
5) 得到最新一步的状态估计
(7)
式中，zk为k时刻的观测向量。
2.2 无迹卡尔曼滤波算法
由于EKF是一种非线性估计问题的线性化方法，其线性化过程会产生一定的误差，且其需要进行复杂的雅克比矩阵求解；而UKF通过UT变换，直接对非线性问题进行处理，不需进行复杂的雅克比矩阵求解，且其精度可以达到二阶泰勒展开的精度。因此，本文同样使用UKF对阻力系数进行辨识，并同EKF进行比较研究，其基本计算步骤为：
1) 初始化滤波参数
(8)
式中：n为增广状态向量的维数；α是很小的一个正数，其取值范围为10-4～1。本文取为0.01；κ的取值为3-n；β的值与x的分布形式有关，对于正态分布的情况，β取2最优。
2) 计算k-1时刻的sigma样本点
(9)
3) 计算k时刻的一步预测值
(10)
4) 计算k时刻的一步预测样本点
(11)
5) 计算P(XZ)k,k-1和P(ZZ)k,k-1
(12)
6) 更新滤波值
(13)
3 基于卡尔曼滤波的阻力系数辨识
由于待辨识参数是阻力系数CD，本文通过增广的方法将其添加到式(1)，以建立包含待辨识参数的状态方程。
利用系数冻结法，设在小区间内弹丸的阻力系数是常数，因此有下式成立：
(14)
则增广状态变量x为
x=[u,θ,x,y,CD]T
(15)
增广之后的状态方程可写为
(16)
ω为系统噪声向量。
通常情况下，弹道跟踪雷达只能得到弹丸的速度和位置信息，因此量测向量取为
z=[u,x,y]T
(17)
量测方程为
z=h(x)+v
(18)
v为观测噪声向量。
4 数值仿真
弹丸的物理参数为：弹重45.5 kg；弹径0.155 m；弹长0.909 m。气象条件为标准气象条件[13]。射击初始值为：初速930 m/s，射角45°。
状态初值为
x0=[930,45°,0,0,CD]Τ
(19)
初始状态误差协方差矩阵为
(20)
系统噪声的协方差矩阵为
(21)
量测噪声的协方差矩阵为
(22)
算例1
假设CD为常数，其真实值为0.35。用CD为0.35计算出来的理论值(速度，位置坐标)加上高斯白噪声v来模拟观测值。对于CD为常数的情况，用EKF和UKF辨识结果如图1～图2所示。从图1和图2可以看出，两种方法辨识的CD最终都收敛在0.35，但UKF具有更快的收敛速度。
算例2
实际情况中，阻力系数CD是马赫数的函数，本算例所采用的数值如表1所示。
根据变化的CD计算出的理论值，加上高斯白噪声v的值，将其用来模拟观测值。将EKF与UKF辨识的马赫数与真实的马赫数之差画成曲线在一起对比，如图3所示。
图1 EKF辨识CD不变
图2 UKF辨识CD不变
Ma0.50.751.01.251.51.75CD0.159 20.1570.299 50.321 40.296 20.275 7Ma2.02.252.52.753.0CD0.2580.239 90.225 70.207 90.191 5
图3 马赫数辨识误差曲线
图3不但表示了CD随马赫数变化的情况，而且可看出UKF辨识的马赫数误差较小，EKF辨识的马赫数误差较大。但是最大的误差不超过10-3，为原马赫数的0.05%，因此认为两种方法辨识的马赫数是准确的。CD虽然是马赫数的函数，但是为了便于比较研究，将真实CD曲线、EKF辨识CD曲线和UKF辨识CD曲线对比，如图4所示。
从图4可以看出，前后两端三者差别不大，现在选取中间20～50 s的时间段，参见图5。
图4 EKF与UKF辨识结果对比(1)
图5 EKF与UKF辨识结果对比(2)
图5中，从上至下分别是真实CD曲线、UKF辨识曲线和EKF辨识曲线。经过比较，UKF辨识曲线的精度较EKF高，两者的差距在0.002 5左右。UKF辨识的误差为0.8%，EKF辨识的误差为1.6%。由此可见，参数辨识的结果具有可行性。
5 结论
1)CD不变的情况下，EKF方法与UKF方法都能快速收敛到真值，且UKF具有较快的收敛速度。
2)CD随马赫数变化的情况下，EKF方法与UKF方法都能有效辨识出阻力系数，且UKF较EKF具有更高的辨识精度。
