内燃机KVMD MHD振动谱图表征与TD 2DPCA编码诊断方法研究

岳应娟 孙钢 蔡艳平 王旭 牟伟杰

摘要： 为了直接对内燃机振动谱图像进行诊断识别，提出一种基于改进变分模态分解（VMD）、MargenauHill（MHD）时频分析与双向二维主成分分析（Twodirectional， Twodimensional PCA，TD2DPCA）的内燃机振动谱图像识别诊断方法。该方法首先针对VMD分解过程中的层数选取问题，提出了一种中心频率筛选的VMD分解层数改进方法（KVMD），然后将内燃机振动信号利用KVMD分解成一组单分量模态信号，并对生成的各个单分量信号进行MHD处理后表征成振动谱图像；在此基础上，对生成的内燃机KVMDMHD振动谱图像采用双向二维主成分分析形成编码矩阵，并采用最近邻分类器（KNNC）对上述编码矩阵直接进行模式识别，以实现内燃机振动谱图像的自动诊断。最后，将该方法应用在气阀机构4种工况下的缸盖表面振动信号诊断实例中，结果表明：该方法不仅改进了传统图像模式识别中的特征参数提取方法，而且能很好地消除时频分布中的交叉干扰项，使各时频分量物理意义明确，能有效诊断出内燃机气阀机构故障，为内燃机振动诊断探索了一条新途径。关键词： 故障诊断； 内燃机； 时频分布； 特征提取； 双向二维主成分分析

中图分类号： TH165+.3； TK428文献标志码：A文章编号：10044523（2017）04068809

DOI：10.16385/j.cnki.issn.10044523.2017.04.021

引言

内燃机缸盖振动信号中包含着丰富的信息，由于其测量的简单方便，分析诊断的不解体和实时性，目前一直是内燃机故障诊断和状态监测的研究前沿和热点。国内外学者针对内燃机的振动响应信号进行了深入的研究，如文献[1]将图像分割理论引入柴油机故障诊断中，提出一种基于时频谱图、图像分割和模糊模式识别的柴油机故障诊断方法；文献[2]提出一种基于局部均值分解边际谱和马氏距离的故障诊断方法；文献[3]将极坐标应用于柴油机燃烧状态的监测，有效地提取了柴油机燃烧特征；文献[4]将高阶累积量与图像纹理特征提取方法相结合，有效提取了柴油机振动信号的故障特征。

根据内燃机的构造和工作机理，气阀与气阀座会发生周期性的冲击作用，若气阀机构有故障， 其故障信息必然会在缸盖振动信号中反映出来[56]。然而内燃机特征信号相互重叠和混淆、特征频率难以确定，还没有形成一个“针对不同故障，采用何种时频分析方法，如何提取特征参数”的共识。究其原因是内燃机振动响应信号十分复杂，既有旋转运动，又有往复运动，且运动部件多，耦合比较严重，具有较强的非线性、非平稳时变等特征[7]。

为有效解决内燃机振动响应信号强耦合、弱故障特征信息提取难题，提出了一种基于改进VMD（KVMD）的MHD时频振动谱图生成，TD2DPCA图像特征参数提取的内燃机故障诊断新方法。KVMDMHD时频分析法有效抑制了MHD分布中的交叉干扰项，并保留了其优良的时频聚集特性，能够准确刻画出内燃机振动信号的时频信息，使各时频分量具有实际物理意义。对生成的KVMDMHD振动谱图，直接采用TD2DPCA提取特征参数的方法，避免了在利用图像分析技术进行特征参数提取时，不同图像特征指标的选择或只是提取图像的单一特征量作为特征参数造成的重要的时频信息遗漏，可对不同的图像自适应地计算图像的特征参数，数据降维效果明显。

最后使用文中方法对内燃机气门间隙的4种工况信号进行了分析和特征参数提取，结合最近邻分类器（KNNC）进行故障诊断分类，并与基于MHD时频分析的二维非负矩阵分解（2DNMF）和双向二维线性判别分析方法（TD2DLDA）特征提取算法进行了对比。

1基于改进变分模态分解的MHD时频分布〖*2〗1.1改进的变分模态分解（KVMD）变分模态分解（Variational Mode Decomposition，VMD）是2014年由Dragomiretskiy等提出的一种新的自适应信号处理方法[8]。對信号进行VMD分解时首先预设分解层数K，信号经过VMD被分解成K个本征模态分量（IMF），每个IMF都可以表示为一个调幅调频uk（t）信号。因此K值选取的恰当与否，直接决定了分解结果的好坏。K值选取过小，对信号的分解不彻底；K值选取过大，会出现过分解现象。经研究发现每个IMF都存在着一个中心频率ωk（t），K值与ωk（t）有着密切的关系，因此本文对通过中心频率对分解层数K进行了优化。

第4期岳应娟，等： 内燃机KVMDMHD振动谱图表征与TD2DPCA编码诊断方法研究振 动 工 程 学 报第30卷KVMD算法的主要步骤为：

Step1初始化K值（K≥2，由于内燃机频带较宽，取K=3）。

Step2对信号进行VMD分解，得到K个IMF分量和每个IMF分量的中心频率ωk（t）。

Step3用前一个IMFk-1分量的中心频率ωk-1（t）比上后一个IMFk分量的中心频率ωk（t），得到一组频率比值λ1，λ2，…，λK-1（λk=ωk+1/ωk，k=1，2，…，K-1）。

Step4设定过分解阈值θ（根据内燃机频带特点取θ=1.1）。当λk>θ时，认为VMD分解不彻底，令K=K+1，重复Step2～Step3。

Step5当λk≤θ时可判定为IMFk和IMFk+1频率混叠，VMD出现了过分解，因此得出结果K=K-1，并输出其分解结果。

1.2基于改进变分模态分解的MHD时频分布

MHD时频分布[9]是一种非平稳信号分析的工具，具有真边缘性、弱支撑性、平移不变性等优良性质。对给定信号x（t）的时频分布p（t，f），Cohen给出一般形式的表达式p（t，f）=∫+∞-∞∫+∞-∞∫+∞-∞xu+τ2x*u-τ2·

τ，ve-j2π（tv+τf-uv）dudtdv（1）式中τ，v是核函数。若≡0，则为WignerVille分布，当=cos（πτv）时，即为MHD分布pMH（t，f）=∫+∞-∞∫+∞-∞∫+∞-∞xu+τ2x*u-τ2·

cos（πτv）e-j2（tv+τf-uv）dudτdv（2）由于双线性核函数的引入，使多个分量在时频平面发生耦合产生了交叉项，MHD时频分布很难将有多个频率成分的信号表示清楚[10]。MHD分布的交叉项是以两个自项构成的矩形对角线顶点，若两个自项分布位于同一频率或同一时间时 ，则自项和交叉项重叠。若对MHD进行加窗处理，就得到了伪魏格纳分布（PMHD）。pPMH（t，f）=∫+∞-∞∫+∞-∞∫+∞-∞h（τ）xu+τ2x*u-τ2·

cos（πτv）e-j2（tv+τf-uv）dudτdv（3）式中h（τ）为窗函数。

信号x（t）的KVMDMHD时频分布定义为KVMDMHDt，f=∑Ki=1∫∞-∞fMHDIMFit，fdfMHDIMFit，fdf（4）KVMDMHD时频分布是利用了线性时频表示满足叠加原理的思想[11]。为消除交叉干扰项，可以将待分析的信号经KVMD分解成一组单分量信号IMF1，IMF2，…，IMFK，先对各个单频率分量信号单独进行MHD分析，这样在频域上就不会产生交叉干扰项，而位于同一频率的时域交叉项会与自项相互叠加，对自项有增强作用，对信号的分析有积极的作用；再将结果线性叠加，这样在保留了MHD时频分布的优良特性的同时，有效地消除了MHD的交叉项的干扰。

为分析该方法的性能，建立一个多分量仿真信号， 设x（t）是由3个原子复合而成，他们的位置分别位于（t1，Ω1）=（28，0.1），（t2，Ω2）=（28，0.4），（t4，Ω4）=（100，0.4），该信号的时域和频域波形如图1所示。

图1仿真信号

Fig.1Simulation Signal图2给出了仿真信号的MHD时频分布图，图中颜色表示信号在对应的时间和频率处的能量幅值大小与图像右侧的颜色标尺对应。仅从时频相平面图中已经很难分辨哪个是自项，那个是交叉项，通过对原始信号进行分析和比对，可以发现在（t2，Ω2）=（100，0.1）为交叉项，其余均为自项与交叉项的叠加。

图2MHD时频分布图

Fig.2Timefrequency distribution of MHD图3PMHD时频分布图

Fig.3Timefrequency distribution of PMHD

对MHD做加窗处理得到x（t）的PMHD分布图如图3所示。从图3可看出，这时交叉项得到抑制，但丢失了自项信号的一些细节，牺牲了时频分辨率，使生成的时频谱图不便于分析和理解。

对x（t）使用文中方法进行分析处理得到KVMDMHD的时频分布图。如图4所示，KVMDMHD时频分析已经去除了交叉项的干扰，使自项分辨的很清楚。

图5所示为2原子和4原子的仿真信号，及其MHD和KVMDMHD时频分布图。从图中可以看出，2原子信号的MHD时频分布在对角构成矩形的另一对角上存在交叉干扰项，图4KVMDMHD时频分布图

Fig.4Timefrequency distribution of KVMDMHD

4原子信号的MHD时频分布图为自项与干扰项的叠加；而两者的KVMDMHD时频分布图均能有效消除MHD分布中的交叉项的干扰。

图5多原子信号时频分布图

Fig.5Timefrequency distribution of multiatoms2TD2DPCA分解

传统的PCA方法需要将二维矩阵数据向量化，样本维数比较大，计算效率低下。2DPCA在特征提取上是直接利用二维投影的方法，数据量少，在提取特征上耗时也更少。但2DPCA仅在图像的行方向上进行运算，忽视了图像列中包含的信息，TD2DPCA[1213]将行和列两种图像信息融合到一个判别分析框架中，识别率得到提高，同时計算复杂度较低。

假设有C类模式：ω1，ω2，…，ωc，总共M个训练样本图像：A1，A2，…，AM，每个大小为m×n。Gt为训练样本总体散度矩阵Gt=1M∑Mi=1Ai-TAi-（5）式中=1M∑Mi=1Ai为训练样本的均值矩阵，可证Gt是n×n的非负定矩阵。

通过线性变换将图像矩阵Ai投影至X上从而获得特征向量Y=AiX（i=1，2，…，k），其中X表示n维单位化的列向量。投影方向X的选取准则是使得投影后的特征向量具有更好的可分性。定义准则函数J（X）=tr（Gt）=XTGtX（6）式中tr（Gt）为Gt的迹。

为了实现投影后得到的特征向量总体分散程度最大，即J（X）最大，需要寻找最优投影向量X。其实，Gt的最大特征值所对应的单位特征向量即为最优投影向量。因Gt为非负定矩阵，则有n个标准正交的特征向量，假定GtXi=λXi，（λ1≥λ2≥…≥λn≥0）（7）为了提高在多类样本情况下的区分性，单一的最优投影方向是不够的，取前d个最大特征值所对应的标准正交的特征向量作为最优投影矩阵P。假设P=[X1，X2，…，Xd]。对图像样本A，利用最优投影矩阵对其进行特征提取，获得相应的特征编码矩阵B，即B=AP。

对第1次提取的特征Bi（i=1，2，…，M）作为训练矩阵进行第2次特征提取，即将BTi作为Ai代入式（5），得到新的散布矩阵t=1M∑Mi=1Bi-Bi-T（8）式中=1M∑Mi=1Bi为首次提取特征后训练集的均值矩阵。

构造与式（6）相似的准则函数，求解t的前h个最大特征值所对应的标准正交的特征向量Z1，Z2，…，Zh，以此作为第2次特征提取的最优投影矩阵Q，则任一图像A经TD2DPCA算法提取的特征矩阵U为U=BT[Z1，Z2，…，Zh]=PTATQ=

[X1，X2，…，Xd]TAT[Z1，Z2，…，Zh]（9）特征矩阵U的维数为h×d，相比于2DPCA第1次提取出的特征维数为m×d，h远小于m，从而进一步压缩特征维数，提高了后续分类效率。

3内燃机故障诊断实例〖*2〗3.1内燃机智能故障诊断流程基于KVMDMHD与TD2DPCA的故障诊断方法对内燃机的故障诊断，共分为以下几个步骤：首先对采集到的内燃机振动信号进行KVMDMHD时频分析得到时频分布图像，然后采用TD2DPCA对时频图像进行分解得到图像的特征参数，最后用分类器对特征参数进行分类，完成对内燃机的故障诊断，具体方法的步骤如图6所示。

图6基于KVMDMHD与TD2DPCA的故障诊断方法的步骤

Fig.6Fault diagnosis method based KVMDMHD and TD2DPCA3.2内燃机实验工况

文中以6135型柴油机为研究对象，实验平台由柴油机、传动轴、电机和控制台4部分组成，如图7所示。取内燃机第2缸盖表面振动信号对内燃机进行故障诊断，采样频率25 kHz，转速为1500 r/min，测试过程中，内燃机空载运行。共设置了4种气门间隙状况，具体情况如表1所示。其中0.06，0.3和0.5 mm分别对应排气阀气门间隙过小、正常与过大，开口表示在气阀上开4 mm（长）×1 mm（宽）的口来模拟严重漏气。实验共采集内燃机气门4种工况（3种故障状态和1种正常状态）下各60种振动信号样本，总计240个。

图7试验平台

Fig.7Experimental platform表14种实验工况设置（单位：mm）

Tab.1Four states of IC engines valve train（Unit：mm）

状态编号进气门排气门10.300.3020.300.0630.300.5040.30开口4×13.3内燃机缸盖振动信號的KVMDPMHD时频分析根据6135柴油机的工作原理可知，引起缸盖振动的原因主要是本缸气体燃烧时产生的爆压、本缸气阀落座撞击以及排气门开启所引起的气流冲击等，其次邻缸的各种振动激励源也会对缸盖的振动产生较大影响[14]。图8所示为进排气阀开闭与曲轴转角的关系。进气门开启的角度在排气上止点前20°附近，关闭的角度在进气下止点后48°附近；排气门开启的角度的在做功下止点前48°附近，关闭的角度在排气上止点后20°附近；柴油机在0°点火。

图8内燃机燃烧和气阀开闭转角图

Fig.8Angular distribution of vibroimpact events from valve train and combustion

对正常工况信号进行VMD分解，不同K值下的中心频率如表2所示。表2不同K值对应的中心频率（单位：Hz）

Tab.2Center frequency corresponding to different K（Unit：Hz）

K=2K=3K=4K=51456130412847767678447744661783—768574644473——81177464———8117

对于正常工况，当K值为4和5时，中心频率出现了比较相近的7464 Hz和8117 Hz，这两个中心频率相近8117/7464=1.08，1.08小于过分解阈值θ，认为出现了过分解现象，因此K应取分解层数4的上一层，即分解层数K=3。图9所示为正常工况信号经KVMD分解后的波形及其功率谱图。

图9正常信号的VMD分解的波形与功率谱

Fig.9Waveform and spectrum of normal signals KVMD decomposition

气阀正常工况下的缸盖信号通过KVMD分解得到的中心频率为1304，4477和7685 Hz，与文献[15]中所述内燃机缸盖振动信号的频域特性“进气门开启和闭合时产生的振动响应相似，其能量主要集中在6～8 kHz；排气门开启和关闭时产生的振动响应相似，其能量主要集中在6.5～8 kHz；燃烧产生的振动能量主要集中在4～6 kHz，燃烧后段产生的振动能量主要集中在0.8～1.7 kHz”相一致。中心频率1304，4477和7685 Hz分别与燃烧后段、燃烧和进排气门开启和关闭的能量相对应，充分验证了过分解阈值θ=1.1的有效性。按照上述方法分别对气门间隙过小，气门间隙过大和气门漏气工况进行分析，大量的实验结果表明KVMD方法对信号的剖分适当，有利于对信号的进一步分析研究。

分别绘制缸盖表面振动信号MHD时频分布图和KVMDMHD时频分布图，如图10和11所示，每幅图中从上至下依次为气门间隙正常、过小、过大和漏气4种工况。图中最上方的曲线为信号的时域波形图，横坐标表示时间，纵坐标表示幅值；左边的曲线为信号的功率谱图，横坐标表示频率，纵坐标表示幅值（将功率谱图顺时针旋转90°看）；中间的图像为时频相平面图，横坐标表示时间，纵坐标表示频率；图中的颜色代表能量幅值的大小，与右边颜色标尺图对应。

图10振动信号的MHD时频图

Fig.10MHD timefrequency image of vibration signal图11振动信号的KVMDMHD时频图

Fig.11KVMDMHD timefrequency image of vibration signal

从图中可以发现缸盖表面的振动信号具有时变和非平稳的特性，随气门间隙的变化各个工况的时频相平面图呈现出较大差异，冲击分量信息在时频相平面图上出现和消失时间不同，幅值大小不同，并且它们频率组成更不相同。对比图10和11可得，图10用MHD方法对缸盖振动信号进行时频分析时频域上存在严重的交叉干扰项，只能分辨出较大冲击位于曲轴转角的位置，无法表示该位置存在的具体频率，造成频域信息丢失，增加了故障诊断难度。图11中KVMDMHD方法可以有效地抑制MHD方法中的交叉干扰项，清晰地分辨出各较大振幅处存在的频率分量，具有良好的时频聚集特性，更有利于后续特征提取与分类。

从能量的分布的角度来看：图11中可以看出，图（a）气门间隙正常时缸盖振动信号的能量主要集中在7～8.5 kHz之间的频带；图（b），（c），（d）当气门间隙处于故障状态时，主要能量会集中在9～12 kHz高频区，相比于正常工况，主要能量分布有向高频移动的趋势。

从燃烧做功的角度来看：气缸内的混合可燃气体做功与否或是否充分燃烧，其特征信息必然会在曲轴转角0°附近体现。图11（a）中，气阀间隙正常时，内燃机正常工作，缸内气体燃烧正常，其对应的冲击分量十分明显。但这一振动分量在图11（b），（c），（d）中很不明显，这说明气门间隙异常（过大或过小）对柴油机的燃烧影响比较大。因为排气阀气门间隙过小或漏气，就会引起气门密封不严，产生漏气；过大，则将使气门迟开、早关，排气时间缩短，影响混合气体的更新，影响正常燃烧。

从振动分量分布的角度来看：4种工况的进气阀都正常，所以进气阀落座产生的冲击分量在4幅图中曲轴转角-132°附近位置均得以体现。图11（a）中排气阀处于正常状态，所以其位置对应于在曲轴转角-340°附近和频率为7.8 kHz附近；图11（b）中，排气阀气门间隙过小，因此冲击分量在图中表现的不是很明显；图11（c）和图11（d）中排气阀处于气门间隙过大和漏气状态，因此频率区别于正常工况的7～8.5 kHz，迁移至高频部分9～12 kHz。曲轴转角为132°和340°附近时，排气阀和进气阀先后开启，由于气阀开启时引起的冲击相比于气阀关闭或是燃烧引起的冲击要小的多，因此在时频分布图中体现的并不是很明显。

3.4KVMDMHD时频谱的TD2DPCA特征提取

取采集到的240个信号作为研究对象并分别绘制KVMDMHD时频相平面图，相应得到240个300×400像素点的时频矩阵。由于得到时频矩阵维数高，计算量大，不利于进行特征参数的提取，对图像进行预处理，将其转化为灰度图像。

KVMDMHD时频相平面图的局部非负矩阵特征参数提取流程如下：

Step 1从4类工况时频分布图中每一类随机选取30幅共120幅，组成TD2DPCA样本集T；

Step 2对样本集T进行TD2DPCA特征提取，最优投影矩阵P300×d和Q400×h。d和h分别表示两次提取的特征维数，它的取值对特征提取结果和后续的识别精度有较大影响；

Step 3将所有时频相平面图向矩阵P和Q投影，可得其对应得编码系数矩阵H，H的维数为h×d，每一个编码系数矩阵H代表了它所对应的时频相平面图。

图12给出的是特征维数h×d=5×5时，KVMDMHD时频相平面图训练集对应的特征系数，图中每个像素所对应的色柱值严格与样本的编码系数值保持一致，文章篇幅有限，每种工况下选取5图12TD2DPCA提取的测试集特征系数

Fig.12Test set parameters for TD2DPCA

个样本的编码系数矩阵进行显示。图中每一行代表一种内燃机工况，从上到下依次为气门间隙正常、过小、过大和漏气。可以看出TD2DPCA对数据进行了非常有效的降维，将300×400维数据压缩到5×5维，大大降低了识别复杂度和计算量。从图中可以看出，同种工况样本的编码矩阵系数较为相似，不同工况间编码矩阵系数区别较大，这有利于后续参数的分类识别。

3.5故障识别

在对内燃机气门间隙工况进行分类时，选择最近邻分类器（KNNC）作为内燃机运行工况判别的智能学习机器。从4类工况中每一类中随机选出30个编码矩阵H共120个，组成训练样本集合。然后用剩余的120个系数向量作为测试集合进行分类测试，重复以上实验10次取平均值。用识别正确率为指标来评价文中方法的性能。为进行对比分析，分别采用双向二维主成分分析（TD2DPCA）、双向二维线性判别分析方法（TD2DLDA）和二维非负矩阵分解（2DNMF）算法对生成的MHD和KVMDMHD时频分布图像进行特征提取并分类。由于在上述特征提取方法对时频分布图进行特征提取过程中涉及特征维数对识别准确率的影响，为增强3种方法的对比性，保持编码矩阵维数的一致性，令2次提取的特征维数d×h=r×r=[2×2，3×3，…，10×10]，识别率准确率结果如图13所示。

图13三种特征提取方法的识别率

Fig.13Recognition rate of three feature extraction methods从图13（a）中可以看出，使用3种分类方法对MHD时频分布图进行特征提取，其中TD2DLDA的识别率相对较低，在不同的特征维数下识别准确率平均值低于85%；TD2DPCA和2DNMF的识别准确率相差不大，各个特征维度的识别准确率都在90%以上，但在试验过程中，由于2DNMF在特征提取过程中要将所有训练数据矩阵按行和列进行拼合，初始分解矩阵维数较大，迭代过程效率低，计算耗时较长。图13（b）中，3种分类方法对KVMDMHD时频分布图进行特征提取，识别率与图13（a）中相比，3种方法均有较大提高。这是由于采用KVMDMHD对内燃机振动信号进行时频分析时，生成的时频分布图时频聚集性好，各个工况间的差异更明显，更利于分类器的分类。TD2DPCA在特征矩阵维度为5×5和更高维度时，识别准确率高达100%。对比图13（a）和（b），表明采用基于KVMDMHD与TD2DPCA的故障诊断方法适用于内燃机气门间隙的故障诊断，并具有較高的诊断精度。

4结论

1） 优化了VMD分解中K值，增强了分解的自适应性，将其与MHD时频分析法相结合，提出了KVMDMHD时频分布，该方法有效抑制了MHD分布中存在的交叉干扰项，具有很高的时频分辨率。用该方法对不同气门间隙工况进行分析，各工况的时频分布特征明显，时频分量物理意义明确。