基于Allan方差与ARMA模型分析提高陀螺仪测量精度的方法

技术领域

本发明涉及一种精确建模惯性导航系统用陀螺仪随机误差的方 法，具体是指基于Allan方差与ARMA模型分析提高陀螺仪测量精度 的方法。

背景技术

陀螺仪作为惯性导航系统中角运动传感器，其测量精度直接决定 了导航系统的性能，目前得到广泛应用的陀螺仪包括激光陀螺仪、光 纤陀螺仪、机械转子陀螺仪和微机电陀螺仪。陀螺仪的输出数据中不 可避免的包含了测量误差，测量误差可分为确定性误差和随机性误差 两类；其中确定性误差包括零偏误差、刻度因子误差、安装误差等， 可通过高精度转台进行误差标定，依据误差模型进行补偿；而随机误 差则包括量化噪声误差、角度随机游走误差和速率随机游走误差等， 通常无法根据实际物理模型进行建模，只能根据数据的统计规律建立 统计误差模型，再对各种误差成分进行评估。实际应用中，在补偿了 陀螺仪确定性误差的基础上，随机误差成为了影响陀螺仪乃至惯性导 航系统性能的主要因素。因此，建立陀螺仪随机误差的时域误差模型 并对其进行补偿，是提高陀螺仪性能并增强其工程实用性的重要环节。

Allan方差方法最初是由美国国家标准局的David Allan提出的。它 是一种时域分析方法，可以利用它来确定数据噪声中随机过程的特性， 有助于识别数据中特定噪声项的来源。噪声的Allan方差与功率谱密度 之间存在某种定量关系，利用这个关系可以根据信号在时域上的输出 数据分析出各误差源的类型，适用于对精密仪器的随机噪声研究，是 一种强有力的时域分析工具。因此，利用Allan方差方法对陀螺仪输出 的随机误差数据进行分析，可得到对应于该类陀螺仪随机误差中各种 误差源的类型。

ARMA模型（Auto-Regressive and Moving Average Model）是由美 国统计学家G.E.P.Box和英国统计学家G.M.Jenkins在二十世纪七十年 代提出的时序分析模型，是研究时间序列的重要方法，由自回归模型 （简称AR模型）与滑动平均模型（简称MA模型）为基础“混合”构 成。ARMA模型描述了一组随机时间序列中前后数据的相关性以及随 机白噪声对该组数据的影响。因此，可将陀螺仪的随机误差表示成 ARMA模型的形式，从而实现对随机误差的预测和补偿。

本发明结合上述两种方法的优点，首先利用Allan方差方法辨识出 陀螺仪中随机误差的误差源，根据误差源的时域表达形式进一步确定 该陀螺随机误差对应的ARMA模型的阶数，然后再采用ARMA模型 参数估计的格林函数法模型中各项系数，从而得到针对该陀螺仪的随 机误差在时间域下的误差模型。本发明能够有效辨识出陀螺仪随机误 差中的噪声源，进而有针对性的建立随机误差的时间序列误差模型， 为提高陀螺仪性能设计及其随机误差的补偿提供有利依据。

发明内容

本发明的目的是提供一种基于Allan方差与ARMA模型分析提高 陀螺仪测量精度的方法，该方法能够有效辨识出陀螺仪随机误差中的 噪声源，从而有针对性的建立随机误差的时间序列误差模型，为提高 陀螺仪性能设计及其随机误差的补偿提供有利依据。

本发明的技术解决方案为：基于Allan方差与ARMA模型分析提 高陀螺仪测量精度的方法，该方法包括如下步骤：

（1）采集陀螺仪在静态环境下输出的原始角速率随机误差数据： 将陀螺仪安装在三轴速率转台上，调整三轴速率转台使得陀螺仪的测 量轴与当地东向重合，保持三轴速率转台静止，从而保证陀螺仪的输 入角速度为0，记录10小时陀螺仪输出的原始角速率数据为 x⁽⁰⁾(n)(n＝1,2,…,N)，N为数据量；

（2）陀螺仪原始角速率随机误差数据野值剔除：首先采用宽度为 5的滑动窗口提取出陀螺仪输出的原始角速率数据x⁽⁰⁾(n)中的中位数序 列x⁽¹⁾(n)，然后再采用宽度为3的滑动窗口提取出序列x⁽¹⁾(n)中的中位 数序列x⁽²⁾(n)，再由获得的序列x⁽²⁾(n)构造海宁平滑滤波器，得到对原 始角速率数据的平滑估计结果x⁽³⁾(n)，最后分析序列x⁽⁰⁾(n)-x⁽³⁾(n)，看 是否有数据出现|x⁽⁰⁾(n)-x⁽³⁾(n)|＞ε，ε取值为x⁽⁰⁾(n)标准差的三倍，如 果有，则用一内插值代替x⁽⁰⁾(n)，从而得到野值剔除后的陀螺仪角速率 数据x(n)；

（3）采用Allan方差方法分析陀螺仪角速率数据x(n)，辨识出陀 螺仪角速率数据x(n)中含有的误差源：首先按照分组样本容量值从小 到大的顺序依次计算各分组样本容量对应的Allan方差及Allan标准 差，然后将Allan标准差及其对应的分组样本容量值绘制在双对数曲线 坐标系下，再根据绘制曲线的斜率范围进一步辨识出陀螺仪随机误差 数据中的各种噪声源；

（4）根据辨识出的噪声源及其对应的时间序列表达式，确定陀螺 仪随机误差所对应的等效ARMA模型中自回归部分的阶数p和滑动平 均部分的阶数q，由此得到阶数确定而系数待定的陀螺仪随机误差等效 时间序列ARMA(p,q)模型；

（5）采用ARMA模型参数估计的格林函数法计算步骤（4）获得 的陀螺仪随机误差等效时间序列ARMA(p,q)模型中的待定系数：首先 建立陀螺仪随机误差等效时间序列ARMA(p,q)模型对应的长自回归 AR(P)模型，根据陀螺仪角速率数据x(n)中的随机噪声数据，采用线 性最小二乘法估计长自回归AR(P)模型的参数，然后求出长自回归 AR(P)模型对应的格林函数，最后再利用长自回归AR(P)模型格林函数 与陀螺仪随机误差等效时间序列ARMA(p,q)模型系数的对应关系，计 算得到陀螺仪随机误差等效时间序列ARMA(p,q)模型的各项系数，从 而获得陀螺仪随机误差等效时间序列ARMA(p,q)模型；

（6）根据步骤（5）获得的陀螺仪随机误差等效时间序列ARMA(p,q) 模型，得到陀螺仪输出数据中各种随机误差项随时间变化的规律，从 而根据该模型对陀螺的输出数据进行误差补偿，以提高陀螺仪的测量 精度，同时能够为提高陀螺仪性能设计及其随机误差的补偿提供依据。

本发明的方法针对惯性导航系统中使用的多种类型的陀螺仪随机 误差的自身特点，提出一种辨识随机误差中所含的噪声源，进而有针 对性的建立随机误差时间序列模型的通用方法，该方法首先采集陀螺 仪在静态环境下输出的角速率随机误差数据，并对数据进行野值剔除 处理，然后利用Allan方差方法对不含野值的角速率数据进行分析，从 而辨识出陀螺仪随机误差中含有的噪声源，在此基础上，根据各种噪 声源对应的时间序列表达式，利用等效ARMA模型重新对角速率随机 误差数据进行时间序列建模，得到陀螺仪随机误差在时间域下的误差 模型。

本发明中，上述步骤中所述的陀螺仪包括激光陀螺仪、光纤陀螺 仪、机械转子陀螺仪和微机电陀螺仪；所述的陀螺仪随机误差包括量 化噪声误差、角度随机游走误差和速率随机游走误差。

本发明的原理是：首先采集陀螺仪在静态环境下输出的角速率随 机误差数据，并对数据进行野值剔除处理，然后利用Allan方差方法对 角速率数据进行分析，从而辨识出陀螺仪随机误差中含有的噪声源， 在此基础上，根据各种噪声源对应的时间序列表达式，利用等效ARMA 模型重新对角速率数据进行时间序列建模，得到陀螺仪随机误差在时 间域下的误差模型。本发明能够有效辨识出各种陀螺仪随机误差中的 噪声源，从而有针对性的建立随机误差的时间序列误差模型，提高了 模型的精度。

本发明与现有技术相比的优点在于：本方法是建立多种不同类型 陀螺仪随机误差模型的通用方法，将Allan方差方法中随机误差源的辨 识能力与ARMA模型的时域误差建模能力相结合，从而有效辨识出各 种陀螺仪随机误差中的噪声源，进而有针对性的建立随机误差的时间 序列误差模型，提高了模型精度。

附图说明

下面结合附图和具体实施方式对本发明做进一步详细说明。

图1为本发明的基于Allan方差与ARMA模型的陀螺仪随机误差 建模方法的流程图；

图2a为本发明的第2步中值序列提取示意图一；

图2b为本发明的第2步中值序列提取示意图二；

图3为本发明的第3步典型Allan方差图像；

图4为本发明的第4步ARMA模型参数估计的格林函数法。

具体实施方式

本发明基于Allan方差与ARMA模型分析以提高陀螺仪测量精度 的方法，具体实施步骤如下：

1、采集陀螺仪在静态环境下输出的角速率数据：将陀螺仪安装在 三轴速率转台上，调整转台使得陀螺仪的测量轴与当地东向重合，保 证转台静止，从而保证陀螺仪的角速率输入为零，设定陀螺仪输出数 据频率为f_s Hz，记录10小时陀螺仪输出的原始数据为 x⁽⁰⁾(n)(n＝1,2,…,N)，N为数据总量；

2、剔除原始数据x⁽⁰⁾(n)中的野值：首先采用宽度为5的滑动窗口 提取原始数据x⁽⁰⁾(n)的中位数序列x⁽¹⁾(n)，然后再采用宽度为3的滑动 窗口提取出序列x⁽¹⁾(n)中的中位数序列x⁽²⁾(n)，再由x⁽²⁾(n)构造海宁平 滑滤波器得到对原始角速率数据的平滑结果x⁽³⁾(n)，最后分析序列 x⁽⁰⁾(n)-x⁽³⁾(n)，看是否有数据出现|x⁽⁰⁾(n)-x⁽³⁾(n)|＞ε，ε为设定阈值，ε 取值为x⁽⁰⁾(n)标准差的三倍，如果有，则用一内插值代替x⁽⁰⁾(n)，得到 野值剔除后的陀螺仪角速率数据x(n)，具体计算步骤如下：

首先按照式（1）提取出原始数据x⁽⁰⁾(n)的中位数序列x⁽¹⁾(n)，

x⁽¹⁾(n)＝Media[x⁽⁰⁾(n-2),x⁽⁰⁾(n-1),x⁽⁰⁾(n),x⁽⁰⁾(n+1),x⁽⁰⁾(n+2)]  n＝3,4,…,N-2

                                                          （1）

同时设定x⁽¹⁾(1)＝x⁽¹⁾(2)＝x⁽¹⁾(N-1)＝x⁽¹⁾(N)＝0；然后按照式（2）提取出 序列x⁽¹⁾(n)的中位数序列，即原始数据x⁽⁰⁾(n)的二次中位数序列x⁽²⁾(n)，

x⁽²⁾(n)＝Media[x⁽¹⁾(n-1),x⁽¹⁾(n),x⁽¹⁾(n+1)]  n＝4,5,…,N-3        （2）

并设定x⁽²⁾(1)＝x⁽²⁾(2)＝x⁽²⁾(3)＝x⁽²⁾(N-2)＝x⁽²⁾(N-1)＝x⁽²⁾(N)＝0；在此基 础上，按照式（3）所述海宁平滑滤波器得到对原始数据的平滑估计结 果x⁽³⁾(n)，

$x^{\left(3\right)}\left(n\right)=\frac{1}{4}{x}^{\left(2\right)}(n-1)+\frac{1}{2}{x}^{\left(2\right)}\left(n\right)+\frac{1}{4}{x}^{\left(2\right)}(n+1),n=\mathrm{5,6},...,N-4---\left(3\right)$

同时

x⁽³⁾(1)＝x⁽³⁾(2)＝x⁽³⁾(3)＝x⁽³⁾(4)＝x⁽³⁾(N-3)＝x⁽³⁾(N-2)＝x⁽³⁾(N-1)＝x⁽³⁾(N)＝0 ；此后，按照式（4）计算得到判定阈值ε，

$ϵ=10\sqrt{\frac{1}{N-8}\underset{n=5}{\overset{N-4}{\Sigma }}{[x^{\left(0\right)}\left(n\right)-x^{\left(3\right)}\left(n\right)]}^2}---\left(4\right)$

最后，按照式（5）计算得到野值剔除后的陀螺仪角速率数据x(n)，

$x\left(n\right)=\left(\begin{array}{cc}{x}^{\left(0\right)}\left(n\right)& \mathrm{if}|x^{\left(0\right)}\left(n\right)-x^{\left(3\right)}\left(n\right)|<ϵ\\ \frac{1}{2}[x^{\left(0\right)}(n-1)+x^{\left(0\right)}(n+1)]& \mathrm{if}|x^{\left(0\right)}\left(n\right)-x^{\left(3\right)}\left(n\right)|>ϵ\end{array}\right)n=\mathrm{6,7},...N-5---\left(5\right)$

并按照x(n)＝x(n+5)n＝1,2,…,M重新标注x(n)的序号，其中M＝N-10。

3、采用Allan方差方法分析角速率数据x(n)并进一步辨识出数据 中含有的误差成分：首先按照分组样本容量值从小到大的顺序依次计 算各分组样本容量对应的Allan方差及Allan标准差，其次将Allan标 准差及其对应的分组样本容量值绘制在双对数（log-log）曲线坐标系 下，再根据绘制曲线的斜率范围辨识出陀螺仪随机误差数据中的各种 噪声源，具体计算步骤如下：

首先，按照每组m个数据对角速率数据x(n)进行分组，得到 K＝Mm组数据，计算每组数据的平均值为：

${\bar{x}}_k\left(\tau \right)=\frac{1}{m}\underset{i=1}{\overset{m}{\Sigma }}x[(k-1)m+i],k=\mathrm{1,2},...,K---\left(5\right)$

其中，τ＝mf_s表示相关时间。然后，按照（6）式和（7）式所示计算 角速率数据x(n)对应于相关时间τ的Allan方差和标准差。

${\sigma }_{\mathrm{Allan}}^2\left(\tau \right)=\frac{1}{2(K-1)}\underset{k=1}{\overset{K-1}{\Sigma }}{[{\bar{x}}_{k+1}\left(\tau \right)-{\bar{x}}_k\left(\tau \right)]}^2---\left(6\right)$

${\sigma }_{\mathrm{Allan}}\left(\tau \right)=\sqrt{{\sigma }_{\mathrm{Allan}}^2\left(\tau \right)}---\left(7\right)$

按照m＝1,2,…,m₀的顺序依次计算每一个值对应的Allan标准差的大小， 等同于计算每一个相关时间τ对应的Allan标准差的大小，其中 m₀＝int[M/7]。由此，得到对应于每一个τ_m＝m/f_s的Allan标准差序列 σ_Allan(τ_m)。

在此基础上，按照式（8）和式（9）计算τ_m和σ_Allan(τ_m)对应的对 数值，

${{\tau }_m}^{\prime }=\lg \left(\frac{{\tau }_m}{3600}\right)---\left(8\right)$

σ_Allan'(τ_m)＝lg[σ_Allan(τ_m)]                       （9）

然后以τ_m'为横轴，σ_Allan'(τ_m)为纵轴，将计算得到的数据绘制在双对数 曲线上，观察曲线的斜率，可根据斜率的大小来辨识陀螺仪随机误差 数据中的各种噪声源。

Allan标准差对数曲线的斜率与噪声源的对应关系如下：a)斜率-1 对应了量化噪声；b)斜率-1/2对应了角度随机游走噪声；c)斜率+1/2对 应了速率随机游走噪声。

4、根据辨识出的噪声源及其对应的时间序列表达式，确定陀螺仪 随机误差所对应的等效ARMA模型中自回归部分和滑动平均部分的阶 数，分别为p和q，从而得到阶数确定而系数待定的陀螺仪随机误差等 效时间序列模型ARMA(p,q)。各种噪声源与其时间序列表达式的对应 关系如下：

a)量化噪声

x_Q(n)＝w_Q(n)-w_Q(n-1)＝(1-B)w_Q(n)            （10）

其中，x_Q(n)为量化噪声，w_Q(n)为形成量化噪声的白噪声，B表示后 向移位算子，量化噪声实际对应了一阶滑动平均过程ARMA(0,1)。

b)角度随机游走

x_AW(n)＝w_AW(n)                       （11）

其中，x_AW(n)表示角度随机游走噪声，w_AW(n)为形成角度随机游走的 白噪声。

c)速率随机游走

$x_{\mathrm{RW}}\left(n\right)=x_{\mathrm{RW}}(n-1)+w_{\mathrm{RW}}\left(n\right)=\frac{1}{1-B}{w}_{\mathrm{RW}}\left(n\right)---\left(12\right)$

其中，x_RW(n)表示速率随机游走噪声，w_RW(n)为形成速率随机游走的 白噪声，B为后向移位算子，速率随机游走实际对应了一阶自回归过 程ARMA(1,0)。

由此，根据ARMA模型等效叠加理论，可根据步骤3中判断的 Allan标准差曲线的斜率做出如下判断：

（1）若斜率仅为-1，则陀螺仪原始角速率数据对应了ARMA(0,1)过程；

（2）若斜率仅为-1/2，则陀螺仪原始角速率数据对应了白噪声过程， 无需建模，直接计算x(n)方差即可；

（3）若斜率仅为+1/2，则陀螺仪原始角速率数据对应了ARMA(1,0)过 程；

（4）若斜率范围包括-1和-1/2，陀螺仪原始角速率数据对应了 ARMA(0,2)过程；

（5）若斜率范围包括-1和+1/2，陀螺仪原始角速率数据对应了 ARMA(1,1)过程；

（6）若斜率范围包括-1/2和+1/2，陀螺仪原始角速率数据对应了 ARMA(1,1)过程；

（7）若斜率范围包括-1，-1/2和+1/2，陀螺仪原始角速率数据对应了 ARMA(1,2)过程；

5、采用ARMA模型参数估计的格林函数法计算步骤4得到的 ARMA(p,q)模型中的待定系数；参数待估计的ARMA(p,q)模型的具体 形式为：

$\left(\begin{array}{c}x\left(n\right)=-a_{1}x(n-1)-a_{2}x(n-2)-...-a_{p}x(n-p)\\ +w\left(n\right)+b_{1}w(n-1)+b_{2}w(n-2)+...+b_{q}w(n-q)\end{array},---\left(13\right)\right)$

其中x(n)为步骤2得到的陀螺仪角速率随机误差数据，w(n)为白噪声 序列，系数$\beta ={[{\beta }_1^T,{\beta }_2^T]}^T={[a_1,a_2,...,a_p,b_1,b_2,...,b_q]}^T$为待估计参数。具体步 骤如下：

（1）建立ARMA(p,q)模型对应的长自回归模型AR(P)，其中阶数 P＝(lgM)^1+λ，其中0≤λ≤1；AR(P)模型的参数为θ＝(θ₁,θ₂,…,θ_P)^T；基于 数据x(n)，采用显性最小二乘法，按照（14）式对模型参数进行估计；

$\theta ={({\theta }_1,{\theta }_2,...,{\theta }_P)}^T={\left(X_P^T{X}_P\right)}^{-1}{X}_P^T{Y}_P---\left(14\right)$

其中，$X_P=\left(\begin{array}{cccc}x\left(P\right)& x(P-1)& ...& x\left(1\right)\\ x(P+1)& x\left(P\right)& ...& x\left(2\right)\\ .& .& & .\\ .& .& & .\\ .& .& & .\\ x(M-1)& x(M-2)& ...& x(M-P)\end{array}\right),$ Y_P＝[x(P+1),x(P+2),…,x(N)]^T；

（2）由AR(P)模型的参数θ＝(θ₁,θ₂,…,θ_P)^T，利用（15）式求出AR(P) 模型对应的格林函数；

$\left(\begin{array}{cc}{G}_0^{\mathrm{AR}}=1& \\ G_k^{\mathrm{AR}}=\underset{j=1}{\overset{P}{\Sigma }}{\theta }_j{G}_{k-j}& k=\mathrm{1,2},...(p+q)\end{array}\right)---\left(15\right)$

（3）利用AR(P)模型格林函数与ARMA(p,q)模型系数的关系，求 解β；具体分两步求解，首先按照式（16）求出ARMA(p,q)模型中AR 部分系数β₁＝[a₁,a₂,…,a_p]^T，

${\beta }_1={[a_1,a_2,...,a_p]}^T={\left(G_1^T{G}_1\right)}^{-1}{G}_1^T{G}_2---\left(16\right)$

其中，$G_1=\left(\begin{array}{cccc}{G}_{q-1}& G_{q-2}& ...& G_0\\ G_q& G_{q-1}& ...& G_1\\ .& .& & .\\ .& .& & .\\ .& .& & .\\ G_{q+p-1}& G_{q+p-2}& ...& G_p\end{array}\right),$G₂＝[G_q,G_q+1,…,G_q+p]^T；然后再根据 计算得到的β₁，进一步求出MA部分系数β₂＝[b₁,b₂,…,b_q]^T，此时分为 两种情况，①当p≥q时，可按照（17）式进行计算，

β₂＝[b₁,b₂,…,b_q]^T＝G₃β₁'-G₄                          （17）

其中，β₁'＝[a₁,a₂,…,a_q]^T，G₄＝[G₁,G₂,…,G_q]^T；

②当p＜q时，可按照（18）式进行计算，

β₂＝[b₁,b₂,…,b_q]^T＝G₅β₁-G₆                       （18）

其中，G₆＝[G₁,G₂,…,G_q]^T。至此得到对 应了陀螺仪随机误差的时间序列误差模型ARMA(p,q)。

（6）根据步骤（5）获得的陀螺仪随机误差等效时间序列ARMA(p,q) 模型，得到陀螺仪输出数据中各种随机误差项随时间变化的规律，从 而根据该模型对陀螺的输出数据进行误差补偿，以提高陀螺仪的测量 精度，同时能够为提高陀螺仪性能设计及其随机误差的补偿提供依据。

本发明的上述实施例并不是对本发明保护范围的限定，本发明的实 施方式不限于此，凡此种种根据本发明的上述内容，按照本领域的普 通技术知识和惯用手段，在不脱离本发明上述基本技术思想前提下， 对本发明上述结构做出的其它多种形式的修改、替换或变更，均应落 在本发明的保护范围之内。