CN114459302A - 一种适用于高旋弹丸的滚转角速率测量方法 - Google Patents

Abstract

本发明提出了一种适用于高旋弹丸的滚转角速率测量方法，通过Oy_b与Ox_b和Oz_b构成右手坐标系，通过右手坐标系建立微机械MEMS陀螺仪随机误差模型，进一步通过给定的陀螺仪的零偏以及输出量建立状态方程以及量测方程，采用卡尔曼滤波器进行数据融合，实时解算出角速率输出；再通过定义选择向量，即可得到融合之后的角速率值，最后将角速率值带入发明中的公式即可得到弹丸的角速率，实现高旋弹丸滚转角速率的高精度测量。

Description

一种适用于高旋弹丸的滚转角速率测量方法

技术领域

本发明涉及军工制作技术领域，具体为一种适用于高旋弹丸的滚转角速率测量方法。

背景技术

精确制导弹药能够显著提高火力打击的精度和效率，成为了目前战场上重要的打击手段之一，也成为了各国军方兵器装备发展的方向。制导弹药实现精确打击的关键在于如何精确地获取弹体的实时运动信息。弹丸极高的自旋率使精确制导弹药的制导问题复杂化。因此，需要对弹丸的滚转角速率进行高精度测量。

目前，用于姿态测量的传感器主要有：太阳方位传感器、连续波雷达、红外传感器、全球导航卫星系统(GNSS)、磁强计、加速度计、陀螺仪等。但是由于弹载工作环境的恶劣性，适合的传感器只有陀螺仪。同时，陀螺仪由于具有低成本，小体积，高输出率等优点被广泛应用于各种场合。然而，当弹丸高速飞行时，具有大量程陀螺仪的误差将随着时间的推移迅速累积，并导致较大的姿态误差。因此，对传统结构的MEMS陀螺仪来说，精确测量高自旋弹丸的滚转角速率是一项亟需解决的技术问题。

发明内容

为解决上述技术问题，本发明提出了一种适用于高旋弹丸的滚转角速率测量方法，能够有效实现高旋弹丸滚转角速率的高精度测量。

一种适用于高旋弹丸的滚转角速率测量方法，其特征在于，包括以下步骤：

步骤1：选取地理坐标系作为导航坐标系，导航坐标系为n系，惯性测量单元IMU固定于弹丸中心，弹体的质心为载体坐标系的原点O，Ox_b轴指向炮弹的前进方向并与炮弹的纵轴保持一致，Oz_b轴正交于Ox_b并指向上且位于炮弹纵向对称面内，Oy_b与Ox_b和Oz_b构成右手坐标系；

步骤2：建立微机械MEMS陀螺仪随机误差模型：

微机械MEMS陀螺仪随机误差模型表示为

其中，b为陀螺仪的零偏；n_a为角度随机游走且服从高斯分布，表示为

E[n_a(t)]＝0,E[n_a(t)n^T _a(t+τ)]＝R_kδ(τ) (2)

δ(τ)是一个delta相关的噪声过程，表示为

E[n_r(t)]＝0,E[n_r(t)n^T _r(t+τ)]＝Q_kδ(τ) (4)

两个白噪声过程不相关，表示为

E[n_a(t)n^T _r(t+τ)]＝0 (5)

步骤3：通过给定的陀螺仪的零偏、输出量建立状态方程以及量测方程，采用卡尔曼滤波器进行数据融合，实时解算出角速率输出；

步骤4：通过定义选择向量，得到融合之后的角速率值；

步骤5：将步骤4得到的角速率值带入下列方程，得到弹丸的角速率，

其中，ω_ie是地球的自转角速率，ω_en是导航坐标系因载体运动而在地球上旋转产生的角速率，ω_nb为姿态速率。

作为本发明的进一步改进：在步骤1中，所述Ox_b轴上安装有三个独立的MEMS陀螺仪。

作为本发明的进一步改进：在步骤2中，所述微机械MEMS陀螺仪随机误差模型使用ARIMA模型，表示为

其中，

Φ(B)为自回归系数，Θ(B)为移动平滑系数，B表示后移运算符,p自回归模型的阶数，q称为移动平均模型的阶数，d称为差分阶数。作为本发明的进一步改进：所述ARIMA模型的阶数越高，即p和q的值越大，陀螺仪系统输出数据的特性就越好，随着阶数的增加，滤波计算的维数随之增加，因此，选择ARIMA(2,1,1)模型对时间序列进行建模，ARIMA(2,1,1)中的参数如下所示：

因此，基于ARIMA模型的MEMS陀螺仪随机误差模型可以表示为：

其中，Y(t)＝[y₁,y₂,y₃]^T，B(t)＝[b₁,b₂,b₃]^T，N_a(t)＝[n_a1,n_a2,n_a3]^T，

N_r(t)＝[n_r1,n_r2,n_r3]^T；Y(t)是一个三维矢量，由三个陀螺仪的测量值组成，B(t)是放置在X_b轴上的陀螺仪的漂移，N_a(t)和N_r(t)分别为ARW和RRW的噪声矢量。作为本发明的进一步改进：在步骤3中，所述卡尔曼滤波器为自适应滤波，在弹丸处在低动态状态下，MEMS陀螺仪输出数据的融合采用标准卡尔曼滤波器，当弹丸处在高动态状态下，基于ARIMA模型的MEMS陀螺仪随机误差模型，MEMS陀螺仪输出数据的融合采用Sage-Husa自适应卡尔曼滤波器，状态向量表示为：

X＝[b₁,b₂,b₃,x_k,x_k-1,x_k-2,ω]^T, (12)

ω是真实角速率，由n_ω驱动，

E[n_ω(t)n^T _ω(t+τ)]＝q_ωδ(τ), (14)

其中，n_ω是delta相关噪声过程，q_ω是n_ω的方差，根据方程(11)，状态方程可以表示为：

其中，

W(t)＝[n_r1,n_r2,n_r3,a_k,a_k-1,n_ω]^T (16)

将状态方程进行离散化处理，可以表示为

X_k＝Φ_k/k-1X_k-1+Γ_k/k-1W_k-1 (17)

其中，Φ_k/k-1是系统转移矩阵，Γ_k/k-1是系统噪声分布矩阵，W_k-1是系统噪声向量，在方程(17),Φ_k/k-1为7×7的矩阵，Γ_k/k-1为7×6的矩阵，分别表示为

其中，T_s代表采样时间间隔，

将MEMS陀螺仪的输出作为量测量，表示为

Z_k＝[y₁,y₂,y₃]^T (20)

量测方程可以表示为

Z_k＝H_kX_k+V_k (21)

其中，H_k是连接状态向量和量测向量之间的转移矩阵，V_k表示量测噪声矢量，分别表示为

V_k＝[n_a1,n_a2,n_a3]^T (23)

在方程(17)和(21),W_k-1和V_k分别为高斯白噪声向量序列，即：

W_k-1～WN(q_k,Q_k),V_k-1～WN(r_k,R_k)，当两噪声之间互不相关，噪声均值为0时或噪声值不为0时，即有：

作为本发明的进一步改进：所述Sage-Husa自适应卡尔曼滤波器的自适应滤波融合算法中的滤波增益由下式给出：

其中，α_k∈(0,1]是自适应因子，K_k是SHAKF算法的增益。

作为本发明的进一步改进：所述Sage-Husa自适应卡尔曼滤波器的自适应滤波融合算法中的自适应因子采用三段差分函数构造，表示为

其中，1≤c₀≤1.5，3≤c₁≤4.5，是两个常数，预测状态误差模型

的学习统计量表示为

其中，tr(·)代表矩阵的秩

与现有技术相比，本发明的有益效果为：

本发明通过Oy_b与Ox_b和Oz_b构成右手坐标系，通过右手坐标系建立微机械MEMS陀螺仪随机误差模型，进一步通过给定的陀螺仪的零偏以及输出量建立状态方程以及量测方程，采用卡尔曼滤波器进行数据融合，实时解算出角速率输出；再通过定义选择向量，即可得到融合之后的角速率值，最后将角速率值带入

即可得到弹丸的角速率，实现高旋弹丸滚转角速率的高精度测量。

附图说明

图1是弹丸的几何结构图；

图2是陀螺仪的分布图；

图3是Sage–Husa自适应卡尔曼滤波(SHAKF)算法流程图。

具体实施方式

下面结合附图与具体实施方式对本发明作进一步详细描述：

一种适用于高旋弹丸的滚转角速率测量方法，包括以下步骤：

步骤1：选取地理坐标系作为导航坐标系，n系，惯性测量单元IMU固定于弹丸中心，弹体的质心为载体坐标系的原点O，Ox_b轴指向炮弹的前进方向并与炮弹的纵轴保持一致，Oz_b轴正交于Ox_b并指向上且位于炮弹纵向对称面内，Oy_b与Ox_b和Oz_b构成右手坐标系，如图1所示，为弹丸的几何结构图；

步骤2：建立微机械(MEMS)陀螺仪随机误差模型，表示为

其中，b为陀螺仪的零偏；n_a为角度随机游走且服从高斯分布，表示为

E[n_a(t)]＝0,E[n_a(t)n^T _a(t+τ)]＝R_kδ(τ) (2)

δ(τ)是一个delta相关的噪声过程，表示为

E[n_r(t)]＝0,E[n_r(t)n^T _r(t+τ)]＝Q_kδ(τ) (4)

两个白噪声过程不相关，表示为

E[n_a(t)n^T _r(t+τ)]＝0 (5)

步骤3：通过给定的陀螺仪的零偏以及输出量建立状态方程以及量测方程，采用卡尔曼滤波器进行数据融合，实时解算出角速率输出；

步骤4：通过定义选择向量，即可得到融合之后的角速率值；

步骤5：将步骤4得到的角速率值带入下列方程，即可得到弹丸的角速率；

其中，ω_ie是地球的自转角速率，ω_en是导航坐标系因载体运动而在地球上旋转产生的角速率，ω_nb为姿态速率。

在Ox_b轴上安装三个独立的MEMS陀螺仪，如图2所示，为陀螺仪的分布图。

MEMS陀螺仪的随机误差模型使用ARIMA模型，可以表示为

其中，

Φ(B)为自回归系数，Θ(B)为移动平滑系数。B表示后移运算符,p自回归模型的阶数，q称为移动平均模型的阶数，d称为差分阶数。假设白噪声a_k独立且同分布，且均值为0，方差为

在弹丸的实际飞行过程中，陀螺仪的输出数据不满足平稳序列。而平稳数据可以通过一阶差分得到的。ARIMA模型的阶数越高，即p和q的值越大，陀螺仪系统输出数据的特性就越好。但是，随着阶数的增加，滤波计算的维数将大大增加。通过对五组实验数据的综合分析，选择ARIMA(2,1,1)模型对时间序列进行建模。ARIMA(2,1,1)中的参数如下所示：

θ₁＝-0.6754。

因此，基于ARIMA模型的MEMS陀螺仪随机误差模型可以表示为：

其中，Y(t)＝[y₁,y₂,y₃]^T，B(t)＝[b₁,b₂,b₃]^T，N_a(t)＝[n_a1,n_a2,n_a3]^T，N_r(t)＝[n_r1,n_r2,n_r3]^T。Y(t)是一个三维矢量，由三个陀螺仪的测量值组成。B(t)是放置在X_b轴上的陀螺仪的漂移。N_a(t)和N_r(t)分别为ARW和RRW的噪声矢量。

步骤3中的卡尔曼滤波器采用自适应滤波方式。在低动态状态下，MEMS陀螺仪输出数据的融合采用标准卡尔曼滤波器。标准卡尔曼滤波器获得最优估计的前提是需要准确地获知随机动力系统的结构参数和统计噪声参数。然而，当弹丸处在高动态的飞行环境时，MEMS陀螺仪的过程噪声往往是时变的或未知的。根据修正后的随机误差模型，采用Sage-Husa自适应卡尔曼滤波器，用于MEMS陀螺仪输出数据的融合，该滤波器能够在弹丸高动态飞行条件下抑制时变过程和测量噪声。状态向量表示为：

X＝[b₁,b₂,b₃,x_k,x_k-1,x_k-2,ω]^T, (12)

ω是真实角速率，由n_ω驱动。

E[n_ω(t)n^T _ω(t+τ)]＝q_ωδ(τ), (14)

其中，n_ω是delta相关噪声过程，q_ω是n_ω的方差。根据方程(11)，状态方程可以表示为：

其中，

W(t)＝[n_r1,n_r2,n_r3,a_k,a_k-1,n_ω]^T (16)

将状态方程进行离散化处理，可以表示为

X_k＝Φ_k/k-1X_k-1+Γ_k/k-1W_k-1 (17)

其中，Φ_k/k-1是系统转移矩阵，Γ_k/k-1是系统噪声分布矩阵，W_k-1是系统噪声向量。在方程(17),Φ_k/k-1为7×7的矩阵，Γ_k/k-1为7×6的矩阵，分别表示为

其中，T_s代表采样时间间隔。

将MEMS陀螺仪的输出作为量测量，表示为

Z_k＝[y₁,y₂,y₃]^T (20)

量测方程可以表示为

Z_k＝H_kX_k+V_k (21)

其中，H_k是连接状态向量和量测向量之间的转移矩阵，V_k表示量测噪声矢量，分别表示为

V_k＝[n_a1,n_a2,n_a3]^T (23)

在方程(17)和(21),W_k-1和V_k分别为高斯白噪声向量序列，即：

W_k-1～WN(q_k,Q_k),V_k-1～WN(r_k,R_k)。假设两噪声之间互不相关，噪声均值不一定为0，即有：

自适应滤波融合算法带有过程噪声时变估计器，时变噪声估计器由下式给出：

Sage-Husa自适应卡尔曼滤波器的自适应滤波融合算法中的滤波增益由下式给出：

该算法的流程图为图3。

其中，α_k∈(0,1]是自适应因子，K_k是SHAKF算法的增益。自适应因子有助于平衡观测状态和预测状态。

Sage-Husa自适应卡尔曼滤波器的自适应滤波融合算法中的自适应因子采用三段差分函数构造，表示为

其中，1≤c₀≤1.5，3≤c₁≤4.5，是两个常数。预测状态误差模型

的学习统计量表示为

其中，t_r(·)代表矩阵的秩。

本发明通过Oy_b与Ox_b和Oz_b构成右手坐标系，通过右手坐标系建立微机械MEMS陀螺仪随机误差模型，进一步通过给定的陀螺仪的零偏以及输出量建立状态方程以及量测方程，采用卡尔曼滤波器进行数据融合，实时解算出角速率输出；再通过定义选择向量，即可得到融合之后的角速率值，最后将角速率值带入

即可得到弹丸的角速率，从而实现高旋弹丸滚转角速率的高精度测量。

以上所述，仅是本发明的较佳实施例而已，并非是对本发明作任何其他形式的限制，而依据本发明的技术实质所作的任何修改或等同变化，仍属于本发明所要求保护的范围。

Claims (7)

1.一种适用于高旋弹丸的滚转角速率测量方法，其特征在于，包括以下步骤：

步骤1：选取地理坐标系作为导航坐标系，导航坐标系为n系，惯性测量单元IMU固定于弹丸中心，弹体的质心为载体坐标系的原点O，Ox_b轴指向炮弹的前进方向并与炮弹的纵轴保持一致，Oz_b轴正交于Ox_b并指向上且位于炮弹纵向对称面内，Oy_b与Ox_b和Oz_b构成右手坐标系；

步骤2：建立微机械MEMS陀螺仪随机误差模型：

微机械MEMS陀螺仪随机误差模型表示为

其中，b为陀螺仪的零偏；n_a为角度随机游走且服从高斯分布，表示为

E[n_a(t)]＝0,E[n_a(t)n^T _a(t+τ)]＝R_kδ(τ) (2)

δ(τ)是一个delta相关的噪声过程，表示为

E[n_r(t)]＝0,E[n_r(t)n^T _r(t+τ)]＝Q_kδ(τ) (4)

两个白噪声过程不相关，表示为

E[n_a(t)n^T _r(t+τ)]＝0 (5)

步骤3：通过给定的陀螺仪的零偏、输出量建立状态方程以及量测方程，采用卡尔曼滤波器进行数据融合，实时解算出角速率输出；

步骤4：通过定义选择向量，得到融合之后的角速率值；

步骤5：将步骤4得到的角速率值带入下列方程，得到弹丸的角速率，

其中，ω_ie是地球的自转角速率，ω_en是导航坐标系因载体运动而在地球上旋转产生的角速率，ω_nb为姿态速率。

2.根据权利要求1所述的一种适用于高旋弹丸的滚转角速率测量方法，其特征在于：在步骤1中，所述Ox_b轴上安装有三个独立的MEMS陀螺仪。

3.根据权利要求1所述的一种适用于高旋弹丸的滚转角速率测量方法，其特征在于：在步骤2中，所述微机械MEMS陀螺仪随机误差模型使用ARIMA模型，表示为

Φ(B)▽^dx_k＝Θ(B)a_k, (8)

其中，▽^d＝(1-B)^d,Φ(B)为自回归系数，Θ(B)为移动平滑系数，B表示后移运算符,p自回归模型的阶数，q称为移动平均模型的阶数，d称为差分阶数。

4.根据权利要求3所述的一种适用于高旋弹丸的滚转角速率测量方法，其特征在于：所述ARIMA模型的阶数越高，即p和q的值越大，陀螺仪系统输出数据的特性就越好，随着阶数的增加，滤波计算的维数随之增加，因此，选择ARIMA(2,1,1)模型对时间序列进行建模，ARIMA(2,1,1)中的参数如下所示：

θ₁＝-0.6754

因此，基于ARIMA模型的MEMS陀螺仪随机误差模型可以表示为：

其中，Y(t)＝[y₁,y₂,y₃]^T，B(t)＝[b₁,b₂,b₃]^T，N_a(t)＝[n_a1,n_a2,n_a3]^T，N_r(t)＝[n_r1,n_r2,n_r3]^T；Y(t)是一个三维矢量，由三个陀螺仪的测量值组成，B(t)是放置在X_b轴上的陀螺仪的漂移，N_a(t)和N_r(t)分别为ARW和RRW的噪声矢量。

5.根据权利要求1所述的一种适用于高旋弹丸的滚转角速率测量方法，其特征在于：在步骤3中，所述卡尔曼滤波器为自适应滤波，在弹丸处在低动态状态下，MEMS陀螺仪输出数据的融合采用标准卡尔曼滤波器，当弹丸处在高动态状态下，基于ARIMA模型的MEMS陀螺仪随机误差模型，MEMS陀螺仪输出数据的融合采用Sage-Husa自适应卡尔曼滤波器，状态向量表示为：

X＝[b₁,b₂,b₃,x_k,x_k-1,x_k-2,ω]^T, (12)

ω是真实角速率，由n_ω驱动，

E[n_ω(t)n^T _ω(t+τ)]＝q_ωδ(τ), (14)

其中，n_ω是delta相关噪声过程，q_ω是n_ω的方差，根据方程(11)，状态方程可以表示为：

其中，

W(t)＝[n_r1,n_r2,n_r3,a_k,a_k-1,n_ω]^T (16)

将状态方程进行离散化处理，可以表示为

X_k＝Φ_k/k-1X_k-1+Γ_k/k-1W_k-1 (17)

其中，Φ_k/k-1是系统转移矩阵，Γ_k/k-1是系统噪声分布矩阵，W_k-1是系统噪声向量，在方程(17),Φ_k/k-1为7×7的矩阵，Γ_k/k-1为7×6的矩阵，分别表示为

其中，T_s代表采样时间间隔，

将MEMS陀螺仪的输出作为量测量，表示为

Z_k＝[y₁,y₂,y₃]^T (20)

量测方程可以表示为

Z_k＝H_kX_k+V_k (21)

其中，H_k是连接状态向量和量测向量之间的转移矩阵，V_k表示量测噪声矢量，分别表示为

V_k＝[n_a1,n_a2,n_a3]^T (23)

在方程(17)和(21),W_k-1和V_k分别为高斯白噪声向量序列，即：

W_k-1～WN(q_k，Q_k),W_k-1～WN(r_k，R_k)，当两噪声之间互不相关，噪声均值为0时或噪声值不为0时，即有：

6.根据权利要求5所述的一种适用于高旋弹丸的滚转角速率测量方法，其特征在于：所述Sage-Husa自适应卡尔曼滤波器的自适应滤波融合算法中的滤波增益由下式给出：

其中，α_k∈(0,1]是自适应因子，K_k是SHAKF算法的增益。

7.根据权利要求6所述的一种适用于高旋弹丸的滚转角速率测量方法，其特征在于：所述Sage-Husa自适应卡尔曼滤波器的自适应滤波融合算法中的自适应因子采用三段差分函数构造，表示为

其中，1≤c₀≤1.5，3≤c₁≤4.5，是两个常数，预测状态误差模型

的学习统计量表示为

其中，tr(·)代表矩阵的秩。

Images