技术领域
[0001] 本
发明属于
传感器信号处理领域，具体涉及融合卡尔曼和小波的MEMS陀螺仪自适应抗野值去噪方法。
背景技术
[0002] 微
电子技术生产的MEMS(Microelectro Mechanical Systems)传感器是一种新型的微型惯性传感器，在低成本的捷联惯导系统中被广泛应用。使用MEMS陀螺仪测量数据前，为了减少其随机漂移误差，提高传感器
精度，需对其进行去噪处理。
[0003] 对MEMS陀螺仪信号去噪，一般采用卡尔曼滤波或小波分析，卡尔曼滤波计算简单易处理，但实际过程中由于各种因素的影响，会导致
滤波器发散，往往通过增加当前量测值的权重抑制滤波器的发散，且信号中的野值等干扰数据也会影响滤波器的精度，影响去噪效果。小波分析常用于
信号处理。通过对信号的分解，对其高频部分
阈值化处理，再对其低频和处理后的高频信号进行小波重构得到去噪信号。当大量的噪声都处于高频分量中，小波去噪对信号中的噪声有良好的抑制效果。但MEMS陀螺仪信号中的
角度随机游走和零偏不
稳定性常常处于低频分量中。
发明内容
[0004] 本发明的技术问题是MEMS陀螺仪测量数据往往存在随机漂移误差，需对其进行去噪处理，而现有的卡尔曼滤波方法的精度受信号中的野值等干扰数据的影响，现有的小波去噪只对信号的高频分量进行阈值处理，忽略了MEMS陀螺仪信号中常常处于低频分量的角度随机游走和零偏不稳定性。
[0005] 本发明的目的是为解决上述问题，提供一种融合卡尔曼和小波的MEMS陀螺仪自适应抗野值去噪方法，对陀螺仪信号进行自适应抗野值卡尔曼滤波之后，继而在确定小波分解
水平的前提下，对低频分量和高频分量同时进行小波阈值处理，融合卡尔曼滤波和小波去噪的优点，更好地抑制陀螺仪信号噪声，提高陀螺仪的测量精度。
[0006] 本发明的技术方案是融合卡尔曼和小波的MEMS陀螺仪自适应抗野值去噪方法，用于对MEMS陀螺仪信号去噪，包括以下步骤，
[0007] 步骤1：建立MEMS陀螺仪信号的时间序列ARMA模型；
[0008] 步骤2：采用卡尔曼滤波对陀螺仪信号去噪；
[0009] 步骤3：利用小波分析，分别对卡尔曼滤波后的陀螺仪信号的低频分量、高频分量进行阈值处理；
[0010] 步骤4：对阈值处理后的陀螺仪信号的高频和低频信号采用matlab
软件中waverec函数进行小波重构。
[0011] 进一步的，步骤2中，所述卡尔曼滤波为自适应抗野值卡尔曼滤波。
[0012] 进一步的，步骤1中，所述建立陀螺仪信号的时间序列模型，先提取陀螺仪原始信号中的趋势项和周期项，再对其残差建立时间序列ARMA模型，采用最终预报误差FPE准则对模型定阶，求得陀螺随机漂移误差模型
[0013] mk+1＝amk+ξk+1
[0014] 式中mk,mk+1分别表示tk,tk+1时刻的陀螺仪误差，a为自回归系数，ξk+1为tk+1时刻的白噪声。
[0015] 进一步的，步骤2包括以下子步骤：
[0016] 步骤2.1：根据步骤1的陀螺随机漂移误差模型建立卡尔曼滤波器的
状态方程和量测方程；
[0017] 步骤2.2：进行自适应抗野值卡尔曼滤波器初始化，在滤波开始时刻t0，利用陀螺随机漂移误差模型估计的陀螺仪误差m，初始化滤波器的状态向量 和方差矩阵P0；
[0018] 步骤2.3：进行自适应抗野值卡尔曼滤波，在完成tk时刻的一次滤波更新之后，计算前k+1时刻信息序列ε(ε1、ε2、…εk+1)的方差 其中εk+1为tk+1时刻的信息，计算||εk+1|-|E[ε]||并与 相比较，其中，β为常数，β∈(0,1)；
[0019] 步骤2.3.1：如果 则认为此数据是干扰数据，修正tk+1时刻卡尔曼滤波增益Kk+1，将αKk+1赋值给Kk+1，α为常数，α∈(0,1)，重新计算tk+1时刻的系统状态更新 和协方差更新Pk+1|k+1以及滤波器输出 进入下一时刻滤波循环；
[0020] 步骤2.3.2：如果 则修正tk时刻系统状态向量的一步预测值 将 赋值给 γ∈(0,1)，并重新计算tk+1时刻的信息εk+1，继
而重新计算tk+1时刻的系统状态更新 和协方差更新Pk+1|k+1以及滤波器输出 进入下一时刻滤波循环；
[0021] 步骤2.3.3：如果 则进入下一时刻的滤波过程循环。
[0022] 进一步的，步骤3包括以下子步骤：
[0023] 步骤3.1：选择小波函数db2；
[0024] 步骤3.2：确定小波分解水平：对滤波后输出的陀螺仪信号采用
选定的db2小波进行小波分解，计算每层高频分量的峰值比，即
[0025]
[0026] 式中j为分解层数，Nj为第j层高频分量的个数，max(|Wj|)为位于第j层的高频分量绝对值的最大值，|Wj,i|为位于第j层i第个高频分量的绝对值，Jj为第j层的峰值比，[0027] 当Jj≤ρ且Jj+1＞ρ时，确定分解层数L＝j，其中ρ为区分某层高频分量中是否含有噪声信息的一个阈值，ρ∈(0,1)；
[0028] 步骤3.3：对L层的低频信号和每一层的高频信号进行阈值处理，采用Visushrink阈值 其中 |c1|为小波分解第一层的小波系数的绝对值，M为信号的长度，采用软阈值函数
[0029]
[0030] 式中r为每一个小波系数。
[0031] 优选地，α＝0.5，β＝0.7，γ＝0.5。
[0033] 1)本发明采用的是自适应抗野值卡尔曼滤波，相较于经典卡尔曼滤波而言，添加了抗野值环节，避免了信号中野值对滤波的影响，提高了滤波精度；
[0034] 2)本发明在小波去噪过程中在保证信号不失真的情况下，对低频信号也进行了阈值处理。相较于传统的小波去噪而言，更好的抑制了低频噪声；
[0035] 3)本发明在自适应抗野值卡尔曼滤波之后，再将信号进行小波处理，相比以往的卡尔曼滤波和小波融合的方法，本
算法更简单有效；
[0036] 4)本发明为MEMS陀螺仪信号提供一种融合卡尔曼滤波和小波的自适应抗野值去噪方案，可以更有效的提高传感器的精度，减少误差。
附图说明
[0037] 下面结合附图和
实施例对本发明作进一步说明。
[0038] 图1为本发明方法的流程示意图；
[0039] 图2为实施例的陀螺仪x轴去噪效果分析图；
[0040] 图3为实施例的陀螺仪y轴去噪效果分析图；
[0041] 图4为实施例的陀螺仪z轴去噪效果分析图；
具体实施方式
[0042] 如图1示，融合卡尔曼和小波的MEMS陀螺仪自适应抗野值去噪方法，具体包括以下步骤，
[0043] 步骤1：建立MEMS陀螺仪信号的时间序列ARMA模型：
[0044] 提取MEMS陀螺仪原始信号中的趋势项和周期项，再对其残差建立时间序列ARMA模型，采用最终预报误差FPE准则对模型定阶，求得MEMS陀螺随机漂移误差模型为AR(1)模型，表示为：
[0045] mk+1＝amk+ξk+1
[0046] 式中mk,mk+1分别表示tk时刻和tk+1时刻的陀螺仪误差，a为自回归系数，ξk+1为tk+1时刻的白噪声；
[0047] 步骤2：采用卡尔曼滤波对陀螺仪信号去噪；
[0048] 步骤2.1：根据步骤1所得的AR(1)模型建立卡尔曼滤波器的状态方程和量测方程，具体为：
[0049]
[0050] 其中，A＝a，B＝1，H＝1，Vk+1为AR(1)模型的估计误差，其方差为R，Xk+1、Xk分别为tk+1、tk时刻的状态量，Zk+1为tk+1时刻的量测量，ηk为系统过程噪声，其方差为Q，其值为步骤1所述ξk+1白噪声的方差，Q和R互不相关，将陀螺仪误差m作为滤波器输入；
[0051] 步骤2.2：自适应抗野值卡尔曼滤波器初始化，在滤波开始时刻t0，利用ARMA模型估计的陀螺仪信号误差m，初始化滤波器的状态向量 和方差矩阵P0，具体为：
[0052]
[0053] 其中E[·]表示均值，(·)T为矩阵的转置；
[0054] 步骤2.3：进行自适应抗野值卡尔曼滤波；
[0055] 步骤3：确定小波基db2；
[0056] 确定小波分解水平：对滤波后输出的陀螺仪信号采用选定的db2小波进行小波分解，计算每层高频分量的峰值比，即 其中，j为分解层数，Nj为第j层高频分量的个数，max(|Wj|)为位于第j层的高频分量的绝对值的最大值，|Wj,i|为位于第j层第i个高频分量的绝对值，Jj为第j层的峰值比；当Jj≤ρ且Jj+1＞ρ时，确定分解层数L＝j，其中ρ为区分某层高频分量中是否含有噪声信息的一个阈值，ρ∈(0,1)；
[0057] 对L层的低频信号和每一层的高频信号进行阈值处理，采用Visushrink阈值其中 |c1|为小波分解第一层的小波系数的绝对值，M为信号的长度。实施例中采用软阈值函数, 其中，r为每一个小
波系数。
[0058] 步骤4：对处理后的高频和低频信号采用matlab软件中waverec函数进行小波重构。
[0059] 步骤2中，自适应抗野值卡尔曼滤波的滤波过程具体为：
[0060] 第一步，状态一步预测：
[0061]
[0062] 第二步，协方差阵一步预测：
[0063] Pk+1|k＝APk|kAT+BQBT
[0064] 第三步，计算滤波增益矩阵：
[0065] Kk+1＝Pk+1|kHT[HPk+1|kHT+R]-1
[0066] 第四步，计算信息：
[0067]
[0068] 第五步，状态更新：
[0069]
[0070] 第六步，协方差更新：
[0071] Pk+1|k+1＝[In-Kk+1|kH]Pk+1|k
[0072] 第七步，k+1时刻的滤波器输出：
[0073]
[0074] 式中 为tk时刻系统状态向量的一步预测值， 分别为tk，tk+1时刻的系统状态向量的时间更新。Pk+1|k为tk时刻的系统协方差矩阵的一步预测值，εk+1为tk+1时刻的信息，Pk|k，Pk+1|k+1分别为tk，tk+1时刻的系统协方差更新，Kk+1为tk+1时刻的卡尔曼滤波增益，In为n阶单位矩阵，n为滤波器状态向量的维数，[·]-1为矩阵的逆矩阵， 为k+1时刻的滤波器输出；
[0075] 其中的自适应抗野值环节：在完成tk时刻的一次滤波更新之后，计算前k+1时刻信息序列ε(ε1、ε2、…εk+1)的方差Sk+1：
[0076]
[0077] 并比较：
[0078]
[0079] 式中的|·|表示为绝对值，若满足式(1)，则认为此数据是干扰数据，将αKk+1赋值给Kk+1，α∈(0,1)实施例取α＝0.5，返回第五步重新进行状态更新；
[0080] 若不满足(1)式，继续比较是否满足:
[0081]
[0082] 实施例中取β＝0.7，若满足式(2)，将 赋值给 γ∈(0,1)，实施例中取γ＝0.5,返回第四步重新计算信息；若既不满足式(1)又不满足式(2)，则进入下一时刻的滤波过程循环。
[0083] 通过实验加以验证，实验内容如下：
[0084] 1.本试验数据来源于某型MEMS陀螺仪静
基座实测数据，静置时间为1小时，
采样频率200HZ。
[0085] 2.根据步骤1，陀螺仪x、y、z轴的自回归系数a分别为0.0092、0.0299、0.0080，其对应的白噪声的方差分别为3.8595×10-7(rad2/s2)、4.6831×10-7(rad2/s2)、3.6487×10-7(rad2/s2)。
[0086] 3.陀螺仪x、y、z轴三轴的卡尔曼滤波器初始条件分别设为：
[0087] P0x＝3.8598×10-7(rad2/s2)
[0088] P0y＝4.6873×10-7(rad2/s2)
[0089] P0z＝2.4301×10-7(rad2/s2)
[0090] 其中 分别为x、y、z轴滤波器的初始状态向量，P0x、P0y、P0z分别为x、y、z轴滤波器的初始协方差。
[0091] 4.根据步骤3，针对此陀螺仪信号，当陀螺仪误差信号的峰值比，Jj≤ρ,Jj+1＞ρ,ρ＝0.05时，确定x、y、z三轴数据的分解层数分别为Lx＝Ly＝Lz＝13。
[0092] 5.分别对陀螺仪x、y、z三轴数据进行经典卡尔曼滤波、Visushrink去噪、本方案去噪，并通过allan方差估计陀螺仪噪声系数进行方法对比。
[0093] 通过计算机仿真，陀螺仪x、y、z轴的去噪效果对比图如图2、图3、图4所示，陀螺仪x、y、z轴的噪声系数对比分别如表1、表2、表3所示。
[0094] 表1陀螺仪x轴噪声系数对比分析表
[0095]
[0096] 表2陀螺仪y轴噪声系数对比分析表
[0097]
[0098] 表3陀螺仪z轴噪声系数对比分析表
[0099]
[0100] 从表1-3中可以看出使用本方案的方法量化噪声系数、角度随机游走和零偏不稳定性均小于卡尔曼滤波和小波去噪，尤其陀螺仪x、y、z轴零偏不稳定性在本方案下比卡尔曼滤波分别提高了31.0％、29.3％、30.5％、比传统小波分别提高了2.4％、12.1％、12.4％。
[0101] 陀螺仪x、y、z轴去噪后方差对比如表4所示。从表4中可以看出，本方案去噪后方差均小于卡尔曼滤波和小波去噪。
[0102] 表4陀螺仪x、y、z轴方差对比分析表
[0103]
[0104] 上述实施例仅是本发明的优选实施方式，应当指出：对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和等同替换，这些对本发明
权利要求进行改进和等同替换后的技术方案，均落入本发明的保护范围。
