一种减振器总成异响检测方法与流程

1.本发明属于减振器异响检测技术领域，更具体地，涉及一种减振器总成异响检测方法。


背景技术：

2.对于汽车，尤其是轿车来说，减振器总成异响通常是客户抱怨的重点问题之一。各大汽车厂商投入了大量的人力、物力、财力对减振器总成异响进行检测，以防止不良产品流向市场。目前，对车辆减振器总成异响的检测有两种方法。
3.第一种方法为主观评价法，依靠有经验的人员对减振器总成噪声进行检查，判断有无异响。这种方式主观性强(对同一减振器总成噪声，不同的人有不同的判断结果)，工作效率低下，且人员长时间工作时易引起听觉疲劳从而导致错检、漏检等问题，不利于产品的自动化生产。
4.第二种方法为频谱分析法，采集减振器总成台架试验噪声信号，对噪声总体进行频谱分析，依据频谱结果判断减振器总成有无异响。这种方式有很大的局限性，由于减振器总成台架试验噪声不是稳态的，只能对总体噪声信号进行时频分析，得到的结果不直观，无法对有无异响提供有力的判据，致使该方法不实用。


技术实现要素：

5.针对现有技术的以上缺陷或改进需求，本发明提出了一种减振器总成异响检测方法，该方法对减振器总成异响识别率高，工作效率高，能对减振器总成是否存在异响做出快速、有效的判断，并能对存在异响的减振器总成作出故障诊断，判断是何种原因导致的异响。
6.为实现上述目的，本发明提供了一种减振器总成异响检测方法，包括：
7.采集待测减振器总成台架试验噪声信号，并对待测减振器总成台架试验噪声信号进行预处理；
8.对预处理后的数据进行频谱分析得到能量谱数据，对能量谱数据进行滤波；
9.对滤波后的数据求对数，进行离散余弦变换得到待测减振器总成噪声信号的mfcc；
10.由待测减振器总成噪声信号的mfcc生成待测减振器总成异响特征距离矩阵；
11.将待测减振器总成异响特征距离矩阵与减振器总成异响检测模板集进行匹配，生成待测减振器总成异响特征匹配矩阵；
12.基于待测减振器总成异响特征匹配矩阵对待测减振器总成异响特征进行判定。
13.在一些可选的实施方案中，所述对待测减振器总成台架试验噪声信号进行预处理，包括：
14.对待测减振器总成台架试验噪声信号进行分帧，将n个连续采样点集合成一帧，将第i帧数据记为x0(i)，n为帧长，n为2的指数次幂；
15.对各帧数据x0(i)进行加窗处理，得到加窗后的第i帧数据x(i)。
16.在一些可选的实施方案中，所述对预处理后的数据进行频谱分析得到能量谱数据，包括：
17.对加窗后的各帧数据x(i)进行fft变换，得到加窗后的第i帧数据x(i)中第k条谱线的频谱x(i,k)；
18.对第i帧数据x(i)中第k条谱线的频谱x(i,k)取模，得到第i帧数据x(i)中第k条谱线的频谱x(i,k)的幅值谱a(i,k)；
19.对第i帧数据x(i)中第k条谱线的频谱x(i,k)的幅值谱a(i,k)取平方，得到第i帧数据x(i)中第k条谱线的频谱x(i,k)的能量谱e(i,k)。
20.在一些可选的实施方案中，所述对能量谱数据进行滤波，包括：
21.根据待测减振器总成台架试验噪声信号的采样频率，得到待测减振器总成台架试验噪声信号的频率上限，待测减振器总成台架试验噪声信号的频率下限设置为与人耳所能感知的噪声的频率下限相等，由频率上限得到mel三角滤波器的最大mel值，由频率下限得到mel三角滤波器的最小mel值；
22.由最大mel值、最小mel值、采样频率和数据帧长度计算mel三角滤波器组各滤波器的中心频率，由各滤波器的中心频率得到各滤波器的传递函数；
23.每一帧频谱数据通过mel滤波器后的频谱能量f(i,m)为能量谱e(i,k)与mel三角滤波器组传递函数的乘积和，1≤m≤m，m为滤波器个数。
24.在一些可选的实施方案中，所述对滤波后的数据求对数，进行离散余弦变换得到待测减振器总成噪声信号的mfcc，包括：
25.由得到待测减振器总成噪声各帧数据的mfcc，t＝1,2
……
m，m，c(i,t)为第i帧数据的mfcc的第t个系数，m为三角滤波器序号，表示第m个三角滤波器。
26.在一些可选的实施方案中，所述减振器总成异响检测模板集的生成包括：
27.采集正常减振器总成和存在异响的减振器总成的噪声信号，对进行噪声信号采集的减振器总成进行编号，并记录故障状态；
28.对采集的各噪声信号进行预处理，对预处理后的各噪声信号进行频谱分析得到能量谱数据，对能量谱数据进行滤波；
29.对滤波后的数据求对数，进行离散余弦变换得到各噪声信号的mfcc，组成r个mfcc矩阵，每个矩阵的每一列代表相应减振器总成某帧数据的m个mel倒谱系数，每个矩阵的列数不尽相同，列数表示相应减振器总成噪声的帧数，记为col_num，col_num(n)表示编号为n的减振器总成的mfcc矩阵的列数，r表示采集的减振器总成的噪声信号的个数；
30.将r个mfcc矩阵组合成一个三维矩阵，记为m_template，m_template(n,i2,j2)表示编号为n的减振器总成mfcc矩阵中的第i2行，第j2列元素，m_template为减振器总成异响检测模板集。
31.在一些可选的实施方案中，所述由待测减振器总成噪声信号的mfcc生成待测减振器总成异响特征距离矩阵，包括：
32.将待测减振器总成噪声各帧数据的mfcc记为m_test；
33.由采用余弦相似度来表征两帧mfcc间的距离度量参数，其中，similarity_0为一三维矩阵，1≤n≤r，1≤i2≤t，1≤j2≤col_num(n)，similarity_0(n,i2,j2)表示待测减振器总成第i2帧mfcc与检测模板集中编号为n的减振器总成第j2帧mfcc间的余弦相似度，t表示待测减振器总成的帧数量，m_test(t,i2)表示待测减振器总成噪声第i2帧数据的第t个mel倒谱系数，m_template(n,t,j2)表示检测模板集中编号为n的减振器总成第j2帧数据的第t个mel倒谱系数；
34.由similarity＝e-similarity_0
对similarity_0进行变换得到similarity，以改变similarity_0的趋势，使similarity越小，两向量越贴近，最后由similarity构成待测减振器总成异响特征距离矩阵。
35.在一些可选的实施方案中，所述生成待测减振器总成异响特征匹配矩阵，包括：
36.根据待测减振器总成mfcc与检测模板集中编号为n的减振器总成mfcc间的距离矩阵similarity(n,:,:)，计算待测减振器总成异响特征匹配矩阵m_match(n,:,:)，其中，m_match(n,1,1)＝similarity(n,1,1)，m_match(n,i2,j2)＝min{m_match(n,i2-1,j2)+similarity(n,i2,j2)，m_match(n,i2-1,j2-1)+similarity(n,i2,j2)，m_match(n,i2,j2-1)+similarity(n,i2,j2)}。
37.在一些可选的实施方案中，所述基于待测减振器总成异响特征匹配矩阵对待测减振器总成异响特征进行判定，包括：
38.取出待测减振器总成与检测模板集中编号为n的减振器总成异响特征匹配矩阵中的右下角元素m_match(n,t,col_num(n))，将此元素记为distance(n)，distance(n)表示待测减振器总成与检测模板集中编号为n的减振器总成的总匹配距离；
39.遍历n，1≤n≤r，求取所有元素m_match(n,t,col_num(n))，得到数组distance的所有值；
40.取数组distance的最小值所对应的序号，记为match_no，其表示待测减振器总成与检测模板集中编号为match_no的减振器总成异响特征最为匹配，与match_no对应的检测模板集中的减振器总成故障状态即为待测减振器总成检测结果。
41.总体而言，通过本发明所构思的以上技术方案与现有技术相比，能够取得下列有益效果：
42.本发明提供了一种基于mfcc和dtw算法的减振器总成异响检测方法。相比传统的减振器总成异响检测方法，本发明中的方法由于考虑了人耳的听觉特性，同时对噪声进行了细致深入的分析，得到其mfcc参数，更加深刻地描述其音色特征。运用dtw算法计算待测减振器总成噪声的mfcc序列与模板集中各减振器总成的mfcc序列的距离，通过最小距离将待测减振器总成与模板集中的减振器总成进行匹配，从而判断是否存在异响，存在何种异响。该方法对减振器总成异响识别率高，工作效率高，能对减振器总成是否存在异响做出快
速、有效的判断，有利于实现减振器总成的自动化检测。该方法能对存在异响的减振器总成作出快速的故障诊断，判断是何种原因导致的异响，方便减振器总成的返修，避免部品的浪费。
附图说明
43.图1是本发明实施例提供的一种减振器总成异响检测方法流程图；
44.图2是本发明实施例提供的一种mel刻度和频率(hz)的关系曲线；
45.图3是本发明实施例提供的一种名义mel三角滤波器组图形；
46.图4是本发明实施例提供的一种实际运算mel三角滤波器组图形。
具体实施方式
47.为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。
48.本发明提供了一种基于mel频率倒谱系数(mel frequency cepstrum coefficient，mfcc)和动态时间规整(dynamic time warping，dtw)算法的减振器总成异响检测方法。相比传统的减振器总成异响检测方法，本发明中的方法由于考虑了人耳的听觉特性，同时对噪声进行了细致深入的分析，得到其mfcc参数，更加深刻地描述其音色特征。运用dtw算法计算待测减振器总成噪声的mfcc序列与模板集中各减振器总成的mfcc序列的距离，通过最小距离将待测减振器总成与模板集中的减振器总成进行匹配，从而判断是否存在异响，存在何种异响。该方法对减振器总成异响识别率高，工作效率高，能对减振器总成是否存在异响做出快速、有效的判断，并能对存在异响的减振器总成作出故障诊断，判断是何种原因导致的异响。
49.如图1所示是本发明实施例提供的一种减振器总成异响检测方法的流程图，在图1所示的方法中，包括以下步骤：
50.s110：采集减振器总成台架试验噪声信号；
51.采集减振器总成台架试验的噪声信号{s0(q)}(1≤q≤n
s0
,n
s0
为信号{s0(q)}的长度)。本实施例中，采样频率为8000hz，采样时长为减振器总成在典型恶劣路面(如减振器总成异响经常发生的路面)上的一段位移行程。采集到的数据以wave格式保存。
52.在本实施例中，上述数据采集环境、条件应一致；须在安静的环境下进行采集，最好在消声室进行采集。
53.s120：读取减振器总成台架试验噪声信号，并进行预处理；
54.在本实施例中，预处理过程包括分帧和加窗，其中，分帧的目的是保证噪声分析有足够高的时间分辨率，分帧可以通过以下方式实现：
55.将n个连续采样点集合成一个数据分析单位，称为帧，n即为帧长。n须为2的指数次幂，即n＝2z(z为一自然数)，以方便后续对数据进行频谱分析。本实施例中n取值为256，减振器总成噪声的采样频率为8000hz，则一帧对应的时间长度为256/8000
×
1000＝32ms。分帧后，第i帧数据记为x0(i)。
56.其中，加窗(施加窗函数)的目的是减少频谱泄露，本实施例中，窗函数选为汉宁窗(hanning)，其表达式为：
57.w(p)＝0.5*(1-cos(2*π*p/p))，0≤p≤p
ꢀꢀꢀꢀꢀꢀꢀꢀ
(1)
58.其中，p，p为整数，p为自变量，w(p)为因变量，为窗函数的幅度，p为常数(窗函数的总长度l＝p+1，l须等于帧长n)。
59.对各帧数据进行加窗，公式如下：
60.x(i)＝x0(i).*w(p)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(2)
61.其中，x(i)为加窗后的第i帧数据，.*为数量乘积，即数列中各对应元素相乘。
62.s130：对预处理后的数据进行频谱分析；
63.对各帧数据进行fft变换，得到各帧的频谱：
64.x(i,k)＝fft[x(i)]
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(3)
[0065]
其中，x(i)表示预处理后的第i帧数据；x(i,k)表示第i帧频谱数据中第k条谱线的频谱。
[0066]
其中，频谱为复数，对其取模，得到幅值谱：
[0067]
a(i,k)＝|x(i,k)|
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(4)
[0068]
a(i,k)表示第i帧频谱数据中第k条谱线的幅值。
[0069]
对a(i,k)取平方，得到数据的能量谱：
[0070]
e(i,k)＝|x(i,k)|2ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(5)
[0071]
e(i,k)表示第i帧频谱数据中第k条谱线的能量。
[0072]
s140：对能量谱数据进行滤波；
[0073]
在本实施例中，步骤s140包括两小步s141和s142。
[0074]
其中，s141：计算mel三角滤波器组数据；
[0075]
由于人耳对声音的感知程度，与频率(hz)不成线性关系，但在mel刻度上，人对音调的主观感知与mel刻度则为线性关系。
[0076]
mel刻度和频率(hz)的关系如下：
[0077]
mel＝2595*log
10
(1+f0/700)
ꢀꢀꢀꢀꢀꢀꢀ
(6)
[0078]
其中，f0为频率(hz)，mel为mel值。
[0079]
mel刻度和频率(hz)的关系曲线如图2所示。
[0080]
由式(6)及图2可以看出，mel刻度在低频(hz)部分的分辨率高，在高频(hz)的分辨率低，与人耳的听觉特性是相符的，mel刻度模拟了人耳的听觉特性。同时，结合人耳的掩蔽效应，设计了mel三角滤波器组，这些滤波器在低频段较密集，在高频段较稀疏。
[0081]
第m个三角滤波器的传递函数为：
[0082][0083]
其中，1≤m≤m，m为滤波器个数；f(m)为第m个滤波器中心频率的名义序号(因f(m)
不一定为整数，故称为名义序号)，f(m-1)为第m-1个滤波器中心频率的名义序号，f(m+1)为第m+1个滤波器中心频率的名义序号，hm(k)为第m个滤波器中第k条谱线的幅值。
[0084]
mel三角滤波器组各滤波器的端点频率及中心频率名义序号f(m)可依照如下方式计算：
[0085]
对于本实施例中减振器总成噪声采样频率8000hz，依据香农采样定理，其分析频率上限为8000/2＝4000hz，依据式(6)，将4000hz代入其中，求出mel值，记为mel_max。
[0086]
选定待测减振器总成台架试验噪声的频率下限，其设置为与人耳所能感知的噪声的频率下限相等，如20hz，代入式(6)，求出mel值，记为mel_min。
[0087]
运用下式(8)计算mel三角滤波器组各滤波器的中心频率(端点频率)名义序号f(j)：
[0088][0089]
其中，0≤j≤(m+1)，m为滤波器个数；n为数据帧的长度；fs为减振器总成噪声采样频率。本实施例中，n＝256，fs＝8000hz。
[0090]
则第一个滤波器的中心频率名义序号为f(1)，左端点频率名义序号为f(0)，右端点频率名义序号为f(2)；第j个滤波器的中心频率名义序号为f(j)，左端点频率名义序号为f(j-1)，右端点频率名义序号为f(j+1)；最后一个滤波器m的中心频率名义序号为f(m)，左端点频率名义序号为f(m-1)，右端点频率名义序号为f(m+1)。
[0091]
图3为n＝256，fs＝8000hz，m＝24的mel三角滤波器组图形。图3实际上为名义mel三角滤波器组图形，实际运算时的滤波器组图形是对名义滤波器组图形在能量谱e(i,k)各频谱处采样得到的。图4为n＝256，fs＝8000hz，m＝24的实际运算mel三角滤波器组图形。
[0092]
s142：对频谱数据进行mel滤波mel处理；
[0093]
每一帧频谱数据通过mel滤波器后的频谱能量为能量谱e(i,k)与mel三角滤波器组传递函数hm(k)的乘积和：
[0094]
f(i,m0＝∑ke(i,k)*hm(k),1≤m≤m
ꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(9)
[0095]
其中，m为滤波器个数，典型值为24。
[0096]
至此，每一帧的数据长度由n个变为m个，数据维数大幅缩减。
[0097]
s150：对滤波后的数据求对数，进行离散余弦变换；
[0098][0099]
其中，m为滤波器个数，c即为mfcc(mel倒谱系数)，c(i,t)为第i帧数据的mfcc的第t个系数(t为整数，为离散余弦变换变量，作用如同fft中的谱线数k)，m为三角滤波器序号，即第m个三角滤波器。
[0100]
s160：生成减振器总成异响检测模板集；
[0101]
按照步骤s110，采集正常减振器总成和存在异响的减振器总成的噪声信号，对正
常减振器总成，可采集1～2个总成的噪声信号，对存在异响的减振器总成，由于其故障原因各不相同(如减振器阀系故障、减振器分液阀力矩不足、减振器压缩阀力矩不足、减振器油液黏度不合规格等原因导致的异响)，其异响特征也不尽相同，对每种故障的减振器总成，可采集1～2个总成的噪声信号。对上述进行噪声信号采集的减振器总成进行编号，并建立数组fault_state记录其故障状态，如fault_state(n)＝fault_stype表示编号为n的减振器总成的故障类型为fault_stype。其中，fault_stype为故障代码变量，可直接命名为数字，如0代表“正常”，1代表“减振器压缩阀力矩不足
”……
设共采集r个减振器总成的噪声信号，则可得到r段音频。
[0102]
按照步骤s120至步骤s150，对这r段音频分别提取其mfcc特征值，可得到r个矩阵，每个矩阵的每一列代表相应减振器总成某帧数据的m个mel倒谱系数，即每个矩阵的行数均是m，每个矩阵的列数不尽相同，列数即相应减振器总成噪声的帧数，记为col_num，col_num(n)表示编号为n的减振器总成的mfcc特征值矩阵的列数，将这r个矩阵组合成一个三维矩阵，记为m_template，m_template(n,i2,j2)表示编号为n的减振器总成mfcc特征值矩阵中的第i2行，第j2列元素，m_template即为减振器总成异响检测模板集。
[0103]
s170：生成待测减振器总成异响特征距离矩阵；
[0104]
依据步骤s110至步骤s150，可以得到待测减振器总成噪声各帧数据的mfcc，如得到t帧mfcc，形成m行t列矩阵，记为m_test；
[0105]
采用余弦相似度来表征两帧mfcc参数间的距离度量参数。余弦相似度是描述向量间靠近程度的参数。在本实施例中，运用余弦相似度描述两帧mfcc参数的贴近程度，即各帧音色(mfcc参数可理解为音色)的贴近程度，计算公式如下：
[0106][0107]
其中，similarity_0为一三维矩阵，1≤n≤r，1≤i2≤t，1≤j2≤col_num(n)，similarity_0(n,i2,j2)表示待测减振器总成第i2帧mfcc参数与检测模板集中编号为n的减振器总成第j2帧mfcc参数间的余弦相似度，其值范围为[-1,1]，其值越大，表示两向量越贴近；其值越小，表示两向量越不相似。
[0108]
由于下述运用dtw(动态时间规整)进行匹配计算时，其算法要求两向量越贴近，其距离越小，这个要求与变量similarity_0不一致(变量similarity_0越大，表示两向量越贴近)。因此，需要对变量similarity_0进行变换。变换公式如式(12)。
[0109]
similarity＝e-similarity_0
ꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(12)
[0110]
式(12)的作用改变similarity_0的趋势(使similarity越小，两向量越贴近)，e为自然对数函数的底数。similarity即为待测减振器总成异响特征距离矩阵。
[0111]
s180：生成待测减振器总成异响特征匹配矩阵；
[0112]
运用dtw(动态时间规整)进行匹配计算。dtw是dynamic time warping的缩写，用来计算两个时间序列之间的相似性。在时间序列分析中，有时两段时间序列的长度可能并不相等，传统的欧几里得距离无法计算两个时间序列之间的距离(或相似性)，而dtw通过把时间序列进行延伸和缩短，能有效地衡量两个时间序列性之间的相似性。
[0113]
根据待测减振器总成mfcc参数与检测模板集中编号为n的减振器总成mfcc参数间
的距离矩阵，记为similarity(n,:,:)，计算其异响特征匹配矩阵，记为m_match(n,:,:)。计算方法如下：
[0114]
首先进行初始化，令m_match(n,1,1)＝similarity(n,1,1)，然后按式(13)进行迭代计算：
[0115]
m_match(n,i2,j2)＝min{m_match(n,i2-1,j2)+similarity(n,i2,j2)，
[0116]
m_match(n,i2-1,j2-1)+similarity(n,i2,j2)，
[0117]
m_match(n,i2,j2-1)+similarity(n,i2,j2)}
ꢀꢀꢀꢀꢀ
(13)
[0118]
式(13)中m_match(n,i2,j2)的计算方法是有顺序的，可依次计算其第一列、第二列、
……
第col_num(n)列。也可依次计算其第一行、第二行、
……
第t行。
[0119]
式(13)中，只对i2-1≥1，j2-1≥1的数值进行取小计算。如求m_match(n,1,2)时，因i2＝1，j2＝2，i2-1＝0，故只需取m_match(n,1,2)＝min{m_match(n,i2,j2-1)}＝min{m_match(n,1,1)}。因为m_match(n,i2-1,j2)，m_match(n,i2-1,j2-1)这两个元素不存在。
[0120]
式(13)的含义是待测减振器总成某帧mfcc参数mfcc_a可能与检测模板集中编号为n的减振器总成某帧mfcc参数mfcc_b一一对应，也可能mfcc_a对应编号为n的减振器总成多帧mfcc参数，也可能mfcc_b对应待测减振器总成多帧mfcc参数。帧之间的匹配必须保持连续，不能出现跨帧匹配的情况。
[0121]
依据式(13)求出待测减振器总成与检测模板集中所有减振器总成的异响特征匹配矩阵，即得到异响特征匹配矩阵m_match。
[0122]
s190：对待测减振器总成异响特征进行判定。
[0123]
取出待测减振器总成与检测模板集中编号为n的减振器总成异响特征匹配矩阵中的右下角元素，即元素m_match(n,t,col_num(n))，将此元素记为distance(n)，表示待测减振器总成与检测模板集中编号为n的减振器总成的总匹配距离；
[0124]
遍历n，1≤n≤r，求取所有元素m_match(n,t,col_num(n))，得到数组distance的所有值；
[0125]
取数组distance的最小值所对应的序号，记为match_no，则：
[0126]
match_no＝argminn(distance(n))
ꢀꢀꢀꢀꢀꢀꢀꢀ
(14)
[0127]
即待测减振器总成与检测模板集中编号为match_no的减振器总成异响特征最为匹配，fault_state(match_no)即为待测减振器总成检测结果，即若fault_state(match_no)表示正常，则待测减振器总成为正常减振器；若fault_state(match_no)表示减振器压缩阀力矩不足，则待测减振器总成有相应故障，其余结果类推。
[0128]
需要指出，根据实施的需要，可将本技术中描述的各个步骤/部件拆分为更多步骤/部件，也可将两个或多个步骤/部件或者步骤/部件的部分操作组合成新的步骤/部件，以实现本发明的目的。
[0129]
本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。