航空矢量重力测量中光纤陀螺随机漂移误差实时补偿方法

王 峥，李建成,2

1. 武汉大学测绘学院，湖北 武汉 430079； 2. 武汉大学地球空间环境与大地测量教育部重点实验室，湖北 武汉 430079

航空矢量重力测量中光纤陀螺随机漂移误差实时补偿方法

王 峥1，李建成1,2

1. 武汉大学测绘学院，湖北 武汉 430079； 2. 武汉大学地球空间环境与大地测量教育部重点实验室，湖北 武汉 430079

光纤陀螺随机漂移误差是影响航空矢量重力测量系统姿态解算精度的关键因素。建立模型并在输出中对其补偿是抑制该项误差的有效方法。针对传统ARMA模型只能对平稳随机漂移误差建模，且模型无法满足实时滤波需求的问题，本文引入适用于非平稳随机漂移误差的ARIMA模型，同时给出详细的建模过程，并提出采用实时平均算法消除原始采样序列中常值分量的思路，实现了随机漂移误差的实时Kalman滤波估计。基于本文所提出的模型和实时滤波算法，对光纤陀螺实测数据进行分析，结果表明处理后信号中随机漂移误差的方差减小了46.5%。Allan方差分析结果表明，滤波后角度随机游走系数和角速率随机游走系数分别降低了约50%和40%。本文的结果说明ARIMA模型能够准确描述陀螺的非平稳随机漂移误差。基于实时平均算法的Kalman滤波可实现随机漂移误差的在线估计，有望提高航空矢量重力测量系统的姿态解算精度。

光纤陀螺；ARIMA模型；Kalman滤波；Allan方差

航空重力测量是以飞机为载体快速地获取分布均匀、精度良好、大面积区域或局部重力场信息的有效方法[1]。目前，航空标量重力测量可提供精度为2 mGal(1 Gal=10-2m/s2)，分辨率为2 km的重力场信息[2-3]。航空矢量重力测量技术与航空标量重力测量技术相比，不仅能获取重力扰动矢量的垂直分量(重力异常)，而且也能获取重力扰动矢量的水平分量(垂线偏差)。由于受重力加速度的放大作用，重力扰动矢量的水平分量受姿态误差影响较大，1″的水平姿态误差将会引起约4.75 mGal的水平分量误差[1]。文献[4—5]通过大量试验发现陀螺漂移误差是影响水平姿态误差的主要因素，如果能够在航空矢量重力测量解算过程中估计出陀螺漂移误差并对其补偿，可显著提高重力扰动矢量水平分量估算精度。目前航空矢量重力测量主要是捷联式航空重力测量，葡萄牙波尔图大学研究表明光纤陀螺可用于航空矢量重力测量[1,6-7]。光纤陀螺漂移主要包含固定漂移和随机漂移[8-9]。固定漂移不随时间变化，可通过实验室标定加以补偿。随机漂移是一个十分复杂的随机、时变过程，不易消除，是影响航空矢量重力测量水平分量精度的关键因素。因此，针对光纤陀螺随机漂移开展研究对于提高航空矢量重力测量系统精度具有重要意义。

国内外研究人员通过各种方法对光纤陀螺随机漂移误差进行了建模和滤波分析。文献[10—13]中用自回归(autoregressive，AR)模型或自回归滑动平均(autoregressive moving average，ARMA)模型建立了光纤陀螺(fiber optic gyroscope，FOG)的随机误差模型，并用卡尔曼滤波器对随机漂移误差进行了滤波，得到了较好的结果。但是，由于ARMA模型要求信号必须为平稳、正态和零均值的时间序列，因此必须首先对陀螺信号进行平稳、正态和零均值处理，才能建立该模型。所以这种建模和滤波方法只适用于对陀螺信号进行事后处理与分析，而不满足捷联式航空矢量重力测量系统初始对准过程的实时性需求。为解决此问题，文献[14—15]通过将陀螺原始采样值中的常值分量作为新的状态变量引入Kalman滤波系统状态方程中，实现了随机漂移的在线估计。但该方法的应用前提是陀螺随机漂移序列满足平稳性假设条件。在实际应用中，陀螺随机漂移误差可能是非平稳的，由该方法形成的Kalman滤波可观性矩阵将具有奇异性，系统不完全可观，无法实现陀螺原始采样值序列中常值分量的分离。因此，该方法同样不适用于陀螺非平稳随机漂移误差的实时在线滤波。

鉴于此，本文基于某型光纤陀螺仪实测数据，详细探讨了光纤陀螺非平稳随机漂移误差的建模过程，引入求和自回归滑动平均模型(autoregressive integrated moving average，ARIMA)。该模型可直接对非平稳随机漂移误差进行建模。在此模型的基础上，本文通过实时平均算法补偿掉陀螺原始采样值序列中由地球自转和零漂引起的常值分量，实现了陀螺随机漂移误差的实时Kalman滤波估计，并对滤波后的数据进行了Allan方差分析。试验结果表明，基于该ARIMA模型的实时Kalman滤波方法能够有效抑制陀螺仪随机漂移误差，将有效提高航空矢量重力测量系统的姿态解算精度。

1 ARIMA原理及建模

1.1 ARIMA模型结构

(1)

式中，B为延迟算子；Φ(B)=1-φ1B-φ2B-…-φpBp，为平稳可逆ARMA(p,q)模型的自回归系数多项式；Θ(B)=1-θ1B-θ2B-…-θqBq，为平稳可逆ARMA(p,q)模型的移动平滑系数多项式。

根据差分运算性质[16-17]

(2)

将式(1)和式(2)联合，可求出非平稳序列xt表达式为

(3)

由式(3)可知，ARIMA模型的实质就是差分运算与ARMA模型的组合。

1.2 ARIMA模型建模步骤

(1) 时间序列预处理(序列平稳性、正态性检验，趋势项提取等)；

(2) 根据时间序列模型识别规则(表1)建立相应模型；

(3) 模型参数估计；

(4) 模型适用性检验，主要诊断残差序列是否为白噪声。

表1 ARMA模型的识别

2 基于ARIMA模型的光纤陀螺非平稳随机漂移误差建模

2.1 数据预处理

数据预处理主要包括原始信号采集、随机漂移信号提取和序列平稳性分析等步骤。

2.1.1 光纤陀螺随机漂移信号的采集

本文以高精度组合导航系统SPAN-FSAS系统中的IMU-FSAS为研究对象，以fs=200 Hz的采样率，采集了3 h的原始陀螺漂移数据，结果如图1所示。

图1 某光纤陀螺仪原始漂移信号Fig.1 Original drift signals for a certain kind of FOG

2.1.2 随机漂移信号提取

光纤陀螺的原始信号中包含常值分量和随机分量。图2为去掉常值分量之后的随机漂移信号。

2.1.3 序列平稳性检验

图2 光纤陀螺仪随机漂移信号Fig.2 Random drift for a certain kind of FOG

表2 陀螺随机漂移数据平稳性检验

图3 一阶差分陀螺随机漂移Fig.3 First-order differenced random drift of FOG

2.2 模型参数识别及优化

2.2.1 ARIMA(p,d,q)模型识别

模型识别的方法很多，使用较广泛的是Box-Jenkins方法。此方法根据时间序列的自相关系数(autocorrelation function，ACF)和偏自相关系数(partial autocorrelation function，PACF)的截尾、拖尾性来识别模型[17-18]。一阶差分后陀螺随机漂移信号的自相关系数和偏自相关系数计算结果如图4所示。

图4 一阶差分后序列的自相关系数与偏自相关系数Fig.4 ACF and PACF of first-order differenced random drift of FOG

由图4可知，自相关系数在滞后阶数为1时显著不为零；在滞后阶数大于1时基本处于±2σ置信带内。偏自相关系数随延迟阶数k的增大逐渐衰减到零，表现出显著的拖尾性。根据表1中模型的定阶原则，模型应该在滑动平均MA和自回归滑动平均ARMA中选取。

2.2.2 ARIMA(p,d,q)模型优化

同一时间序列可由AR、MA和ARMA模型中的一种或多种进行拟合。最优拟合模型可通过BIC信息准则确定[19-20]。利用BIC准则对一阶差分后的陀螺漂移数据进行定阶，各拟合模型的BIC值如表3所示。

由表3中各模型的BIC值可知，一阶差分后陀螺随机漂移序列的最优模型为ARMA(1,2)。但在实际情况中ARMA模型相当于一个最小实现的线性系统，模型中的自回归阶数应大于或等于滑动平均阶数(p≥q)， 故表3中满足的模型只

有ARMA(1,1)、ARMA(2,1)两种。比较二者的BIC值，选取值较小的ARMA(1,1)模型作为一阶差分后陀螺随机漂移误差模型。

表3 光纤陀螺信号各模型参数及BIC值

Tab.3 Parameters and BIC values for various models of FOG random drift

模型MA(1)MA(2)AR(1，1)ARMA(1，2)ARMA(2，1)θ10．983650．954020．985010．004690．98511θ2—0．03048—0．96267—ϕ1——0．03112-0．984760．03114ϕ2————0．00276BIC值-7．60454-7．60454-7．60225-7．62278-7．59974

2.3 模型参数估计及适用性检验

2.3.1 ARIMA模型参数估计

本文采用最小二乘估计法对差分后平稳模型ARMA(1,1)参数进行估计[21]，各模型对应的自回归系数φi(i=1,2)和滑动平均系数θi(i=1,2)的估值如表3所示。将ARMA(1,1)模型所对应系数代入式(1)得一阶差分后陀螺随机漂移数据模型为

(4)

结合式(3)和式(4)，求出光纤陀螺随机漂移数据的ARIMA(1,1,1)模型为

ωk=1.031 12ωk-1-0.031 12ωk-2+εk- 0.985 01εk-1ε～N(0,0.000 44)

(5)

2.3.2 模型的适用性检验

本文利用QLB统计量对ARIMA(1,1,1)模型的残差序列进行白噪声检验[17]，结果如表4所示。由表4可知统计量QLB的P值显著大于显著性水平α，故残差序列为白噪声，拟合模型ARIMA(1,1,1)显著有效。

表4 ARIMA(1,1,1)模型适用性检验(α=0.05)

3 基于ARIMA模型的陀螺随机漂移误差实时Kalman滤波

Kalman滤波可实现线性最小均方误差估计意义下随机信号的最优线性滤波[22]。在建立了陀螺随机漂移ARIMA(1,1,1)模型后，可采用Kalman滤波对陀螺随机漂移误差进行补偿。根据式(5)建立系统的状态方程

Xk=Φk,k-1Xk-1+Γk,k-1Wk-1

(6)

由于式(6)是针对零均值非平稳陀螺随机漂移误差的系统状态方程，要实现陀螺随机漂移误差的在线Kalman滤波估计，需实时补偿掉陀螺原始采样值序列中由地球自转和常值零漂引起的常值部分。本文根据实时平均算法来分离陀螺原始采样值中的常值分量，采用当前时刻及其前9个时刻采样点的平均值作为当前时刻的常值估计值。

(7)

(8)

根据系统状态方程(6)和量测方程(8)，可得系统的可观测性矩阵为[23-24]

(9)

由可观性矩阵O的秩rank(O)=2，可知该随机系统是完全可观的，由Kalman滤波能对状态实现最优估计。将实测光纤陀螺原始采样数据作为卡尔曼滤波器输入对陀螺随机漂移误差进行实时在线Kalman滤波估计，滤波处理结果如图5所示。

由表5可知，滤波前后数据均值基本一致，其大小约为0.003 54 °/s，因为Kalman滤波属于线性无偏最优滤波。滤波后数据方差较滤波前减小46.5%，说明基于ARIMA模型的实时卡尔曼滤波算法能实现对非平稳随机漂移误差的在线补偿。

表5 陀螺随机漂移数据滤波前后统计量比较

Tab.5 Comparison of the statistical results before and after the filtering for FOG random drift

统计量滤波前滤波后均值/(°/s)0．00353860．0035412方差/((°/s)2)4．3738e⁃0042．3400e⁃004

4 Allan方差分析

由图7可知，光纤陀螺随机漂移误差通过实时Kalman滤波补偿后的Allan方差曲线较原始Allan方差曲线幅值有了明显下降，特别是在高频段下降了一个量级。由图7中τ-σA(τ)双对数曲线斜率特性可知，该型号光纤陀螺仪的随机漂移误差主要由角度随机游走和角速率随机游走组成，中间有稍弱的零偏稳定性表现出来，具体误差系数分析结果如表6所示。

表6 滤波前后噪声源误差系数比较

Tab.6 Comparison of error coefficients of noise before and after the Kalman filtering

参数误差系数滤波前滤波后角度随机游走/(°/h1/2)N0．10340．0507零偏不稳定性/(°/h)B0．07970．0531角速率随机游走/(°/h3/2)K0．47100．2779

由表6可知，经Kalman滤波补偿后，光纤陀螺随机漂移误差中的角度随机游走系数N、零偏不稳定性系数B和角速率随机游走系数K都较滤波前显著减小，尤其是角度随机游走系数N和角速率随机游走系数K分别降低了约50%和40%。

图5 基于ARIMA模型的Kalman滤波结果Fig.5 Results of Kalman filtering based on the ARIMA model

图的幅频特性

5 结 论

本文从航空矢量重力测量系统对惯性器件的实际需求出发，针对传统ARMA模型只能对平稳随机漂移误差建模，且模型无法满足实时滤波需求的问题，引入了适用于非平稳随机漂移误差的ARIMA模型，并给出详细的建模过程。在此模型的基础上，基于实时平均算法的Kalman滤波器对陀螺随机漂移误差进行了实时在线滤波与补偿，并对滤波结果进行Allan方差分析。主要结论如下：

图7 光纤陀螺随机漂移数据Kalman滤波前后Allan方差曲线Fig.7 Allan-variance results of FOG random drift with and without Kalman filtering

(1) 根据光纤陀螺随机漂移误差的统计特性建立了ARIMA(1,1,1)模型，该模型能够较准确地描述陀螺的非平稳随机漂移误差。

(2) 实时平均算法可有效补偿陀螺原始采样序列中由地球自转和零漂引起的常值分量，实现随机漂移误差的实时Kalman滤波与补偿。

(3) 实测数据处理表明，基于本文所提出的方法滤波后随机漂移的主要误差较滤波前明显减小，信号中的高频噪声得到有效抑制，噪声水平降低了46.5%，角度随机游走系数和角速率随机游走系数分别降低了50%和40%。

在实际应用中，可通过约当标准型或可控标准型等方法将光纤陀螺随机漂移误差的ARIMA模型转换为等效的连续状态方程表达式，并扩充到航空矢量重力测量系统的误差方程中，通过Kalman滤波实时补偿来提高航空矢量重力测量系统的姿态解算精度。

[1] 蔡劭琨. 航空重力矢量及误差分离方法研究[D]. 长沙: 国防科学技术大学, 2014. CAI Shaokun. The Research about Airborne Vector Gravimeter and Methods of Errors Separation[D]. Changsha: National University of Defense Technology, 2014.

[2] FORSBERG R, OLESEN A V. Airborne Gravity Field Determination[M]∥XU Guochang. Sciences of Geodesy-I. Berlin Heidelberg: Springer, 2010: 83-104.

[3] LI X. Examination of Two Major Approximations Used in the Scalar Airborne Gravimetric System—A Case Study Based on the LCR System[J]. Journal of Geodetic Science, 2013, 3(1): 32-39.

[4] GLEASON D M. Gravity Vector Estimation from Integrated GPS/Strapdown IMU Data[J]. Navigation, 1992, 39(2): 237-253.

[5] JEKELI C. Airborne Vector Gravimetry Using Precise, Position-aided Inertial Measurement Units[J]. Bulletin Géodésique, 1994, 69(1): 1-11.

[6] DEURLOO R A,MARTIN J, BASTOS M L, et al. Comparison of the Performance of Two Types of Inertial Systems for Strapdown Airborne Gravimetry[C]∥Proceedings of the IEEE/ION Position Location and Navigation Symposium. 2012.

[7] AYRES-SAMPAIO D, DEURLOO R, BOS M, et al. A Comparison Between Three IMUs for Strapdown Airborne Gravimetry[J]. Surveys in Geophysics, 2015, 36(4): 571-586.

[8] 毛奔, 林玉荣. 惯性器件测试与建模[M]. 哈尔滨: 哈尔滨工程大学出版社, 2008. MAO Ben, LIN Yurong. Testing and Modelling of Inertial Device[M]. Harbin: Harbin Engineering University Press, 2008.

[9] 严恭敏, 李四海, 秦永元. 惯性仪器测试与数据分析[M]. 北京: 国防工业出版社, 2012. YAN Gongmin, LI Sihai, QIN Yongyuan. Testing and Data Analysis of Inertial Instruments[M]. Beijing: National Defense Industry Press, 2012.

[10] ORAVETZ A S, SANDBERG H J. Stationary and Nonstationary Characteristics of Gyro Drift Rate[J]. AIAA Journal, 1970, 8(10): 1766-1772.

[11] CHEN Xiyuan. Modeling Random Gyro Drift by Time Series Neural Networks and by Traditional Method[C]∥Proceedings of the 2003 International Conference on Neural Networks and Signal Processing. Nanjing, China: IEEE, 2003, 1: 810-813.

[12] 吴富梅, 杨元喜. 基于高阶AR模型的陀螺随机漂移模型[J]. 测绘学报, 2007, 36(4): 389-394. WU Fumei, YANG Yuanxi. Gyroscope Random Drift Model Based on the Higher-order AR Model[J]. Acta Geodaetica et Cartographica Sinica, 2007, 36(4): 389-394.

[13] 白俊卿, 张科, 卫育新. 光纤陀螺随机漂移建模与分析[J]. 中国惯性技术学报, 2012, 20(5): 621-624. BAI Junqing, ZHANG Ke, WEI Yuxin. Modeling and Analysis of Fiber Optic Gyroscope Random Drifts[J]. Journal of Chinese Inertial Technology, 2012, 20(5): 621-624.

[14] 王立冬, 张春熹. 高精度光纤陀螺信号的在线建模与滤波[J]. 光电工程, 2007, 34(1): 1-3, 58. WANG Lidong, ZHANG Chunxi. On-line Modeling and Filter of High-precise FOG Signal[J]. Opto-Electronic Engineering, 2007, 34(1): 1-3, 58.

[15] 李家垒, 许化龙, 何婧. 光纤陀螺随机漂移的实时滤波方法研究[J]. 宇航学报, 2010, 31(12): 2717-2721. LI Jialei, XU Hualong, HE Jing. Real-time Filtering Methods of Random Drift of Fiber Optic Gyroscope[J]. Journal of Astronautics, 2010, 31(12): 2717-2721.

[16] BOX G E P,JENKINS G M.Time Series Analysis: Forecasting and Control[M]. San Francisco: Holden-day Press, 1970.

[17] 杨叔子, 吴雅, 轩建平. 时间序列分析的工程应用[M]. 2版. 武汉: 华中科技大学出版社, 2007. YANG Shuzi, WU Ya, XUAN Jianping. Time Series Analysis in Engineering Application[M]. 2nd ed. Wuhan: Huazhong University of Science and Technology Press, 2007.

[18] GOODWIN G C,PAYNE R L.Dynamic System Identification: Experiment Design and Data Analysis[M]. New York: Academic Press, 1977.

[19] AKAIKE H. Fitting Autoregressive Models for Prediction[J]. Annals of the Institute of Statistical Mathematics, 1969, 21(1): 243-247.

[20] AKAIKE H. A Bayesian Extension of the Minimum AIC Procedure of Autoregressive Model Fitting[J]. Biometrika, 1979, 66(2): 237-242.

[21] ANSLEY G F. An Algorithm for the Exact Likelihood of a Mixed Autoregressive-moving Average Process[J]. Biometrika, 1979, 66(1): 59-65.

[22] KALMAN R E. A New Approach to Linear Filtering and Prediction Problems[J]. Journal of Basic Engineering, 1960, 82(1): 35-45.

[23] TSIEN H S. Engineering Cybernetics[M]. New York: McGraw-Hill, 1954.

[24] CANON M D, CULLUM C D JR, POLAK E. Theory of Optimal Control and Mathematical Programming[M]. New York: McGraw-Hill, 1970.

[25] IEEE Std. 952-1997 IEEE Standard Specification Format Guide and Test Procedure for Single-axis Interferometric Fiber Optic Gyros[S]. [S.l.]: IEEE, 1998.

(责任编辑：宋启凡)

Research on the Real-time Compensation of the Fiber Optic Gyroscope Random Drift in Airborne Vector Gravimetry

WANG Zheng1,LI Jiancheng1,2

1. School of Geodesy and Geomatics, Wuhan University, Wuhan 430079, China; 2. Key Laboratory of Geospace Environment and Geodesy, Ministry of Education, Wuhan 430079, China

Random drift error of fiber optic gyroscope is the crucial factor that influences the calculation accuracy of the attitude of airborne vector gravimetry. Modeling and compensating it can restrain this type of error significantly. Given the problem that traditional ARMA model can be only applied in the case of stable random drift, which cannot meet the need of real-time filtering, an ARIMA model (autoregressive integrated moving average) which is suitable for non-stable random drift is introduced along with the detailed procedure in this paper. The algorithm that can eliminate the constant component of original sampling sequence with real-time averaging method is also proposed as well as the real-time Kalman filtering estimation of the random drift. With the methods proposed above, the variance of random drift can be reduced by 46.5%. The analysis of Allan variance suggests that the coefficients of random drift for angle and angular speed have decreased about 50% and 40%, respectively. The results showed that non-stable random drift can be accurately characterized by ARIMA model and that online estimation of random drift can be realized by real-time average algorithm, indicating the potential to improve the calculation accuracy of the attitude of airborne vector gravimetry.

fiber optic gyroscope (FOG); ARIMA model; Kalman filter; Allan variance

The National Natural Science Foundation of China (Nos. 41210006；41504016); The National Basic Research Program of China (973 Program) (No. 2013CB733301); Basic Research Foundation Program of Surveying by Key Laboratory of Geospace Environment and Geodesy, Ministry of Education (No. 14-02-01)

WANG Zheng (1986—), female, PhD candidate, majors in satellite geodesy, airborne gravimetry.

王峥，李建成.航空矢量重力测量中光纤陀螺随机漂移误差实时补偿方法[J].测绘学报，2017,46(2)：144-150.

10.11947/j.AGCS.2017.20160174. WANG Zheng,LI Jiancheng.Research on the Real-time Compensation of the Fiber Optic Gyroscope Random Drift in Airborne Vector Gravimetry[J]. Acta Geodaetica et Cartographica Sinica,2017,46(2):144-150. DOI:10.11947/j.AGCS.2017.20160174.

P224

A

1001-1595(2017)02-0144-07

国家自然科学基金(41210006；41504016)；国家973计划(2013CB733301)；地球空间环境与大地测量教育部重点实验室测绘基础研究基金(14-02-01)

2016-04-18

王峥(1986—)，女，博士生，研究方向为卫星大地测量、航空重力测量。

E-mail： zhengwang@whu.edu.cn

修回日期： 2016-12-16