CN108490349B - 基于Mel频率倒谱系数的电机异音检测方法 - Google Patents

Abstract

本发明公开了基于Mel频率倒谱系数的电机异音检测方法，包括设定采样频率和采样时长t；对音频信号进行分帧加窗处理，设置每帧的帧长L和相邻两帧的重叠长度(帧移)M，将音频信号分为N帧信号；使用基于Mel频率倒谱系数提取音频特征；计算出当前音频信号每一帧的mfcc参数，得到每个mfcc参数随帧数(时间)的变化曲线等步骤；本发明具有诸多优势，进行电机异音识别，提高检测效率并且保证产品的出厂质量，进而提高企业的整体的生产效率，降低企业的制造成本，同时了保护了工人的身体健康，可以有效地解决电机音频信号非稳态的问题，有效检测非稳态电机异音故障，识别准确率高。

Description

基于Mel频率倒谱系数的电机异音检测方法

技术领域

本发明涉及电机故障检测领域，特别是电机异音的检测方法。

背景技术

我国是洗衣机、家用空调、电冰箱、电风扇等各类小型家电用电动机的主要产地，仅家用空调电机的年产量就逾数十亿台之多。

在小型电机生产线上，产品下线之前普遍采用人工听音的方法分辨良、次品，即工人在隔音房内，用耳朵依次听取电机运行时的声音，通过工人的个人经验判断电机是否存在故障。

由于需要人做主观判断，长期以来一直难以由自动化装置替代。并且评价者的个人经验，无法建立统一的评价标准，不同的评价者可能产生不同的结论。而大批量生产过程中该工序不仅要消耗大量劳动成本，而且重复、单调的听音工作极易引起人员疲劳，容易出现误判，若个别不良品混入整批成品中，会给工厂带来严重经济损失，甚至严重影响产品声誉。

CN201510266743.8(公布号CN104992714A)披露了一种电机异音的检测方法，包括以下步骤：1、将电机处于空载状态下进行音频采集；2、将所采集到的电机的时域音频信号经傅里叶变换转换为频域波形；3、若电机的正常频域范围的最高值外存在波形，则认为该电机存在异音；若电机的正常频域范围的最高值外不存在波形，则表明此电机不存在异音。这种电机异音的检测方法的缺点在于：1、音频信号属于一种准稳态信号，即短时稳定；傅里叶变换是对稳态信号处理的常用手段，只能对稳态的信号进行信号特征提取；但是电机异音故障中，存在大量瞬时非稳态信号的故障样例；对于这些非稳态的信号，傅里叶变换无能为力。2、通过判断在指定的最高值之外是否存在波形来判断异音，不具备自适应能力，通用性不高：因为对于不同型号的电机产品就需要技术人员重新设定阈值；并且通过传感器采集音频信号时，传感器距离音源的远近也会影响阈值设定。3、无论是设置最大值还是比较波形图，都需要专业人员来判断，无法实现全自动识别。

发明内容

本发明的目的在于提供一种能够识非稳态信号的故障样例，通用性高，且故障识别过程中无需人员参与的电机异音检测方法。

基于Mel频率倒谱系数的电机异音检测方法，包括以下步骤：

步骤1、设定采样频率和采样时长t，将电机处于空载状态下进行音频信号采集x(n)，x(n)＝x₁(n)*x₂(n)；音频信号x(n)是一种振动信号，音频信号是由于振动源(激励信号，不妨设为x1(n))作用于腔体产生的，在信号处理中，只要知道腔体的单位冲击响应函数(不妨设为x2(n))，通过x1(n)*x2(n)就能得到输出信号x(n)，其中“*”表示卷积操作。对于电机异音检测来说，合格电机的激励信号与故障电机的激励信号不一样，但是腔体的冲击响应都是一样的(腔体的冲击响应只和腔体的大小、材料、形状等有关，对于同一种电机，这些自然都是一样的)。但是通常来说，无法直接区分合格电机的x(n)和故障电机的x(n)，所以需要一种方法，能够将卷积操作变成加法操作，所以就之后的“倒谱”操作，它就是一种解卷积方法，将卷积变成加法。

步骤2、对音频信号进行分帧加窗处理，设置每帧的帧长L和相邻两帧的重叠长度(帧移)M，将音频信号分为N帧信号；帧长L以每帧内的信号能看做稳态信号为宜，从而避免整段音频信号非稳态、时变的影响；窗函数选择汉宁窗。

步骤3、使用基于Mel频率倒谱系数提取音频特征：

步骤3-1、对每一帧信号进行FFT，将时域数据转变为频域数据：

X(i，k)＝FFT[x(i，n)]

步骤3-2、对每一帧FFT后数据的X(i，k)计算谱线能量：E(i，k)＝[X(i，k)]²

步骤3-3、计算每帧信号通过Mel滤波器组之后，在各个滤波器中的能量：

步骤3-4、计算DCT倒谱；

对每一帧信号，对每个Mel滤波器的能量求对数后再进行离散余弦变换

步骤4、计算出当前音频信号每一帧的mfcc参数，得到每个mfcc参数随帧数(时间)的变化曲线；

步骤5、选取N个合格的电机音频样本，重复步骤2-4，对于每个电机音频样本，得到以帧数为横坐标每个MFCC参数随帧数变化曲线图。对于每个电机音频样本，提取每个mfcc参数随帧数变化曲线的平均值作为一组特征，训练一个基于支持向量数据描述的学习模型modelA。然后对每个电机音频样本，提取每个mfcc参数随帧数变化曲线的峰峰值作为一组特征，训练另一个基于支持向量数据描述的学习模型modelB。

步骤6、进行电机异音诊断时，重复步骤1-4获得电机音频的音频样本特征，提取每个mfcc参数随帧数变化曲线的平均值作为特征，输入学习模型modelA，判断是否为故障样本；然后提取每个mfcc参数随帧数变化曲线的峰峰值作为一组特征，输入学习模型modelB，判断是否为故障样本。只有当两个学习模型都判断为合格时，才认定该样本为合格样本，否者都判断为故障样本。

进一步，步骤3-4中的mfcc参数为静态特征的mfcc参数，计算公式为：

S(i，m)是由步骤3-3中求出的Mel滤波器能量；m是指第m个Mel滤波器(共M个)；i是指第i帧；n是DCT后的谱线。

或者，步骤3-4中可以加上mfcc参数的动态特征，

计算公式为

其中，i表示帧号，n表示DCT后的谱线，也即MFCC的参数序号。然后将一阶微分参数放入标准MFCC参数之后。标准静态mfcc参数一共12维，加上12维动态mfcc参数，每一帧信号可以提取24维参数。由于首尾两帧无法计算动态mfcc参数，故在计算m帧音频的mfcc参数时，理论上应该得到m×24大小的mfcc参数矩阵，其中第i行代表第i帧音频的mfcc参数，第j列表示第j个mfcc参数随时间变化趋势，但是由于第1、2、m-1、m这四帧无法计算动态mfcc参数，而实际中音频的帧数一般比较多，所以在最后得到的MFCC参数中，去掉这四帧，剩下一个(m-4)×24大小的矩阵。

进一步，步骤2中，帧长L＝1S，重叠长度M＝0.5s。

进一步，步骤3-3中的滤波器如下：在人耳能够感受的音频范围内，设置M个带通滤波器Hm(k)，0≤m＜M，每个滤波器都具有三角形滤波特性，中心频率为f(m)；在Mel刻度上，M个带通滤波器均为等带宽；带通滤波器的传递函数为

其中

f_l表示滤波器组频率范围内的最低频率，f_h为滤波器组频率范围的最高频率，N为FFT的长度，f_s为采样频率，F_mel＝2059log₁₀(1+f/700)，

是F_mel的逆函数：

本发明的优点在于：

1、可以辅助工人进行电机异音识别，提高检测效率并且保证产品的出厂质量，进而提高企业的整体的生产效率，降低企业的制造成本，同时了保护了工人的身体健康。

2、可以有效地解决电机音频信号非稳态的问题，有效检测非稳态电机异音故障，识别准确率高。

3、本文通过使用支持向量机来自动判断电机音频是否存在异音故障，使用过程不需要人工干预。设置支持向量机的罚函数、核函数和异常样本比例三个参数，对不同电机类型使用不同电机特征进行训练，就可以得到适用于对应电机的支持向量机模型。并且三个参数对不同电机都具有通用性，无需更改。使用机器学习的方法，也不需要为对不同电机设定不同的阈值。支持向量机会根据训练样本设定好针对该型号电机的判别函数。这种方法的特点是以生产线上极易获得的正常声音样本为基础建立判别函数。并且，随着样本数量增加、对特征分布的覆盖范围扩大，判别的准确性也随之增加。

4、Mel滤波器的带宽是随着频率的升高而逐渐增加的，而人耳对音频的分辨率也是随着频率的升高而逐渐降低的，使用MFCC能够很好的模拟人耳对音频的感知。MFCC提取的特征可以很好地辨识出电机音频信号中的异音故障成分。

附图说明

图1是Mel频率与实际频率的关系曲线。

图2是汉宁窗分帧加窗的示意图。

图3是Mel组在频率上的分布图。

图4是MFCC特征参数提取流程图。

图5是同态处理中实现解卷积流程图。

图6是16个Mel滤波器的频率响应曲线图。

图7是MFCC识别电机的合格信号和故障信号的曲线图。

具体实施方式

电机异音检测方法，包括以下步骤：

步骤1、设定采样频率和采样时长t，将电机处于空载状态下进行音频信号采集x(n)，x(n)＝x₁(n)*x₂(n)；

步骤2、设置帧长L和相邻两帧的重叠长度M，使用汉宁窗对音频信号进行分帧加窗处理，如图2所示，将音频信号分为N帧信号；帧长L以每帧内的信号能看做稳态信号为宜，从而避免整段音频信号非稳态、时变的影响；

步骤3、使用基于Mel频率倒谱系数提取音频特征：

步骤4、计算出当前音频信号的第i帧的mfcc参数；

步骤5、选取N个合格的电机音频样本，重复步骤2-4，对于每个电机音频样本，得到以帧数为横坐标每个MFCC参数随帧数变化曲线图。对于每个电机音频样本，提取每个mfcc参数随帧数变化曲线的平均值作为一组特征，训练一个基于支持向量数据描述的学习模型modelA。然后对每个电机音频样本，提取每个mfcc参数随帧数变化曲线的峰峰值作为一组特征，训练另一个基于支持向量数据描述的学习模型modelB。

步骤6、进行电机异音诊断时，重复步骤1-4获得电机音频的音频样本特征，提取每个mfcc参数随帧数变化曲线的平均值作为特征，输入学习模型modelA，判断是否为故障样本；然后提取每个mfcc参数随帧数变化曲线的峰峰值作为一组特征，输入学习模型modelB，判断是否为故障样本。只有当两个学习模型都判断为合格时，才认定该样本为合格样本，否者都判断为故障样本。

Mel频率倒谱系数(Mel Frequency Cepstrum Coefficient，MFCC)是基于人的听觉机理，即根据人的听觉实验结果来分析音频信号的频谱，以期望能够获得更好的音频特征。MFCC分析依据的听觉机理有两个。

第一，人的主观感知频域的划分不是线性的，根据Stevens和Volkman(1940)的工作，可以得出实际音频信号的频率与人耳感知的频率满足式(3-1)，曲线如图1所示：

F_mel＝2059log₁₀(1+f/700) (3-1)

其中，f是视频信号的实际频率，单位Hz。F_mel是人耳的主观感知频率，单位是美尔(Mel)。

第二，人耳基底膜具有与频率分析器相似的作用。具体来说，在人耳频率的响应范围(20～22050Hz)可以分为成25个频率群，如表1所示。

表1

人耳的基底膜可以划分成许多很小的部分，每一部分都对应一个频率群。大脑对同一个频率群内的信号是叠加在一起评价的。但是这些滤波器在频率坐标轴上却不是统一分布的，在低频区域有很多的滤波器，他们分布比较密集，但在高频区域，滤波器的数目就变得比较少，分布很稀疏。

在人耳能够感受的音频范围内，设置M个带通滤波器H_m(k)，0≤m＜M。每个滤波器都具有三角形滤波特性，其中心频率为f(m)。在Mel刻度上，这些滤波器都是等带宽的。每个带通滤波器的传递函数如式(3-2)所示。

f(m)可以用下面的方法加以定义：

式中，f_l表示滤波器组频率范围内的最低频率，f_h为滤波器组频率范围的最高频率。N为FFT的长度，f_s为采样频率。

是F_mel的逆函数：

由这M个滤波器组构成了MFCC特征中关键的Mel滤波器组。图3是Mel组在频率上的分布图(使用三角窗)。由图3可以看出，这些三角滤波器的中心频率在低频分布密集，在高频部分分布稀疏。但是转换到Mel频谱上则是均匀分布的。这个特征与人耳的感知实验结果想匹配。

选择使用三角带通滤波器的主要目的有两个：

1)对频谱进行平滑化，并消除谐波的作用，突显原先语音的共振峰。(因此一段语音的音调或音高，是不会呈现在MFCC参数内，换句话说，以MFCC为特征的语音辨识系统，并不会受到输入语音的音调不同而有所影响)

2)降低运算量。

MFCC特征参数提取流程图如图4所示：

1)音频信号x(n)首先经过预处理，即分帧加窗。

2)对每一帧信号进行FFT，将时域数据转变为频域数据：

X(i，k)＝FFT[x(i，n)] (3-5)

3)对每一帧FFT后数据的X(i，k)计算谱线能量：

E(i，k)＝[X(i，k)]² (3-6)

4)计算每帧信号通过Mel滤波器组之后，在各个滤波器中的能量：

在频域中的操作就相当于用3)中的谱线能量与Mel滤波器的频响H_m(k)，就可以得到该滤波器下具有的能量。

5)计算DCT倒谱

倒谱是同态处理中常用的一种技术。音频信号x(n)一般都可以表示为一个激励信号x₁(n)和物体单位冲击响应函数x₂(n)的卷积x(n)＝x₁(n)*x₂(n)，符号“*”表示卷积。在音频识别中求取激励信号x1(n)非常重要。因为对于某种特定的物体，它的单位冲击响应函数相差不大。不同音频在听觉上的差异主要是x₁(n)不同造成的。

由卷积结果求取参与卷积的各个信号，称为解卷积。在音频信号中处理解卷积常用的有两种：一种是线性预测分析；另一种是同态处理。在MFCC中用的就是同态处理技术，它实现了将卷积关系变成求和关系。

求倒谱就是同态处理中实现解卷积的关键技术，系统流程图如图5所示：

其中FT[]表示傅里叶变换。FT^-1表示傅里叶的逆变换。可以看出音频信号x(n)＝x₁(n)*x₂(n)经过该系统后变成

由此就可以得到需要求激励信号和物体单位冲击响应信号。

信号也是时域序列，但是它所处在的离散时间域显然不同于x(n)所在的离散时间域，所以把它称为复倒频谱域。

称为x(n)的复倒频谱域，简称为复倒谱。如果在傅里叶变换时只考虑实数部分，那么得到的

称为x(n)的倒频谱，简称倒谱。

离散余弦变换(Discrete Cosine Transform，DCT)具有信号频谱分量丰富、能量集中，且不需要对语音相位进行估算等优点，能在较低的运算复杂度下取得较好的增强效果。

设x(n)是N个有限值的一维实数信号序列，n＝0,1,…,N-1，DCT的完备正交归一函数是：

式子中，n＝0，1，…，N-1，k＝0，1，…，N-1，α(k)的定义为：

X(k)可以用系数矩阵C来表示：

在计算MFCC使用求DCT倒谱的原理和上述求FFT倒谱的原理类似，把Mel滤波器的能量求对数后在进行离散余弦变换(DCT)：

式子中，S(i，m)是由4)中求出的Mel滤波器能量；m是指第m个Mel滤波器(共M个)；i是指第i帧；n是DCT后的谱线。这样就计算出了第i帧的mfcc参数。

使用DCT的原因是因为Mel滤波器的子带有重叠，所以相邻子带能量有一定冗余。可以采用Karhunen-Loeve变换去相关并选择主成分进行降维，也就是主成分分析。这个变换比较复杂，DCT变换是K-L变换的一个很好近似。通过DCT变换对Mel滤波器的子带能量特征进行去相关并降维。

在MFCC应用过程中通过大量的实验发现，MFCC的各个分量对音频识别率的贡献是不一样的。在音频采集过程中，中间部分的MFCC分量相对低阶和高阶MFCC分量来说，不易受到噪声影响，具有很好的鲁棒性^[24]。现在有的方法是使用半升正弦函数。对MFCC系数中某些谱线增强，加权函数如式(3-13)：

其中，Q表示DCT的阶数，一般都取Q＝12。然后进行归一化权重：

w′(n)＝w(n)/max(w) n＝1，2，...，Q (3-14)计算经过加权增强以后的MFCC参数：

mfcc′(i，n)＝mfcc(i，n)*w′(n) (3-15)

6)标准的倒谱参数MFCC只反映了音频的静态特性，因为它将每一帧语音信号认为是孤立的，提取本帧MFCC参数时只考虑本帧的音频特征。但是实际当中，音频信号是连续变化的，每帧音频都不是孤立的。音频信号的动态特征可以用这些静态特征的差分来描述，也即声学特征在相邻帧间的变化情况。实验证明：把动、静态特征结合起来才能有效提高系统的识别性能。一阶差分参数的计算可以采用下面的公式：

其中

或者取其他固定值。在本文中，C取1/3。

Mel频率倒谱系数在电机音频特征提取中的应用：

本文选择Mel滤波器组的最高频率为

最低频率f_l＝0。本文选择的采样频率为48000Hz，即f_h＝24000。Mel滤波器的个数p＝16，即在f_l和f_h之间一共有16个Mel滤波器。Mel滤波器设计时选用汉明窗函数。得到16个Mel滤波器的频率响应曲线如图6所示：

由上文可知，MFCC是将Mel滤波器的能量经过余弦变换以后得到的系数，相当于将Mel滤波器能量变换到余弦空间。本文选择余弦变换的阶数为12，即余弦变换矩阵C是一个12×16的矩阵：

其中，j＝1，2，3，...，12，k＝0，1，2，3，...，p-1。p表示Mel滤波器的个数，本文中p＝16。

本文除了选择标准的MFCC参数之外，为了描述音频的动态特征，还计算一阶差分MFCC参数，计算公式如下：

其中i表示帧号，n表示DCT后的谱线，也即MFCC的参数序号。然后将一阶微分参数放入标准MFCC参数之后。标准MFCC参数一共12维，加上12维动态MFCC参数，每一帧信号可以提取24维动态参数。由于首尾两帧无法计算动态MFCC参数，故在计算m帧音频的MFCC参数时，理论上应该得到m×24大小的MFCC参数矩阵。其中第i行代表第i帧音频的MFCC参数。第j列表示第j个MFCC参数随时间变化趋势。由于第0、1、m-1、m这四帧无法计算动态MFCC参数，而实际中音频的帧数一般比较多，所以在最后得到的MFCC参数中，去掉这四帧，剩下一个(m-4)×24大小的矩阵。

本文设定信号每帧时长L为1s，两帧之间的帧移M＝0.5s。选取每个电机音频样本的前60帧，计算每帧信号的MFCC参数。以帧数为横坐标，绘制每个MFCC参数随帧数变化曲线图。因为每一帧都是按照时间顺序依次获取得到，故每个MFCC参数随帧数变化的曲线实质上就是每个MFCC参数随时间变化的曲线。

MFCC参数曲线能够有效区分合格信号和故障信号，如图7所示。

本说明书实施例所述的内容仅仅是对发明构思的实现形式的列举，本发明的保护范围不应当被视为仅限于实施例所陈述的具体形式，本发明的保护范围也及于本领域技术人员根据本发明构思所能够想到的等同技术手段。

Claims (3)

1.基于Mel频率倒谱系数的电机异音检测方法，包括以下步骤：

步骤1、设定采样频率和采样时长t，将电机处于空载状态下进行音频信号采集x(n)，x(n)＝x₁(n)*x₂(n)，x(n)为音频信号，x₁(n)为激励信号，x₂(n)为物体单位冲击响应函数；

步骤2、对音频信号进行分帧加窗处理，设置每帧的帧长L和相邻两帧的重叠长度M₁，将音频信号分为N帧信号；

步骤3、使用基于Mel频率倒谱系数提取音频特征：

步骤3-1、对每一帧信号进行FFT，将时域数据转变为频域数据：

X(i，k)＝FFT[x(i，n)]，i是指第i帧，k是频率，单位Hz；n是DCT后的谱线；

步骤3-2、对每一帧FFT后数据的X(i，k)计算谱线能量：E(i，k)＝[X(i，k)]²；

步骤3-3、计算每帧信号通过Mel滤波器组之后，在各个滤波器中的能量：

S(i，m)是Mel滤波器能量，H_m(k)为频响，i是指第i帧，m是指第m个Mel滤波器，M是指M个Mel滤波器；

步骤3-4、计算DCT倒谱；

对每一帧信号，对每个Mel滤波器的能量求对数后再进行离散余弦变换

S(i，m)是由步骤3-3中求出的Mel滤波器能量；m是指第m个Mel滤波器，共M个；i是指第i帧；n是DCT后的谱线；

步骤4、计算出当前音频信号每一帧的mfcc参数，得到每个mfcc参数随帧数的变化曲线；

步骤5、选取N个合格的电机音频样本，重复步骤2-4，对于每个电机音频样本，得到以帧数为横坐标每个MFCC参数随帧数变化曲线图；

步骤6、进行电机异音诊断时，重复步骤1-4获得电机音频的音频样本特征，提取每个mfcc参数随帧数变化曲线的平均值作为特征，输入学习模型modelA，判断是否为故障样本；然后提取每个mfcc参数随帧数变化曲线的峰峰值作为一组特征，输入学习模型modelB，判断是否为故障样本；只有当两个学习模型都判断为合格时，认定该样本为合格样本，否者都判断为故障样本。

2.如权利要求1所述的基于Mel频率倒谱系数的电机异音检测方法，其特征在于：步骤3-4中加入mfcc参数的动态特征，

计算公式为

其中，i表示帧号，n表示DCT后的谱线，k是频率，单位Hz，然后将一阶微分参数放入标准MFCC参数之后；

标准静态mfcc参数一共12维，加上12维动态mfcc参数，每一帧信号提取24维参数；

在计算m帧音频的mfcc参数时，得到(m-4)×24大小的参数矩阵，其中第i行代表第i帧音频的mfcc参数，第j列表示第j个mfcc参数随时间变化趋势，参数矩阵不包含第1帧、第2帧、第m-1帧和第m帧的动态mfcc参数。

3.如权利要求1-2之一所述的基于Mel频率倒谱系数的电机异音检测方法，其特征在于：步骤2中，帧长L＝1S，重叠长度M₁＝0.5s。

Images