基于经验模分解的陀螺信号消噪

甘 雨，隋立芬

信息工程大学测绘学院，河南郑州450052

基于经验模分解的陀螺信号消噪

甘 雨，隋立芬

信息工程大学测绘学院，河南郑州450052

陀螺随机漂移是影响惯性导航精度的重要因素。小波消噪方法对异常噪声效果不明显，且对小波基和分解尺度等因素依赖性较强。提出陀螺信号经验模分解（EMD）消噪方法，将信号进行经验模分解得到一个本征模态函数（IMF）组，先基于2σ准则处理异常噪声IMF分量，再利用相关系数确定高频噪声IMF分量个数，将噪声分量去除以实现陀螺信号消噪。详细对比小波方法与EMD方法，利用交叠式Allan方差分析两者的消噪效果，通过惯导算例进一步验证EMD方法的实效性。结果表明，相比小波方法，EMD消噪法能剔除异常噪声，可以更有效地抑制陀螺漂移。

陀螺随机漂移；小波；经验模分解；消噪

1 引 言

陀螺随机漂移是影响惯性导航精度的重要因素。抑制陀螺漂移的方法主要有两种［1］：① 建立漂移模型，使用Kalman滤波等方法进行补偿；②对陀螺输出信号进行消噪处理。由于随机漂移往往表现为弱非线性、非平稳、慢时变，且易受到外部环境等多种不确定因素的影响［2］，无法建立其准确的系统模型，故需要采用陀螺信号消噪的方法。

目前对陀螺信号消噪主要采用小波方法。小波具有优良的多分辨率分析特性，小波消噪不需要系统的误差模型，因此被广泛用于陀螺信号的消噪处理中［1－4］。然而，小波分解虽然也能实现对非平稳信号的滤波，但其实质是带通滤波器，限制了滤波的精确性，且小波变换中小波基一经选定，整个信号分析过程中就只能使用这一个小波基，即小波变换是非适应性的［5］。有研究表明［6］：对于如陀螺信号一类的非平稳信号，当异常噪声淹没了有用信号时，采用小波消噪效果也不甚理想，对这类信号的消噪目前还未发现较好的方法。

文献［7］提出一种分析非平稳、非线性信号的自适应分解方法：经验模分解方法（empirical mode decomposition，EMD），它将复杂的信号分解成若干个按频率高低排列的本征模态函数（intrinsic mode function，IMF），每个IMF是一个单分量信号。该方法与小波分析的区别在于它不需要事先选定基函数，而是根据信号本身的特性自适应地产生合适的模态函数，这些模态函数能很好地反映信号在任何时间局部的频率特征。该方法已经用于机械振动信号分析、SAR影像滤波、气象以及GPS信号处理等领域［8－11］。

将EMD方法引入到陀螺信号的处理中，给出陀螺EMD消噪的方法，按照2σ准则剔除振幅偏大的异常噪声，利用相关系数确定随机高频噪声IMF分量的个数。与小波消噪依赖小波基、分解尺度、阈值估计方法不同，EMD消噪过程完全依赖陀螺信号本身特性。用交叠式Allan方差对比分析本文方法与小波消噪方法，通过惯性导航算例进一步验证方法的实效性。

2 EMD消噪方法

2.1 EMD基本原理

EMD分解方法认为任何待分解信号都由一组固有振动模式构成，并据此将信号分解为若干本征模态函数IMF的和。这些本征模态函数既可以是线性的，也可以是非线性的；既可以是平稳的，也可以是非平稳的［11］。分解得到的IMF分量满足：①零点数目与极值点数目相同或至多相差1；②函数由局部极大值点构成的包络线和由局部极小值构成的包络线的均值为零。分解过程通过一个称为“筛选”的步骤来完成。“筛选”过程可以表示为［7］：

（1）分别由原始信号x（t）的极小值点与极大值点用三次样条插值得到x（t）的上下包络线，计算上下包络的均值m1，进而计算x（t）和m1的差值h1＝x（t）－m1，判断h1是否满足IMF的两个条件。若满足，则h1为x（t）的第一个分量imf1。

（2）若不满足，则将h1作为新的信号继续步骤（1），得到h11，判断h11是否满足条件。若不满足则继续（1）的步骤，直到重复k次后h1k满足IMF的条件，则有imf1＝h1k，求出原始信号与imf1的差值r1＝x（t）－imf1。

（3）将r1作为新的“原始”信号重复上述步骤，直到提取出第2个，第3个，直至第n个IMF分量imf2，imf3，…，imfn，有r2＝r1－imf2，…，rn＝rn－1－imfn，当IMF分量imfn或余项rn小于预先设定的值，或者余项rn已经成为单调函数时，整个筛选分解过程结束，本文以余项为单调函数作为分解的终止条件。

经过上述步骤后，原始信号x（t）可分解为n个IMF分量和1个余项的和

从上述的经验模态分解方法可以看出，越早分解出来的IMF频率越高［7］，第一个分解出来的代表原信号的最高频率成份，各个IMF的频率几乎是按2的负幂次方的形式递减［6，12］。对于混有随机噪声的信号，其先分解出的IMF分量通常对应于信号的高频噪声［8］。若IMF组去除了先分解的几个IMF，把其余的IMF组合起来形成一个信号，可以削弱信号的噪声。可以看出，用EMD对陀螺信号消噪的关键在于从信号分解出的IMF组中辨识出噪声IMF分量。

2.2 处理异常噪声

由于外界环境的异常变化，陀螺信号会在局部时段内受到振幅较大的异常噪声的干扰，导致信号输出失真，对陀螺信号的处理首先要消除这类噪声的影响，而小波对这种噪声的消噪效果不佳。经验模分解将信号按照其固有振动模式分成IMF组［13］，而异常噪声的频率和振动模式一般既不同于随机高频噪声，也不同于真实信号，对信号进行分解之后，异常噪声一般被分解到其中一个或几个IMF之中，只要能将相应的IMF进行识别并加以处理，就能减弱异常噪声的影响。

利用误差理论中制定极限误差的2σ准则，可以对陀螺信号的异常噪声IMF分量进行识别。设原始信号为x（t），求得其采样标准差为σ，计算IMF各个分量imfi的最大振幅Ai，对各分量按如下判断准则进行处理

异常噪声作用的时间相对比较短，σ通常反映陀螺信号的整体观测精度，若某个imfi的最大振幅大于2σ，则可认为它是异常噪声分量。IMF分量的频率各不相同，且如前面所述按2的负幂次方递减，因而异常噪声通常会被分解到少数IMF分量且是中高频分量，非常便于进行识别和剔除。

判别和处理异常噪声是陀螺EMD消噪方法的第一步，除非陀螺的观测环境非常好，否则都应该在对信号进行EMD分解后进行异常噪声的辨识和剔除。

图1给出某次试验中陀螺（标称陀螺漂移为1°／h）静基座下X轴8 000个历元的输出信号，在1 500～2 000、4 100～4 300、5 600～6 900历元处受到异常噪声的干扰。对该信号进行EMD分解，得到10个IMF分量，图2显示为前6个IMF以及最终的余项r10。可直观看出，1 500～2 000、4 100～4 300历元的异常噪声主要分解在imf3内，5 600～6 900历元的分量分解到imf2和imf3内。

计算陀螺信号标准差为σ＝34.02，则限差2σ＝68.04，各IMF分量的最大振幅分别为A1＝10.66，A2＝216.38，A3＝87.51，A4＝39.85，A5＝42.96，A6＝12.28。按照式（2），imf2和 imf3的振幅超限，应该将这两个分量置零，这与图形所示的结果一致。

图1 陀螺原始信号Fig.1 Original gyro signal

图2 陀螺信号EMD分解的IMF分量Fig.2 IMFs from EMD of gyro signal

2.3 确定噪声IMF分量个数

剔除异常噪声后，需要消除陀螺信号中的高频噪声。设前m级IMF分量（已经过异常噪声的处理）为高频噪声分量，则消噪过程可表示为

若选取的m偏小，滤波后的信号x′（t）仍会受到噪声干扰，若选取的m偏大，可能会把有用的信号滤去。为避免此问题，有文献借鉴小波阈值消噪的思想提出对每个IMF分量也进行阈值处理［14］。然而，每个IMF分量都是具有特定物理意义的信号成分，反映原始信号某一方面的内在特征［7，13］，利用阈值消噪的方法可能会破坏信号特征的完整性，且阈值消噪法存在确定阈值估计方法的难点，其实用性有待进一步研究。EMD消噪的关键还是在于确定噪声分量的个数m。

长度为N的信号序列x（t）和y（t）的采样相关系数可定义为

相关系数反映两个信号之间的相似程度或依赖程度。设前m个IMF分量之和为若这m个分量都是噪声IMF分量，不含原信号的有用成分，则按照式（3）滤波后的信号x′m（t）与原信号x（t）之间应该具有较强的依赖程度。因此可考虑由相关系数确定噪声分量的个数，计算

从m＝1开始计算，当ρxx′m＞c（c为常量，可取经验值0.75～0.8）时，m＝m＋1，计算下一个ρxx′m，直到某个ρxx′m＜c，停止计算，说明第m个分量imfm为信号的有用成分，不是噪声IMF，取M＝m－1作为噪声分量个数。最终的EMD消噪结果为

陀螺信号的EMD消噪方法归纳为：

（1）对陀螺信号进行EMD分解，得到IMF分量；

（2）按2σ准则处理异常噪声IMF；

（3）根据相关系数确定噪声IMF个数，原信号减去噪声IMF完成消噪过程。

以上的EMD过程无需事先的参数设置，只与陀螺信号本身有关。陀螺信号的EMD分解是按照信号的固有振动模式进行的，采样标准差σ和相关系数ρxx′m都由陀螺信号及其分解的IMF分量来计算，整个消噪过程不受外在因素影响。需要注意的是，常量c的取值及2σ准则的确定是EMD消噪中的重要环节，如果选取不合适可能会影响最终的消噪结果，因此实用中需要结合实际的数据类型进行分析确定。

确定图2中噪声分量个数，从m＝4开始，相关系数均小于0.7，则取噪声分量个数为3，EMD消噪的结果如图3所示。

图3 EMD消噪信号Fig.3 EMD de－noised signal

3 消噪方法对比

3.1 消噪信号直接比较

为对比小波消噪法与本文EMD消噪法，对图1中的陀螺原始信号进行小波阈值消噪，采用软阈值函数，阈值估计方法为SUREShrink阈值［3，15］。分别取不同小波基与分解尺度进行消噪试验，消噪的效果总体一致，图4、图5分别为db8小波9尺度消噪结果和Haar小波7尺度消噪结果。

图4 db8小波消噪信号Fig.4 db8wavelet de－noised signal

图5 Haar小波消噪信号Fig.5 Haar wavelet de－noised signal

比较图1、图3、图4和图5，可以看出：

（1）EMD利用异常噪声与有用信号和高频噪声不同的内在振动模态将异常噪声分离，因此EMD消噪法既能消去随机高频噪声，又能剔除异常噪声；

（2）小波消噪方法能有效地去除高频噪声，但是不能抵制异常噪声的干扰；

（3）小波分析的本质决定了它不能分离出异常噪声，通过选取不同的小波基和分解尺度也无法克服这一缺陷。

3.2 交叠式Allan方差对比分析

采用交叠式Allan方差进一步对比两种方法的消噪效果。Allan方差适合于分析非平稳的随机信号，是IEEE推荐的陀螺随机误差分析方法［16］，其中的交叠式Allan方差在相同的置信水平下比普通Allan方差分析方法具有更大的置信区间，是对Allan方差的改进［17］，其具体定义见文献［17—18］。

信号的Allan方差可表示为σ2（τ）（τ＝nτ0，τ0为采样间隔），则Allan标准差为σ（τ）～τ双对数曲线图可以描述陀螺的各种随机误差成分，不同的成分具有不同的斜率特性，这些成分包括量化噪声，角度随机游走，零偏不稳定性，速率随机游走和速率斜坡，可由Allan方差拟合得到，具体计算方法和它们在σ（τ）～τ图中的斜率特性可参考文献［18］。

分别求取陀螺原始信号、小波消噪信号、EMD消噪信号的交叠式Allan标准差σ（τ），σ（τ）～τ双对数曲线如图6所示，其中db8和Haar小波的结果接近，只显示db8小波结果。根据计算的Allan方差拟合得到的各种误差成分的值如表1所示。

图6 σ（τ）～τ曲线对比Fig.6 Comparison ofσ（τ）～τcurve

表1 陀螺误差成分对比Tab.1 Comparison of error components of gyro

分析图6和表1的结果可知：

（1）小波对陀螺信号具有一定的消噪作用，但是由于受到异常噪声的干扰，小波消噪后的陀螺信号中仍然含有较大的误差成分；

（2）EMD消噪方法有效地削弱了陀螺各种误差成分；

（3）在τ＝0.2～0.3s处，EMD及小波消噪后信号的Allan标准差与原始信号相比都没有明显改善，是因为受到残余低频有色噪声的影响。一般可通过建立时间序列模型的方法削弱陀螺信号的有色噪声成分。

4 惯导解算与分析

所用数据为一组静态惯性测量单元（inertial measurement unit，IMU）数据，标称陀螺漂移和加速度计偏置分别为1°／h和10－4g，IMU采样频率100Hz，取历元282 930.83～283 232.82s（GPST，对应的周数为1 467）的数据进行试验。利用282 930.83～283 112.83s的数据进行Kalman滤波精对准得到初始姿态角，初始航向角、俯仰角和翻滚角分别为8.916°、－0.013°、1.094°，分别采取3种方案对283 112.83～283 232.82s的IMU数据进行惯导解算，得到各历元的速度值。各方案的加速度计数据均为原始输出，但陀螺数据的处理方法不同。由于静止状态下惯导的速度真值实际上是“零”，因此，计算的速度值其实就是速度误差。采用如下3种方案进行解算：

方案1 使用IMU输出的各轴陀螺原始数据解算；

方案2 使用db8小波消噪后的各轴陀螺数据进行解算；

方案3 使用本文EMD方法消噪后的各轴陀螺数据进行解算。

3种方案的速度误差见图7～图9，东向和北向的误差结果类似，这里给出的是东向的结果，RMS的比较如表2所示。

图7 方案1东向速度误差Fig.7 Veerror of scheme 1

图8 方案2东向速度误差Fig.8 Veerror of scheme 2

图9 方案3东向速度误差Fig.9 Veerror of scheme 3

表2 3种方案RMS比较Tab.2 Comparison of RMS for three schemes（m／s）

分析解算结果，可知：

（1）受到加速度计误差和残余陀螺漂移的影响，单独惯导解算的误差随着历元数增加而不断积累，3种方案的误差均呈现扩大趋势，需要引入其他的导航信息（如卫星导航）才能予以修正。

（2）方案2计算的速度误差明显小于方案1的结果，说明小波具有削弱陀螺漂移、减小惯导误差的能力。

（3）用EMD方法不仅削弱了陀螺漂移中随机高频噪声的影响，而且消除了异常噪声的干扰，相比小波又进一步减小了惯导误差。

5 结束语

小波消噪方法能在一定程度上抑制陀螺漂移，但是它有一系列固有缺陷：不能消除异常噪声；分解精度受测不准原理影响；对小波基、分解尺度和阈值估计方法等依赖性太大，需要繁琐的调试才能达到好的效果。EMD消噪方法根据陀螺信号本身特性进行分解得到IMF，无需任何事先的参数设置，消噪过程中所需的阈值及相关系数由原始陀螺信号及其IMF分量计算得到，因此EMD方法的消噪过程只与陀螺信号本身有关。利用EMD方法对陀螺信号消噪，既能够抵制异常噪声干扰，又能削弱高频噪声，减少陀螺各项误差成分，提高惯导解算的精度。

［1］ HUO Ju，WANG Shijing，YANG Ming，et al.Noise Processing of FOG Signal Based on Wavelet Threshold－value［J］.Journal of Chinese Inertial Technology，2008，16（3）：343－347.（霍炬，王石静，杨明，等.基于小波变换阈值法处理光纤陀螺信号噪声［J］.中国惯性技术学报，2008，16（3）：343－347.）

［2］ WAN Yanhui，QING Yongyuan.Application of Wavelet Analysis in Gyro Signal Filtering［J］.Piezoelectrics and Acoustooptics，2005，27（4）：455－457.（万彦辉，秦永元.小波分析在陀螺信号滤波中的研究［J］.压电与声光，2005，27（4）：455－457.）

［3］ WU Fumei，YANG Yuanxi.GPS／INS Integrated Navigation by Adaptive Filtering Based on Wavelet Threshold De－noising［J］.Acta Geodaetica et Cartographica Sinica，2007，36（2）：124－128.（吴富梅，杨元喜.基于小波阈值消噪自适应滤波的GPS／INS组合导航［J］.测绘学报，2007，36（2）：124－128.）

［4］ TANG Wei，LI Shixin，LIU Luyuan，et al.Select of Wavelet Basis in Gyro Signal Processing［J］.Journal of Chinese Inertial Technology，2002，10（5）：28－30.（汤巍，李士心，刘鲁源，等.关于陀螺信号处理中小波基选取的研究［J］.中国惯性技术学报，2002，10（5）：28－30.）

［5］ LIU Bin，JIANG Jinshui，YU Weikai.EMD De－noising Method and Its Application in the Processing for Rolling Mill Signals［J］.Acta Metrologica Sinica，2009，30（1）：73－77.（刘彬，蒋金水，于伟凯.EMD相关度消噪及其在轧机信号处理中的应用［J］.计量学报，2009，30（1）：73－77.）

［6］ JIANG Li，LI Changyun.A Study of Wavelet Threshold Filtering Based on Empirical Mode Decomposition［J］.Signal Processing，2005，21（6）：659－662.（江力，李长云.基于经验模分解的小波阈值滤波方法研究［J］.信号处理，2005，21（6）：659－662.）

［7］ HUANG N E，SHEN Z，LONG S R，et al.The Empirical Mode Decomposition and the Hilbert Spectrum for Nonlinear and Nonstationary Time Series Analysis［C］∥Proceedings of Royal Society.London：［s.n.］，1998：903－993.

［8］ DAI Wujiao，DING Xiaoli，ZHU Jianjun，et al.EMD Filter Method and Its Application in GPS Multipath［J］.Acta Geodaetica et Cartographica Sinica，2006，35（11）：321－327.（戴吾蛟，丁晓利，朱建军，等.基于经验模式分解的滤波消噪法及其在GPS多路径效应中的应用［J］.测绘学报，2006，35（11）：321－327.）

［9］ KATHLEEN T.Stratospheric and Tropospheric Signals Extracted Using the Empirical Mode Decomposition Method［D］.Washington：University of Washington，2003.

［10］ WANG Jian，GAO Jingxiang，WANG Jinling.GPS Baseline Solution Based on Empirical Mode Decomposition［J］.Acta Geodaetica et Cartographica Sinica，2008，37（1）：10－14.（王坚，高井祥，王金岭.基于经验模态分解的GPS基线解算模型［J］.测绘学报，2008，37（1）：10－14.）

［11］ CHEN Jun，XU Youli.Structural Damage Identification Based on Response of Forced Vibration［J］.Journal of Vibration，Measurement and Diagnosis，2005，25（2）：101－104.（陈隽，徐幼麟.经验模分解在信号趋势项提取中的应用［J］.振动、测试与诊断，2005，25（2）：101－104.）

［12］ FLANDRIN P，RILLING G，GONCALVES P.Empirical Mode Decomposition as a Filter Bank［J］.IEEE Signal Processing Letters，2004，11（2）：112－114.

［13］ HUANG N E，WU Zhaohua.A Review on Hilbert－Huang Transform：Method and Its Applications to Geophysical Studies［J］.Reviews of Geophysics，2008，46：1－23.

［14］ BOUDRAA A O，CEXUS J C，SAIDI Z.EMD－based Signal Noise Reduction［J］.International Journal of Signal Processing，2004，1：33－36.

［15］ DONOHO D L，JOHNSTONE I M.Adapting to Unknown Smoothness via Wavelet Shrinkage［J］.Journal of the American Statistical Association，1995，90（12）：1200－1224.

［16］ IEEE STD 952－1997.IEEE Standard Specification Format Guide and Test Procedure for Single－Axis Interferometric Fiber Optic Gyros［S］.［S.l.］：IEEE Standard Board，1997.

［17］ IEEE 1139.Definitions of Physical Quantities for Fundamental Frequency and Time Metrology－Random Instabilities［S］.［S.l.］：IEEE Standard Board，2008.

［18］ LI Xiaoying，HU Min，ZHANG Peng，et al.Applying Overlapping Allan Variance Theory to Better Stochastic Modeling of Microgyro［J］.Journal of Northwestern Polytechnical University，2007，25（2）：225－229.（李晓莹，胡敏，张鹏，等.交叠式Allan方差在微机械陀螺随机误差辨识中的应用［J］.西北工业大学学报，2007，25（2）：225－229.）

De－noising Method for Gyro Signal Based on EMD

GAN Yu，SUI Lifen
Institute of Surveying and Mapping，Information Engineering University，Zhengzhou 450052，China

Gyro random drift is a remarkable factor that can affect the precision of inertial navigation system（INS）.Wavelet de－noising method is poor in coping with exceptional noise，and it depends greatly on the selection of wavelet base and decomposition scale.Empirical mode decomposition（EMD）de－noising method for gyro signal is presented.The signal is decomposed into an intrinsic mode function（IMF）group.Based on this group，IMFs of exceptional noise are first disposed by2σcriterion and then the number of IMFs of high frequency noise is determined by correlation coefficient.The de－noising process is finally done by removing the noisy IMFs.Detailed comparison between EMD method and wavelet method is given.Overlapping Allan variance is used to analyze the effect of the two methods，and the applicable ability of EMD method is tested through an INS calculation.It is shown that EMD method outperforms wavelet method in removing exceptional noise and is more efficient in weakening random drift.

gyro random drift；wavelet；empirical mode decomposition；de－noising

GAN Yu（1988—），male，postgraduate，majors in dynamic geodetic data processing.

1001－1595（2011）06－0745－06

P228

A

国家自然科学基金（40974010）；信息工程大学测绘学院硕士学位论文创新与创优基金（S201101）

宋启凡）

2010－12－06

2010－12－30

甘雨（1988—），男，硕士生，主要从事动态大地测量数据处理研究。

E－mail：ganyu099＠163.com