基于极小化矢量距离准则的有限模型滤波方法

摘要

本发明克服了在惯性系统的惯性器件的随机模型不能精确获知、或者不能用单一线性模型描述时的惯性器件随机误差降噪时经典卡尔曼滤波算法不能应用的困难，提出一种基于极小化矢量距离准则的有限模型滤波方法，步骤一、采集数据并对数据进行预处理；步骤二、对步骤一中预处理后的时间序列数据建模并整定模型参数；步骤三、对每个系统模型设计卡尔曼滤波器并应用基于最小化矢量距离的有限模型算法对模型进行在线实时切换。

说明书

技术领域

本发明涉及一种有限模型滤波方法，属于自适应滤波以及陀螺随机误差降噪滤波领域。

背景技术

惯性导航系统具有完全自主性、全天候、抗外接干扰等优点，作为惯性导航系统的核心 器件——陀螺仪在整个惯性导航系统中起着十分重要的作用。惯性导航系统的精度、成本主 要决定于惯性仪表（陀螺仪和加速度计）的精度和成本，尤其是陀螺仪的漂移对惯性导航系 统位置误差增长的影响是时间的三次方函数。对于MEMS陀螺，影响其精度的一个重要因素就 是MEMS陀螺的随机误差。随机误差具有不确定性，不能像确定性误差那样通过简单的方法进 行补偿，它是MEMS陀螺的主要误差，因此通过减小随机误差从而提高MEMS陀螺测量精度便 尤为重要。

卡尔曼滤波理论自1960年由卡尔曼提出后，经过50多年的发展，针对不同的应用背景， 卡尔曼滤波理论已经在不同的工程领域得到了相应的理论推广和应用。卡尔曼滤波是一种时 域滤波方法，采用状态空间方法描述系统，从与被提取信号有关的量测量中通过算法估计出 所需信号。其中被估计信号是由白噪声激励引起的随机响应，激励源与响应之间的传递结构 (系统方程)已知，量测量与被估计量之间的函数关系(量测方程)也已知。估计过程利用了如 下信息：系统方程、量测方程、白噪声激励的统计特性、量测误差的统计特性。

卡尔曼滤波算法是MEMS陀螺随机误差降噪的常用算法。但由于应用卡尔曼滤波算法需要 预先知道精确的系统模型和噪声的统计特性才能获得最优的滤波估计，这种要求在实际的应 用过程中限制了算法的应用。当系统模型无法精确获知、或者不能用单一的线性模型对系统 描述时，经典的卡尔曼滤波算法就不再适用。

发明内容

本发明的目的是为了克服在惯性系统的惯性器件的随机模型不能精确获知、或者不能用 单一线性模型描述时的惯性器件随机误差降噪时经典卡尔曼滤波算法不能应用的困难，提出 基于极小化矢量距离准则的有限模型滤波方法。

本发明的目的是通过下述技术方案实现的。

该基于极小化矢量距离准则的有限模型滤波方法，包括如下步骤：

步骤一、采集数据并对数据进行预处理；首先进行异常值剔除，然后进行零均值处理， 其次为去除趋势项部分，最后是提取周期项形成时间序列数据；

步骤二、对步骤一中预处理后的时间序列数据建模并整定模型参数。

2.1建立系统模型：时间序列线性模型的参数估计是指在辨识得到模型类别和阶数的基础 上，求出模型中的自回归系数和滑动平均系数的数值；

2.2检验模型的适用性：时间序列线性模型的适用性检验是指用样本的数据检验按上面 的方法确定的模型是否适用，即检验残差序列是否为白噪声序列；

步骤三、对每个系统模型设计卡尔曼滤波器并应用基于最小化矢量距离的有限模型算法 对模型进行在线实时切换：

3.1对每一个建立的系统模型设计卡尔曼滤波器，在辨识出的系统模型的基础上给定系 统噪声和观测噪声阵；

3.2对每一个系统模型根据基于极小化矢量距离准则拟合评价函数来评定系统模型的准 确度；

3.3根据上述给出的评价函数将待选择的模型中准确度最高的系统模型的滤波结果作为 最终的结果输出，

步骤3.3中通过精确度结合贝叶斯概率将所有系统模型滤波结果的融合结果作为结果输 出。

本发明的基本原理是：对于惯性器件的随机误差，通过不同的模型机构、或不同的机理 可以建立不同的系统模型，在所建立的不同模型中每个模型都有其擅长处理的动态特性，在 惯性器件实时测量过程中由于其系统的复杂性，不可能用单一的线性系统模型对系统完成精 确描述；因此，我们建立一个模型池，其中包括有限个不同的系统模型，在采集陀螺动态输 出的数据时，根据实际数据与模型池中备选模型的匹配情况，动态地选择出与实际数据最为 匹配的模型用来进行实时在线滤波估计。借助本发明提出的滤波机制，不需要预先知道陀螺 随机误差的精确系统模型，其本质是利用实际系统量测到的输出数据，用模型池中的有限个 模型来动态逼近实际系统，达到在不知系统真实模型情况下进行最优滤波估计的效果。为了 在滤波过程中选择需激活的模型，需要对基于模型池中各个模型对应的子滤波器的估计输出 进行评估、决择，本发明给出了一种对模型选择的准则，以此来完成对模型的选择，并且给 出几种常用的系统模型建立的方法。

本发明的有益效果：

本发明相对于传统的卡尔曼滤波算法，降低了对模型精确性的要求，可以针对不同的系 统动态特性建立多个存在较大差别的线性系统模型，且可以对类似的模型选取不同的模型参 数来提高系统的鲁棒性。且在MEMS惯性器件的算法验证中看出，该算法由于其计算简单且具 有良好的实时性，利于工程应用和实现。

附图说明

图1为本发明提出的有限模型滤波算法的系统框图；

图2为本发明MEMS静态信号的降噪测试结果图；

图3为本发明MEMS常值角速率输入信号的降噪测试结果图；

图4为本发明MEMS高动态低幅值信号的降噪测试结果图；

图5为本发明MEMS中动态高幅值信号的降噪测试结果图。

具体实施方式

以下所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说， 在不脱离本发明原理的前提下，还可以做出若干改进，或者对其中部分技术特征进行等同替 换，这些改进和替换也应视为本发明的保护范围。

步骤一、通过器件的采集数据，对数据进行预处理（采用经典的数据与处理的方法）。

在对陀螺随机漂移数据进行模型辨识之前，需要进行数据预处理，使之成为零均值、平 稳、正态的时序信号，之后才能对信号建立数学模型。

在数据预处理中，第一步是异常值剔除。剔除方法可以采用莱以特准则。莱以特准则 又称<3σ>准则，根据误差理论，

P{|x-μ|≤3σ}≈99.7%    （1）

该式说明，误差ε＝|x-μ|≤3σ的概率约为99.7%，如果将ε＞3σ的值舍去，犯“弃真” 错误的概率最大为0.3%。

数据预处理的第二步是零均值处理。对于有限长的时间序列x_k(k＝1,2,…N)计算其均值。 陀螺数据求出其平均值后，将陀螺每一时刻的数据都减去平均值，即可得到零均值处理后的 数据。

第三步是去除趋势项部分。实际工程中测得的陀螺漂移数据序列往往为非平稳随机序列。 对此，应去除其中的有规律部分即趋势项，实现数据的平稳化处理。一个非平稳时间序列y_k可认为是一个确定性的趋势项A_k和一个均值为零的高度相关平稳时间序列x_k的线性叠加， 即：

y_k＝A_k+x_k    (2)

式中，k＝1,2,…N。

其中，趋势项A_k一般可表示为时间k的多项式，即：

A_k＝a₀+a₁k+a₂^k2+…+a_mk^m    (3)

式中，a₀,a₁,…a_m为多项式的系数。

对含趋势项的非平稳序列，首先应进行趋向性检验，然后再进行趋势项的提取。

趋势项检验中，可以采用下式作为检测量。

$u=\frac{A+\frac{1}{2}-\frac{K(K-1)}{4}}{\sqrt{\frac{{2K}^3+{3K}^2-5K}{72}}}---\left(4\right)$

该检测量理论上渐进服从N(0,1)分布。

对序列通过对上式的计算，得到统计量u，在一定的显著性水平下即可判断是否可接受 序列的趋向性存在。

若存在趋势项，可以应用最小二乘法求取趋势项的多项式系数。补偿后进行趋势项检测， 若仍存在趋势项则升高多项式阶次，重新求估计值多项式系数，直到检验结果无趋向为止， 从而求得趋势项的系数。

数据预处理的第四步是提取周期项。陀螺漂移序列经过趋势项检验并提取趋势向后，如 果漂移数据中仍存在着隐含周期项，则漂移数据从传统意义上说还不能视为平稳，需要进一 步识别并提取序列中的隐含周期项。周期图分析方法不仅能识别出隐含周期项的周期，而且 还能识别出幅值和相位，因此在陀螺漂移的处理中采用周期图分析方法。

步骤二、对预处理后的时间序列数据建模（模型为自回归滑动平均模型），并整定模型参 数。经预处理后，如果检验得到的序列为平稳、正态、零均值的平稳随机时间序列，就可以 对其建立时间序列模型。

数据建模的第一步是建立系统模型。时间序列线性模型的参数估计，是指在辨识得到模 型类别和阶数的基础上，求出模型中的自回归系数和滑动平均系数的数值。

将AR(n)模型：x_k+a₁x_k-1+…+a_nx_k-n＝ε_k改写成(设μ＝0)

x_k＝φ₁x_k-1+…+φ_nx_k-n+ε_k    (5)

已知x₁,x₂,…,x_N，求φ₁,φ₂,…,φ_n的LS估计。

表示成矩阵形式：Y＝Xβ+e，式中：

参数β的最小二乘估计为：$\hat{\beta }={\left(X^{T}X\right)}^{-1}{X}^{T}Y---\left(7\right)$

利用最小二乘法可得出AR模型参数估计的具体表达式为：

本发明采用的是时序理论估计法中的参数比较估计的长自回归模型法。其基本思想如下： 基于观测时序建立起来的AR模型、MA模型、ARMA模型均是等价的数学模型，因此，由这些 模型确定的等价系统的传递函数实际上是相等的。所以，可先估计出AR模型，再根据传递函 数相等的关系估计出ARMA模型的参数φ_i与θ_j。

建模的第二部是检验模型的适用性。时间序列线性模型的适用性检验，是指用样本的数 据检验按上面的方法确定的模型是否适用。在估计出模型参数后进行实用性检验，其实质是 检验残差序列是否为白噪声序列。检验方法主要有χ²检验和F检验。本发明采用的方法是χ²检验。

假设检验是对应于一定的显著性水平α下，对残差序列{w_k}能否视为白噪声序列作出判 断。对应于一点的显著性水平，由χ²分布表可查得自由度为K的χ²变量值界限值若 是满足条件：

$Q\le {\chi }_{K:\alpha }^2---\left(9\right)$

则在显著性水平α下，{w_k}是白噪声序列这一假设成立，从而认为所确定的模型是适用 的。如果则应作出相反的判定，从而认为所确定的模型不适用。在此情形下，则 需要改变模型的类别或阶数，重新进行参数估计，并对模型进行适用性检验。

步骤三、对每个系统模型设计卡尔曼滤波器并应用基于最小化矢量距离的有限模型算法 对模型进行在线实时切换：

第一步是：对每一个建立的系统模型设计卡尔曼滤波器，主要是在辨识出的系统模型的 基础上给定系统噪声和观测噪声阵。

通用的线性系统模型如下：

$\left(\begin{array}{c}{x}_k={\phi }_{k,k-1}{x}_{k-1}+{\mathrm{\Gamma \omega }}_K\\ z_k=H_k{x}_k+v_k\end{array}\right)---\left(10\right)$

假设现在模型池中有N个系统模型，且最优逼近线性模型γ∈M，其中 M□{M₁，M₂，…M_N}。M即为建立的系统模型池，模型池中的模型并不要求必须有相似的形 式，可以有较大的差别。不过模型池中的每个模型必须包含需要估计的公共状态。

第二步是：对于系统中的每一个系统，基于极小化矢量距离准则拟合评价函数评定系统 模型的准确度。对于每一个模型M_i，在第k步卡尔曼滤波估计器给出状态的估计输出相对应的我们可以定义滤波器的观测估计

$z_k^{\left(i\right)}=H_k{\hat{x}}_k^{\left(i\right)}---\left(11\right)$

假设存在一个函数d_i(k)，定义如下：

$d_i\left(k\right)={\left|\right|f_{z_k}-f_{{\hat{z}}_k^{\left(i\right)}}\left|\right|}^2---\left(12\right)$

其中是量测序列{z_k}的函数，是滤波器观测估计序列的函数。d_i(k)函数可 以用来刻画模型M_i距离最优线性逼近模型的矢量距离。d_i(k)的值越小说明待选择模型M_i距 离最优模型越接近。

对于第k步模型的选择，可以采用以下准则：

$i_k=\underset{1\le i\le N}{\arg \min }\left|\right|d_i\left(k\right)\left|\right|---\left(13\right)$

下面就给出d_i(k)的一种形式。

首先给出定义的形式。设j＝1,2,…,2Δ+1为一个选定的时域小视窗。定义z_k(j)□z_k-j为选定时域小视窗内的观则数据。给出的二次多项式形式：

$f_{z_k}\left(j\right)={\alpha }_0+{\alpha }_{1}j+{\alpha }_2{j}^2---\left(14\right)$

其中α₀,α₁和α₂是待定参数，可用最小二乘的形式给出。即极小化下式：

$\underset{j=1}{\overset{2\Delta +1}{\Sigma }}{(z_k\left(j\right)-({\alpha }_0+{\alpha }_{1}j+{\alpha }_2{j}^2))}^2---\left(15\right)$

根据上述形式可给出如下等式：

$\left(\begin{array}{ccc}1& 1& 1\\ 1& 2& 4\\ .& .& .\\ .& .& .\\ .& .& .\\ 1& 2\Delta +1& {(2\Delta +1)}^2\end{array}\right)\left(\begin{array}{c}{\alpha }_0\\ {\alpha }_1\\ {\alpha }_2\end{array}\right)=\left(\begin{array}{c}{z}_k\left(1\right)\\ z_k\left(2\right)\\ .\\ .\\ .\\ z_k(2\Delta +1)\end{array}\right)---\left(16\right)$

定义：

根据最小二乘拟合算法可得(α₀ α₁ α₂)^T＝(M^TM)^-1M^TZ_k，然后可得

$f_{z_k}=M{\left(M^{T}M\right)}^{-1}{M}^T{Z}_k---\left(19\right)$

类似的可以给出的二次项逼近表达式：

$f_{{\hat{z}}_k}\left(j\right)={\hat{\alpha }}_0+{\hat{\alpha }}_{1}j+{\hat{\alpha }}_2{j}^2---\left(20\right)$

定义：

可得：

$f_{{\hat{z}}_k}=M{\left(M^{T}M\right)}^{-1}{M}^T{\hat{Z}}_k---\left(22\right)$

通过上述表述方法的变换可得：

$d_i\left(k\right)={\left|\right|f_{z_k}-f_{{\hat{z}}_k}\left|\right|}^2={\left|M{\left(M^{T}M\right)}^{-1}{M}^T(Z_k-{\hat{Z}}_k)\right|}^2---\left(23\right)$

第三步是：根据上述给出的评价函数，将待选择的模型中准确度最高的系统模型的滤波 结果作为最终的结果输出，亦可以通过精确度结合贝叶斯概率将所有系统模型滤波结果的融 合结果作为结果输出。

下面以两模型情况作为例子，来给出有限模型滤波算法的滤波过程：滤波过程大致可以 分为两部分：

a)滤波初始化

$P_0^{\left(1\right)}=P_0$

$P_0^{\left(2\right)}=P_0---\left(24\right)$

${\hat{x}}_0={\bar{x}}_0$

其中P₀是一个正定矩阵是滤波初值的估计方差阵。是初值的数学期望值，即 （如果没有x₀的先验信息也可以设置）。

b)滤波估计

第二步是根据每一个线性模型和设计好的滤波器并行滤波估计给出每一个滤波器的状态 估计。

${\hat{x}}_{k,k-1}^{\left(1\right)}={\phi }_{k,k-1}^{\left(1\right)}{\hat{x}}_{k-1}$

$P_{k,k-1}^{\left(1\right)}={\phi }_{k,k-1}^{\left(1\right)}{P}_{k-1}{\phi }_{k,k-1}^{\left(1\right)T}+{\Gamma }_{k-1}{Q}_{k-1}{\Gamma }_{k-1}^T$

$K_k^{\left(1\right)}=P_{k,k-1}^{\left(1\right)}{H}_k^T{(H_k{P}_{k,k-1}^{\left(1\right)}{H}_k^{\left(T\right)}+R_k)}^{-1}$(25)

$P_k^{\left(1\right)}=(I-K_k^{\left(1\right)}{H}_k)P_{k,k-1}^{\left(1\right)}$

${\hat{x}}_k^{\left(1\right)}={\hat{x}}_{k,k-1}^{\left(1\right)}+K_k^{\left(1\right)}(z_k-H_k{\hat{x}}_{k,k-1}^{\left(1\right)})$

${\hat{z}}_k^{\left(1\right)}=H_k{\hat{x}}_k^{\left(1\right)}$

${\hat{x}}_{k,k-1}^{\left(2\right)}={\phi }_{k,k-1}^{\left(2\right)}{\hat{x}}_{k-1}$

$P_{k,k-1}^{\left(2\right)}={\phi }_{k,k-1}^{\left(2\right)}{P}_{k-1}{\phi }_{k,k-1}^{\left(2\right)T}+{\Gamma }_{k-1}{Q}_{k-1}{\Gamma }_{k-1}^T$

$K_k^{\left(2\right)}=P_{k,k-1}^{\left(2\right)}{H}_k^T{(H_k{P}_{k,k-1}^{\left(2\right)}{H}_k^{\left(T\right)}+R_k)}^{-1}$(26)

$P_k^{\left(2\right)}=(I-K_k^{\left(2\right)}{H}_k)P_{k,k-1}^{\left(2\right)}$

${\hat{x}}_k^{\left(2\right)}={\hat{x}}_{k,k-1}^{\left(2\right)}+K_k^{\left(2\right)}(z_k-H_k{\hat{x}}_{k,k-1}^{\left(2\right)})$

${\hat{z}}_k^{\left(2\right)}=H_k{\hat{x}}_k^{\left(2\right)}$

然后通过不同的滤波器输出的观则估计计算函数d_i(k)：

$d_i\left(k\right)={\left|\right|f_{z_k}-f_{{\hat{z}}_k}\left|\right|}^2(i=\mathrm{1,2})---\left(27\right)$

然后通过判断d_i(k)的函数值来完成模型的选择。即，如果d₁(k)≤d₂(k)则i_k＝1，并且

$P_k=P_k^{\left(1\right)}$(28)

${\hat{x}}_k={\hat{x}}_k^{\left(1\right)}$

如果d₁(k)＞d₂(k)则i_k＝2，并且

$P_k=P_k^{\left(2\right)}$(29)

${\hat{x}}_k={\hat{x}}_k^{\left(2\right)}$

下面的内容给出上述有限模型算法在MEMS去噪中的半实物实验验证。在算法验证试验 中，选取的是MEMS130陀螺作为实验惯性器件，测试温度为25°C，采样时间为T＝5ms，以 陀螺输出的角速率信号为观测数据。建立的模型为两个模型，其中一个是ARMA(2,2)模型， 另外一个是基于陀螺的物理学动态特性建立的模型。

测试的内容主要分为四个部分：第一个部分是静态测试，即当陀螺无角速率输入时，角 速率输出信号的降噪；第二个部分是常值角速率输入测试，即陀螺的有常值角速率输入时， 输出信号的降噪滤波；第三部分是高动态低幅值测试，即输入角速率为周期信号（信号周期 短、幅值低，且频率和幅值都在MEMS130的测量范围内），测试滤波算法的实时特性；第四 部分是中动态高幅值的测试，即输入角速率为另一周期信号（信号周期中、幅值中，且频率 和幅值都在MEMS130的测量范围内），测试滤波算法的幅值特性。