光纤陀螺随机漂移的补偿方法研究

汪徐胜,张瑞民,信 东,李 冰

摘 要:在研究高精度光纤陀螺时,尤其是在捷联惯性导航系统中,随机误差是光纤陀螺误差中不可忽视的部分,对光纤陀螺随机误差的补偿就显得非常必要。这里基于对光纤陀螺随机漂移建模的方法,首先采用ARIMA方法对光纤陀螺仪随机漂移进行建模;然后采用强跟踪卡尔曼滤波器进行滤波补偿,并利用实测数据进行了实验验证。实验结果证明,这种方法能够较好地补偿光纤陀螺的随机漂移。

关键词:光纤陀螺;捷联惯性导航系统;随机误差;卡尔曼滤波器;自适应滤波器

中图分类号:TP274文献标识码:A

文章编号:1004-373X(2009)12-095-04

Research of Compensation Method for Stochastic Shifting of Fiber Optic Gyros

WANG Xusheng1,ZHANG Ruimin1,XIN Dong1,3,LI Bing2

(1.The Second Artillery Officer School,Qingzhou,262500,China;

2.The Second Artillery Officer Stationed in the Forth Space Flight Academe,Xi′an,710025,China;

3.The Second Artillery Engineering College,Xi′an,710025,China)

Abstract:When doing research in optical fiber gyros,especially in strapdown inertial navigation system,stochastic error is an assignable part of the errors of optical fiber gyros.It is very necessary to compensate stochastic error of optical fiber gyros.Basing on the first method,the article builds a model for the stochastic shifting,then adopts strong tracing Kalman filter to filter and compensate.At last,the methods are testified by experimental data.The results prove that the method is effective.

Keywords:optical fiber gyros;strapdown inertial navigation system;stochastic error;Kalman filter;self-adapting filter

0 引 言

在惯性技术中,通常将惯性敏感元件的误差模型分为静态误差模型、动态误差模型和随机误差模型三类[1,2]。静态误差模型和动态误差模型可以比较容易地用代数方程来表示,因此比较容易补偿。随机误差是指误差中随机变化的部分,要用统计规律描述,误差要通过滤波方法来减小。研究低精度光纤陀螺时,由于陀螺的随机漂移较大,严重影响了陀螺的测试精度,必须进行随机误差补偿。光纤陀螺应用于惯性导航中,由于工作时间比较长,积分结果将会对导弹造成不可容忍的误差。因此,对光纤陀螺的漂移进行补偿具有重要的意义。

人们对这一部分的研究也相对比较成熟[3-5],但在建立随机误差模型时,大都局限于采用简单的自回归模型,影响了误差模型的精度;同时,滤波大都采用普通kalman滤波方法,由于这种滤波方法的局限性(噪声模型必须准确已知),容易导致滤波的不稳定甚至发散,本文采用能够跟踪信号变化的强跟踪kalman滤波方法[6],有效地克服了这一问题。

1 时间序列建模方法研究

1.1 建模前的数据检验

一般,光纤陀螺的漂移数据可以近似看作时间序列来处理,而常用的时间序列建模方法[7]有自回归模型法(AR)、滑动模型法(MA)、自回归滑动模型法(ARMA)。通常ARMA模型更能够准确地反映数据的随机特性。无论在用那种方法进行建模,首先时间序列必须是平稳的。但实际情况是,由于受各种干扰的影响,所测得的时间序列一般不是平稳的,这就需要进行数据的检验和平稳化处理。主要包括平稳性检验、正态性检验以及零均值检验[8,9],在建模前,要对数据进行这3种检验。由于零均值检验和正态性检验在参考文献中有详细的介绍,文中仅对平稳性检验进行介绍。

时间序列的平稳性是建模的重要前提,因此需要对测试数据进行平稳性检验。这里采用一种工程中常用的逆序检验法,其具体过程为:

设样本序列x1,x2,…,xN足够长,把样本序列分成k个子序列,即N=lM,M是一个较大的正整数;l也为正整数。

逆序检验法就是检验各子列均值的差异性,把各子列均值μj构成一个序列μ1,μ2,…,μl。当i>j时,每出现一次μi>μj,定义Aj为μj的一个逆序数。对于某一μj,定义μj的逆序数Aj为μi>μj(i>j)出现的次数,所有逆序数Aj的总和为序列的逆序总数,其表达式为:

A=∑lj=1Aj

由于{xk}是随机过程的一个样本,所以均值μj也应该是随机的。当{xk}是平稳序列时,μ1后面的l-1个随机数μj大于和小于μ1的概率是相等的,故μ1的逆序数A1的理论平均值E[A1]=(l-1)/2;同理,E[A2]=(l-2)/2,…,E[Ak]=(l-k)/2(k=1,2,…,l-1),从而逆序数的理论平均值:

E[A]=∑l-1j=1EAj=l(l-1)/4

逆序总数的理论方差为:

σ2A=l(2l2+3l-5)/72

构造统计量:

u=(A+0.5-E[A])/σA

这样构造的统计量满足标准正态分布,因此当显著性水平为0.05时,如果|u|≤1.96,则可确定μj间无显著差异,即序列是平稳的。

1.2 ARIMA时间序列建模方法

在许多实际问题中,所考虑的时间序列{xk}并不是平稳的,它不能用ARMA模型来表示。但是,如果将{xk} 进行有限次差分后,得到的时间序列{x′k} 是平稳的,{x′k}可以用ARMA模型来描述。原时间序列{xk}可以用ARIMA(n,d,m)(基于序列平稳化的ARMA建模)表示。平稳序列的ARMA(n,m)模型表达式为:

Φ(B)xk=Θ(B)ak(1)

式中:Φ(B)=1-∑ni=1ΦiBi;Θ(B)=1-∑mj=1θjBj;n是AR部分的阶次;m是MA部分的阶次;B是后移算子,Bxk=xk-1;ak是残差,满足ak～(0,σ2a)。

如果序列不是平稳的,则要对序列进行d阶差分使序列平稳化。差分是将时间序列逐项递减而得到一个新的序列,其目的是为了消除原始序列的趋势项,将非平稳序列化为平稳序列。然后建立差分后序列的ARMA(n,m)模型,原序列则可以表示为ARIMA(n,d,m)模型:

Φ(B)(1-B)dxk=Θ(B)ak(2)

通过以上分析,建立ARIMA(n,d,m)是建立在ARMA(n,m)模型基础之上的。在建模前要进行数据平稳化的检验,如果随机序列是平稳的,取d=0。ARIMA(n,d,m)模型简化为ARMA(n,m)。

1.3 建立ARMA模型

1.3.1 拟合模型

按照ARMA模型参数估计中采用的方法和理论,模型参数估计方法分为3种[2]。由于时序理论估计法算比较简单,而且可以保证较高的精度,所以得到比较广泛的应用。在此采用的也是这种参数估计方法。

ARMA(n,m)模型[7]可以表示为:

xk=∑ni=1Φixk-1-∑mj=1θjak-j+ak(3)

对自回归模型参数Φi和滑动回归模型参数θj可以采用长自回归模型法。这种方法的思路是:基于观测序列建立起来的模型AR(n),MA(m),ARMA(n,m)都是等价系统的数学模型,具有相同的传递函数。基于这个原则,可以先估计出AR模型,然后利用传递函数的等价性,估计出ARMA的参数自回归模型参数Φi和滑动回归模型参数θj。下面就长自回归模型法[3,5]做简单介绍。

随机序列{xk}可以用模型AR(p)表示,也可以用ARMA(n,m)表示,AR(p)描述的等价系统传递函数为:

1/Φp(B)=1/(1-∑pi=1IiBi)(4)

式中:B为后移算子,Ii是逆函数,其取值等于AR(p)模型的参数Φi(i=1,2,…,p)。用ARMA(n,m)描述的等价系统传递函数为:

θm(B)/Φn(B)=(1-∑mj=1θjBj)/(1-∑ni=1ΦiBi)(5)

由于不同形式的传递函数描述的系统是等价的,所以上述两式应该相等。把两式展开,得到等式:

(1-I1B-I2B2-…-IpBp)(1-θ1B-θ2B2-

…-θmBm)=1-Φ1B-Φ2B2-…-ΦnBn

所以,对应B的相同阶次的系数应该相等。由此得到:

Φ1=θ1+I1

Φ2=θ2-θ1I1+I2

Φ3=θ3-θ2I1-θ1I2+I3

…

Φn=-θmIn-m-θm-1In-m-1-…-θ1In-1+In

0=-θmIk-m-θm-1Ik-m-1-…-θ1Ik-1+Ik,

k>n(6)

用最后一个方程,当k取n+1,n+2,…,n+m的时候,可以用以下矩阵求出θj(j=1,2,…,m):

In+1In+2驣n+m=InIn-1In-2…In-m+1In+1InIn-1…In-m-2

In+m-1In+m-2In+m-3…Inθ1θ2螃萴(7)

在求出θj(j=1,2,…,m)以后,利用式(1)可以求出Φi(i=1,2,…,n)。关于Φi(i=1,2,…,n)的线性方程为:

Φ1Φ2螃祅=θ1θ2螃萵+100…0-θ110…0

-θn-θn-1-θn-2…0I1I2驣n(8)

求解上述线性方程,便得到Φi(i=1,2,…,n)的值,这就是长自回归模型法的计算原理。长自回归模型法在估计Φi和θj的过程中都是解线性方程组。它将ARMA(n,m)模型参数估计的非线性问题转化为线性问题来解决,大大简化了问题。

2 干涉式光纤陀螺仪随机漂移的建模

这里以一组VG951型光纤陀螺仪实测数据为例,建立光纤陀螺仪的随机漂移模型,这组数据采样间隔为1 s,总体样本长度为3 000,数据中已减去了地球常值的影响。图1为实测数据曲线。

图1 光纤陀螺仪的实测数据

由图1可以看出序列均值不为零,且有下降的趋势。用逆序法检验平稳性,取l=6,得到u=-2.630 1,因为|u|>1.96,所以序列不是平稳的。其平均值为3.846 2,漂移值为2.081 6,因此建模前要对序列进行预处理。

首先对原序列进行一阶差分,得到一新的序列,检验新序列的统计特性,得到如下数据:

平稳性检验:u=0.751 5;

正态性检验:ξ=0.014 3;γ=3.004 1;

经过检验,处理后的数据基本上符合平稳性和正态性的检验,可以对其建立ARMA(n,m)模型。

分析处理后的序列,其可以用ARMA(n,m)模型表示。建模采用BIC准则(见参考文献[5])确定模型的阶次,用长自回归模型法辨识模型的系数。经计算ARMA(3,1)模型为适用模型,建模过程为:

首先确定AR(3)模型的系数,模型可以表示为(其中v(k)为白噪声):

y(k)=a1y(k-1)+a2y(k-2)+a3y(k-3)+

(1-a1-a2-a3)+v(k)

求解关于模型参数的矩阵方程,得系数值:

123=-0.760 4-0.521 8-0.257 5

从而得到ARMA(2,1)的系数:

θ1=0.493 5;Φ1=-0.267 0;Φ2=-0.146 6

代入差分方程,得到ARIMA模型的系数:

Φ=0.7330.120 40.146 6;θ=[0.499 2]

从而,原随机序列模型的表达式为:

xk=0.733xk-1+0.120 4xk-2+0.146 6xk-3-

ak+0.493 5ak-1(9)

3 滤波算法

强跟踪卡尔曼滤波是一种适合计算机计算的递推滤波方法,它成功地采用了状态空间概念,把信号过程视为白噪声作用下一个线性系统的输出。这种方法不需要存储过去的观测数据,当新的数据被观测后,只要根据新的数据与前一个时刻的估计量,再借助于信号过程本身的状态转移方程,按照递推公式,即可计算出该时刻的估计量。

光纤陀螺输入/输出数据的卡尔曼滤波器的设计流程,如图2所示。

图2 卡尔曼滤波的总体流程

离散系统的卡尔曼最优滤波的递推方程为:

(k|k)=(k|k-1)+

K(k)[Z(k)-H(k)(k|k-1)]

(k|k-1)=Φ(k,k-1)(k-1|k-1)

K(k)=P(k|k-1)HT(k)•

[H(k)P(k|k-1)HT(k)+Rk]-1

P(k|k-1)=Φ(k,k-1)P(k-1|k-1)ΦT•

(k,k-1)+Γ(k,k-1)Qk-1ΓT(k,k-1)

P(k|k)=[I-K(k)H(k)]P(k|k-1)•

[I-K(k)H(k)]T+K(k)RkKT(k)

对第2部分随机序列建模的结果,可以直接运用卡尔曼滤波方法减小随机漂移,但是这种方法不能跟踪变化的信号,在仿真中给出的是一种能够跟踪变化信号的滤波方法。

4 仿真验证

此前建立的ARIMA模型可作为卡尔曼滤波器的状态方程,建立观测方程得到随机漂移的观测数据。目前,得到的数据只有光纤陀螺输出值。考虑到光纤陀螺的输出由常值漂移、角速度敏感量、随机漂移和噪声组成,可以假设相邻两时刻的角加速度变化率相等。光纤陀螺输入/输出的关系式可以表示为:

w=(1/ks)(y-y0-x-ε)

对上式进行二阶微分可以得到:

d2wdt2=1ks(-d2ydt2+d2xdt2+d2εdt2)=0

差分后得到结果:

yk-2yk-1+yk-2=xk-2xk-1+xk-2+

εk-2εk-1+εk-2

设:

νk=εk-2εk-1+εk-2

令Zk=yk-2yk-1+yk-2作为随机噪声的观测量,可以写出观测方程:

Zk=xk-2xk-1+xk-2+νk

利用上式作为观测方程,与前面所得到的模型方程一起构成典型卡尔曼滤波器。

令Xk=[xk,xk-1,xk-2],将所建随机模型(1)和观测方程写成矩阵形式为:

Xk=ΦXk-1+GWkZk=HXk+Vk(10)

式中,Φ系统状态转移矩阵,其值为:

Φ=0.7330.120 40.146 6100010

G是系统噪声阵,可表示为:

G=1-0.493 50000000

H是系统观测阵,表示为:H=[1-21]。

式中:Wk是系统噪声;Vk是观测噪声,两者是互不相关的零均值高斯白噪声,在滤波时取为标量。图3为补偿后的光纤陀螺随机漂移数据。

图3 补偿后的随机漂移

经补偿,序列的随机漂移为0.745 5 °/h。可以看出,依据漂移的模型进行滤波能够有效地抑制随机漂移。这里的关键问题是得到光纤陀螺的随机漂移模型。但在一些实际问题中,光纤陀螺的漂移模型是时变的,事先建立的漂移模型不一定能够准确反映当前的陀螺漂移。因此,在实时补偿时,会存在一定的危险。为了达到较好的补偿效果,需要在线建立光纤陀螺的随机漂移模型,这样才能达到较好的动态补偿效果。分析影响光纤陀螺随机漂移的各种因素,采用自适应建模的方法,建立光纤陀螺的自适应补偿模型,根据检测的环境变量和采样数据的统计特性自动修改模型参数,能够有效地克服光纤陀螺漂移模型的时变性,可以达到较理想的效果,其中采用横向自适应滤波器是一种较好的解决方案。

5 结 语

静止时光纤陀螺输出的是一随机序列,在分析了光纤陀螺随机漂移特性以及时间序列的随机建模的基础上,利用实测数据建立了光纤陀螺的ARIMA漂移模型。然后利用强跟踪卡尔曼滤波进行随机漂移的最优估计,用以对输出进行补偿。对检测数据的处理结果证明,依据ARIMA模型对光纤陀螺的随机漂移进行补偿具有较好的效果,大大降低了光纤陀螺的随机漂移。采用的ARIMA建模方法以及强跟踪Kalman滤波方法,不管是在精度上,还是在算法的简单性和可行性上,都比以前的随机建模及滤波方法有一定的提高,这也是本文的创新点所在。

参考文献

[1]王珍熙.捷联式惯性导航系统惯性元件的设置与可靠性[J].中国惯性技术学报,1996,4(1):61-65.

[2]郭秀中.惯导系统陀螺仪理论[M].北京:国防工业出版社,1996.

[3]王新龙.光纤陀螺随机误差模型分析[J].北京航空航天大学学报,2006,32(7):769-772.

[4]吉训生,王寿荣.MEMS陀螺仪随机漂移误差研究[J].宇航学报,2006,27(4):74-76.

[5]李爱华.闭环光纤陀螺的建模方法[J].电光与控制,2006,13(4):43-45.

[6]刘铭,周东华.残差归一化的强跟踪滤波器及其应用[J].中国电机工程学报,2005,25(2):74-78.

[7]王志贤.最优估计与系统辨识[M].西安:西北工业大学出版社,2003.

[8]秦永元.卡尔曼滤波与组合导航原理[M].西安:西北工业大学出版社,2004.

[9]陆光华.随机信号处理[M].西安:西安电子科技大学出版社,2003.