CN112797967B - 一种基于图优化的mems陀螺仪随机漂移误差补偿方法 - Google Patents
一种基于图优化的mems陀螺仪随机漂移误差补偿方法 Download PDFInfo
- Publication number
- CN112797967B CN112797967B CN202110134703.3A CN202110134703A CN112797967B CN 112797967 B CN112797967 B CN 112797967B CN 202110134703 A CN202110134703 A CN 202110134703A CN 112797967 B CN112797967 B CN 112797967B
- Authority
- CN
- China
- Prior art keywords
- data
- equation
- follows
- model
- value
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 69
- 238000005457 optimization Methods 0.000 title claims abstract description 36
- 238000005259 measurement Methods 0.000 claims abstract description 32
- 241001123248 Arma Species 0.000 claims abstract description 23
- 230000003068 static effect Effects 0.000 claims abstract description 14
- 238000010606 normalization Methods 0.000 claims abstract description 13
- 238000007781 pre-processing Methods 0.000 claims abstract description 12
- 238000004364 calculation method Methods 0.000 claims abstract description 11
- 230000006870 function Effects 0.000 claims description 50
- 239000011159 matrix material Substances 0.000 claims description 35
- 239000013598 vector Substances 0.000 claims description 32
- 238000012360 testing method Methods 0.000 claims description 24
- 238000004458 analytical method Methods 0.000 claims description 9
- 238000005311 autocorrelation function Methods 0.000 claims description 9
- 230000008569 process Effects 0.000 claims description 9
- 238000005070 sampling Methods 0.000 claims description 7
- 230000002159 abnormal effect Effects 0.000 claims description 6
- 238000007689 inspection Methods 0.000 claims description 6
- 238000012545 processing Methods 0.000 claims description 6
- 238000010998 test method Methods 0.000 claims description 6
- 230000007704 transition Effects 0.000 claims description 6
- 238000009826 distribution Methods 0.000 claims description 5
- 238000010276 construction Methods 0.000 claims description 3
- 238000010219 correlation analysis Methods 0.000 claims description 3
- 238000013016 damping Methods 0.000 claims description 3
- 230000002441 reversible effect Effects 0.000 claims description 3
- 238000000528 statistical test Methods 0.000 claims description 3
- 230000009466 transformation Effects 0.000 claims description 3
- 238000012163 sequencing technique Methods 0.000 claims description 2
- 238000013112 stability test Methods 0.000 claims description 2
- 230000000694 effects Effects 0.000 description 4
- 230000006872 improvement Effects 0.000 description 3
- 238000013528 artificial neural network Methods 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 238000012300 Sequence Analysis Methods 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 238000003556 assay Methods 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000002349 favourable effect Effects 0.000 description 1
- 238000001914 filtration Methods 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 230000001737 promoting effect Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000012731 temporal analysis Methods 0.000 description 1
- 238000000700 time series analysis Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
- G01C19/00—Gyroscopes; Turn-sensitive devices using vibrating masses; Turn-sensitive devices without moving masses; Measuring angular rate using gyroscopic effects
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
- G01C25/00—Manufacturing, calibrating, cleaning, or repairing instruments or devices referred to in the other groups of this subclass
- G01C25/005—Manufacturing, calibrating, cleaning, or repairing instruments or devices referred to in the other groups of this subclass initial alignment, calibration or starting-up of inertial devices
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Manufacturing & Machinery (AREA)
- Complex Calculations (AREA)
Abstract
本发明公开了一种基于图优化的MEMS陀螺仪随机漂移误差补偿方法，该方法步骤如下：使用单轴速率转台采集MEMS陀螺仪1800±10s静态数据，对得到的原始数据进行预处理操作，包括去除奇异值、常值项、趋势项，以及数据平稳性和正态性检验；对预处理后的数据进行自相关特性和偏相关特性的分析，利用AIC、BIC准则对时间序列的ARMA模型进行定阶后，使用矩估计方法对ARMA模型进行参数识别并建模；建立线性离散的图优化状态模型和量测模型，并使用Levenberg‑Marquart的图优化计算方法对MEMS陀螺仪随机漂移误差进行计算并进行补偿。本发明具有高效准确、计算量小的优点。
Description
技术领域
本发明涉及惯性传感器的输出误差补偿技术，特别是一种基于图优化的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σ准则)，该准则表示如下：
其中，xi为采集上位机软件的角速率，为采集静态角速率的均值，σ为采集角速率的方差标准差值，/>和σ定义如下：
其中，x1,x2,…xn为静态角速率采样点数值，N为采样点数量；
如果采集的陀螺仪数据满足式(1)，则表明此时陀螺仪的数据异常，丢弃该值，并利用上一个采样时刻采集的角速率数据代替此刻的角速率数据；
(2.3)去除常值项
以MEMS陀螺z轴数据为例，对MEMS陀螺z轴数据去均值操作，可以扣除常数误差项以及敏感的地球自转角速度，得到z轴陀螺仪的随机噪声；
(2.4)趋势项去除
受外界因素的影响，实际输出的数据可能存在随时间变化的趋势项；趋势项的存在会导致随机误差序列出现不平稳情况，因此需对趋势性进行处理，采用最小二乘法来拟合趋势项的参数；
设样本的m阶多项式函数为：
其中n为样本数量，ak是待定系数(k＝0,…m，m<n)，拟合的标准是使待拟合函数y(x)在xi处的函数值yi,i＝1,2,...n与f(xi)之间的距离δi的平方和最小；
其中δi＝|f(xi)-yi|；
记
J为最优函数，为使J达到最小，只需使极值的必要条件则关于a1,...,am的线性方程组为
记矩阵矩阵A＝[a0,…,am]T，矩阵Y＝[y1,…,yn]T；
则式(5)可表示为
RTRA＝RTY (6)
{1,x,…,xm}线性无关，RTR可逆，则可得系数方程的最优解为A＝(RTR)-1RTY；
得到预处理后的陀螺仪随机误差数据后，需要对数据进行统计特性的检验，判断数据是否满足时间序列建模的条件，若数据不满足建模条件，则需对数据进行处理以满足建模要求；
(2.5)数据平稳性检验
采用基于Spearman相关系数的Daniel检验法对数据进行平稳性检验；Spearman相关系数是一种秩相关系数；将x1,x2,…,xn顺序排序，顺序统计量为x(1),x(2),…,x(n)；如果xi＝x(k)，那么称K是xi在x1,x2,…,xn中的秩，记为Ri(i＝1,2,…n)，称Ri为第i个秩统计量，R1,R2,…Rn则总称为秩统计量；若记Rt＝R(Xt)为时间序列样本X1,X2,…Xn的秩，那么该样本的Spearman相关系数qs为
构造统计量T
做出如下假设检验：H0：时间序列Xt平稳，H1：序列Xt非平稳；则根据Daniel检验法：给定显著水平α，tα/2(n)为服从自由度为n的t分布的随机变量X满足概率p(X＞tα/2(n))＝α/2的值；当|T|＞tα/2(n-2)时，那么拒绝H0，时间序列不平稳；如果|T|≤tα/2(n-2)，则接受H0，认为序列Xt是平稳序列；
(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)自相关特性分析
自相关是指时间序列中每个值xt,xt-1,…xt-k之间的相关关系，序列的自相关程度可由自相关系数度量；
对于一个平稳序列xt的样本x1,x2,…xn，常用样本均值来估计随机序列均值，即
随机序列的自协方差函数γk定义如下：
随机序列的自相关函数ρk定义如下：
式中，Kd为滞后阶数；
(3.2)偏相关特性分析
偏相关是指时间序列中xt与xt-k之间的条件相关关系，序列的偏相关程度可由偏相关系数来度量，偏相关系数定义如下：
式中，
(3.3)AIC定阶准则
AIC在提高模型拟合度的基础上引入了惩罚项的概念，而惩罚项的引入使得模型参数尽可能少，这将有助于降低过拟合的可能性；ARMA(p,q)的AIC的定阶准则为：选取随机数p,q使得
AIC＝nlnσ2+2(p+q+1) (15)
达到最小值，其中n为样本大小，σ2为拟合残差方差；
(3.4)BIC准则
BIC的惩罚项比AIC的大，当系统的序列样本数目过多时，BIC的惩罚项可有效防止模型的过拟合现象；BIC准则具体算法如下：
选取p,q使得
BIC＝nlnσ2+(p+q+1)lnn (16)
达到最小值；
ARMA模型阶次越高，计算将越复杂，且实时性会变差，因此，在实际应用中，p和q一般不超过3；并且考虑到系统的可建模性，一般要求自回归阶数p大于滑动平均阶数(q)；
(3.5)参数识别
采用矩估计的方法估计模型参数；采用先求列向量再求列向量(θ1,θ2,…θp)T的方法，主要分为以下三步；
自回归参数可以通过求解Yule-Walker方程得到，而Yule-Walker方程则是由AR(p)模型的自相关函数组构成；
AR(p)模型的自相关函数可以写为如下形式：
且ρk满足ρ0＝1,ρk＝ρ-k，式(3.20)中取k＝1,2…,p，则可以得到Yule-Walker方程：
则可求得自回归系数为
然后，令随机变量那么Yt的自协方差函数γk,E(Yt+kYt)为
把Yt近似看做MA(q)序列，那么
Yt≈Xt-θ1Xt-1-…-θqXt-q (21)
MA(q)序列的自协方差满足
其中，为Xt的方差，式(22)可以改写为如下方程：
θ＝(θ1,θ2,…θq)T与的值可由线性迭代算法求出；即给定θ＝(θ1,θ2,…θq)T与θ＝0,/>的一组初值(比如θ＝0,/>)带入式(22)右边，可得到等式左边的值θ(1)＝(θ1(1),θ2(1),…θq(1))T与/>此为θ与/>的一步迭代值，再将它们带入到式(22)的右边得到θ(2)与/>此为θ与/>的二步迭代值；不断迭代，直到某步的θ(m)与/>与前一步的θ(m-1)与/>相差不大为止；此时θ(m)即为所求滑动平均参数值。
进一步地，步骤3中建立线性离散的图优化状态模型和量测模型，并使用Levenberg-Marquart的图优化计算方法对MEMS陀螺仪随机漂移误差进行计算并进行补偿，具体如下：
(4.1)线性离散系统的状态方程描述如下：
Xk＝Φk,k-1Xk-1+Γk-1Wk-1 (24)
其中，Xk是tk时刻的系统状态；Φk,k-1是系统从第tk-1时刻转移到第tk时刻的状态转移矩阵；Γk-1是系统tk-1时刻的过程噪声驱动矩阵；Wk-1描述系统噪声矩阵；
线性离散系统的量测方程描述如下：
Zk＝HkXk+Vk (25)
其中，Zk表示tk时刻的量测值；Hk表示系统在tk时刻的量测矩阵；Vk表示系统在tk时刻的噪声矩阵；
同时，系统噪声和量测噪声需满足如下条件：
其中，Qk是系统噪声方差矩阵；Rk是系统的量测噪声方差矩阵；Qk、Rk均具有正定性；
状态和量测模型的误差函数定义如下：
其中，表示根据真实状态向量和基于状态变换方程的预测状态向量计算出的状态预测误差；/>表示所获得的量测向量和实际量测向量之间的量测误差；
可以通过最小化公式(27-28)中列出的误差来获得状态向量Xk的最佳状态估计；定义总的代价函数L(Xk)如下：
在式(29)中，状态向量Xk的估计被用来最小化代价函数L(Xk)，过去时刻的状态向量也可以加入到方程(31)中，并且可以将它们视为通过最小化代价函数共同估算的变量；当将更多的状态向量添加到公式(29)时，新的代价函数编写如下：
其中Ns表示优化中使用的状态向量的数量；
(4.2)图优化模型的Levenberg-Marquart计算方法
方程(30)中定义了代价函数；接下来的步骤是求解方程并找到最佳估计；在这里，使用Levenberg-Marquart算法来求解优化方程；在更一般的形式中，代价函数可以重写为如下形式；
L(X)＝[F(X)ΩF(X)T] (31)
其中，X代表待估计的参数，Ω代表方差阵的逆矩阵的方差；F(X)代表代价函数的误差函数；以式(31)为例来说明Levenberg-Marquart算法的求解过程，可以分为以下三个步骤：
第一步，以一阶泰勒级数展开L(X)，假设有一个较好的初值那么误差函数可以写成如下形式：
其中，Jk是误差函数在初值附近的雅克比矩阵；
将式(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自回归滑动平均序列模型。设Xt是平稳序列，并且满足下列模型：
其中Xt是阶数为p或q的自回归滑动平均序列，该模型反映着t时刻状态值和过去q个时刻白噪声以及p个过去时刻的状态值之间的关系。并且由该模型看出，状态值Xt与不仅和前q个时刻的白噪声存在相关性，还和p个过去时刻的状态值存在相关性。
(1.2)时间序列数据的采集与预处理。首先采集陀螺仪静止数据，在实际测试环境中，信号受到干扰或者震动会使陀螺仪输出中存在奇异值，奇异值的存在会影响随机误差建模精度和后续惯导解算效果。因此，需要剔除异常数据。为了得到随机噪声序列，还需对陀螺仪数据进行去均值处理。通过去均值处理可以扣除数据中常数误差项以及敏感的地球自转角速度，得到陀螺仪的随机噪声。此外，由于外界因素的影响，实际输出的数据可能存在随时间变化的趋势项。趋势项的存在会导致随机误差序列出现不平稳情况。因此需对趋势性进行处理。
(1.2.1)MEMS数据采集。利用单轴速率转台采集常温情况下陀螺的输出角速率。将陀螺仪置于单轴转台上，给定外部转台的输入为零，利用上位机软件保存采集的数据，得到陀螺仪三轴1800s的输出信号。
(1.2.2)去除奇异值。在实际测试环境中，信号受到干扰或者震动会使陀螺仪输出中存在奇异值，奇异值的存在会影响随机误差建模精度和后续惯导解算效果。因此，需要剔除异常数据。在静止情况下，常用的剔除奇异值的方法为拉布达准则(3σ准则)，该准则表示如下：
其中，xi为采集上位机软件的角速率，为采集静态角速率的均值，σ为采集角速率的方差标准差值，/>和σ定义如下：
如果采集的陀螺仪数据满足式(2)，则表明此时陀螺仪的数据异常，丢弃该值，并利用上一个采样时刻采集的角速率数据代替此刻的角速率数据。
(1.2.3)去除常值项。以MEMS陀螺z轴数据为例。对MEMS陀螺z轴数据去均值操作，可以扣除常数误差项以及敏感的地球自转角速度，得到z轴陀螺仪的随机噪声。
(1.2.4)去除趋势项。受外界因素的影响，实际输出的数据可能存在随时间变化的趋势项。趋势项的存在会导致随机误差序列出现不平稳情况。因此需对趋势性进行处理。采用最小二乘法来拟合趋势项的参数。设样本的m阶多项式函数为：
其中n为样本数量，ak是待定系数(k＝0,…m，m<n)，拟合的标准是使yi,i＝1,2,...n与f(xi)之间的距离δi的平方和最小。
记
为使J达到最小，只需使极值的必要条件则关于a1,...,am的线性方程组为
记A＝[a0,…,am]T,Y＝[y1,…,yn]T。
则式(6)可表示为
RTRA＝RTY (7)
{1,x,…,xm}线性无关，RTR可逆，则可得系数方程的最优解为A＝(RTR)-1RTY。
得到预处理后的陀螺仪随机误差数据后，需要对数据进行统计特性的检验，判断数据是否满足时间序列建模的条件，若数据不满足建模条件，则需对数据进行处理以满足建模要求。
(1.2.5)对数据平稳性进行检验。采用基于Spearman相关系数的Daniel检验法对数据进行平稳性检验。Spearman相关系数是一种秩相关系数。设x1,x2,…,xn是从一元总体样本抽取的容量为n的样本，该样本顺序统计量为x(1),x(2),…,x(n)。如果xi＝x(k)，那么称K是xi在该样本中的秩，记为Ri(i＝1,2,…n)，称Ri为第i个秩统计量。R1,R2,...Rn则总称为秩统计量。记Rt＝R(Xt)为时间序列样本X1,X2,…Xn的秩，那么该样本的Spearman相关系数qs为
构造统计量
做出如下假设检验：H0：时间序列Xt平稳，H1：序列Xt非平稳。
则Daniel检验法为：给定显著水平α，当|T|＞tα/2(n-2)时，那么拒绝H0，时间序列不平稳。如果|T|≤tα/2(n-2)，则接受H0，认为序列Xt是平稳序列。
(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)进行自相关特性分析。自相关是指时间序列中每个值xt,xt-1,…xt-k之间的相关关系，序列的自相关程度可由自相关系数度量。对于一个平稳序列xt的样本x1,x2,…xn，常用样本均值来估计随机序列均值，即
随机序列的自协方差函数定义如下：
随机序列的自相关函数定义如下：
式中，Kd为滞后阶数。
(2.1.2)进行偏相关特性分析。偏相关是指时间序列中xt与xt-k之间的条件相关关系，序列的偏相关程度可由偏相关系数来度量，偏相关系数定义如下：
/>
式中，
(2.2)利用AIC和BIC准则确定数据的自回归、滑动平均的阶数。
(2.2.1)使用AIC定阶准则。AIC在提高模型拟合度的基础上引入了惩罚项的概念，而惩罚项的引入使得模型参数尽可能少，这将有助于降低过拟合的可能性。ARMA(p,q)的AIC的定阶准则为：选取p,q使得
AIC＝nlnσ2+2(p+q+1) (16)
达到最小值，其中n为样本大小，σ2为拟合残差方差。
(2.2.2)使用BIC准则。BIC的惩罚项比AIC的大，当系统的序列样本数目过多时，BIC的惩罚项可有效防止模型的过拟合现象。BIC准则具体算法如下：选取p,q使得
BIC＝nlnσ2+(p+q+1)lnn (17)
达到最小值。
ARMA模型阶次越高，计算将越复杂，且实时性会变差，因此，在实际应用中，p和q一般不超过3。并且考虑到系统的可建模性，一般要求自回归阶数(p)大于滑动平均阶数(q)。
(2.3)使用P-W法提出的一种适合于工程应用的ARMA(2n,2n-1)模型并进行参数识别。
采用矩估计的方法估计模型参数。采用先求再求(θ1,θ2,…θp)T的方法，主要分为以下三步。自回归参数可以通过求解Yule-Walker方程得到，而Yule-Walker方程则是由AR(p)模型的自相关函数组构成。
AR(p)模型的自相关函数可以写为如下形式：
且ρk满足ρ0＝1,ρk＝ρ-k，取k＝1,2…,p，则可以得到Yule-Walker方程：
则可求得自回归系数为：
然后，令那么Yt的自协方差函数为：
把Yt近似看做MA(q)序列，那么
Yt≈εt-θ1εt-1-…-θqεt-q (22)
MA(q)序列的自协方差满足
其中，为εt的方差，改写为如下方程：
θ＝(θ1,θ2,…θq)T与的值可由线性迭代算法求出。即给定θ＝(θ1,θ2,…θq)T与θ＝0,/>的一组初值(比如θ＝0,/>)带入式右边，可得到等式左边的值θ(1)＝(θ1(1),θ2(1),…θq(1))T与/>此为θ与/>的一步迭代值，再将它们带入到式的右边得到θ(2)与/>此为θ与/>的二步迭代值。不断迭代，直到某步的θ(m)与/>与前一步的θ(m-1)与/>相差不大为止。此时θ(m)即为所求滑动平均参数值。
步骤3：建立线性离散的图优化状态模型和量测模型，并使用Levenberg-Marquart算法进行计算。
(3.1)建立图优化状态模型和量测模型。线性离散系统的状态方程描述如下：
Xk＝Φk,k-1Xk-1+Γk-1Wk-1 (25)
其中，Xk是tk时刻的系统状态；Φk,k-1是系统从第tk-1时刻转移到第tk时刻的状态转移矩阵；Γk-1是系统tk-1时刻的过程噪声驱动矩阵；Wk-1描述系统噪声矩阵。
线性离散系统的量测方程描述如下：
Zk＝HkXk+Vk (26)
其中，Zk表示tk时刻的量测值；Hk表示系统在tk时刻的量测矩阵；Vk表示系统在tk时刻的噪声矩阵。
同时，系统噪声和量测噪声需满足如下条件：
其中，Qk是系统噪声方差矩阵；Rk是系统的量测噪声方差矩阵；Qk、Rk均具有正定性。状态和量测模型的误差函数定义如下：
其中，表示根据真实状态向量和基于状态变换方程的预测状态向量计算出的状态预测误差。/>表示所获得的量测向量和实际量测向量之间的量测误差。可以通过最小化公式(28-29)中列出的误差来获得状态向量Xk的最佳状态估计。定义总的代价函数L(Xk)如下：/>
在式(30)中，状态向量Xk的估计被用来最小化代价函数L(Xk)，过去时刻的状态向量也可以加入到方程(30)中，并且可以将它们视为通过最小化代价函数共同估算的变量。当将更多的状态向量添加到公式(30)时，新的代价函数编写如下：
其中N表示优化中使用的状态向量的数量。
(3.2)使用图优化模型的Levenberg-Marquart计算方法。方程(31)中定义了代价函数。接下来的步骤是求解方程并找到最佳估计。在这里，使用Levenberg-Marquart算法来求解优化方程。在更一般的形式中，代价函数可以重写为如下形式。
L(X)＝[F(X)ΩF(X)T] (32)
其中，X代表待估计的参数，Ω代表方差阵的逆矩阵的方差。F(X)代表代价函数的误差函数。以式(32)为例来说明Levenberg-Marquart算法的求解过程，可以分为以下三个步骤：
第一步，以一阶泰勒级数展开L(X)，假设有一个较好的初值那么误差函数可以写成如下形式：
其中，Jk是误差函数在初值附近的雅克比矩阵。
将式(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陀螺仪随机漂移误差补偿值。
Claims (1)
1.一种基于图优化的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σ准则，该准则表示如下：
其中，xi为采集上位机软件的角速率，为采集静态角速率的均值，σ为采集角速率的方差标准差值，/>和σ定义如下：
其中，x1,x2,…xn为静态角速率采样点数值，N为采样点数量；
如果采集的陀螺仪数据满足式(1)，则表明此时陀螺仪的数据异常，丢弃该值，并利用上一个采样时刻采集的角速率数据代替此刻的角速率数据；
(2.3)去除常值项
对MEMS陀螺z轴数据去均值操作，扣除常数误差项以及敏感的地球自转角速度，得到z轴陀螺仪的随机噪声；
(2.4)趋势项去除
受外界因素的影响，实际输出的数据可能存在随时间变化的趋势项；趋势项的存在会导致随机误差序列出现不平稳情况，因此需对趋势性进行处理，采用最小二乘法来拟合趋势项的参数；
设样本的m阶多项式函数为：
其中n为样本数量，ak是待定系数，k＝0,…m，m<n，拟合的标准是使待拟合函数y(x)在xi处的函数值yi,i＝1,2,...n与f(xi)之间的距离δi的平方和最小；
其中δi＝|f(xi)-yi|；
记
J为最优函数，为使J达到最小，只需使极值的必要条件则关于a1,...,am的线性方程组为
记矩阵矩阵A＝[a0,…,am]T，矩阵Y＝[y1,…,yn]T；
则式(5)表示为
RTRA＝RTY (6)
{1,x,…,xm}线性无关，RTR可逆，则得系数方程的最优解为A＝(RTR)-1RTY；
得到预处理后的陀螺仪随机误差数据后，需要对数据进行统计特性的检验，判断数据是否满足时间序列建模的条件，若数据不满足建模条件，则需对数据进行处理以满足建模要求；
(2.5)数据平稳性检验
采用基于Spearman相关系数的Daniel检验法对数据进行平稳性检验；Spearman相关系数是一种秩相关系数；将x1,x2,...,xn顺序排序，顺序统计量为x(1),x(2),...,x(n)；如果xi＝x(k)，那么称K是xi在x1,x2,...,xn中的秩，记为Ri，i＝1,2,…n，称Ri为第i个秩统计量，R1,R2,...Rn则总称为秩统计量；若记Rt＝R(Xt)为时间序列样本X1,X2,…Xn的秩，那么该样本的Spearman相关系数qs为
构造统计量T
做出如下假设检验：H0：时间序列Xt平稳，H1：序列Xt非平稳；则根据Daniel检验法：给定显著水平α，tα/2(n)为服从自由度为n的t分布的随机变量X满足概率p(X＞tα/2(n))＝α/2的值；当|T|＞tα/2(n-2)时，那么拒绝H0，时间序列不平稳；如果|T|≤tα/2(n-2)，则接受H0，认为序列Xt是平稳序列；
(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)自相关特性分析
自相关是指时间序列中每个值xt,xt-1,…xt-k之间的相关关系，序列的自相关程度由自相关系数度量；
对于一个平稳序列xt的样本x1,x2,…xn，常用样本均值来估计随机序列均值，即
随机序列的自协方差函数γk定义如下：
随机序列的自相关函数ρk定义如下：
式中，Kd为滞后阶数；
(3.2)偏相关特性分析
偏相关是指时间序列中xt与xt-k之间的条件相关关系，序列的偏相关程度由偏相关系数来度量，偏相关系数定义如下：
式中，j＝1,2,…k-1；
(3.3)AIC定阶准则
AIC在提高模型拟合度的基础上引入了惩罚项的概念，而惩罚项的引入使得模型参数尽可能少；ARMA(p,q)的AIC的定阶准则为：
选取随机数p,q使得
AIC＝n lnσ2+2(p+q+1) (15)
达到最小值，其中n为样本大小，σ2为拟合残差方差；
(3.4)BIC准则
BIC准则具体算法如下：
选取p,q使得
BIC＝n lnσ2+(p+q+1)ln n (16)
达到最小值；
(3.5)参数识别
采用矩估计的方法估计模型参数；
采用先求列向量再求列向量(θ1,θ2,…θp)T的方法，分为以下三步；
自回归参数可以通过求解Yule-Walker方程得到，而Yule-Walker方程则是由AR(p)模型的自相关函数组构成；
AR(p)模型的自相关函数写为如下形式：
且ρk满足ρ0＝1,ρk＝ρ-k，式中取k＝1,2…,p，则得到Yule-Walker方程：
求得自回归系数为
然后，令随机变量那么Yt的自协方差函数γk,E(Yt+kYt)为
把Yt近似看做MA(q)序列，那么
Yt≈Xt-θ1Xt-1-…-θqXt-q (21)
MA(q)序列的自协方差满足
其中，为Xt的方差，式(22)可以改写为如下方程：
θ＝(θ1,θ2,…θq)T与的值由线性迭代算法求出；即给定θ＝(θ1,θ2,…θq)T与/>的一组初值带入式(22)右边，得到等式左边的值θ(1)＝(θ1(1),θ2(1),…θq(1))T与/>此为θ与/>的一步迭代值，再将它们带入到式(22)的右边得到θ(2)与/>此为θ与/>的二步迭代值；不断迭代，直到某步的θ(m)与/>与前一步的θ(m-1)与/>相差不大于设定值为止；此时θ(m)即为所求滑动平均参数值；
步骤3中建立线性离散的图优化状态模型和量测模型，并使用Levenberg-Marquart的图优化计算方法对MEMS陀螺仪随机漂移误差进行计算并进行补偿，具体如下：
(4.1)线性离散系统的状态方程描述如下：
Xk＝Φk,k-1Xk-1+Γk-1Wk-1 (24)
其中，Xk是tk时刻的系统状态；Φk,k-1是系统从第tk-1时刻转移到第tk时刻的状态转移矩阵；Γk-1是系统tk-1时刻的过程噪声驱动矩阵；Wk-1描述系统噪声矩阵；
线性离散系统的量测方程描述如下：
Zk＝HkXk+Vk (25)
其中，Zk表示tk时刻的量测值；Hk表示系统在tk时刻的量测矩阵；Vk表示系统在tk时刻的噪声矩阵；
同时，系统噪声和量测噪声需满足如下条件：
其中，Qk是系统噪声方差矩阵；Rk是系统的量测噪声方差矩阵；Qk、Rk均具有正定性；
状态和量测模型的误差函数定义如下：
其中，表示根据真实状态向量和基于状态变换方程的预测状态向量计算出的状态预测误差；/>表示所获得的量测向量和实际量测向量之间的量测误差；
通过最小化公式(27)(28)中列出的误差来获得状态向量Xk的最佳状态估计；
定义总的代价函数L(Xk)如下：
在式(29)中，状态向量Xk的估计被用来最小化代价函数L(Xk)，过去时刻的状态向量也可以加入到方程(31)中，并且将它们视为通过最小化代价函数共同估算的变量；当将更多的状态向量添加到公式(29)时，新的代价函数编写如下：
其中Ns表示优化中使用的状态向量的数量；
(4.2)图优化模型的Levenberg-Marquart计算方法
方程(30)中定义了代价函数，接下来的步骤是求解方程并找到最佳估计；
使用Levenberg-Marquart算法来求解优化方程；
代价函数重写为如下形式；
L(X)＝[F(X)ΩF(X)T] (31)
其中，X代表待估计的参数，Ω代表方差阵的逆矩阵的方差；F(X)代表代价函数的误差函数；
针对式(31)，Levenberg-Marquart算法的求解过程，分为以下三个步骤：
第一步，以一阶泰勒级数展开L(X)，假设有一个初值那么误差函数写成如下形式：
其中，Jk是误差函数在初值附近的雅克比矩阵；
将式(32)代入式(31)，得
为了使式(33)达到最小，对式(33)求ΔX的一阶导数并令其为0，则得到下列等式
AΔX＝-B (34)
在Levenberg-Marquardt方法中，将阻尼因子添加到方程中估算增量；新的公式写为：
(A+λI)ΔX＝-B (35)
通过求解方程式(35)中的线性系统，获得增量ΔX，此时被更新为
当ΔX达到预定义的阈值或迭代计数达到设置值时，则停止；否则，返回第一步继续。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110134703.3A CN112797967B (zh) | 2021-01-31 | 2021-01-31 | 一种基于图优化的mems陀螺仪随机漂移误差补偿方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110134703.3A CN112797967B (zh) | 2021-01-31 | 2021-01-31 | 一种基于图优化的mems陀螺仪随机漂移误差补偿方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN112797967A CN112797967A (zh) | 2021-05-14 |
| CN112797967B true CN112797967B (zh) | 2024-03-22 |
Family
ID=75813354
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110134703.3A Active CN112797967B (zh) | 2021-01-31 | 2021-01-31 | 一种基于图优化的mems陀螺仪随机漂移误差补偿方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN112797967B (zh) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US11808780B1 (en) | 2022-05-20 | 2023-11-07 | Honeywell International Inc. | Inertial sensor error modeling and compensation, and system for lifetime inertial sensor calibration and navigation enhancement |
| CN117664117B (zh) * | 2024-01-31 | 2024-04-23 | 西安晟昕科技股份有限公司 | 一种光纤陀螺的漂移数据分析及优化补偿方法 |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105043384A (zh) * | 2015-04-30 | 2015-11-11 | 南京林业大学 | 一种基于鲁棒Kalman滤波的陀螺随机噪声ARMA模型建模方法 |
| CN105787265A (zh) * | 2016-02-23 | 2016-07-20 | 东南大学 | 基于综合集成赋权法的原子自旋陀螺随机误差建模方法 |
| CN109993113A (zh) * | 2019-03-29 | 2019-07-09 | 东北大学 | 一种基于rgb-d和imu信息融合的位姿估计方法 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US9269012B2 (en) * | 2013-08-22 | 2016-02-23 | Amazon Technologies, Inc. | Multi-tracker object tracking |
-
2021
- 2021-01-31 CN CN202110134703.3A patent/CN112797967B/zh active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105043384A (zh) * | 2015-04-30 | 2015-11-11 | 南京林业大学 | 一种基于鲁棒Kalman滤波的陀螺随机噪声ARMA模型建模方法 |
| CN105787265A (zh) * | 2016-02-23 | 2016-07-20 | 东南大学 | 基于综合集成赋权法的原子自旋陀螺随机误差建模方法 |
| CN109993113A (zh) * | 2019-03-29 | 2019-07-09 | 东北大学 | 一种基于rgb-d和imu信息融合的位姿估计方法 |
Non-Patent Citations (2)
| Title |
|---|
| ARMA模型的MEMS陀螺随机误差卡尔曼补偿方法;孙伟;测绘科学;第42卷(第12期);正文第52-56页 * |
| 基于因子图的组合导航方法及其可行性研究;朱晓晗;电光与控制;第26卷(第4期);正文第66-69页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN112797967A (zh) | 2021-05-14 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN112797967B (zh) | 一种基于图优化的mems陀螺仪随机漂移误差补偿方法 | |
| CN111985093A (zh) | 一种带噪声估计器的自适应无迹卡尔曼滤波状态估计方法 | |
| CN101183266A (zh) | 使用粒子滤波器估算移动机器人姿态的方法、设备和介质 | |
| CN108734218B (zh) | 一种多传感器系统的信息融合方法和装置 | |
| CN114912551B (zh) | 面向桥梁变形监测的gnss和加速度计实时融合方法 | |
| CN111561930A (zh) | 一种车载mems陀螺仪随机漂移误差的抑制方法 | |
| CN110677140B (zh) | 一种含未知输入和非高斯量测噪声的随机系统滤波器 | |
| CN115235527A (zh) | 传感器外参标定方法、装置以及电子设备 | |
| CN115727871A (zh) | 一种轨迹质量检测方法、装置、电子设备和存储介质 | |
| CN112747773A (zh) | 基于Allan方差和随机多项式提高陀螺仪精度的方法 | |
| CN111750962B (zh) | 一种基于滤波的物体重量高精度估计方法 | |
| CN114279311A (zh) | 一种基于惯性的gnss变形监测方法与系统 | |
| CN113932815A (zh) | 稳健性优化的Kalman滤波方法、装置、电子设备和存储介质 | |
| CN102679984A (zh) | 基于极小化矢量距离准则的有限模型滤波方法 | |
| CN111735478A (zh) | 基于lstm的行人实时导航零速检测方法 | |
| CN112069592B (zh) | 一种航天器外弹道跟踪测速数据特征点识别方法 | |
| CN117808052B (zh) | 基于真空环境机械臂负载自适应方法及系统 | |
| JP7120361B2 (ja) | 補正装置、姿勢算出装置、補正方法、及び、プログラム | |
| CN113588059B (zh) | 动态称重设备参数配置方法和系统 | |
| CN113076826B (zh) | 一种传感器的滤波方法和装置 | |
| CN112978532B (zh) | 电梯状态的检测方法、装置及存储介质 | |
| CN107976206B (zh) | 一种基于信息熵的mems陀螺性能评价方法 | |
| CN117609737B (zh) | 一种惯性导航系统健康状态预测方法、系统、设备及介质 | |
| Yi et al. | Random Errors Compensation method of Autoregressive Adaptive Kalman Filter based on MEMS Gyroscope | |
| CN117406759A (zh) | 一种管道机器人爬行校准方法和系统 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| CB03 | Change of inventor or designer information |
Inventor after: Chen Shuai Inventor after: Wang Bohao Inventor after: Liu Qingxiu Inventor after: Zhang Mengjun Inventor after: Zhang Min Inventor after: Cheng Yu Inventor after: Jiang Changhui Inventor after: Zhang Kun Inventor before: Wang Bohao Inventor before: Liu Qingxiu Inventor before: Zhang Mengjun Inventor before: Zhang Min Inventor before: Cheng Yu Inventor before: Chen Shuai Inventor before: Jiang Changhui Inventor before: Zhang Kun |
|
| CB03 | Change of inventor or designer information | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |