基于奇异谱分解的微机械加速度计振动噪声抑制方法

伍宗伟， 姚敏立， 马红光， 马帮立， 田方浩

(1. 第二炮兵工程大学 空间工程系，西安 710025;2. 郧阳师专物理与电子工程系，十堰 442700)

低成本姿态估计系统广泛应用于导航控制、光电和惯性稳定平台、机器人，以及无人机等系统中[1-5]，而微机械(Micro-Electro Mechanical Systems, MEMS)惯性器件的应用有效地降低了由于高精度惯性器件所带来的昂贵成本。微机械惯性器件的成本较低，但是其性能较差，不仅测量噪声与零偏变化较大，且振动也对其造成较大的影响。因此，在实际实用中需要采用数据融合方法来克服使用单一器件带来的较大误差。文献[3-5]代表了工程上常用的低成本姿态估计方法，即利用陀螺的短时精度和加速度计的长时精度通过滤波算法融合得到较好的姿态估计，去除陀螺的零偏以及加速度计较大的高频测量噪声，算法简单，但是需要克服两个问题：①MEMS惯性器件受到振动噪声影响很大，虽然卡尔曼滤波算法一定程度上可以认为是低通滤波器，但是仍然无法去除叠加在有用姿态信号上的低频振动噪声；②由于加速度计依靠重力矢量测量得到姿态角，因此易受载体的机动影响。因此，在融合滤波器之前需要采用预处理方法来克服这些误差和干扰。第二个问题已经有较多应对方法，如判断开关以及自适应滤波方法[4-5]，因此如何处理低频振动噪声是提高姿态估计精度的关键。

针对MEMS惯性器件的噪声预处理，主要方法有小波分析以及经验模态分析(Empirical Mode Decomposition，EMD)两种。文献[6]代表了基于小波阈值方法去除高频噪声的方法，并且利用机动性判断来确定小波分解层数，一定程度上克服了小波阈值去噪分解层数难以确定的问题。但是，该方法并没有能够有效地去除与有用姿态信号频域叠加的低频振动噪声，且在一定意义上与后续的融合算法作用重复，达不到提高姿态融合精度的目的，反而增加了运算量。文献[7]指出了小波分解的不足，即合适的小波基难以选取以及小波分解层数难以确定，提出了采用EMD的方法抑制陀螺的随机漂移，EMD是完全数据驱动的一种自适应数据分解方法，但是对于MEMS惯性器件来说，噪声较大且噪声的不确定性很强，分解出来的虚假模式很多，并且各个固有模态之间信息混叠，且姿态变化与振动噪声并没有显著特性，很难提取出真实的振动模式和有用姿态信号，同时，EMD自身存在端点效应等问题，限制了该方法在MEMS惯性传感器噪声处理中的应用。另外，文献[8]提出结合小波和EMD来进行去噪，不仅算法复杂，而且没有从原理上克服二者的缺陷。

本文以车载平台的水平姿态估计为研究对象，提出采用奇异谱分析(Singular Spectrum Analysis, SSA)的方法抑制对姿态融合影响较大的振动噪声，该方法能够有效地提取姿态变化的趋势项，分离与振动噪声有关的谐振项。跑车实验的结果表明，SSA算法能够有效地的分离有用信号和振动噪声，提高低成本姿态估计的精度。

1 振动噪声分析

1.1 姿态估计原理

分析噪声对低成本姿态估计影响的前提是对姿态估计原理与融合滤波器的分析。首先简要介绍姿态估计原理，陀螺积分可以得到姿态角，而加速度计的重力测量分量也能够得到姿态角，定义坐标系为：地理坐标系取为东北天坐标系，车体坐标系的x,y,z轴分别指向车右侧方向、车头方向以及车体正上方。在无机动加速度干扰的情况下，加速度计测量姿态角的公式为[3]：

(1)

式中g为重力加速度值，fx、fy分别为两个水平放置的加速度计测量值，n为测量噪声。对于微机械加速度来说，受到振动噪声影响较大。而陀螺的积分能够提供较为精确的姿态角，但是由于零偏的存在使得姿态角存在长期误差积累，因此将二者融合即可得到不受振动干扰和时间累积误差的精确姿态角。通常选择Kalman滤波器作为融合滤波器，将陀螺零偏扩增为状态变量进行估计，选择合适的噪声阵，能够有效抑制加速度计的测量噪声n。由Kalman滤波器的性质可知，选择较大的测量噪声阵R，可以抑制加速度计的测量噪声n，相当于一个低通滤波器，但是如果选择的R阵过大，则测量值对系统状态方程的修正量过小，会造成滤波器的扰动，图1给出了一组在实验中的估计结果，R阵中的测量噪声方差分别选择为5×10-3与5×10-1的比较结果(实际的测量噪声方差为5×10-4，但是在Kalman滤波中一般适当放大该值以取得较好的滤除高频误差效果)。

图1 测量噪声方程阵对Kalman滤波器的影响

从图中可以看出，较大的噪声阵能够有助于抑制振动噪声，但是会对Kalman滤波器本身带来扰动，实际的应用中选择的R阵要兼顾精度和抗扰动，因此，低频的振动噪声一般难以通过Kalman滤波器得到较好的抑制。

1.2 振动噪声

车辆的振动噪声主要来自于三个方面的激励：①路面不平度；②车辆发动机；③车轮/轮胎(内部激励，主要由于车轮不平衡以及车轮冲击和轮胎质量不均匀导致)[9]。因此，为了研究载体振动噪声，分别采样载体不同运动状态下的传感器测量数据进行频谱分析。MEMS陀螺与加速度计均来源于星网宇达公司的微机械惯性测量单元XW-ADU5220，分析结果如下图2所示。

从频谱分析中我们可看出：载体的振动噪声对MEMS陀螺影响较小，变化较为明显是第五种情况下，即土路的颠簸不平对姿态角的影响较大，从而对陀螺的测量幅值产生较大的影响；而加速度计则明显受到载体振动噪声影响，发动机的振动噪声以及路面不平度引起的噪声在全频段均体现出来，静止且发动机振动情况下，有明显的振动峰值谱以及谐振谱，而车辆在低速平坦路面行驶时同样有频率大致在16 Hz和40 Hz附近的振动噪声，且随着车速以及路面颠簸度增强，该16 Hz附近的振动有加强趋势，低频段的振动噪声也同样加强。

图2 载体不同运动状态下的测量数据频谱分析

图3 低通之后的测量数据频谱分析

首先想到的是通过低通滤波器滤除高频噪声。通过实际的行车实验中参考系统提供的高精度姿态测量，可知车体的姿态变化一般低于5 Hz，因此首先设计低通滤波器滤除大部分的高频噪声。图3显示了通过一个5 Hz的低通滤波器之后的加速度计的测量值与实际通过高精度惯导计算出的加速度计值的比较，选择的是慢速平坦的测量环境。

通过上图我们看到，经过低通滤波之后加速度计仍然含有较大的低频噪声，这个噪声分量主要是由于车体振动引起的。通过大量的跑车实验，发现该低频振动噪声的最大强度在1 Hz左右，实验中采用的是全顺面包车，测量设备均固连于车顶，对于不同的载体与安装位置，均可进行简单的频谱分析得到该主要低频振动噪声的集中频点。

2 奇异谱分析

奇异谱分析的目的是将原始的时间序列分解为一些序列之和，使得这些序列能够分离成为趋势项、周期项以及其他的非周期项，或者噪声，然后选择合适的序列重构，从而达到处理信号的目的。实现的手段类似于经验模态分解或者小波分解，方法类似于主成分分析。由于该方法本质上是一种利用协同延迟矩阵进行主成分分析的方法，因此可以看作是一种从能量域分析信号的方法。该方法不同于小波分析需要所分析信号的先验知识，含有经验模态分解数据驱动的优点，同时又很适合短数据的趋势提取与谐波分量和噪声分量的分离，已经成功的应用于气象与地理学、电信以及金融等领域[10]。将加速度计的测量分为三个方面：①有用的姿态变化信号，即趋势项；②振动引起的谐波分量；③测量噪声。通过奇异谱分解可以分别对上述的三种分量进行基于数据的自适应提取，其基本算法分为两个互补的步骤：分解和重构。具体的实现算法为[11]：

(1) 插值构造

插值构造实际上是一种将单维数据映射到多维的方法。假设我们有一个一维时间序列构成的信号矢量X=[x1,…,xN]T，我们将其扩张为多维序列信号矩阵H=[X1,…,Xk]，其中的每一个Xi=[xi,…,xi+L-1]T，K=N-L-1，都是一个所谓延迟信号矢量，这里的时间窗口长度L是一个整数，其取值范围为2

(2)

从上面的构成过程看，实际上构成的H阵为一个Hankel阵，且构成过程就是一个扩维过程，扩张的维数由窗口长度L确定，SSA根据这个轨线矩阵H来确定主要分量，后续的分析实际上就是找到H的主分量以及每一个主分量对应的子空间，然后选择合适的信号子空间，对轨线矩阵进行投影。

(2) 奇异值分解

计算HHT的奇异值，实际上和PCA类似，求取其特征值λ1,λ2,…λL。由于HHT为对称矩阵，则有HHT=VΛVT，这里的Λ为以λi为对角元素的对角阵，且λi满足λ1≥λ2≥…≥λL>0，令每一个λi对应的分量为Hi，我们得到

(3)

(3) 分类

将上述分解得到的一些初等矩阵H进行分组，令I1={i1,…,ip}是一个I1分组的下标i1,…,ip集合，则组I就定义为HI+H…Hip，总共分为s组，则原来的H阵就变为：HI=HI1+…+HIs。对于给定的组，其占用的权重为:

一般情况下，都是直接分为两类，一类为噪声，一类为趋势项，即s=2。

(4) 对角平均

所谓对角平均即是将上述给定的分组I中的每一个初等矩阵Hi对应的一维信号恢复出来，也就是所谓的Hankel化。假定给定的任意矩阵为Z，其元素为zij，首先将其变换为Hankel矩阵形式Z，然后从Hankel矩阵求序列，也就是说第k个数值是将所有满足i+j=k+2的HZ的元素求平均即可。这里的H算子其定义如下：对于任意的L×K矩阵Z=(zij)，L≤K，对于i+j=c，N=L+K-1，Z的元素可以通过下面的算法进行计算:

(4)

通过SSA算法可以有效地分离出趋势项和谐振项，谐振项可以认为是由于振动噪声引起的，而趋势项则反映了姿态的实际变化，由于Kalman滤波器融合陀螺积分和加速度计测量结果得到姿态估计，因此，这里的加速度计利用的是其长期精度，故取该趋势项足以满足Kalman滤波的要求，同时还能够去除叠加在有用姿态信号上的低频振动噪声。

注意这里的趋势项并不一定是SSA分析得到的分离信号中的第一个谱分量，由于SSA得到的信号序列是按照能量大小进行排列的，对于MEMS加速度计来说，很多情况下都是振动噪声的能量要远大于有用的姿态变化能量，因此，这里需要采用一个合适的方法区分出振动项和趋势项，一个简单的准则就是过零点检测法，由于姿态变化是慢变的趋势项，因此过零点是最少的。较为复杂和精确的方法是采用FFT进行频域分离，在零频处频谱最大值即为趋势项。

另外，在SSA实现过程中，仅有一个参数需要设定，即SSA算法中的时间窗口L，根据轨线矩阵H的构成来看，每一列都是L个采样点，因此，该参数应该与需要SSA分析的信号中的主要周期项的周期大致可比拟才能得到较好的谐振项与趋势项的分离。根据实验中测量的主要振动项大致为1 Hz，而采样频率为100 Hz，故L设定为100。

3 跑车实验数据分析

实验设备包括星网宇达公司的微机械惯性测量单元XW-IMU5220，对比的参考系统为XW-ADU7612。其中XW-IMU5220包含有3个MEMS陀螺和3个MEMS加速度计，陀螺的零偏<0.08°/s，零偏稳定性为<0.05°/s，而加速度计的零偏<0.005g，零偏稳定性为<0.001g。XW-ADU7612能够以100Hz的数据更新速率输出被测对象精确的姿态角(航向、俯仰和横滚角)，精度可达0.1°。如图4所示，实验设备安装于全顺车顶。

图4 实验载体及设备安装照片

实验的行车环境分别为图1中的后两种运动环境，即城市公路以及乡村土路，行车速度大致为30km/h。

(1)城市公路

首先分析城市公路的数据处理情况。图5给出了车辆近似直行情况下加速度计的测量值与高精度参考得到的数据对比结果。从图5中可以看出，加速度计的测量值存在较大的低频振动噪声项，但是其趋势项与参考值大致能够保持一致，因此，如果能够较好的分离出趋势项与低频振动噪声项，则可以在最后的融合估计中对振动噪声项进行抑制。下面我们通过SSA对横滚角测量数据进行详细分析，分离出前6个能量依次递减的谱分量如图6所示。

从图6中可以明显看出趋势项与分离的振动项，对比图6和图5可以看出，SSA方法在提取趋势项方面是很有效的，另外各个谐振项也能够分离出来，由于低通之后的测量数据低频振动较强，因此分离结果中的2、3，以及4、5分量大体一致。通过简单的过零点检测方法即可提取出趋势项，趋势项与参考值对比如图7所示，同时给出主振动项，对比图7中的提取结果与图5原始的测量数据，更进一步说明了SSA方法在趋势项提取以及振动项分离方面的优势。

图5 加速度计的测量数据

图6 横滚角加速度计SSA分离的信号分量

图7 横滚角加速度计通过SSA提取的趋势项与振动项

图8 俯仰角加速度计通过SSA提取的趋势项与振动项

同理，我们给出俯仰角加速度计的SSA处理结果，如图8所示。

利用上述提取出的趋势项进行Kalman滤波融合得到姿态估计，具体的计算过程见参考文献[3,5]，测量方差选择为5×10-2，得到的姿态估计结果对比如图9所示，图中给出了采用SSA抑制振动噪声预处理之后的姿态估计误差与直接进行滤波融合的估计误差结果，表明了该方法能够提高低成本姿态估计的精度。

图9 姿态估计结果对比

(2) 乡村土路

下面给出在乡村土路的一组测量数据处理结果。由于乡村土路上颠簸环境，因此加速度计测量结果中振动噪声更大，因此，采用SSA预处理的效果将体现的更为明显。图10给出了横滚加速度计原始测量值、SSA预处理之后的测量值以及参考值的对比。

图10 乡村土路SSA-ICA预处理结果

利用SSA方法预处理之后得到的姿态估计误差与未处理的估计误差对比如下图11所示，这里由于是在颠簸土路，因此测量方差选择为10×10-2。

图11 乡村土路姿态估计结果对比

从两组不同路况下的实验结果来看，提出的振动噪声抑制方法均能够有效地提高姿态估计精度。

4 结 论

为了提高低成本MEMS惯性传感器估计姿态的精度，抑制载体振动噪声对MEMS加速度计带来的影响，本文提出了采用SSA方法来去除惯性传感器中的低频振动噪声。利用SSA对趋势项以及周期项能够有效分离的优势，抑制低频振动噪声，实验结果表明，提出的传感器噪声预处理方法能够有效地抑制低频振动噪声，提高通过姿态融合算法得到的姿态角的估计精度。同时，该方法有利于提高微机械惯性器件的工程实用性。

参 考 文 献

[1]Maeder U, Morari M. Attitude estimation for vehicles with partial inertial measurement [J]. IEEE Transactions on Vehicular Technology, 2011, 60(4): 1496-1504.

[2]Lai Y C, Jan S S. Attitude estimation based on fusion of gyroscopes and single antenna GPS for small UAVs under the influence of vibration [J]. GPS Solutions, 2011, 15(1): 67-77.

[3]Wu Z W, Yao M L, Ma H G, et al. Low-cost attitude estimation with MIMU and two-antenna GPS for Satcom-on-the-move[J]. GPS Solutions, 2013, 17(1): 75-87.

[4]Suh Y S. Orientation estimation using a quaternion-based indirect Kalman filter with adaptive estimation of external acceleration [J]. IEEE Transactions on Instrument and Measurement, 2010, 59(12): 3296-3305.

[5]伍宗伟, 姚敏立, 马红光, 等. 移动卫星通信低成本多传感器融合姿态估计方法 [J]. 西安交通大学学报, 2012, 46(12): 55-61.

WU Zong-wei, YAO Min-li, MA Hong-guang, et al. Low-cost attitude estimation method using multi-sensor fusion for mobile satellite communication [J]. Journal of Xi’an Jiaotong University, 2012, 46(12): 55-61.

[6]沈晓卫, 姚敏立, 常瑞花. 模糊逻辑提升小波在惯性传感器去噪中的应用 [J]. 振动与冲击, 2011, 30(12): 238-242.

SHEN Xiao-wei, YAO Min-li, CHANG Rui-hua. Application of fuzzy logic and lifted wavelet transformation in de-noising of inertial sensors [J]. Journal of vibration and shock, 2011, 30(12): 238-242.

[7]甘 雨，隋立芬. 基于EMD-小波随机消噪模型的GPS /INS 组合导航[J]. 测绘学报, 2011, 40(6): 745-750.

GAN Yu, SUI Li-fen. De-noising method for gyro signal based on EMD [J]. Acta Geodaetica et Cartographica Sinica, 2011, 40(6): 745-750.

[8]王 坚, 李增科, 高井祥, 等.基于EMD-小波随机消噪模型的GPS/INS组合导航 [J]. 东南大学学报: 自然科学版, 2012, 42(3): 406-412.

WANG Jian, LI Zeng-ke, GAO Jing-xiang, et al. EMD-wavelet based stochastic error reducing model for GPS /INS integrated navigation [J]. Journal of Southeast University: Natural Science Edition, 2012, 42(3): 406-412.

[9]Mitschke M, Wallentowitz H. 汽车动力学 [M]. 第四版. 陈荫三，余 强，译. 北京：清华大学出版社，2009.

[10]Harris T J, Yuan H. Filtering and frequency interpretations of singular spectrum analysis [J]. Physica D, 2010, 239: 1958-1967.

[11]Ma H G, Jiang Q B, Liu Z Q, et al. A novel blind source separation method for single-channel signal [J]. Signal Processing, 2010, 90: 3232-3241.