一种基于图优化的MEMS陀螺仪随机漂移误差补偿方法

技术领域

本发明涉及惯性传感器的输出误差补偿技术，特别是一种基于图优化的MEMS陀螺仪随机漂移误差补偿方法。

背景技术

虽然MEMS陀螺仪具有体积小、成本低等的优点，但是由于MEMS陀螺仪的精度较低，输出的信号中会含有随机噪声，会对惯导系统的正常运行产生影响。对其随机噪声补偿方法进行探索改进，可以指导惯性传感器性能改进的正确方向；有利于推进MEMS惯性器件制造工艺的改进，为MEMS惯性器件的精度提升提供支持。

现有的补偿方法中，利用多结构的LSTM网络对陀螺仪随机误差进行训练并预测，确实能降低陀螺仪随机误差，然而神经网络存在过拟合、局部最优、实时性低等问题。另一种常用的手段是小波阈值去噪法，但是其阈值难以确定且具有较大的计算量，不少学者在选择阈值函数和降低计算量方面做了许多研究，但是补偿效果并没有得到明显提升。

发明内容

本发明目的在于提供一种高效、精确的基于图优化的MEMS陀螺仪随机漂移误差补偿方法。

实现本发明目的的技术解决方案为：一种基于图优化的MEMS陀螺仪随机漂移误差补偿方法，包括以下步骤：

步骤1：使用单轴速率转台采集MEMS陀螺仪1800±10s静态数据，对得到的原始数据进行预处理操作，包括去除奇异值、常值项、趋势项，以及数据平稳性和正态性检验；

步骤2：对步骤1中预处理后的数据进行自相关特性和偏相关特性的分析，利用AIC、BIC准则对时间序列的ARMA模型进行定阶后，使用矩估计方法对ARMA模型进行参数识别并建模；

步骤3：建立线性离散的图优化状态模型和量测模型，并使用Levenberg-Marquart的图优化计算方法对MEMS陀螺仪随机漂移误差进行计算并进行补偿。

本发明与现有技术相比，其显著优点为：(1)相较于传统的MEMS陀螺仪随机误差的建模方法，时间序列分析方法的ARMA模型克服了神经网络存在过拟合、局部最优、实时性低等问题，克服了小波阈值去噪法超大计算量且阈值函数难以选择的问题；(2)使用图优化方法进行分析相较于常用的卡尔曼滤波方法，包含除当前状态外的过去状态向量，可以得到更加平滑的处理结果。

附图说明

图1为基于图优化的MEMS陀螺仪随机漂移误差补偿方法的流程图。

图2为MEMS陀螺仪时间序列数据处理步骤示意图。

具体实施方式

本发明基于图优化的MEMS陀螺仪随机漂移误差补偿方法，包括以下步骤：

步骤1：使用单轴速率转台采集MEMS陀螺仪1800±10s静态数据，对得到的原始数据进行预处理操作，包括去除奇异值、常值项、趋势项，以及数据平稳性和正态性检验；

步骤2：对步骤1中预处理后的数据进行自相关特性和偏相关特性的分析，利用AIC、BIC准则对时间序列的ARMA模型进行定阶后，使用矩估计方法对ARMA模型进行参数识别并建模；

步骤3：建立线性离散的图优化状态模型和量测模型，并使用Levenberg-Marquart的图优化计算方法对MEMS陀螺仪随机漂移误差进行计算并进行补偿。

进一步地，步骤1中所述的使用单轴速率转台采集MEMS陀螺仪1800±10s静态数据，对得到的原始数据进行预处理操作，具体如下：

(2.1)MEMS数据采集

利用单轴速率转台采集常温情况下陀螺的输出角速率；将陀螺仪置于单轴转台上，给定外部转台的输入为零，利用上位机软件保存采集的数据，得到陀螺仪三轴1800±10s的输出信号；

(2.2)去除奇异值

在实际测试环境中，信号受到干扰或者震动会使陀螺仪输出中存在奇异值，奇异值的存在会影响随机误差建模精度和后续惯导解算效果；因此，需要剔除异常数据；在静止情况下，常用的剔除奇异值的方法为拉布达准则(3σ准则)，该准则表示如下：

其中，x

其中，x

如果采集的陀螺仪数据满足式(1)，则表明此时陀螺仪的数据异常，丢弃该值，并利用上一个采样时刻采集的角速率数据代替此刻的角速率数据；

(2.3)去除常值项

以MEMS陀螺z轴数据为例，对MEMS陀螺z轴数据去均值操作，可以扣除常数误差项以及敏感的地球自转角速度，得到z轴陀螺仪的随机噪声；

(2.4)趋势项去除

受外界因素的影响，实际输出的数据可能存在随时间变化的趋势项；趋势项的存在会导致随机误差序列出现不平稳情况，因此需对趋势性进行处理，采用最小二乘法来拟合趋势项的参数；

设样本的m阶多项式函数为：

其中n为样本数量，a

其中δ

记

J为最优函数，为使J达到最小，只需使极值的必要条件

记矩阵

则式(5)可表示为

R

{1,x,…,x

得到预处理后的陀螺仪随机误差数据后，需要对数据进行统计特性的检验，判断数据是否满足时间序列建模的条件，若数据不满足建模条件，则需对数据进行处理以满足建模要求；

(2.5)数据平稳性检验

采用基于Spearman相关系数的Daniel检验法对数据进行平稳性检验；Spearman相关系数是一种秩相关系数；将x

构造统计量T

做出如下假设检验：H

(2.6)数据正态性检验

采用直方图法与Jarque-Bera检验相结合，进行样本正态性的判断；直方图检验是先绘制出样本数据的概率直方图，并将其与理想的正态分布图进行比较；而Jaque-Bera检验法是基于偏度系数和峰度系数的统计检验法，其统计量JB为：

其中，N为样本数量，S为偏度系数，K为峰度系数，后二个字符分别定义如下；

设置显著水平α＝0.05，当|JB|＜5.99，可证明样本满足正态性。

进一步地，步骤2中所述的对步骤1中预处理后的数据进行自相关特性和偏相关特性的分析，利用AIC、BIC准则对时间序列的ARMA模型进行定阶后，使用矩估计方法对ARMA模型进行参数识别并建模，具体如下：

首先要对数据进行自相关和偏相关分析；其次，利用AIC和BIC准则确定数据的自回归、滑动平均的阶数；确定数据的阶数后，对模型的待定参数进行估计；最后，需要对残差数据进行白噪声检验，检验通过则表明模型建立准确；

(3.1)自相关特性分析

自相关是指时间序列中每个值x

对于一个平稳序列x

随机序列的自协方差函数γ

随机序列的自相关函数ρ

式中，K

(3.2)偏相关特性分析

偏相关是指时间序列中x

式中，

(3.3)AIC定阶准则

AIC在提高模型拟合度的基础上引入了惩罚项的概念，而惩罚项的引入使得模型参数尽可能少，这将有助于降低过拟合的可能性；ARMA(p,q)的AIC的定阶准则为：选取随机数p,q使得

AIC＝nlnσ

达到最小值，其中n为样本大小，σ

(3.4)BIC准则

BIC的惩罚项比AIC的大，当系统的序列样本数目过多时，BIC的惩罚项可有效防止模型的过拟合现象；BIC准则具体算法如下：

选取p,q使得

BIC＝nlnσ

达到最小值；

ARMA模型阶次越高，计算将越复杂，且实时性会变差，因此，在实际应用中，p和q一般不超过3；并且考虑到系统的可建模性，一般要求自回归阶数p大于滑动平均阶数(q)；

(3.5)参数识别

采用矩估计的方法估计模型参数；采用先求列向量

自回归参数可以通过求解Yule-Walker方程得到，而Yule-Walker方程则是由AR(p)模型的自相关函数组构成；

AR(p)模型的自相关函数可以写为如下形式：

且ρ

则可求得自回归系数为

然后，令随机变量

把Y

Y

MA(q)序列的自协方差满足

其中，

θ＝(θ

进一步地，步骤3中建立线性离散的图优化状态模型和量测模型，并使用Levenberg-Marquart的图优化计算方法对MEMS陀螺仪随机漂移误差进行计算并进行补偿，具体如下：

(4.1)线性离散系统的状态方程描述如下：

X

其中，X

线性离散系统的量测方程描述如下：

Z

其中，Z

同时，系统噪声和量测噪声需满足如下条件：

其中，Q

状态和量测模型的误差函数定义如下：

其中，

可以通过最小化公式(27-28)中列出的误差来获得状态向量X

在式(29)中，状态向量X

其中N

(4.2)图优化模型的Levenberg-Marquart计算方法

方程(30)中定义了代价函数；接下来的步骤是求解方程并找到最佳估计；在这里，使用Levenberg-Marquart算法来求解优化方程；在更一般的形式中，代价函数可以重写为如下形式；

L(X)＝[F(X)ΩF(X)

其中，X代表待估计的参数，Ω代表方差阵的逆矩阵的方差；F(X)代表代价函数的误差函数；以式(31)为例来说明Levenberg-Marquart算法的求解过程，可以分为以下三个步骤：

第一步，以一阶泰勒级数展开L(X)，假设有一个较好的初值

其中，J

将式(32)代入式(31)，可得

为了使式(33)达到最小，对式(33)求ΔX的一阶导数并令其为0，则可以得到下列等式

AΔX＝-B (34)

在这里，在Levenberg-Marquardt方法中，将阻尼因子添加到方程中，可以更精确地估算增量；新的公式写为：

(A+λI)ΔX＝-B (35)

通过求解方程式(35)中的线性系统，可以获得增量ΔX，此时

当ΔX达到预定义的阈值或迭代计数达到设置值时，则停止；否则，返回第一步继续。

下面结合附图及具体实施例对本发明作进一步详细说明。

实施例

结合图1，一种基于时间序列分析方法ARMA模型的图优化方法，包括以下步骤：

步骤1：使用高精度单轴速率转台采集MEMS陀螺仪1800±10s静态数据；并通过去除奇异值，常值项，趋势项和数据平稳性和正态性检验的数据预处理。结合图2，具体的实施方法如下：

(1.1)根据下面的公式建立ARMA自回归滑动平均序列模型。设X

其中X

(1.2)时间序列数据的采集与预处理。首先采集陀螺仪静止数据，在实际测试环境中，信号受到干扰或者震动会使陀螺仪输出中存在奇异值，奇异值的存在会影响随机误差建模精度和后续惯导解算效果。因此，需要剔除异常数据。为了得到随机噪声序列，还需对陀螺仪数据进行去均值处理。通过去均值处理可以扣除数据中常数误差项以及敏感的地球自转角速度，得到陀螺仪的随机噪声。此外，由于外界因素的影响，实际输出的数据可能存在随时间变化的趋势项。趋势项的存在会导致随机误差序列出现不平稳情况。因此需对趋势性进行处理。

(1.2.1)MEMS数据采集。利用单轴速率转台采集常温情况下陀螺的输出角速率。将陀螺仪置于单轴转台上，给定外部转台的输入为零，利用上位机软件保存采集的数据，得到陀螺仪三轴1800s的输出信号。

(1.2.2)去除奇异值。在实际测试环境中，信号受到干扰或者震动会使陀螺仪输出中存在奇异值，奇异值的存在会影响随机误差建模精度和后续惯导解算效果。因此，需要剔除异常数据。在静止情况下，常用的剔除奇异值的方法为拉布达准则(3σ准则)，该准则表示如下：

其中，x

如果采集的陀螺仪数据满足式(2)，则表明此时陀螺仪的数据异常，丢弃该值，并利用上一个采样时刻采集的角速率数据代替此刻的角速率数据。

(1.2.3)去除常值项。以MEMS陀螺z轴数据为例。对MEMS陀螺z轴数据去均值操作，可以扣除常数误差项以及敏感的地球自转角速度，得到z轴陀螺仪的随机噪声。

(1.2.4)去除趋势项。受外界因素的影响，实际输出的数据可能存在随时间变化的趋势项。趋势项的存在会导致随机误差序列出现不平稳情况。因此需对趋势性进行处理。采用最小二乘法来拟合趋势项的参数。设样本的m阶多项式函数为：

其中n为样本数量，a

记

为使J达到最小，只需使极值的必要条件

记

则式(6)可表示为

R

{1,x,…,x

得到预处理后的陀螺仪随机误差数据后，需要对数据进行统计特性的检验，判断数据是否满足时间序列建模的条件，若数据不满足建模条件，则需对数据进行处理以满足建模要求。

(1.2.5)对数据平稳性进行检验。采用基于Spearman相关系数的Daniel检验法对数据进行平稳性检验。Spearman相关系数是一种秩相关系数。设x

构造统计量

做出如下假设检验：H

则Daniel检验法为：给定显著水平α，当|T|＞t

(1.2.6)对数据正态性进行检验。正态性检验包括直方图法，Jarque-Bera检验等。为避免仅使用一种检验法而产生误判，可以采用直方图法与Jarque-Bera检验相结合，进行样本正态性的判断。直方图检验是先绘制出样本数据的概率直方图，并将其与理想的正态分布图进行比较，Jaque-Bera检验法是基于偏度系数和峰度系数的统计检验法，其统计量为：

其中，S为偏度系数，K为峰度系数，分别定义如下。

设置显著水平α＝0.05，当|JB|＜5.99，可证明样本满足正态性。

步骤2：对数据进行自相关特性和偏相关特性的分析，使用AIC,BIC准则对ARMA模型进行定阶；使用矩估计方法对ARMA模型进行参数识别。

(2.1)首先要对数据进行自相关和偏相关分析。确定数据的阶数后，对模型的待定参数进行估计。最后，需要对残差数据进行白噪声检验，检验通过则表明模型建立准确。

(2.1.1)进行自相关特性分析。自相关是指时间序列中每个值x

随机序列的自协方差函数定义如下：

随机序列的自相关函数定义如下：

式中，K

(2.1.2)进行偏相关特性分析。偏相关是指时间序列中x

式中，

(2.2)利用AIC和BIC准则确定数据的自回归、滑动平均的阶数。

(2.2.1)使用AIC定阶准则。AIC在提高模型拟合度的基础上引入了惩罚项的概念，而惩罚项的引入使得模型参数尽可能少，这将有助于降低过拟合的可能性。ARMA(p,q)的AIC的定阶准则为：选取p,q使得

AIC＝nlnσ

达到最小值，其中n为样本大小，σ

(2.2.2)使用BIC准则。BIC的惩罚项比AIC的大，当系统的序列样本数目过多时，BIC的惩罚项可有效防止模型的过拟合现象。BIC准则具体算法如下：选取p,q使得

BIC＝nlnσ

达到最小值。

ARMA模型阶次越高，计算将越复杂，且实时性会变差，因此，在实际应用中，p和q一般不超过3。并且考虑到系统的可建模性，一般要求自回归阶数(p)大于滑动平均阶数(q)。

(2.3)使用P-W法提出的一种适合于工程应用的ARMA(2n,2n-1)模型并进行参数识别。

采用矩估计的方法估计模型参数。采用先求

AR(p)模型的自相关函数可以写为如下形式：

且ρ

则可求得自回归系数为：

然后，令

把Y

Y

MA(q)序列的自协方差满足

其中，

θ＝(θ

步骤3：建立线性离散的图优化状态模型和量测模型，并使用Levenberg-Marquart算法进行计算。

(3.1)建立图优化状态模型和量测模型。线性离散系统的状态方程描述如下：

X

其中，X

线性离散系统的量测方程描述如下：

Z

其中，Z

同时，系统噪声和量测噪声需满足如下条件：

其中，Q

其中，

在式(30)中，状态向量X

其中N表示优化中使用的状态向量的数量。

(3.2)使用图优化模型的Levenberg-Marquart计算方法。方程(31)中定义了代价函数。接下来的步骤是求解方程并找到最佳估计。在这里，使用Levenberg-Marquart算法来求解优化方程。在更一般的形式中，代价函数可以重写为如下形式。

L(X)＝[F(X)ΩF(X)

其中，X代表待估计的参数，Ω代表方差阵的逆矩阵的方差。F(X)代表代价函数的误差函数。以式(32)为例来说明Levenberg-Marquart算法的求解过程，可以分为以下三个步骤：

第一步，以一阶泰勒级数展开L(X)，假设有一个较好的初值

其中，J

将式(33)代入式(32)，可得

为了使式(34)达到最小，对式(34)求ΔX的一阶导数并令其为0，则可以得到下列等式

AΔX＝-B (35)

在这里，在Levenberg-Marquardt方法中，将阻尼因子添加到方程中，可以更精确地估算增量。新的公式写为：

(A+λI)ΔX＝-B (36)

通过求解方程式(37)中的线性系统，可以获得增量ΔX，此时

当ΔX达到预定义的阈值或迭代计数达到设置值时，则停止。否则，返回第一步继续。

(4.3)将(4.2)中得到的图优化结果进行1σ统计，得到的统计结果即可作为MEMS陀螺仪随机漂移误差补偿值。