基于Hausdorff距离的光纤陀螺阈值滤波算法

曹 健 杨志强 刘晨晨

（长安大学地质工程与测绘学院，陕西西安710054）

采矿机械的导航定位与姿态控制是实现采矿自动化的关键技术之一，惯性导航以导航信息的实时性与自主性使其成为一种在复杂采矿工作面可行的采矿机械定位与姿态控制技术。光纤陀螺是一种基于Sagnac效应的全固态角速度传感器，是惯性导航系统（Inertial navigation system，INS）的核心器件[1]，但其随机漂移是制约INS精度的重要因素。为了抑制光纤陀螺的随机漂移，传统的滤波方法主要有两种，一种是进行时间序列分析，建立自回归滑动平均模型（Auto-Regressive Moving Average Model，ARMA），并设计Kalman滤波器进行补偿[2]。而Kalman滤波需要完整的信号和噪声的统计特性，其初始参数确认较复杂，同时在实际情况中，光纤陀螺会受到测试系统噪声及外界环境等各种因素干扰，使输出信号表现出非平稳、弱非线性的特征，因此难以建立准确的数学模型。第二种是使用小波变换的降噪方法[3-5]。小波具有良好的时频局部性和多分辨分析能力，但小波消噪性能非常依赖于基函数的选择和分解层数的预先设定，这些先验信息的确定使去噪方案失去了自适应性，无法保证时变信号小波去噪的效果。

Huang等[6]针对非平稳、非线性信号提出了一种新型的信号处理方法——经验模态分解。EMD是一种不需要任何信号先验信息的自适应数据处理方法，能够根据信号的本身特征，自适应地将信号分解为若干个由高频到低频的本征模态函数及一个余项。

在利用EMD滤波时，很多文献都是通过经验判断或者多次测试，直接把分解后的前若干个IMF分量当作噪声进行剔除，这类方法过于主观，缺乏严谨的理论框架支撑，且很难恰当地排除噪声干扰。将EMD与时间序列分析相结合，可以从全频率角度对随机误差进行分析建模[7]，但该方法需要对所有IMF分别进行建模与去噪处理，过程繁琐，计算量较大。当采用局部重构方法时，需要对提取的IMF进行物理解释，以确定IMF是纯噪声、纯信号，还是两者都包含[8]。基于相关系数的EMD部分重构方法（EMD-Correlation Coefficient，EMD-Cor），使用相关系数来区分相关和不相关的IMF[9]。对于不同信噪比的噪声信号，在某些情况下，由于各IMF与原始信号的相关性太强或太弱，使得该方法非常不稳定。基于连续均方误差准则（EMD-Consecutive Mean Square Error，EMD-CMSE）的滤波方法[10]对IMF在能量第一次发生重大变化处开始进行部分重构，取得了不错的效果，但是在某些情况下EMD-CMSE准则可能会被限制在局部极小值中。Ali Komaty等[11]提出表示数据分布形状的概率密度函数，可以反映两种信号的差异。根据输入信号与各IMF的PDF之间的相似性来选择相关IMF。这种相似性值越大，说明所比较的信号越具有相似的特征，因此重构信号中应包含相关的IMF。而基于信息理论度量（Kullback-Leibler Divergence，KLD）的概率相似性测量方法对于分布的形式（高斯、均匀、双峰等）非常敏感，在某些情况下可能会造成误判。

基于以上问题，本研究提出了一种基于几何度量的概率相似性测量方法——Hausdorff距离[12]与硬阈值相结合的EMD信号滤波算法。

1 滤波方法原理

EMD将任意随时间t变化的信号x()t分解为一组IMF，每个IMF具有不同的频率尺度[6]。信号x(t)可分解为L个从高频到低频排列的IMF，以及一个余项rL(t)，即：

式中，IMFi(t)为第i阶IMF分量。

根据经验模态分解原理，对于混有随机噪声的信号，低阶IMF分量通常对应于信号的高频噪声[13]。若除去若干个低阶IMF，把剩余的相关IMF与余项重构为一个新信号，可以削弱随机噪声的影响。相关IMF的选择需要根据一个给定的准则，该准则能够识别携带了输入信号信息的IMF分量。

假设y(t)为无噪的纯信号，n(t)为高斯白噪声，则光纤陀螺输出信号x(t)模型可表示为

去噪的目标是找到x(t)的估计y͂(t)，来提高信号的信噪比，即：

式中，kth为局部重构起始分量的阶数，使用该式的关键是如何确定参数kth。

由于PDF包含其对应IMF的完整信息，因此可以使用PDF相似性量度来识别带有y(t)特征的IMF。对原始信号进行EMD分解之后，使用核密度函数估计IMF的PDF。所以就要寻找这样一种度量方法，利用PDF来区分相关和不相关的IMF，也就确定了参数kth。

1.1 Hausdorff距离

Hausdorff距离是描述两个点集之间距离的方法，该方法可作为两个点集之间的相似性量度。在本研究中，HD不是用作两个几何形状间的相似性量度，而是用在两个一维PDF之间。HD相较于其他距离度量方法的优点是对异常值非常敏感，同时顾及目标的整体形状（或边界），而不是将两个空间目标的距离简单地表达为两个点之间的距离。

令d(a,b)=‖a-b‖为元素a、b之间的欧氏距离。元素a到有限集合B={b1,…,bN}之间的距离定义为

两个有限集A={a1,…,aM}和B={b1,…,bN} 之间的单向HD可以定义为

类似地，两个有限集B={b1,…,bN}和A={a1,…,aM}之间的单向HD可定义为

单向HD具有非对称性，这是极大极小函数的特性。绝对HD为2个单向HD的最大值，即：

如果A、B两组点是相似的（距离较近），但B中存在一点b，有远离A的任意一点，那么该距离就是由这个点决定的。这一距离可以表明IMF的分布曲线有多少尖锐和突出的部分。

1.2 IMF的选择及阈值滤波

为了选择相关的IMF，首先将信号x(t)分解为若干个IMF，之后估计x(t)及各IMF的PDF，分别记为PDF[x(t)]与PDF[IMFi(t)]。将估计的分布看作度量空间（欧式空间）的子集，用HD测量相似程度。采用（5）式、（6）式和（7）式计算x(t)与各IMF之间PDF的HD值，记为HD(i):

重构所选择的起始IMF分量是HD在第一个局部最大值之后开始减小的第一个IMF。记这个指数为kHD，即：

得到kHD后，则可确定信号的局部直接重构信号y1(t)：

若直接重构，则重构信号中仍然含有部分噪声。所以在重构之前需要对相关IMF进一步去噪，方法可以借鉴小波阈值去噪。阈值去噪方法有硬阈值和软阈值[14]，两种阈值后的各阶IMF分别见式（11）、式（12）。

式中，kHD≤i≤L；Ti表示第i个IMF的阈值，Ti的选择可采用小波阈值估计方法[15-16]。

软阈值的优点是具有良好的连续性和平滑性，但它是对所有系数进行压缩，如果阈值估计产生误差，该误差会扩散到所有系数上，导致有偏的结果。而硬阈值只对低于阈值的部分进行处理，对高于阈值的部分予以保留，可以免受阈值估计误差的影响，故而本研究采用硬阈值滤波方法。

由于使用了硬阈值的消噪策略，随着i的增加，IMF系数的大部分幅值都高于阈值，阈值的作用会逐渐降低。所以无需判别噪声信号混合的IMF与有用信号为主的高阶IMF的边界，不加区分地对所有需要被重构的模态进行阈值滤波。所以最终的重构信号为

综上所述，本研究提出的基于Hausdorff距离的EMD阈值滤波（EMD-Hausdorff Distance&Threshold Filtering，EMD-HD&TF）算法的具体实现步骤如下：

（1）对原始信号进行EMD分解，得到若干个IMF分量。

（2）使用核密度函数，估计原始信号及各IMF的PDF。

（3）使用HD方法比较IMF分量与原始信号PDF的相似性，确定重构的起始分量——第kHD个IMF。

（4）计算各需要重构的IMF的噪声阈值Ti，按照式（13）重构输出信号。

2 试验与结果分析

2.1 数据仿真试验

为验证硬阈值的优越性，对测试信号（以Doppler信号为例）加入不同大小的白噪声，输入信噪比为-6～6 dB，分别使用硬阈值与软阈值对所有IMF进行滤波并重构，结果如图1。硬阈值的去噪效果在各信噪比下都优于软阈值，故本研究选择HD筛选准则与硬阈值滤波相结合。

为测试本研究滤波方法的有效性，对4个代表性的测试信号（Doppler、Heavysine、Bumps和Blocks）分别加入以上噪声，使用EMD-Cor、EMD-CMSE、EMDHD（EMD-Hausdorff Distance）及 EMD-HD&TF 4种方法进行去噪，对去噪结果进行比较研究。

以输入信噪比为-4 dB的Doppler含噪信号为例，其经过EMD分解后得到的10个IMF见图2。各阶IMF的HD取值见图3，由此可得kHD=6。在EMDHD&TF方法中，对IMF6～10分别进行阈值滤波，结果如图4所示。

最终按式（13）将滤波后的IMF与余项重构。从图5中可以看出，相比于其他几种方法，由EMDHD&TF方法得到的信号更能保留Doppler标准信号的特征，初步判断为最佳的去噪结果。表1给了的4种滤波结果的信噪比（Signal to Noise Ratio，SNR）及均方根误差（Root Mean Square Error，RMSE）。

从表1可知，4种方法都对白噪声具有一定的去噪能力，但EMD-Cor方法效果相对较差。两种基于HD准则的EMD局部重构方法的去噪性能明显优于另外两种方法，在此基础上使用硬阈值方法进行滤波可以进一步提升去噪效果。

统计4种信号在不同输入信噪比下滤波后的信噪比及均方根误差，如图6与图7所示。由此可知，基于EMD的滤波方法对于不同类型的信号都表现出不同程度的去噪效果，具有一定的适用性。其中EMD-Cor方法在输入信噪比较低时性能较差，但会随着输入信噪比的提高而改善，逐渐达到与其他方法相近的水平。但在多次试验中，经常会有重构后的信号与原信号相似度过高或者过低的现象发生，导致无法进行筛选，所以该方法稳定性较差。EMD-HD&TF在EMD-HD基础上加入硬阈值，虽然滤波效果在本研究中只是略有提升，但在绝大多数情况下都能得到最佳的滤波效果。上述去噪效果不仅支持了EMD-HD和阈值滤波混合方法的有效性，而且证明了IMF中所携带的信号信息可以通过它们的PDF来表现，同时HD比EMD-Cor和EMD-CMSE法更有效地量化了信号与相关IMF之间的相互作用，进而使用HD进行筛选更为有效。

2.2 陀螺滤波试验

在常温（25℃）下，对某型号陀螺漂移信号进行了试验测试。采样周期为100 ms，为了避免启动过程中的温度变化，选取中间3 h的数据，样本总数为108 000个。对该组数据分别采用文中的4种方法进行处理。其中EMD-Cor、EMD-CMSE、与两种HD方法重构的起始IMF分量分别为第5阶、第7阶和第9阶。为了验证降噪效果，采用了惯性器件误差分析中最常用的艾伦（Allan）方差。Allan方差的双对数曲线图可以识别陀螺的各种随机误差成分，包括量化噪声（Q）、角度随机游走（N）、零偏不稳定性（B）、角速率随机游走（K）和速率斜坡（R），由试验所得的Allan方差图如图8所示。不同的误差具有不同的斜率特性，可由Allan方差拟合得到[17-18]，拟合得到的各种误差成分如表2所示。

由图8和表2可知，4种方法对角速率随机游走（K）和速率斜坡（R）均有良好的抑制效果。但经EMD-Cor和EMD-CMSE方法处理后，量化噪声（Q），角度随机游走（N）和零偏不稳定性（B）不但没有减少，甚至还有小幅度增加。基于HD方法得到的重构IMF较为恰当，其抑制5种噪声的效果优于前两者。在此基础上加入阈值滤波的EMD-HD&TF方法又能进一步滤除噪声，使各项噪声得到了进一步降低，取得了最佳滤波效果。

?

3 结 论

（1）本研究提出了一种新的EMD局部重构与阈值相结合的EMD-HD&TF方法，该方法首先以各阶IMF和原始信号PDF的HD为依据对IMF进行筛选，减少了IMF的重构个数；然后采用硬阈值处理方法，对信号的噪声进行了进一步滤除。该方法的优点是对异常值敏感、复杂度低，在滤波应用中较为实用。

（2）为了验证该方法的有效性，在不同的信噪比情况下进行了数值模拟。与EMD-Cor、EMD-CMSE和单纯的EMD-HD方法相比，该方法具有更好的去噪性能。最后，将该方法应用于对实测陀螺漂移信号进行试验。从Allan方差曲线可以看出，EMD-HD&TF方法大大降低了各项随机误差。因此，本研究提出的光纤陀螺消噪算法在采矿机械导航定位与姿态控制的自动化方面具有广泛的应用前景。

?