CN103162678B - 批量mems陀螺信息融合方法 - Google Patents

Abstract

批量MEMS陀螺信息融合方法，属于MEMS陀螺信息整合技术领域。本发明是为了解决采用现有针对单一的陀螺信息进行处理的方法对多传感器的信息进行融合，运算量大及输出结果精度低的问题。它将批量MEMS陀螺均分成N组陀螺阵列；采用初级融合模块对每组陀螺阵列中所有MEMS陀螺采集获取的数据进行基于支持度的初级融合，获得每组陀螺阵列的融合数据；再通过时间序列模型获得对应的状态空间方程；再进行多级序贯式滤波，获得第N级融合状态估计值，经提取获得批量MEMS陀螺k时刻融合后的角速度。本发明用于批量MEMS陀螺信息的融合。

Description

批量MEMS陀螺信息融合方法

技术领域

本发明涉及批量MEMS陀螺信息融合方法，属于MEMS陀螺信息整合技术领域。

背景技术

目前，随着微机电系统(Micro Electronic Mechanical System，MEMS)的发展，MEMS陀螺技术开始不断成熟起来，MEMS陀螺可以看作一种微型角速度传感器。现有MEMS陀螺的种类很丰富，均具有体积小、重量轻、功耗低、成本低等特点。目前，为提高MEMS陀螺精度通过以下两种方式：

一方面是从硬件上，对陀螺的设计及加工工艺加以改进；二是从软件上，对陀螺进行滤波处理，滤除多余的噪声，它采用时间序列与Kalman滤波相结合的方法，对陀螺的数据进行滤波处理。该方法主要针对单一的陀螺信息进行处理，面对多传感器的信息融合，仍旧采用这种普通Kalman滤波处理信息时，存在运算量过大的问题，并且，获得的MEMS陀螺的输出结果精度不理想。

发明内容

本发明目的是为了解决采用现有针对单一的陀螺信息进行处理的方法对多传感器的信息进行融合，运算量大及输出结果精度低的问题，提供了一种批量MEMS陀螺信息融合方法。

本发明所述批量MEMS陀螺信息融合方法，它包括以下步骤：

步骤一：将批量MEMS陀螺均分成N组陀螺阵列，每组陀螺阵列包括n个MEMS陀螺，n取值为4、5、6或7，N为大于2的正整数；

步骤二：采用初级融合模块对每组陀螺阵列中所有MEMS陀螺采集获取的数据进行基于支持度的初级融合，获得每组陀螺阵列的融合数据Y(k)，其中k表示时刻值；

步骤三：对每组陀螺阵列输出的融合数据Y(k)进行时间序列建模获得时间序列模型，并通过时间序列模型获得对应的状态空间方程；

步骤四：采用多级序贯式滤波模块对N组陀螺阵列的融合数据Y(k)和N个状态空间方程顺次进行多级序贯式滤波，其中每个陀螺阵列的融合数据Y(k)和状态空间方程作为多级序贯式滤波模块中对应的子滤波器的输入，获得第N级融合状态估计值

步骤五：将第N级融合状态估计值输入角速度提取模块，经提取获得批量MEMS陀螺k时刻融合后的角速度ω(k)。

所述步骤二中获得每组陀螺阵列的融合数据Y(k)的具体方法为：

根据鲁棒理论中的尺度估计，定义相互支持度函数s_ij(k)和每组陀螺阵列中第i个MEMS陀螺的自身支持度函数s_i(k)为：

$s_{\mathrm{ij}}\left(k\right)=\left\{\begin{array}{cc}1-\sin [\frac{|Z_i\left(k\right)-Z_j\left(k\right)|}{\mathrm{\βr}\left(k\right)}\·\frac{\mathrm{\π}}{2}],& |Z_i\left(k\right)-Z_j\left(k\right)|\≤\mathrm{\βr}\left(k\right),\\ 0,& \mathrm{other}\end{array}\right.$

相互支持度函数s_ij(k)表示k时刻，每组陀螺阵列中MEMS陀螺i和MEMS陀螺j的相互支持度，其中i，j＝1，2，…，n，并且i≠j，

Z_i(k)为每组陀螺阵列中第i个MEMS陀螺在k时刻的输出值，Z_j(k)为每组陀螺阵列中第j个MEMS陀螺在k时刻的输出值，β为阈值尺度参数，r(k)为相互支持度阈值因子；

$s_i\left(k\right)=\left\{\begin{array}{cc}1-\sin [\frac{|Z_i\left(k\right)-\mathrm{median}\left(Z_i\right)|}{\mathrm{\&beta;}{r}_i}\&CenterDot;\frac{\mathrm{\&pi;}}{2}],& |Z_i\left(k\right)-\mathrm{median}\left(Z_i\right)|<{\mathrm{\&beta;r}}_i,\\ 0,& \mathrm{other}\end{array}\right.$

自身支持度函数s_i(k)表示在k时刻每组陀螺阵列中第i个MEMS陀螺采集获取的数据Z_i对其输出值Z_i(k)的支持程度，r_i为每组陀螺阵列中第i个MEMS陀螺的自身支持度阈值因子；

建立每组陀螺阵列的支持度矩阵D(k)：

$D\left(k\right)=\left[\begin{array}{ccccc}{d}_{11}\left(k\right)& d_{12}\left(k\right)& ...& d_{1n}\left(k\right)& d_{1(n+1)}\left(k\right)\\ d_{21}\left(k\right)& d_{22}\left(k\right)& ...& d_{2n}\left(k\right)& d_{2(n+1)}\left(k\right)\\ ...& ...& ...& ...& ...\\ d_{n1}\left(k\right)& d_{n2}\left(k\right)& ...& d_{\mathrm{nn}}\left(k\right)& d_{n(n+1)}\left(k\right)\\ d_{(n+1)1}\left(k\right)& d_{\left(n+1\right)2}\left(k\right)& ...& d_{(n+1)n}\left(k\right)& 1\end{array}\right],$

其中，

d_ii(k)＝s_i ²(k)，

d_ij(k)＝s_i(k)s_ij(k)s_j(k)，

$d_{i(n+1)}\left(k\right)=\left\{\begin{array}{cc}(1-\sin [\frac{|Z_i\left(k\right)-Y(k-1)|}{{\mathrm{\&beta;r}}_i}\&CenterDot;\frac{\mathrm{\&pi;}}{2}\left]\right)\&times;d_{\mathrm{ii}}\left(k\right),& |Z_i\left(k\right)-Y(k-1)|<{\mathrm{\&beta;r}}_i,\\ 0,& \mathrm{other}\end{array}\right.$

$d_{(n+1)j}\left(k\right)=\left\{\begin{array}{cc}(1-\sin [\frac{|Z_j\left(k\right)-Y(k-1)|}{{\mathrm{\&beta;r}}_j}\&CenterDot;\frac{\mathrm{\&pi;}}{2}\left]\right)\&times;d_{\mathrm{jj}}\left(k\right),& |Z_j\left(k\right)-Y(k-1)|<{\mathrm{\&beta;r}}_j,\\ 0,& \mathrm{other}\end{array}\right.$

其中，当k＝1时，Y(1)＝(Z₁(1)+Z₂(1)+…+Z_n(1))/n；

通过支持度矩阵D(k)的第i行判断获得每组陀螺阵列中第i个MEMS陀螺的权值系数w_i(k)，则权值向量Δ(k)＝[w₁(k)，w₂(k)，…，w_n+1(k)]^T；

设置一组非负数系数a₁(k)，a₂(k)，…，a_n+1(k)，使得：

w_i(k)＝a₁(k)d_i1(k)+d₂(k)d_i2(k)+…+a_n+1(k)d_i(n+1)(k)，

获得系数矩阵A(k)＝[a₁(k)，a₂(k)，…，a_n+1(k)]^T，

即Δ(k)＝D(k)A(k)；

支持度矩阵D(k)为非负对称矩阵，则D(k)的最大模特征值λ＞0，并有：

Δ(k)＝λA(k)；

每组陀螺阵列中第v个MEMS陀螺的权值系数w_v(k)与第l个MEMS陀螺的权值系数w_l(k)的关系为：

$\frac{w_v\left(k\right)}{w_l\left(k\right)}=\frac{a_v\left(k\right)}{a_l\left(k\right)}$ v，l＝1，2，…，n+1，

对每组陀螺阵列中第i个MEMS陀螺的权值系数w_i(k)归一化处理，得到：

$w_i\left(k\right)=\frac{a_i\left(k\right)}{a_1\left(k\right)+a_2\left(k\right)+\&CenterDot;\&CenterDot;\&CenterDot;+a_{n+1}\left(k\right)},$

则获得每组陀螺阵列的融合数据Y(k)：

$Y\left(k\right)=\frac{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{w}_i\left(k\right)Z_i\left(k\right)+w_{n+1}\left(k\right)Y(k-1)}{\underset{i=1}{\overset{n+1}{\mathrm{\&Sigma;}}}{w}_i\left(k\right)},$

其中，当k＝1时，Y(1)＝(Z₁(1)+Z₂(1)+…+Z_n(1))/n。

步骤三中获得的时间序列模型为自回归滑动平均模型ARMA(p，q)，其中p为自回归模型的阶数，q为滑动平均模型的阶数；

自回归滑动平均模型ARMA(p，q)的表达式为：

y(k)＝b₁y(k-1)+…+b_py(k-p)+ε(k)-a₁ε(k-1)-…-a_qε(k-q)，

其中y(k)为时间序列模型在k时刻的响应，b_u为自回归系数，υ＝1，2，…，p，ε(k)为随机扰动项白噪声，a_κ为滑动平均系数，κ＝1，2，…，q；

根据自回归滑动平均模型ARMA(p，q)建立状态空间方程，并且扩维引入真实角速度ω，扩维后的状态方程为：

X(k)＝Φ(k-1)X(k-1)+Γ(k-1)W(k)，

Y(k)＝H(k)X(k)+V(k)，

式中，X(k)为状态变量，Φ(k-1)为一步转移阵，Γ(k-1)为系统噪声驱动阵，W(k)为系统激励噪声序列，Y(k)由每组陀螺阵列的真实角速度ω和陀螺阵列的随机漂移η组成，H(k)为量测阵，V(k)为量测噪声序列。

步骤四中获得第N级融合状态估计值的具体方法为：

多级序贯式滤波模块第h级融合在k时刻获得的状态估计值为：

${\hat{X}}_{f\left(h\right)}\left(k\right)=(1-{\mathrm{\&lambda;}}_{h-1}\left(k\right)){\hat{X}}_{f(h-1)}\left(k\right)+{\mathrm{\&lambda;}}_{h-1}\left(k\right){\hat{X}}_h\left(k\right),$

式中为子滤波器h的状态估计，h＝2，3，…，N；

多级序贯式滤波模块第h级融合在k时刻获得的估计误差方差P_f(h)为：

P_f(h)(k|k)＝[1-λ_h(k)]²P_f(h-1)(k|k)+[λ_h(k)]²P_h(k|k)+2[1-λ_h(k)]λ_h(k)P_f(h-1)，h(k|k)，

其中，P_h为子滤波器h的估计误差方差；λ_h(k)为最优估计系数；

${\mathrm{\&lambda;}}_h\left(k\right)=\frac{{\mathrm{trP}}_{f(h-1)}\left(k\right|k)-{\mathrm{trP}}_{f(h-1),h}\left(k\right|k)}{{\mathrm{trP}}_{f(h-1)}\left(k\right|k)+{\mathrm{trP}}_h\left(k\right|k)-2\mathrm{tr}{P}_{f(h-1),h}\left(k\right|k)},$

$P_{f(h-1),h}\left(k\right|k)=\underset{\mathrm{\&alpha;}=2}{\overset{h-1}{\mathrm{\&Sigma;}}}[1-{\mathrm{\&lambda;}}_{h-2}\left(k\right)]P_{f\left(1\right),h}\left(k\right|k)+\underset{\mathrm{\&alpha;}=2}{\overset{h-1}{\mathrm{\&Sigma;}}}\left[\underset{\mathrm{\&beta;}=\mathrm{\&alpha;}+1}{\overset{h-1}{\mathrm{\&Sigma;}}}\right[1-{\mathrm{\&lambda;}}_{\mathrm{\&beta;}}\left(k\right)\left]\right]{\mathrm{\&lambda;}}_n\left(k\right)P_{\mathrm{\&alpha;},h}\left(k\right|k),$

式中，tr表示矩阵的迹，P_f(i-1)，i为和的互协方差阵，将第一级子滤波器滤波后获得的状态估计值作为第一级融合数据则P_f(1)，h(k|k)＝P_1，h(k|k)，并且P_1，h(k|k)为：

P_1，h(k|k)＝[I-K₁(k)H₁(k)]×[Φ(k-1)P_1，h(k-1|k-1)Φ^T(k-1)+

Γ(k-1)Q(k-1)Γ^T(k-1)]×[I-K_h(k)H_h(k)]^T

式中，K_h(k)为k时刻多级序贯式滤波模块中第h个子滤波器的滤波增益，H_h(k)为k时刻多级序贯式滤波模块中第h个子滤波器的量测阵，Q(k-1)为k-1时刻多级序贯式滤波模块中子滤波器系统噪声序列的方差阵；设置P_1，h(0|0)的初值以获得P_1，h(k|k)，最后获得第N级融合状态估计值

步骤五中获得批量MEMS陀螺k时刻融合后的角速度ω(k)的具体方法为：

$\mathrm{\&omega;}\left(k\right)=\left[\begin{array}{cccc}0& 0& ...& 1\end{array}\right]{\hat{X}}_g=\left[\begin{array}{cccc}0& 0& ...& 1\end{array}\right]\left[\begin{array}{c}\mathrm{\&eta;}\left(k\right)\\ \mathrm{\&eta;}(k-1)\\ .\\ .\\ .\\ \mathrm{\&eta;}(k-p+1)\\ \mathrm{\&omega;}\left(k\right)\end{array}\right].$

本发明的优点：本发明方法能够有效进行大批量MEMS陀螺的信息融合，无需知道陀螺数据的先验知识，能够极大地提高低成本、低精度陀螺的精度，减小其零偏不稳定性。其中陀螺阵列的初级融合是基于支持度的信息融合方法，经过改进后，不但能够提高陀螺的精度，降低陀螺的零偏不稳定性，还对陀螺的软、硬故障具有极强的容错能力。通过故障检测与隔离模块，若多级序贯式滤波的某一级发生故障并被检测到，则可以跳过这一级，顺序执行下一级融合；当故障被修复时，还可以恢复对该级的融合。

本发明方法采用多级分布式信息融合方法，避免了普通Kalman滤波处理多传感器信息时的运算量过大的问题，具有极强的容错能力，避免了MEMS陀螺故障对整个信息融合的污染，具有实际的应用意义。

本发明方法能够对多个相同或近似精度的MEMS陀螺的输出数据进行信息融合，通过对多个低精度的MEMS陀螺的输出结果进行多级融合，得到近似于高精度陀螺的输出结果，且具有极强的容错能力。

附图说明

图1是本发明所述批量MEMS陀螺信息融合方法的流程示意图；

图2是获得每组陀螺阵列的融合数据Y(k)的流程框图；

图3是批量MEMS陀螺信息融合方法的详细流程框图。

具体实施方式

具体实施方式一：下面结合图1至图3说明本实施方式，本实施方式所述批量MEMS陀螺信息融合方法，它包括以下步骤：

步骤一：将批量MEMS陀螺均分成N组陀螺阵列，每组陀螺阵列包括n个MEMS陀螺，n取值为4、5、6或7，N为大于2的正整数；

步骤二：采用初级融合模块对每组陀螺阵列中所有MEMS陀螺采集获取的数据进行基于支持度的初级融合，获得每组陀螺阵列的融合数据Y(k)，其中k表示时刻值；

步骤三：对每组陀螺阵列输出的融合数据Y(k)进行时间序列建模获得时间序列模型，并通过时间序列模型获得对应的状态空间方程；

步骤四：采用多级序贯式滤波模块对N组陀螺阵列的融合数据Y(k)和N个状态空间方程顺次进行多级序贯式滤波，其中每个陀螺阵列的融合数据Y(k)和状态空间方程作为多级序贯式滤波模块中对应的子滤波器的输入，获得第N级融合状态估计值

步骤五：将第N级融合状态估计值输入角速度提取模块，经提取获得批量MEMS陀螺k时刻融合后的角速度ω(k)。

所述子滤波器可以选择普通Kalman滤波、EKF或UKF等。

在多级序贯式滤波模块中设置有故障检测与隔离FDI模块，它可以检测对应的前一级子滤波器故障与否，从而判断该子滤波器是否参与多级滤波。FDI模块可选择状态残差χ²检验法等。

具体实施方式二：下面结合图2说明本实施方式，本实施方式对实施方式一作进一步说明，所述步骤二中获得每组陀螺阵列的融合数据Y(k)的具体方法为：

根据鲁棒理论中的尺度估计，定义相互支持度函数s_ij(k)和每组陀螺阵列中第i个MEMS陀螺的自身支持度函数s_i(k)为：

$s_{\mathrm{ij}}\left(k\right)=\left\{\begin{array}{cc}1-\sin [\frac{|Z_i\left(k\right)-Z_j\left(k\right)|}{\mathrm{\&beta;r}\left(k\right)}\&CenterDot;\frac{\mathrm{\&pi;}}{2}],& |Z_i\left(k\right)-Z_j\left(k\right)|\&le;\mathrm{\&beta;r}\left(k\right),\\ 0,& \mathrm{other}\end{array}\right.$

相互支持度函数s_ij(k)表示k时刻，每组陀螺阵列中MEMS陀螺i和MEMS陀螺j的相互支持度，其中i，j＝1，2，…，n，并且i≠j，

Z_i(k)为每组陀螺阵列中第i个MEMS陀螺在k时刻的输出值，Z_j(k)为每组陀螺阵列中第j个MEMS陀螺在k时刻的输出值，β为阈值尺度参数，r(k)为相互支持度阈值因子；

$s_i\left(k\right)=\left\{\begin{array}{cc}1-\sin [\frac{|Z_i\left(k\right)-\mathrm{median}\left(Z_i\right)|}{\mathrm{\&beta;}{r}_i}\&CenterDot;\frac{\mathrm{\&pi;}}{2}],& |Z_i\left(k\right)-\mathrm{median}\left(Z_i\right)|<{\mathrm{\&beta;r}}_i,\\ 0,& \mathrm{other}\end{array}\right.$

自身支持度函数s_i(k)表示在k时刻每组陀螺阵列中第i个MEMS陀螺采集获取的数据Z_i对其输出值Z_i(k)的支持程度，r_i为每组陀螺阵列中第i个MEMS陀螺的自身支持度阈值因子；

建立每组陀螺阵列的支持度矩阵D(k)：

$D\left(k\right)=\left[\begin{array}{ccccc}{d}_{11}\left(k\right)& d_{12}\left(k\right)& ...& d_{1n}\left(k\right)& d_{1(n+1)}\left(k\right)\\ d_{21}\left(k\right)& d_{22}\left(k\right)& ...& d_{2n}\left(k\right)& d_{2(n+1)}\left(k\right)\\ ...& ...& ...& ...& ...\\ d_{n1}\left(k\right)& d_{n2}\left(k\right)& ...& d_{\mathrm{nn}}\left(k\right)& d_{n(n+1)}\left(k\right)\\ d_{(n+1)1}\left(k\right)& d_{\left(n+1\right)2}\left(k\right)& ...& d_{(n+1)n}\left(k\right)& 1\end{array}\right],$

其中，

d_ii(k)＝s_i ²(k)，

d_ij(k)＝s_i(k)s_ij(k)s_j(k)，

$d_{i(n+1)}\left(k\right)=\left\{\begin{array}{cc}(1-\sin [\frac{|Z_i\left(k\right)-Y(k-1)|}{{\mathrm{\&beta;r}}_i}\&CenterDot;\frac{\mathrm{\&pi;}}{2}\left]\right)\&times;d_{\mathrm{ii}}\left(k\right),& |Z_i\left(k\right)-Y(k-1)|<{\mathrm{\&beta;r}}_i,\\ 0,& \mathrm{other}\end{array}\right.$

$d_{(n+1)j}\left(k\right)=\left\{\begin{array}{cc}(1-\sin [\frac{|Z_j\left(k\right)-Y(k-1)|}{{\mathrm{\&beta;r}}_j}\&CenterDot;\frac{\mathrm{\&pi;}}{2}\left]\right)\&times;d_{\mathrm{jj}}\left(k\right),& |Z_j\left(k\right)-Y(k-1)|<{\mathrm{\&beta;r}}_j,\\ 0,& \mathrm{other}\end{array}\right.$

其中，当k＝1时，Y(1)＝(Z₁(1)+Z₂(1)+…+Z_n(1))/n；

通过支持度矩阵D(k)的第i行判断获得每组陀螺阵列中第i个MEMS陀螺的权值系数w_i(k)，则权值向量Δ(k)＝[w₁(k)，w₂(k)，…，w_n+1(k)]^T；

设置一组非负数系数a₁(k)，a₂(k)，…，a_n+1(k)，使得：

w_i(k)＝a₁(k)d_i1(k)+a₂(k)d_i2(k)+…+a_n+1(k)d_i(n+1)(k)，

获得系数矩阵A(k)＝[a₁(k)，a₂(k)，…，a_n+1(k)]^T，

即Δ(k)＝D(k)A(k)；

支持度矩阵D(k)为非负对称矩阵，则D(k)的最大模特征值λ＞0，并有：

Δ(k)＝λA(k)；

每组陀螺阵列中第v个MEMS陀螺的权值系数w_v(k)与第l个MEMS陀螺的权值系数w_l(k)的关系为：

$\frac{w_v\left(k\right)}{w_l\left(k\right)}=\frac{a_v\left(k\right)}{a_l\left(k\right)}$ v，l＝1，2，…，n+1，

对每组陀螺阵列中第i个MEMS陀螺的权值系数w_i(k)归一化处理，得到：

$w_i\left(k\right)=\frac{a_i\left(k\right)}{a_1\left(k\right)+a_2\left(k\right)+\&CenterDot;\&CenterDot;\&CenterDot;+a_{n+1}\left(k\right)},$

则获得每组陀螺阵列的融合数据Y(k)：

$Y\left(k\right)=\frac{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{w}_i\left(k\right)Z_i\left(k\right)+w_{n+1}\left(k\right)Y(k-1)}{\underset{i=1}{\overset{n+1}{\mathrm{\&Sigma;}}}{w}_i\left(k\right)},$

其中，当k＝1时，Y(1)＝(Z₁(1)+Z₂(1)+…+Z_n(1))/n。

本实施方式中，；r(k)＝c·median|T(k)-median[T(k)]|，式中，median()为求中值运算；c＝1.4826；T(k)＝{Z_i(k)|i＝1，2，…，n}为k时刻每个陀螺阵列中所有MEMS陀螺测量值组成的集合；r_i(k)＝c·median|T_i(t)-median[T_i(t)]|，式中，T_i(t)＝{Z_i(t)|t＝1，2，…，k}为k时刻内陀螺阵列中第i个MEMS陀螺的所有测量值组成的集合；

获得每组陀螺阵列的融合数据Y(k)的过程如图2所示，其可以根据陀螺测量数据之间的差异性大小判断陀螺之间的相互支持程度，并由此推导出每个陀螺对应的权值，进行加权融合。

阈值尺度参数β一般选择为3，它可以调节尺度的大小。

实际中每组陀螺阵列中第i个MEMS陀螺的自身支持度阈值因子r_i的计算应在每一时刻下进行并不考虑权值为零的数据，为防止数据饱和，可引入限定记忆法或在足够长度的数据下计算。易知，s_i(k)∈(0，1]。s_i(k)随陀螺的精度增大而增大，当陀螺发生故障，偏差较大时，s_i(k)就会趋近于零。

具体实施方式三：本实施方式对实施方式二作进一步说明，步骤三中获得的时间序列模型为自回归滑动平均模型ARMA(p，q)，其中p为自回归模型的阶数，q为滑动平均模型的阶数；

自回归滑动平均模型ARMA(p，q)的表达式为：

y(k)＝b₁y(k-1)+…+b_py(k-p)+ε(k)-a₁ε(k-1)-…-a_qε(k-q)，

其中y(k)为时间序列模型在k时刻的响应，b_υ为自回归系数，υ＝1，2，…，p，ε(k)为随机扰动项白噪声，a_κ为滑动平均系数，κ＝1，2，…，q；

根据自回归滑动平均模型ARMA(p，q)建立状态空间方程，并且扩维引入真实角速度ω，扩维后的状态方程为：

X(k)＝Φ(k-1)X(k-1)+Γ(k-1)W(k)，

Y(k)＝H(k)X(k)+V(k)，

式中，X(k)为状态变量，Φ(k-1)为一步转移阵，Γ(k-1)为系统噪声驱动阵，W(k)为系统激励噪声序列，Y(k)由每组陀螺阵列的真实角速度ω和陀螺阵列的随机漂移η组成，H(k)为量测阵，V(k)为量测噪声序列。

以ARMA(2，1)为例，y(k)＝b₁y(k-1)+b₂y(k-2)+ε(k)-a₁ε(k-1)。其状态空间方程为：

$\left[\begin{array}{c}\mathrm{\&eta;}\left(k\right)\\ \mathrm{\&eta;}(k-1)\\ \mathrm{\&omega;}\left(k\right)\end{array}\right]=\left[\begin{array}{ccc}{b}_1& b_2& 0\\ 1& 0& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}\mathrm{\&eta;}(k-1)\\ \mathrm{\&eta;}(k-2)\\ \mathrm{\&omega;}(k-1)\end{array}\right]+\left[\begin{array}{ccc}1& -a_1& 0\\ 0& 0& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}\mathrm{\&epsiv;}\left(k\right)\\ \mathrm{\&epsiv;}(k-1)\\ \mathrm{\&xi;}\left(k\right)\end{array}\right],$

$Z\left(k\right)=\left[\begin{array}{ccc}1& 0& 1\end{array}\right]\left[\begin{array}{c}\mathrm{\&eta;}(k-1)\\ \mathrm{\&eta;}(k-2)\\ \mathrm{\&omega;}(k-1)\end{array}\right]+V\left(k\right).$

其中，ξ(k)为角速度ω(k)的驱动白噪声序列。

具体实施方式四：本实施方式对实施方式三作进一步说明，步骤四中获得第N级融合状态估计值的具体方法为：

多级序贯式滤波模块第h级融合在k时刻获得的状态估计值为：

${\hat{X}}_{f\left(h\right)}\left(k\right)=(1-{\mathrm{\&lambda;}}_{h-1}\left(k\right)){\hat{X}}_{f(h-1)}\left(k\right)+{\mathrm{\&lambda;}}_{h-1}\left(k\right){\hat{X}}_h\left(k\right),$

式中为子滤波器h的状态估计，h＝2，3，…，N；

多级序贯式滤波模块第h级融合在k时刻获得的估计误差方差P_f(h)为：

P_j(h)(k|k)＝[1-λ_h(k)]²P_f(h-1)(k|k)+[λ_h(k)]²P_h(k|k)+2[1-λ_h(k)]λ_h(k)P_f(h-1)，h(k|k)，

其中，P_h为子滤波器h的估计误差方差；λ_h(k)为最优估计系数；

${\mathrm{\&lambda;}}_h\left(k\right)=\frac{{\mathrm{trP}}_{f(h-1)}\left(k\right|k)-{\mathrm{trP}}_{f(h-1),h}\left(k\right|k)}{{\mathrm{trP}}_{f(h-1)}\left(k\right|k)+{\mathrm{trP}}_h\left(k\right|k)-2\mathrm{tr}{P}_{f(h-1),h}\left(k\right|k)},$

$P_{f(h-1),h}\left(k\right|k)=\underset{\mathrm{\&alpha;}=2}{\overset{h-1}{\mathrm{\&Sigma;}}}[1-{\mathrm{\&lambda;}}_{h-2}\left(k\right)]P_{f\left(1\right),h}\left(k\right|k)+\underset{\mathrm{\&alpha;}=2}{\overset{h-1}{\mathrm{\&Sigma;}}}\left[\underset{\mathrm{\&beta;}=\mathrm{\&alpha;}+1}{\overset{h-1}{\mathrm{\&Sigma;}}}\right[1-{\mathrm{\&lambda;}}_{\mathrm{\&beta;}}\left(k\right)\left]\right]{\mathrm{\&lambda;}}_n\left(k\right)P_{\mathrm{\&alpha;},h}\left(k\right|k),$

式中，tr表示矩阵的迹，P_f(i-1)，i为和的互协方差阵，将第一级子滤波器滤波后获得的状态估计值作为第一级融合数据则P_f(1)，h(k|k)＝P_1，h(k|k)，并且P_1，h(k|k)为：

P_1，h(k|k)＝[I-K₁(k)H₁(k)]×[Φ(k-1)P_1，h(k-1|k-1)Φ^T(k-1)+

Γ(k-1)Q(k-1)Γ^T(k-1)]×[I-K_h(k)H_h(k)]^T

式中，K_h(k)为k时刻多级序贯式滤波模块中第h个子滤波器的滤波增益，H_h(k)为k时刻多级序贯式滤波模块中第h个子滤波器的量测阵，Q(k-1)为k-1时刻多级序贯式滤波模块中子滤波器系统噪声序列的方差阵；设置P_1，h(0|0)的初值以获得P_1，h(k|k)，最后获得第N级融合状态估计值

具体实施方式五：下面结合图1至图3说明本实施方式，本实施方式对实施方式四作进一步说明，步骤五中获得批量MEMS陀螺k时刻融合后的角速度ω(k)的具体方法为：

$\mathrm{\&omega;}\left(k\right)=\left[\begin{array}{cccc}0& 0& ...& 1\end{array}\right]{\hat{X}}_g=\left[\begin{array}{cccc}0& 0& ...& 1\end{array}\right]\left[\begin{array}{c}\mathrm{\&eta;}\left(k\right)\\ \mathrm{\&eta;}(k-1)\\ .\\ .\\ .\\ \mathrm{\&eta;}(k-p+1)\\ \mathrm{\&omega;}\left(k\right)\end{array}\right].$

Claims (2)

1.一种批量MEMS陀螺信息融合方法，其特征在于，它包括以下步骤：

步骤一：将批量MEMS陀螺均分成N组陀螺阵列，每组陀螺阵列包括n个MEMS陀螺，n取值为4、5、6或7，N为大于2的正整数；

步骤二：采用初级融合模块对每组陀螺阵列中所有MEMS陀螺采集获取的数据进行基于支持度的初级融合，获得每组陀螺阵列的融合数据Y(k)，其中k表示时刻值；

步骤三：对每组陀螺阵列输出的融合数据Y(k)进行时间序列建模获得时间序列模型，并通过时间序列模型获得对应的状态空间方程；

步骤四：采用多级序贯式滤波模块对N组陀螺阵列的融合数据Y(k)和N个状态空间方程顺次进行多级序贯式滤波，其中每个陀螺阵列的融合数据Y(k)和状态空间方程作为多级序贯式滤波模块中对应的子滤波器的输入，获得第N级融合状态估计值

步骤五：将第N级融合状态估计值输入角速度提取模块，经提取获得批量MEMS陀螺k时刻融合后的角速度ω(k)；

所述步骤二中获得每组陀螺阵列的融合数据Y(k)的具体方法为：

根据鲁棒理论中的尺度估计，定义相互支持度函数s_ij(k)和每组陀螺阵列中第i个MEMS陀螺的自身支持度函数s_i(k)为：

$s_{\mathrm{ij}}\left(k\right)=\left\{\begin{array}{cc}1-\sin [\frac{|Z_i\left(k\right)-Z_j\left(k\right)|}{\mathrm{\&beta;r}\left(k\right)}\&CenterDot;\frac{\mathrm{\&pi;}}{2}],& |Z_i\left(k\right)-Z_j\left(k\right)|<\mathrm{\&beta;r}\left(k\right)\\ 0,& \mathrm{other}\end{array}\right.,$

相互支持度函数s_ij(k)表示k时刻，每组陀螺阵列中MEMS陀螺i和MEMS陀螺j的相互支持度，其中i,j＝1,2,…,n，并且i≠j，

Z_i(k)为每组陀螺阵列中第i个MEMS陀螺在k时刻的输出值，Z_j(k)为每组陀螺阵列中第j个MEMS陀螺在k时刻的输出值，β为阈值尺度参数，r(k)为相互支持度阈值因子；

$s_i\left(k\right)=\left\{\begin{array}{cc}1-\sin [\frac{|Z_i\left(k\right)-\mathrm{median}\left(Z_i\right)|}{\mathrm{\&beta;}{r}_i}\&CenterDot;\frac{\mathrm{\&pi;}}{2}],& |Z_i\left(k\right)-\mathrm{median}\left(Z_i\right)|<\mathrm{\&beta;}{r}_i\\ 0,& \mathrm{other}\end{array}\right.,$

自身支持度函数s_i(k)表示在k时刻每组陀螺阵列中第i个MEMS陀螺采集获取的数据Z_i对其输出值Z_i(k)的支持程度，r_i为每组陀螺阵列中第i个MEMS陀螺的自身支持度阈值因子；

建立每组陀螺阵列的支持度矩阵D(k)：

$D\left(k\right)=\left[\begin{array}{ccccc}{d}_{11}\left(k\right)& d_{12}\left(k\right)& ...& d_{1n}\left(k\right)& d_{1(n+1)}\left(k\right)\\ d_{21}\left(k\right)& d_{22}\left(k\right)& ...& d_{2n}\left(k\right)& d_{2(n+1)}\left(k\right)\\ ...& ...& ...& ...& ...\\ d_{n1}\left(k\right)& d_{n2}\left(k\right)& ...& d_{\mathrm{nn}}\left(k\right)& d_{n(n+1)}\left(k\right)\\ d_{(n+1)1}\left(k\right)& d_{(n+1)2}\left(k\right)& ...& d_{(n+1)n}\left(k\right)& 1\end{array}\right],$

其中，

d_ii(k)＝s_i ²(k)，

d_ij(k)＝s_i(k)s_ij(k)s_j(k)，

$d_{i(n+1)}\left(k\right)=\left\{\begin{array}{cc}(1-\sin [\frac{|Z_i\left(k\right)-Y(k-1)|}{\mathrm{\&beta;}{r}_i}\&CenterDot;\frac{\mathrm{\&pi;}}{2}\left]\right)\&times;d_{\mathrm{ii}}\left(k\right),& |Z_i\left(k\right)-Y(k-1)|<\mathrm{\&beta;}{r}_i\\ 0,& \mathrm{other}\end{array}\right.,$

$d_{(n+1)j}\left(k\right)=\left\{\begin{array}{cc}(1-\sin [\frac{|Z_j\left(k\right)-Y(k-1)|}{\mathrm{\&beta;}{r}_j}\&CenterDot;\frac{\mathrm{\&pi;}}{2}\left]\right)\&times;d_{\mathrm{jj}}\left(k\right),& |Z_j\left(k\right)-Y(k-1)|<\mathrm{\&beta;}{r}_j\\ 0,& \mathrm{other}\end{array}\right.,$

其中，当k＝1时，Y(1)＝(Z₁(1)+Z₂(1)+…+Z_n(1))/n；

通过支持度矩阵D(k)的第i行判断获得每组陀螺阵列中第i个MEMS陀螺的权值系数w_i(k)，则权值向量Δ(k)＝[w₁(k),w₂(k),…,w_n+1(k)]^T；

设置一组非负数系数a₁(k),a₂(k),…,a_n+1(k)，使得：

w_i(k)＝a₁(k)d_i1(k)+a₂(k)d_i2(k)+…+a_n+1(k)d_i(n+1)(k)，

获得系数矩阵A(k)＝[a₁(k),a₂(k),…,a_n+1(k)]^T，

即Δ(k)＝D(k)A(k)；

支持度矩阵D(k)为非负对称矩阵，则D(k)的最大模特征值λ>0，并有：

Δ(k)＝λA(k)；

每组陀螺阵列中第v个MEMS陀螺的权值系数w_v(k)与第l个MEMS陀螺的权值系数w_l(k)的关系为：

$\frac{w_v\left(k\right)}{w_l\left(k\right)}=\frac{a_v\left(k\right)}{a_l\left(k\right)},v,l=\mathrm{1,2},...,n+1,$

对每组陀螺阵列中第i个MEMS陀螺的权值系数w_i(k)归一化处理，得到：

$w_i\left(k\right)=\frac{a_i\left(k\right)}{a_1\left(k\right)+a_2\left(k\right)+...+a_{n+1}\left(k\right)},$

则获得每组陀螺阵列的融合数据Y(k)：

$Y\left(k\right)=\frac{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{w}_i\left(k\right)Z_i\left(k\right)+w_{n+1}\left(k\right)Y(k-1)}{\underset{i=1}{\overset{n+1}{\mathrm{\&Sigma;}}}{w}_i\left(k\right)},$

其中，当k＝1时，Y(1)＝(Z₁(1)+Z₂(1)+…+Z_n(1))/n；

步骤三中获得的时间序列模型为自回归滑动平均模型ARMA(p,q)，其中p为自回归模型的阶数，q为滑动平均模型的阶数；

自回归滑动平均模型ARMA(p,q)的表达式为：

y(k)＝b₁y(k-1)+…+b_py(k-p)+ε(k)-a₁ε(k-1)-…-a_qε(k-q)，

其中y(k)为时间序列模型在k时刻的响应，b_υ为自回归系数，υ＝1,2,…,p，ε(k)为随机扰动项白噪声，a_κ为滑动平均系数，κ＝1,2,…,q；

根据自回归滑动平均模型ARMA(p,q)建立状态空间方程，并且扩维引入真实角速度ω，扩维后的状态方程为：

X(k)＝Φ(k-1)X(k-1)+Γ(k-1)W(k)，

Y(k)＝H(k)X(k)+V(k)，

式中，X(k)为状态变量，Φ(k-1)为一步转移阵，Γ(k-1)为系统噪声驱动阵，W(k)为系统激励噪声序列，Y(k)由每组陀螺阵列的真实角速度ω和陀螺阵列的随机漂移η组成，H(k)为量测阵，V(k)为量测噪声序列；

步骤四中获得第N级融合状态估计值的具体方法为：

多级序贯式滤波模块第h级融合在k时刻获得的状态估计值为：

${\hat{X}}_{f\left(h\right)}\left(k\right)=(1-{\mathrm{\&lambda;}}_{h-1}\left(k\right)){\hat{X}}_{f(h-1)}\left(k\right)+{\mathrm{\&lambda;}}_{h-1}\left(k\right){\hat{X}}_h\left(k\right),$

式中为子滤波器h的状态估计，h＝2,3,…,N；

多级序贯式滤波模块第h级融合在k时刻获得的估计误差方差P_f(h)为：

P_f(h)(k|k)＝[1-λ_h(k)]²P_f(h-1)(k|k)+[λ_h(k)]²P_h(k|k)+2[1-λ_h(k)]λ_h(k)P_f(h-1),h(k|k)，

其中，P_h为子滤波器h的估计误差方差；λ_h(k)为最优估计系数；

${\mathrm{\&lambda;}}_h\left(k\right)=\frac{{\mathrm{trP}}_{f(h-1)}\left(k\right|k)-\mathrm{tr}{P}_{f(h-1),h}\left(k\right|k)}{\mathrm{tr}{P}_{f(h-1)}\left(k\right|k)+\mathrm{tr}{P}_h\left(k\right|k)-2\mathrm{tr}{P}_{f(h-1),h}\left(k\right|k)},$

$P_{f(h-1),h}\left(k\right|k)\underset{\mathrm{\&alpha;}=2}{\overset{h-1}{\mathrm{\&Pi;}}}[1-{\mathrm{\&lambda;}}_{h-2}\left(k\right)]P_{f\left(1\right),h}\left(k\right|k)+\underset{\mathrm{\&alpha;}=2}{\overset{h-1}{\mathrm{\&Sigma;}}}\left[\underset{\mathrm{\&beta;}=\mathrm{\&alpha;}+1}{\overset{h-1}{\mathrm{\&Pi;}}}\right[1-{\mathrm{\&lambda;}}_{\mathrm{\&beta;}}\left(k\right)\left]\right]{\mathrm{\&lambda;}}_h\left(k\right)P_{\mathrm{\&alpha;},h}\left(k\right|k),$

式中，tr表示矩阵的迹，P_f(i-1),i为和的互协方差阵，将第一级子滤波器滤波后获得的状态估计值作为第一级融合数据则P_f(1),h(k|k)＝P_1,h(k|k)，并且P_1,h(k|k)为：

P_1,h(k|k)＝[I-K₁(k)H₁(k)]×[Φ(k-1)P_1,h(k-1|k-1)Φ^T(k-1)+，

Γ(k-1)Q(k-1)Γ^T(k-1)]×[I-K_h(k)H_h(k)]^T

式中，K_h(k)为k时刻多级序贯式滤波模块中第h个子滤波器的滤波增益，H_h(k)为k时刻多级序贯式滤波模块中第h个子滤波器的量测阵，Q(k-1)为k-1时刻多级序贯式滤波模块中子滤波器系统噪声序列的方差阵；设置P_1,h(0|0)的初值以获得P_1,h(k|k)，最后获得第N级融合状态估计值

2.根据权利要求1所述的批量MEMS陀螺信息融合方法，其特征在于，步骤五中获得批量MEMS陀螺k时刻融合后的角速度ω(k)的具体方法为：

$\mathrm{\&omega;}\left(k\right)=\left[\begin{array}{cccc}0& 0& ...& 1\end{array}\right]{\hat{X}}_g=\left[\begin{array}{cccc}0& 0& ...& 1\end{array}\right]\left[\begin{array}{c}\mathrm{\&eta;}\left(k\right)\\ \mathrm{\&eta;}(k-1)\\ .\\ .\\ .\\ \mathrm{\&eta;}(k-p+1)\\ \mathrm{\&omega;}\left(k\right)\end{array}\right].$