杨金显, 陈超, 李志鹏

(河南理工大学电气工程与自动化学院, 河南 焦作 454000)

0 引 言

基于MEMS惯性器件构建低成本的微型惯性测量组合(Miniature Inertial Measurement Unit,MIMU)其体积小、质量轻、价格低、抗震动冲击能力强等优点,越来越受到重视。但是,其系统中存在多个误差源,包括惯性元器件的安装误差、测量误差、初始对准误差以及计算误差,这些误差大大影响了系统信号的准确性与实时性。其中惯性元器件的测量误差对系统误差影响最为明显,而惯导单元测量误差主要由陀螺漂移产生。所以,陀螺仪的测量误差决定了系统的精度[1-2]。为解决陀螺仪测量过程中容易产生漂移的问题,近年来不少学者提出利用小波算法或神经网络理论对陀螺仪进行降噪处理[3-4],以减少随机误差对陀螺漂移的影响。但是小波算法考察的是整个时域过程的频域特征或整个频域过程的时域特征,对于平稳过程,有很好的效果,但对于非平稳过程却有诸多不足。鉴于陀螺仪的数据分析是动态的非平稳性的,其并不适合。文献[5-6]提出在时间序列模型的基础上建立AR模型,对数据进行滤波处理,但是对不能使用低阶AR模型描述的陀螺漂移,文献[5-6]的方法滤波效果并不好。

为了克服小波算法的局限性,本文提出一种基于ARIMA模型的MEMS陀螺仪信号卡尔曼滤波处理算法。相对于AR模型,ARIMA模型运用差分方法去除了时间序列的趋势等因素,并不需要时间序列是均值为0的静态时间序列的条件,它的适用性得到了很大的提高。

1 ARIMA模型的判别

按照时间顺序把随机事件变化发展的过程记录下来就构成了一个时间序列。对时间序列进行观察、研究,找寻它变化发展的规律,预测它将来的走势就是时间序列分析。其中分为平稳序列和非平稳序列,而ARIMA模型适用于差分平稳序列拟合。

ARIMA模型全称为差分自回归移动平均模型。其中,p为自回归AR的自回归项;q为移动平均MA的移动平均项数;d为时间序列成为平稳时所做的差分次数。ARIMA模型的基本思想:将预测对象随时间推移而形成的数据序列视为一个随机序列,用一定的数学模型来近似描述这个序列。

1.1 模型结构

ARIMA(p,d,q)模型差分算子为

Δxt=xt-xt-1=xt-lxt=(1-l)xt

Δ2xt=Δxt-Δxt-1=(1-l)xt-

(1-l)xt-1=(1-l)2xt

Δdxt=(1-l)dxt

(1)

对d阶单整序列xt

wt=Δdxt=(1-l)dxt

(2)

若wt是平稳序列,可对wt建立ARMA(p,q)模型,所得到的模型称为ARIMA(p,d,q),模型形式

wt=φ1wt-1+φ2wt-2+…+φpwt-p+δ+

ut+θ1ut-1+θ2ut-2+…+θqut-q

φ(l)Δdxt=δ+Θ(l)ut

(3)

式中,ut为白噪声序列;δ是常数。由此便可转化为ARMA模型。

1.2 样本自相关系数和偏自相关系数的识别

(1) 平稳AR(p)的自相关系数

xt=φ1xt-1+φ2xt-2+…+φpxt-p+ut

|φi|<1,i=1,2,…,p,E(ut)=0

(4)

xt-kxt=φ1xt-kxt-1+…+φpxt-kxt-p+xt-kxt

(5)

γk=φ1λk-1+φ2λk-2+…+φpλk-p,k>0

(6)

平稳AR(p)的自相关系数是

pk=φ1pk-1+φ2pk-2+…+φppk-p,k>0

(7)

(2)k阶平稳自回归过程AR(k)偏自相关系数

xt=φk1xt-1+φk2xt-2+…+φkkxt-k+ut

(8)

xtxt-j=φk1xt-1xt-j+φk2xt-2xt-j+…+

φkkxt-kxt-j+utxt-j

(9)

γj=φk1γj-1+φk2γj-2+…+φkkγj-k

(10)

两边同除以γ0

pj=φk1pj-1+φk2pj-2+…+φkkpj-k

(11)

对任意j>0都成立。根据ρ0=1和对称性ρj=ρ-j,得到Yule-Walker方程组

p1=φk1+φk2p1+…+φkkpk-1

p2=φk1p1+φk2+…+φkkpk-2

⋮

pk=φk1pk-1+φk2pkk-2+…+φkk

(12)

对于给定的k,p1,p2,…,pk已知,每个方程组最后一个解就是相应的偏自相关系数φ11,φ22,…,φkk。

(3) MA(q)自相关系数

xt=μ+ut+θ1ut-1+θ2ut-2+…+θqut-q

(13)

γk=E(xtxt-k)=

σ2(θk+θ1θk+1+…+θk-qθq)0

0k>q

(14)

pk=rkr0=1k=0

0k>q

(15)

当k>q时,pk=0,xt与xt+k不相关,这种现象称为截尾,可根据自相关系数是否从某一点开始一直为0判断MA(q)模型的阶数q。

(4) MA(q)偏自相关系数。MA(q)模型对应一个AR(∞),通过AR(∞)解决。

1.3 模型参数识别

通过观察数据相关图中自相关函数(ACF)和偏相关函数(PACF)可以得出数据适用于AR、MA还是ARIMA模型,并通过图形形状判断p、d、q。表1和表2为ACF与PACF的理论模式图形判断的标准。

表1 ACF的典型模式

表2 PACF的典型模式

确定p、d、q,步骤:①通过单位根检验法确定xt～I(d)的d;②确定xt～AR(p)中的p;③确定xt～MA(q)中的q。平稳序列自相关函数

pk=cov(xt,xt+k)var(xt)var(xt+k)=

cov(x0,xk)var(x0)var(x0)=rkr0

(16)

其中,p0=1,p-k=pk(对称)。

1.4 建模步骤

(1) 识别。找出适当的p、d、和q值。通过相关图和偏相关图可以解决。

(2) 估计。估计模型周所含自回归和移动平均项的参数。有时可以用最小二乘法,有时候需要用非线性估计方法。

(3) 诊断(检验)。检验计算出来的残差是不是白噪音,是,则接受拟合;不是,则重新再做。

(4) 预测。短期更为可靠。

2 MEMS陀螺仪数据建模

2.1 数据采集与分析

因为MEMS陀螺角速度为0时的陀螺输出信号能够较好地反映其噪声特性[10],所以本文首先以静态漂移信号为例进行建模处理。由于陀螺噪声信号是连续信号,而时序建模的对象是离散的时间序列,这就需要对连续信号进行采样。以x轴的陀螺为例,采样周期为10 ms。即每10 ms采集1次陀螺输出的瞬时数据。进行20 min试验。陀螺漂移原始测量部分信号如图1所示。

图1 原始数据

由图1可知陀螺仪静态信号由常量的零位偏移与噪声信号组成,样本序列是随机时间序列。从文献[10]对陀螺仪信号噪声成分的分析可以知道,陀螺仪的随机误差主要为角度随机游走和零偏不稳定性;角度随机游走是指由角速率随机白噪声积分引起的具有随机游走特性的误差角增量,角度随机游走具有角速率白噪声功率谱,是输出信号的测量白噪声。所以本文首先对陀螺仪原始数据进行时间序列分析建立ARIMA模型,然后使用卡尔曼滤波算法减小陀螺仪随机误差。

2.2 模型阶数确定

由图1可知,原始数据并未有明显的上升与下降趋势,只是在0值附近不断的偏移,所以可以把陀螺仪静态漂移信号看作平稳信号进行处理,ARIMA模型参数中d=0。

观察数据的相关分析图(见图2)。从图2中可以看到,自相关函数序列ACF呈现正弦波形状,这是AR(2)模型的特征。对时间序列,可以初步判定序列属于二阶自回归模型;此外偏相关函数序列PACF呈现正弦波形状,其自相关函数序列只有一个显著不为0,因此判定时间序列适用于二阶移动平均模型MA(1)。结合AR(2)、MA(1),可以判断此时间序列满足ARIMA(2,0,1)。

图2 相关性分析图

图3 ARIMA模型拟合图

图3为ARIMA模型的拟合效果及残差。计算ARIMA模型与AR模型的拟合残差均值、拟合残差方差见表3,可以看出陀螺仪静态漂移信号的模型中,相对于常用的AR模型,本文选用的模型ARIMA(2,0,1)在拟合效果上有了一定的提高,是更为精确的模型。

表3 模型拟合数据比较

3 卡尔曼滤波方程的建立

卡尔曼滤波是一个高效的递归滤波器,它可以实现从一系列的噪声测量中估计动态系统的状态[11-12]。

由以上得到的ARIMA(2,0,1)模型的参数,即自回归系数,由最小二乘法可以算出模型参数为[0.7822;-0.04502;0.05241],则卡尔曼滤波方程的状态空间模型

Xk=AXk-1+BWk-1

Zk=HXk+Vk

(17)

同时,Wk和Vk满足

系统状态为Xk=[t,t-1]T,过程噪声为Vk=[at,0]T,状态空间模型系数A、B为,

A=0.7822-0.04502

10,B=0,052410T

设系统的输出为Zk=Xk,则输出方程中的系数为H=[1,0]T

k,k-1=Ak-1,k-1
k,k=k,k-1+Kg(Zk-Hk,k-1)
Kg=Pk,k-1HT(HPk,k-1HT+R)-1
Pk,k-1=APk-1,k-1AT+BQBT
Pk,k=[I-KgH]Pk,k-1

(19)

图4为卡尔曼滤波处理后的效果图。由图4可以看到滤波后的数据明显比滤波前的数据平滑了,把漂移误差角度有效的控制在0.01°范围内。表4中,σ2表示数据的方差,这个统计量是描述数据离散状态很重要的指标,较小时,表示数据比较集中;较大时,表示数据的取值比较分散。从表4可以看出陀螺仪随机漂移信号在经过卡尔曼滤波后,σ2明显减小。表4中2个模型滤波前后信号方差的对比,ARIMA模型的滤波效果相比于AR模型提高了42%,这些结果表明,选取以ARIMA误差模型的卡尔曼滤波更能有效地抑制微机械陀螺的随机漂移。

图4 卡尔曼滤波处理效果图

数据方差滤波前AR模型滤波ARIMA模型滤波σ2121E⁃0045281E⁃0053068E⁃005

图5 动态非平稳过程拟合图

为验证ARIMA建模在动态非平稳过程条件下的处理效果,实验室条件下,使用Z1C-DW-26H型电镐分别模拟了钻头在空气、砖块、石灰、水泥等不同地质情况下的钻井工作,采集了不同的信号数据。图5是对不同数据进行ARIMA模型拟合的效果图。从图5可以看出,在不同的动态非平稳过程条件下ARIMA模型的拟合效果都很好,说明模型在动态非平稳过程条件下依旧适合数据的拟合。

表5是在动态非平稳过程条件下,陀螺仪随机漂移信号经过滤波处理后的方差对比。从表5可以看出,动态漂移信号在经过卡尔曼滤波后,方差σ2明显减小,说明ARIMA模型对于动态非平稳信号的处理效果在大多数情况下优于其他模型。表6分别在静态和动态条件下,给出了1500个数据处理时2种模型算法的消耗时间,从表6数据可以看出,ARIMA误差模型的滤波算法效率相比于AR模型也有着一定的优势。

表5 动态非平稳过程滤波方差对比

表6 算法效率对比

4 结 论

(1) 提出了一种基于ARIMA模型的MEMS陀螺仪信号卡尔曼滤波处理算法,从算法思路及算法步骤等方面做了详细讲解,并对实测的陀螺仪的随机噪声信号进行降噪处理,处理结果与常用AR模型进行了对比说明。

(2) 该方法选用的误差模型能够较其他的模型更为精确,可以有效增强滤波效果,大大降低陀螺仪的零位偏移,减小陀螺仪的随机误差,这将对于随钻测量系统姿态估计的精度提高有很大帮助。

参考文献:

[1] 李荣, 赵建辉, 李帆. 轴向速度组合的陀螺随钻测量方法研究 [J]. 测井技术, 2012, 36(2): 192-196.

[2] Erdinc T, Said E A, Tayfun A. Quadrature-error Compensation and Corresponding Effects on the Performance of Fully Decoupled MEMS Gyroscopes [J]. Journal of Micro-electromechanical Systems, 2012, 21(3): 656-667.

[3] 赵世峰, 张海, 范耀祖. MEMS陀螺随机漂移多尺度滤波方法 [J]. 中国惯性技术学报, 2007, 15(2): 229-232.

[4] 席绪奇, 姚志成, 何志昆, 等. 基于多项式模型和BP神经网络的光纤陀螺温度补偿 [J]. 计算机应用与软件, 2013, 30(11): 54-56.

[5] 张智永, 范大鹏, 李凯, 等. 微机电陀螺零点漂移数据滤波方法的研究 [J]. 中国惯性技术学报, 2006, 14(4): 67-69.

[6] 侯青剑, 缪栋, 彭云辉. 激光陀螺随机漂移数据建模与滤波 [J]. 中国惯性技术学报, 2005, 13(4): 39-42.

[7] 周鑫, 张锦, 赵研科, 等. GPRS小区流量预测中时序模型的比较研究 [J]. 计算机应用, 2010, 30(4): 884-887.

[8] 张冉, 赵成龙. ARIMA模型在网络流量预测中的应用研究 [J]. 计算机仿真, 2011, 28(2): 171-174.

[9] 杨叔子, 吴雅, 轩建平, 等. 时间序列分析的工程应用 [M]. 武汉: 华中理工大学出版社, 2007: 175-184.

[10] 金光明, 张国良, 陈林鹏, 等. MEMS陀螺仪静态漂移模型与滤波方法研究 [J]. 传感器与微系统, 2007, 26(11): 48-50.

[11] 赵建辉, 王鹏鹏, 李帆, 等. 陀螺测斜仪随钻对准中的可观测性研究 [J]. 测井技术, 2009, 33(3): 279-283.

[12] 杨金显. 基于航弹的动基座传递对准研究 [J]. 弹箭与制导学报, 2009, 29(6): 76-80.