基于量子高斯混合模型的振动信号降噪方法

杨望灿， 张培林， 陈彦龙， 吴定海， 李海平

(1.91404部队,秦皇岛 066004;2.陆军工程大学石家庄校区 七系，石家庄 050003;3.陆军特种作战学院,桂林 541000；4.陆军工程大学石家庄校区 六系，石家庄 050003)

在机械设备运行过程中，由于其工作条件和工(作环境较为恶劣，因此机械设备部件容易出现磨损和发生故障[1]。为了及时了解机械设备的运行状态和及早发现机械设备的故障情况，振动信号监测和分析方法是一种行之有效的技术手段[2]。但是，由于机械设备的运行工况较为复杂，振源较多，采集到的振动信号经常被强背景噪声所淹没，导致机械设备故障状态特征不明显，降低了对机械设备运行状态判断的准确程度。因此，如何对采集的机械振动信号进行降噪处理得到了广大学者的热切关注[3-5]。

由于小波变换具有多分辨率的时频分析特性，因此被用于多种类型信号的降噪分析处理中[6-7]。其中，部分学者通过假设小波分解系数符合某一概率统计模型，结合贝叶斯先验知识，估计噪声小波系数方差，得到小波系数收缩函数，从而实现信号的降噪。文献[8]将双树复小波分解信号得到的高频系数和低频系数分别采用最大后验估计算法和广义形态滤波进行降噪处理，消除可见光近红外光谱噪声。文献[9]应用高斯混合模型(Gaussian Mixture Model，GMM)对小波分解的各类系数建模，通过噪声方差估计对信号进行降噪。其中，GMM能够较好地描述小波系数的分布情况，降噪效果比较理想。但是传统的GMM使用全局统一的阈值函数，算法的局部自适应较差，使算法的应用受到一定的限制[10]。

自从Eldar等[11]提出量子信号处理(Quantum Signal Processing，QSP)理论后，基于量子理论的信号处理技术逐渐得到了广泛地应用[12-14]。文献[15]通过建立量子模型，对声音信号进行量子傅里叶编码，相比传统的声音信号处理方法，该方法在鲁棒性和计算复杂度方面均有改善。文献[16]将图像像素的位置信息和灰度值信息分别用比特量子系统的基态和量子比特描述，实现了量子图像的水印加注。受量子衍生信号处理技术的启发，本文提出了一种基于量子高斯混合模型(Quantum Gaussian Mixture Model，QGMM)的振动信号降噪方法，将量子叠加态理论应用于双树复小波包系数建立的高斯混合模型中，改善高斯混合模型的局部自适应性，使小波包系数自适应非线性收缩，提升机械振动信号的降噪处理效果。

1 高斯混合模型的建立

1.1 双树复小波包变换

双树复小波变换是一种有限冗余的小波变换方法，具有近似平移不变性，改善了传统离散小波变换的Gibbs效应和平移敏感性。由于双树复小波变换没有对高频部分进一步分解，因此为了提高振动信号的频率分辨率，本文采用双树复小波包变换(Dual-Tree Complex Wavelet Packet Transform，DTCWPT)对信号的高频部分和低频部分同时进行双树复小波变换，用两棵并行的小波树实现振动信号的分解。双树复小波包变换如下式所示

(1)

式中：l=(1，2，…，J)为变换尺度因子;J为最大变换尺度。

双树复小波包变换第一层采用传统的非下抽样小波变换分解，然后对小波系数进行奇偶分离，分别采用Mallat算法进行分解和重构，双树复小波包分解的小波系数如下式所示

(2)

(3)

式中:h(n)和g(n)为希尔伯特变换对的滤波器组，满足半采样延迟，即

g(n)≈h(n-0.5)

(4)

信号重构时，根据下式对树a和树b进行联合重构

(5)

1.2 高斯混合模型

高斯混合模型是一种应用较为广泛的数理统计模型，其理论基础为概率统计，具有较强的灵活性.机械振动信号经过双树复小波包变换后，其小波分解系数呈现稀疏和聚集的分布特性，即在零值附近小波系数较为稀疏，但幅值较大；在分布两端小波系数较为聚集，但幅值较小[17-18]。由于高斯混合模型的分布曲线具有“高峰值”和“长拖尾”的特性，所以可以用高斯混合模型来描述机械振动信号小波系数的分布规律。

对含有加性高斯白噪声的信号y(t)(y(t)=x(t)+n(t))进行双树复小波包变换，根据小波变换的线性性质，信号分解得到的小波系数满足以下关系

Y=X+N

(6)

式中:Y=Yr+iYi，X=Xr+iXi，N=Nr+iNi分别为含噪信号，有用信号和噪声信号的双树复小波包系数。

对含噪信号进行双树复小波包分解后，将同一层次、同一子带内复小波系数的实部系数和虚部系数按照一一间隔的方式排列，即Yd(s)=[Yr(1,s)，Yi(1,s)，Yr(2,s)，Yi(2,s)，…，Yr(nk,s,s)]，其中nk,s为第s层第k个子带中小波系数的个数。对小波系数排列Yd(s)采用两状态的高斯混合模型建立模型，其概率密度函数如下式所示

F(Yd)=a1Gaussian(Yd,σ1)+a2Gaussian(Yd,σ2)=

(7)

式中：a1，a2分别为各个单高斯分布的权重;σ1，σ2分别为各个单高斯分布的标准差。

对于高斯混合模型中参数a1，a2，σ1和σ2，采用最大期望(Expectation Maximization，EM)算法计算确定。EM算法首先计算对数似然函数的期望，然后搜索使期望值最大的参数，得到混合高斯模型参数的最优值，具体步骤参见文献[19-20]。

对于高斯白噪声信号，噪声信号的小波系数近似服从均值为0，标准差为σN的高斯分布，其概率密度函数如下所示

(8)

1.3 贝叶斯最大后验估计

对于振动信号的降噪，可以通过从观测信号y(t)中尽可能精确地估计出真实信号x(t)来完成信号的降噪。本文采用贝叶斯最大后验估计(Maximum A Posterior，MAP)算法对建立的高斯混合模型进行求解，实现信号的降噪处理。

(9)

根据贝叶斯理论

(10)

则式(9)变换为

(11)

此式等价于

(12)

令f(X)=lnpX(X)，对式(12)求一阶导数，并使其等于0，得到下式

(13)

将式(7)、式(8)代入式(13)，经过推导可得X的估计式，即小波系数收缩函数，如下所示

(14)

式中：σN为噪声信号的标准差，其鲁棒估计为双树复小波包系数的绝对值的中值[21]，即

σN=median(Yd)/0.674 5

(15)

根据建立的混合高斯模型和MAP估计，可以从含噪信号的小波系数中估计出有用信号的小波系数，然后对估计的小波系数进行重构，实现信号的降噪。但是从式(14)可知，传统的高斯混合模型为全局概率模型，其参数均为全局统一的参数，不具有局部自适应性。当信号中信噪比较低时，传统的高斯混合模型对小波系数的分布不能达到理想的拟合精度，影响信号的降噪效果。

离合词“A了个B”与网络词“A了(嘞)个B”同属概念重组的结果，形式上相同，两者都不能接宾语，但如果我们深究入语法功能、语体运用等方面，可发现它们之间的巨大差别。

2 量子高斯混合模型

传统的高斯混合模型本质上是一种层内模型，其仅考虑了同一层次内小波系数的分布情况，但实际上信号的小波系数还存在着高阶相关性，如尺度间的传递性。对于含噪信号的小波系数，其实质上是有用信号小波系数和噪声信号小波系数的叠加，与量子信号处理理论中量子叠加态原理相似，因此，根据量子叠加态理论和小波系数尺度间的相关性，提出了一种量子高斯混合模型，使高斯混合模型具有局部自适应性，改善信号的降噪效果。

2.1 量子叠加态原理

在量子理论中，量子比特是描述量子世界的基本单位，量子比特有|0〉和|1〉两个基本状态。根据量子叠加态原理，量子比特所表达的状态如下所示

|Ψ〉=a|0〉+b|1〉

(16)

式中，系数a和b称为|0〉和|1〉两个基态的量子概率幅，其模的平方称为量子概率，表示对应基态出现的概率值。量子概率幅需满足归一化条件，即

(17)

含噪信号的小波系数中包含有用信号和噪声信号两种小波系数，如果用基态|0〉表示噪声信号的小波系数;基态|1〉表示有用信号的小波系数。这样量子概率幅a和b就可以表示噪声信号和有用信号小波系数的概率，从而可以分析处理二者之间的关系。

2.2 基于量子叠加态的模型参数估计

由于小波系数具有尺度间的传递性，当父代系数的模较大时，那么其子代系数的模也大，反之，当父代系数的模较小时，其子代系数的模也较小。所以将父代和子代的小波系数取模相乘得

(18)

(19)

(20)

(21)

(22)

将小波系数的量子比特叠加态应用于传统高斯混合模型得到的小波系数收缩函数式(14)可得

(23)

从式(23)可知，量子高斯混合模型通过小波系数的量子叠加态特性，根据信号自身的特点，对传统高斯混合模型中参数进行了局部自适应调整。当尺度s+1下第j个小波系数中出现有用信号的概率大时，则自适应地增大该位置有用信号的标准差，反之，当尺度s+1下第j个小波系数中出现噪声信号的概率大时，则自适应地增大该位置噪声信号的标准差，提高高斯混合模型对小波系数的分布的拟合精度，改善信号的降噪效果。

2.3 基于量子高斯混合模型的降噪流程

基于量子高斯混合模型的信号降噪流程步骤如表1所示。

表1 量子高斯混合模型降噪算法

3 仿真信号分析

为了模拟行星齿轮箱齿轮局部故障，构造了如下的仿真信号

(24)

仿真信号的采样频率为4 096 Hz，采样点数为2 048个。由式(24)中可知，仿真信号由一个正弦谐波分量，一个调幅分量和一个调频分量线性叠加而成。在仿真信号中加入高斯白噪声，双树复小波包采用4层分解，采用信噪比σSNR作为评价指标，信噪比计算公式如下

(25)

在实验过程中，引入传统的高斯混合模型、硬阈值和软阈值的降噪方法作为对比实验，传统的硬阈值、软阈值的阈值计算按下式计算

(26)

式中:nk,s为第s层第k个子带中小波系数的总个数。硬阈值降噪方法将绝对值低于阈值的小波系数置0，保留其他小波系数进行重构信号；软阈值降噪方法将绝对值低于阈值的小波系数置0，将绝对值高于阈值的小波系数减去阈值后进行重构信号。

不同方法对仿真信号的实验结果如图1所示。从图1可知，对于仿真信号，量子高斯混合模型去噪效果优于其他几种方法，降噪后的信噪比最高。硬阈值和软阈值降噪方法，在去除噪声的同时，也消除了部分原始的真实信号，导致信号部分失真，信噪比提升有限。传统的高斯混合模型采用全局统一的模型参数，导致噪声残留较多，降噪后的信噪比也不高。而量子高斯模型将小波系数尺度间的传递性用量子比特的叠加态特性描述，自适应调整模型参数，在保留原始波形特征的基础上，最大限度的消除噪声，取得更优的降噪效果。

4 实测信号分析

实测信号来自如图2所示的行星齿轮箱传动系统实验台，加速度传感器安装在行星齿轮箱的箱体上，采集实验过程中的振动信号，采样频率为20 kHz。该实验通过机械加工的手段模拟了太阳轮齿轮的单个齿齿面的轻微磨损故障。实验过程中，电机转速为1 200 r/min，磁粉制动器提供的负载为1.2 N·m，采样时长为2 s。行星齿轮箱的结构参数如表2所示。根据理论计算，太阳轮齿轮故障的特征频率为55.5 Hz。

图1 不同方法对仿真信号的降噪处理

Fig.1 Denoising process for simulated signal by different methods

图2 行星齿轮箱实验台

齿轮太阳轮行星轮(个数)齿圈齿数1364(3)146

图3为行星齿轮箱实验台采集的太阳轮齿轮故障的原始振动信号的时域波形和包络频谱图，为了便于图形展示，时域波形的时间取0～0.5 s，频谱图中频率取0～500 Hz。从图3可知，由于故障较为轻微，故障齿轮所产生的冲击信号不明显，基本淹没在背景噪声中，在包络谱中也看不到太阳轮故障的特征频率，无法判断行星齿轮箱是否发生故障。

图4～图7为分别采用不同方法对实测太阳轮故障振动信号的降噪结果。实验过程中，双树复小波包分解都为4层分解。从图4可知,采用量子高斯混合模型降噪后，振动信号的时域波形出现了冲击特征，对其进行包络解调，从包络谱图中可以明显地看到太阳轮故障的特征频率55.5 Hz及其2倍频、3倍频和4倍频，符合太阳轮故障的频谱特点，所以可以判断行星齿轮箱中太阳轮轮齿发生了故障。图5为传统高斯混合模型的降噪结果，虽然从包络谱图中能看到太阳轮故障的特征频率55.5 Hz，但是仍有大量的干扰噪声，特征频率淹没在其他干扰频率中，不能准确判断行星齿轮箱中太阳轮是否存在故障。图6和图7为采用小波硬阈值和小波软阈值的降噪结果，噪声被大量去除，但是由于太阳轮齿轮磨损轻微，太阳轮故障特征信号也被当做噪声信号被去除或是削弱，所以包络谱图中也无法得到太阳轮故障的特征频率。

(a) 时域波形

(b) 包络谱

(a) 时域波形

(b) 包络谱

(a) 时域波形

(b) 包络谱

(a) 时域波形

(b) 包络谱

5 结 论

本文提出了一种基于量子高斯混合模型的机械振动信号降噪方法。该方法在传统高斯混合模型的基础上，根据量子信号处理理论，将信号小波系数相邻尺度间的传递特性转化为量子比特叠加态，使高斯混合模型中的参数根据信号特征自适应地调整大小，增强了高斯混合模型的局部自适应性。仿真信号分析表明，相比于其他方法，量子高斯混合模型降噪方法能够有效提高信号的信噪比。实测行星齿轮箱振动信号结果表明，本文所提方法能够抑制信号中的噪声，保留太阳轮齿轮轻微故障的特征信号，准确提取其故障特征频率，有效地判断行星齿轮箱的故障状态。

(a) 时域波形

(b) 包络谱