技术领域
[0001] 本
发明涉及参数辨识技术领域，特别是一种基于扩展卡尔曼滤波的高速旋转弹气动参数辨识方法。
背景技术
[0002] 火炮武器在战场中反应迅速，能够提供强大且持续的火
力支援，是未来部队联合作战的重要组成单元。获取传统高速旋转稳定弹丸准确的气动参数，对于提高火炮射表
精度、减小落点散布、增强打击精度具有重要的意义。获取弹丸气动参数的方法主要有三种：第一种方法通过理论计算得出气动参数，第二种方法采用
风洞吹风法，第三种方法利用弹丸的自由飞行数据对弹丸的气动参数进行离线辨识。其中，理论计算方法虽然简单，但由于模型中的未建模因素和不确定因素导致计算结果存在一定的误差；风洞吹风法作用于弹丸模型，结果较为准确但由于其成本较高，不能够精准地模拟高速旋转等状态，因而该方法也有
缺陷；利用弹丸自由飞行数据辨识弹丸的气动参数，不仅符合实际情况，还能根据辨识结果，及时对弹丸进行调整，从而提高炮弹的打击精度。
[0003] 用于参数辨识的方法通常有递推最小二乘法、递推极大似然法、卡尔曼滤波法等。史继刚基于粒子群初值选取的
牛顿
迭代优化
算法辨识弹丸的零升阻力系数，该方法基于最大似然准则，相比于卡尔曼滤波法实现更加困难；管军等提出一种新的自适应混沌变异粒子群算法来求解该准则下的气动参数最优解，进而得到弹丸的气动参数，但是其在工程上比较难以实现；Rogers等提出了一种基于证据理论的参数估计方法，偏于理论计算；史金光等利用扩展卡尔曼滤波法对弹道修正弹的阻力和升力符合系数进行了辨识，并且由此对后续弹道进行了修正，然而该方法要求较高，难以在实际应用中实现。
发明内容
[0004] 本发明提出了一种基于扩展卡尔曼滤波的高速旋转弹气动参数辨识方法。
[0005] 实现本发明的技术解决方案为：一种基于扩展卡尔曼滤波的高速旋转弹气动参数辨识方法，具体步骤为：
[0006] 步骤1、建立高速旋转弹四
自由度动力学模型；
[0007] 步骤2、根据高速旋转弹四自由度动力学模型，利用扩展卡尔曼滤波辨识气动参数，具体为：
[0008] 步骤2-1、将待辨识参数加入状态变量中形成增广状态向量，得到增广之后的
状态方程和量测方程；
[0009] 步骤2-2、根据实际工程的初始值以及增广之后的状态方程和量测方程进行滤波计算，得到每一时刻的增广状态变量的估计值以及待辨识参数估计值。
[0010] 进一步地，建立的高速旋转弹四自由度动力学模型具体为：
[0011]
[0012]
[0013]
[0014]
[0015]
[0016]
[0017]
[0018]
[0019] 式中，Vx、Vy、Vz分别为弹丸在x、y、z方向的速度，x、y、z为弹丸的
位置坐标；γ为
滚转角，为滚转
角速度；m为弹丸
质量；d为弹丸最大横截面直径；l为弹长；g为重力
加速度；ρ为空气
密度；S为特征面积，通常取为弹丸的最大横截面积；Cx0和Cx2分别为零升阻力系数和诱导阻力系数；αe为动力平衡角；Vr为弹丸相对于空气的速度；C′y为升力系数导数；C″z为
马氏力系数的联合偏导数；C为极
转动惯量；A为赤道转动惯量；m′xz为极阻尼力矩系数导数；Wx、Wz分别为x和z方向的风速；m″y为马格努斯力矩系数的联合偏导数；m′z为静力矩系数导数，Vr为弹丸相对于空气的速度。
[0020] 进一步地，步骤2-1中得到的增广之后的状态方程和量测方程为：
[0021]
[0022] y(t)＝ga[xa(t),u(t)]
[0023] z(k)＝y(k)+Gv(k)
[0024] 式中，F为系统噪声转移矩阵，G为量测噪声转移矩阵，w(t)为系统噪声，v(k)为量测噪声，xa(t)为增广状态变量，u(t)为外部输入，β为待辨识参数，y(t)为连续时间下的观测向量，z(k)为量测向量。
[0025] 进一步地，步骤2-2中得到每一时刻的增广状态变量的估计值以及待辨识参数估计值的具体步骤为：
[0026] 步骤2-2-1、确定初始时刻的增广状态变量估计值及增广状态变量误差协方差矩阵估计值；
[0027] 步骤2-2-2、根据k-1时刻的增广状态变量估计值及增广状态变量误差协方差矩阵估计值得到k时刻的增广状态变量预测值及增广状态变量误差协方差矩阵的预测值；
[0028] 优选地，k时刻的增广状态变量预测值及增广状态变量误差协方差矩阵的预测值具体为：
[0029]
[0030]
[0031] 表示k时刻的增广状态变量的预测值， 为外部输入， 表示连续时间t时刻下的估计的增广状态变量，t(k)表示对应着k时刻的时间， 为k时刻的增广状态变量误差协方差矩阵的预测值，Δt为k-1时刻到k时刻的时间间隔， 为k-1时刻的增广状态变量误差协方差矩阵的估计值，Φa(k)为状态转移矩阵，且 Aa(k)
为增广后的状态函数fa对xa的雅克比矩阵。
[0032] 步骤2-2-3、根据求得的k时刻的增广状态变量误差协方差矩阵预测值以及k时刻的增广状态变量的预测值得到k时刻的增广状态变量误差协方差矩阵估计值以及k时刻的增广状态变量的估计值；
[0033] 优选地，k时刻的增广状态变量误差协方差矩阵估计值以及k时刻的增广状态变量的估计值具体为：
[0034] 增广状态变量误差协方差矩阵估计值：
[0035]
[0036] k时刻的增广状态变量的估计值：
[0037]
[0038] 式中，u(k)为外部输入，I为单位矩阵，z(k)为k时刻的量测向量，Ka(k)为k时刻的滤波增益矩阵，Ca(k)为ga对xa的雅克比矩阵， 增广状态变量误差协方差矩阵预测值，G为噪声转移矩阵。
[0039] 步骤2-2-4、返回步骤2-2-2，直至量测向量完结。
[0040] 在实际运用中，根据实际情况
选定初始条件，进行参数辨识，运用本方法时，要先确定炮弹发射时的气象条件，对应不同型号的高速旋转弹，具有不同的物理参数。确定高速旋转弹发射的初速度、射角、初始滚转角速度，确定观测值和量测方程，选择计算机语言进行编程实现。
[0041] 本发明与
现有技术相比，其显著优点为：本发明采用扩展卡尔曼滤波进行参数辨识，基于实际项目测得的数据，具有实际性能；同时，用扩展卡尔曼滤波辨识的成本较低，精度更高；扩展卡尔曼滤波法适合于各种各样的弹道模型，对于不同控制系统的炮弹还能做出及时调整；本发明根据炮弹前半段的飞行数据辨识出气动参数，然后根据辨识结果对炮弹进行及时调整以提高射击精度，具有实时性；本发明适用于非线性系统，其算法简单、实用性强、精度高，在已经获得观测值的情况下，通过弹丸的运动方程和滤波算法就可以获得相应的气动参数。
[0042] 下面结合
说明书附图对本发明做进一步说明。
附图说明
[0044] 图2为
实施例1中辨识零升阻力系数的结果示意图。
具体实施方式
[0045] 下面结合附图对本发明作进一步详细描述。
[0046] 如图1所示，一种基于扩展卡尔曼滤波的辨识高速旋转弹气动参数滤波方法，包括以下步骤：
[0047] 步骤1、建立高速旋转弹四自由度动力学模型，四自由度动力学模型建立在地面
坐标系下。
[0048] 进一步的实施例中，所述高速旋转弹四自由度动力学模型具体为：
[0049]
[0050]
[0051]
[0052]
[0053]
[0054]
[0055]
[0056]
[0057] 弹丸相对于空气的速度为：
[0058]
[0059] 动力平衡角αe的直接计算公式为：
[0060]
[0061]
[0062]
[0063]
[0064] 直接计算公式中的参数为：
[0065]
[0066]
[0067]
[0068]
[0069] 式中，Vx、Vy、Vz分别为弹丸在x、y、z方向的速度；x、y、z为弹丸的位置坐标；γ为滚转角，为滚转角速度；m为弹丸质量；d为弹丸最大横截面直径；l为弹长；g为
重力加速度；ρ为空气密度；S为特征面积，通常取为弹丸的最大横截面积；Cx0和Cx2分别为零升阻力系数和诱导阻力系数；αe为动力平衡角；Vr为弹丸相对于空气的速度；C′y为升力系数导数；C″z为马氏力系数的联合偏导数；C为极转动惯量；A为赤道转动惯量；m′xz为极阻尼力矩系数导数；Wx、Wz分别为x和z方向的风速；m″y为马格努斯力矩系数的联合偏导数；m′z为静力矩系数导数。
[0070] 步骤2、根据高速旋转弹四自由度动力学模型，利用扩展卡尔曼滤波辨识气动参数，具体为：
[0071] 步骤2-1、将待辨识参数加入状态变量中形成增广状态向量，得到增广之后的状态方程和量测方程；
[0072] 根据步骤1中所列出的动力学模型，有8个状态变量，状态向量为：
[0073]
[0074] 本发明是根据系统测量量间接估计系统状态量，所以把系统建模过程中的未知参数(气动参数)当做状态变量处理，建立增广式状态方程。方程(1-8)可以写成的形式，则状态方程和量测方程为：
[0075]
[0076] 由于在实际工程中，不能得到所有的状态变量，只能得到部分状态变量或者和状态变量有关的变量作为观测值，所以量测方程是状态向量x(t)与观测值z(k)之间的关系方程。式中，x(t)为状态向量，x(t0)为初始状态向量，u(t)为外部输入，y(t)为连续时间下的观测向量，y(k)为离散时间下的观测向量，β为未知的参数向量，z(k)为量测向量，F为系统噪声转移矩阵，G为量测噪声转移矩阵，w(t)为系统噪声，v(k)为量测噪声。
[0077] 设待辨识参数为β，在小范围内可视为常量，因此有如下方程：
[0078]
[0079] 将其加入状态方程之中，将β加入状态变量之中，形成增广状态变量：
[0080]
[0081] 进一步的实施例中，增广之后的状态方程和量测方程为：
[0082]
[0083] y(t)＝ga[xa(t),u(t)]
[0084] z(k)＝y(k)+Gv(k)
[0085] 步骤2-2、根据实际工程的初始值以及增广之后的状态方程和量测方程，得到每一时刻的增广状态变量的估计值以及待辨识参数估计值，具体步骤为：
[0086] 步骤2-2-1、确定初始时刻的增广状态变量估计值及增广状态变量误差协方差矩阵估计值。
[0087] 根据实际工程中可以得到的初始数据，得到初始增广状态变量：
[0088]
[0089] 根据以往的统计数据得到初始状态误差协方差矩阵：
[0090]
[0091] 式中， 为初始增广状态变量估计值。
[0092] 设连续时间条件下t＝t0时刻，对应着离散时间k＝1时刻，t0时刻(初始时刻)的增广状态变量估计值及增广状态变量误差协方差矩阵估计值分别等于 和 即：
[0093]
[0094] 用 表示k-1时刻的估计的增广状态变量，则 表示0(即k＝1)时刻的估计的增广状态变量； 表示k＝1时刻的估计的增广状态变量误差协方差矩阵；
[0095] 步骤2-2-2、从步骤2-2-1确定的初始时刻的增广状态变量估计值及增广状态变量误差协方差矩阵估计值开始，根据k-1时刻的增广状态变量估计值及增广状态变量误差协方差矩阵估计值得到k时刻的增广状态变量预测值及增广状态变量误差协方差矩阵预测值。
[0096] 将增广后的状态方程线性化，用增广后的状态函数fa对k-1时刻的增广状态变量估计值 求雅克比矩阵，得到Aa(k)：
[0097]
[0098] 从而得到状态转移矩阵Φa(k)：
[0099]
[0100] Δt为k-1时刻到k时刻的时间间隔。
[0101] 由k-1时刻的估计值得到k时刻的预测值：
[0102]
[0103]
[0104] 表示k时刻的增广状态变量的预测值， 为外部输入， 表示连续时间t时刻下的估计的增广状态变量，t(k)表示对应着k时刻的时间， 为k时刻的增广状态变量误差协方差矩阵的预测值， 为k-1时刻的增广状态变量误差协方差矩阵的估计值；
[0105] 步骤2-2-3、根据步骤2-2-2求得的k时刻的增广状态变量误差协方差矩阵预测值以及k时刻的增广状态变量的预测值得到k时刻的增广状态变量误差协方差矩阵估计值以及k时刻的增广状态变量的估计值。
[0106] 将增广后的量测方程线性化，用增广后的量测函数ga对k时刻的增广状态变量预测值 求雅克比矩阵，得到Ca(k)：
[0107]
[0108] k时刻的滤波增益矩阵为：
[0109]
[0110] k时刻的观测向量为：
[0111]
[0112] k时刻的增广状态变量误差协方差矩阵的估计值为：
[0113]
[0114] k时刻的增广状态变量的估计值为：
[0115]
[0116] u(k)为外部输入，I为单位矩阵，z(k)为k时刻的量测向量。在实际工程项目中，z(k)为仪器设备能够测得的数据。从k-1时刻到k时刻，滤波得到的最终结果为 和为k时刻的增广状态变量误差协方差矩阵的估计值，用于下一步迭代运算；为求得的k时刻的增广状态变量的估计值， 中含有待辨识参数，k时刻的待辨识参数(气动参数)就可以得到。
[0117] 步骤2-2-4、返回步骤2-2-2，用k时刻的估计值得到k+1时刻的估计值，直到量测向量完结。这样就得到每一时刻的增广状态变量估计值，也就得到了每一时刻的待辨识参数估计值。
[0118] 下面结合实施例对本发明做进一步说明。
[0119] 实施例1
[0120] 下面结合实施例进行更详细的描述。
[0121] 某型号炮弹的参数设置如下表所示：
[0122] 表1高速旋转弹的参数值
[0123]
[0124] 在表2的初始条件下，对零升阻力系数Cx0进行参数辨识：
[0125] 表2高速旋转弹的初始条件
[0126]
[0127] 待辨识参数β为Cx0，则增广状态向量为：
[0128]
[0129] 增广后的状态方程和量测方程为：
[0130]
[0131] y(t)＝ga[xa(t),u(t)]
[0132] z(k)＝y(k)+Gv(k)
[0133] 忽略外部输入u(t)，将F和G视为单位矩阵。
[0134] 实际情况下，只有三个速度和位置坐标可以测得，所以量测值为：
[0135] z＝[Vx Vy Vz x y z]T
[0136] 本例中的 和 分别为：
[0137]
[0138]
[0139] 零升阻力系数Cx0随马赫数变化，马赫数在高速旋转弹运动过程中不断变化。最终辨识Cx0的结果如图2所示，其中EKF为扩展卡尔曼滤波的英文缩写。可以看出，该方法将高速旋转弹运动过程的零升阻力系数辨识了出来。
