靳凯迪，柴洪洲，宿楚涵，向民志

信息工程大学地理空间信息学院，河南 郑州 450001

由于元件构造和环境影响，光纤陀螺(fiber optic gyroscope，FOG)会产生复杂的噪声源，直接影响惯性导航系统初始对准和导航解算整个过程的精度。光纤陀螺噪声主要由系统性噪声与随机噪声组成，系统性噪声可通过试验标校有效补偿，随机噪声由于性质复杂、随机性强、补偿难度大，成为衡量光纤陀螺精度的重要指标[1-3]。因此，对光纤陀螺随机噪声进行数学建模并通过滤波等手段加以补偿，对准确分析光纤惯性导航系统误差和提高导航解算精度等方面具有重要意义。目前，用于光纤陀螺随机误差建模的方法主要包括Allan方差分析、功率谱密度分析和时间序列分析等[4-8]。文献[9]使用自回归滑动平均模型(auto-regressive moving average，ARMA)建模光纤陀螺随机噪声，并给出原始观测序列的Kalman滤波方程，但并未考虑MA模型存在时状态噪声为有色噪声的情况。文献[10]使用高阶AR模型计算ARMA模型的噪声估值，达到白化有色噪声的目的，但其作为一种近似方法并不严密。后续研究中或是将随机噪声建模为AR模型[11]，或是沿用文献[9—10]的思想[12-13]。针对随机噪声滤波模型先验信息不准确的问题，文献[11]使用Sage-Husa自适应滤波(Sage-Husa adaptive Kalman filter，SHAKF)在线估计系统噪声和量测噪声。文献[12]使用动态Allan方差分析法(dynamic Allan variance，DAVAR)估计量测噪声。文献[7]针对自适应滤波参数耦合问题，使用DAVAR和SHAKF分别估计状态噪声和量测噪声。然而，SHAKF易丧失矩阵正定性或非正定性[14-16]，尤其在ARIMA模型收敛初期，此现象更加显著。

为此，本文使用并扩展了Harvey法,将光纤陀螺差分自回归滑动平均模型(auto-regressive integrated moving average，ARIMA)构造为原始观测序列的Kalman滤波方程，实现MA参数存在时的有色噪声白化。针对滤波模型及先验噪声不准确的问题，在使用序贯平差更新ARIMA模型参数的同时，计算状态噪声递推解。针对模型更新参数与状态噪声相互耦合，分析了动态Allan方差估计量测噪声的不准确性，使用变分贝叶斯自适应滤波(variational Bayesian adaptive Kalman filter,VBAKF)修正状态噪声和量测噪声，实现Kalman滤波函数模型与随机模型的实时更新。实测数据表明，该方法可有效降低光纤陀螺随机误差，提高陀螺输出精度。

1 光纤陀螺随机噪声建模

ARIMA模型实质上是差分运算与ARMA模型的融合。ARIMA建模光纤陀螺随机噪声序列的主要内容包括噪声随机序列平稳性检验、ARIMA模型阶数确定、模型参数估计及适用性检验等[17-20]。ARMA模型要求平稳的时间序列，因此，首先要对采样光纤陀螺原始随机噪声进行平稳性和随机性检验，如检验不通过，则需进行差分直至满足平稳随机过程。使用中海达iPos光纤惯导以100 Hz采样频率采集静态数据如图1所示。采用单位根检验(包括ADF检验与KPSS检验)对信号检测通过，因此，ARIMA模型差分阶数为0。设ARIMA(p,0,q)模型如下[21]

图1 FOG随机噪声原始输出信号Fig.1 Original random noise signals of FOG

(1)

式中，xk为历元k时随机噪声值；ai为自回归项系数；bj为滑动平均项系数；εk为未知方差的白噪声；c为常值。

使用前2000个历元数据建立ARIMA模型，计算过程分为模型阶数(p,q)确定和模型系数(ai，bj)计算。(p,q)可以通过时间序列的偏自相关函数图确定，但由偏自相关函数图得到的ARIMA模型阶数往往不准确。因此采用赤池信息准则(akaike information criterion，AIC)与贝叶斯信息准则(Bayesian information criterionm，BIC)确定随机噪声序列ARIMA阶数[6]，计算结果为p=3，q=3。模型阶数确定后，通过最小二乘法计算模型系数值，所得ARIMA模型见表1。

表1 光纤陀螺随机噪声ARIMA(p,0,q)模型参数Tab.1 Parameters for ARIMA(p,0,q) of FOG random nosie

图2为对所建模型采用标准化残差检验结果。由图2可知，所得残差序列符合正态分布。Durbin-Watson统计量是常用的自相关统计量，设残差序列为vk(k=1,2,…,n)，构造DW统计量为[22]

图2 标准化残差检验Fig.2 Standardized residuals test

(2)

式中，ρ为自相关因子。ρ越接近0，等价于DW统计量越接近2，序列越不存在一阶相关性，本次残差序列相关性检验所得结果为1.953 9，认为残差序列符合高斯正态分布。综上，认为所建模型较为准确，通过适用性检验。

2 ARIMA模型有色噪声Kalman滤波

在光纤陀螺随机噪声建模为ARIMA模型后，使用Kalman滤波是进行随机噪声降噪的有效方法。以往处理存在MA系数的ARIMA模型时，往往忽视噪声的相关性或使用噪声估值进行白化。为此，建立基于Harvey法的Kalman滤波模型白化ARIMA有色噪声。

2.1 ARIMA模型滤波方程

获得如式(1)所示的随机噪声ARIMA(p,0,q)模型后，令状态向量Xk=[xkxk-1…xk-p+1c]T，建立Kalman滤波状态方程和量测方程为[9]

(3)

2.2 Harvey法有色噪声滤波方程

当存在MA系数时，式(3)中的状态噪声为有色噪声，此时，Kalman滤波将会导致结果发散或失真[14，21]。传统上，通常将有色噪声建模为一阶AR模型[23-24]，但ARIMA的有色噪声数学模型是明确的，对噪声再进行时间序列建模并不完全准确。在标准Harvey法滤波方程中[25]，针对如下所示ARIMA模型

(4)

令y1,k=xk，有

a1y1,k-1+y2,k-1+εk

(5)

考虑y2,k

a2y1,k-1+y3,k-1+b1εk

(6)

递推可得如下关系

(7)

式中，bi=0(i=q+1,q+2,…,p-1)；m=1,2，…,p-1。构建状态方程为

(8)

式中，F为状态转移矩阵；G为噪声驱动矩阵。可见式(8)仅适用于p>q时，且建模时，没有考虑常数项c。对于本文所述的q=p，且含常数项的光纤陀螺随机噪声ARIMA模型，需要对Harvey方程进行扩展。考虑常数项c，对于式(7)递推至m=p-1时，有

(9)

继续进行递推

(10)

(11)

式中

(12)

式中，εk和Vk均为白噪声且两者不相关，满足随机模型

(13)

扩展后的Harvey方程多出了一维由“常值+白噪声”构成的随机常值状态参数。实现了附带常数项，且q=p的ARIMA模型状态方程的构建。通过Harvey法构造的状态噪声为高斯白噪声，满足标准Kalman滤波条件。

3 自适应Kalman滤波

3.1 ARIMA模型及状态噪声在线更新

光纤陀螺随机噪声复杂且时变，实际应用中使用有限数量或不同时期的数据所得到的ARIMA模型随时间的推移，误差会逐渐增大。因此，有必要依据当前观测数据进行ARIMA参数在线更新。将式(1)改写为

(14)

式中

(15)

使用序贯平差实现模型参数在线更新的主要流程如下

(16)

(17)

3.2 变分贝叶斯自适应滤波

基于Allan方差的滤波器是一种带通滤波器，可直接滤除部分低频系统噪声的特性，传统上，将Allan方差值近似为宽频量测噪声的方差[1]，选择取样间隔为最短采样时间的Allan方差为

(18)

对于上文所述的量测模型，量测值Zk：N(xk,R)，其中N(μ,∑)表示均值为μ和协方差阵为∑的高斯概率密度函数，由式(18)可得

(19)

(20)

可见，当且仅当序列真值X=[X1X2…Xk]为均值时，即量测序列为均值+白噪声时，偏差项等于0，最小时间间隔的Allan方差才与真实噪声方差相等，而对于时变序列如本文所建模的ARIMA模型序列，将会存在上述偏差项。

VBAKF是一种近似的贝叶斯法，它利用已知的模型信息、观测信息和先验信息来获得状态向量和未知参数联合后验密度的近似解。可同时修正Qk与Rk的VBAKF的主要流程[15]如下。

(1) 状态更新过程

(21)

(22)

式中，IW分布的两参数分别为自由度参数和逆尺度矩阵。

(2) 变分量测迭代更新过程。

返回②继续迭代。

4 试验与分析

对图1所示的光纤陀螺随机噪声进行试验分析，原始陀螺随机噪声Allan方差曲线如图3所示，可见该陀螺随机漂移主要由角度随机游走、零偏不稳定性组成。图4为通过序贯平差建立的ARIMA系数在线更新，在试验条件下，光纤陀螺随机噪声序列的ARIMA模型较为稳定，估值很快收敛，因此，在实际测量中，可仅于初始阶段进行模型参数更新，以提升计算效率。

图3 光纤陀螺随机噪声Allan方差Fig.3 Allan-variance of FOG random noise

图4 光纤陀螺随机噪声ARIMA系数更新Fig.4 ARIMA coefficients update for FOG random noise

将动态Allan方差应用于量测噪声自适应估计是一种近似方法。为对比DAVAR与VBAKF对量测噪声的估计能力，分别在原始训练数据中，增加均方差为1×10-3、1.5×10-3、2×10-3(°)/s的高斯白噪声进行解算，量测噪声估计结果如图5所示，应认为原始训练数据与所建ARIMA模型相符性较好，因此，量测噪声与添加噪声大小近似相同，在图5中第1000历元，VBAKF量测噪声估值分别收敛至1.05×10-3、1.44×10-3、1.95×10-3(°)/s，DAVAR量测噪声估值收敛至1.46×10-3、1.84×10-3、2.36×10-3(°)/s。因此，可认为VBAKF较DAVAR对量测噪声有更好的估计性能。使用Harvey法构造Kalman滤波方程，结合VBAKF(以Harvey VBAKF代替)对光纤陀螺随机噪声进行实时滤波处理。同时，分别使用不考虑有色噪声的滤波方程(以KF代替)及不使用VBAKF的Harvey法滤波方程(以Harvey代替)对同一数据进行处理，原始数据与滤波结果如图6所示，数据统计特性见表2。

图5 DAVAR与VBAKF噪声自适应对比Fig.5 Adaptive comparison of DAVAR and VBAKF measurement noise

表2 滤波结果统计特性对比Tab.2 Comparison of statistical results for filtering data (°)/s

图6 滤波输出数据比较Fig.6 Comparison of filtering output data

比较图1、图6可知，总体上，光纤陀螺随机噪声经Kalman滤波后明显减小。3种滤波策略结果的均值与原始数据相同(表2)，文献[5]使用滑动平均消除常值结果与原始数据均值相差在10-5量级，因此认为本文滤波结果的均值没有变化。均方差上，相较原始数据，传统Kalman滤波、Harvey法滤波、Harvey法结合VBAKF滤波分别降低了29%、41%、54%，随机噪声显著减小，可见在保证无偏估计的同时，Harvey法建立的Kalman滤波方程有效消除了ARIMA模型滤波中有色噪声影响，VBAKF实现了对ARIMA模型与状态噪声的解耦，可以有效估计状态噪声与量测噪声，它对光纤陀螺随机噪声抑制效果更加明显。

5 结 论

使用Harvey法构造的光纤陀螺ARIMA滤波模型可以有效减弱有色噪声影响，VBAKF具有良好的收敛性和准确的噪声估计精度。本文将Harvey法构造ARIMA滤波方差引入光纤陀螺随机误差降噪，并针对其中存在的状态噪声耦合及先验信息不准确问题，使用变分贝叶斯自适应卡尔曼滤波实现对状态噪声和量测噪声的同时修正，并且验证了VBAKF对量测噪声有比动态Allan方差更准确的估计精度。该方法在光纤陀螺随机误差建模与补偿中可提供一定的参考。