基于短时傅立叶变换和主分量分析的发动机失火故障诊断方法

摘要

本发明提供一种基于短时傅立叶变换和主分量分析的发动机失火故障诊断方法：首先基于短时傅立叶变换精确获得工作循环内反映各缸工作状态的瞬时转速波形，并在此基础上进行工作循环内波形的特征参数提取，提取反映各缸工作状态的高维特征向量。本发明以提取的高维特征向量为分析对象，利用主分量分析的T2和SPE监测指标对失火故障进行监测，并采用主分量分析的贡献图法对故障缸进行精确定位。

说明书

技术领域

本发明属于发动机故障诊断领域，具体涉及一种基于短时傅立叶变换和主分量分析的发动机失火故障诊断方法。

背景技术

失火故障是多缸发动机的一种常见故障，它不仅降低了发动机的动力性、经济性，而且加剧排放污染。因此对发动机失火故障进行有效监测和诊断十分必要。发动机在稳定工况下工作时，其曲轴平均转速是不变的，当发生失火故障时，故障缸的工作状态发生变化，这会引起曲轴的瞬时转速也发生显著波动。因此利用曲轴的瞬时转速波动能对失火故障进行有效诊断，而准确及时获得曲轴瞬时转速是对发动机进行有效失火故障诊断的关键和基础。目前最为常见的瞬时转速测量手段是采用霍尔式传感器，直接利用内燃机固有的编码盘(如飞轮齿圈)通过计数法进行测量。这种测量方法影响测量精度的因素较多，需要对传感器原始信号进行处理。处理的方法主要包括硬件计数法和软件计数法两种。前者测量成本高昂，后者主要通过软件插值的方法计算瞬时转速，此方法受噪声信号影响较大，精度差，当齿圈存在缺齿时(用于指示气缸压缩上止点)，瞬时转速计算结果会出现严重畸变。

发明内容

本发明的目的在于提供一种基于短时傅立叶变换和主分量分析的发动机失火故障诊断方法。

为了实现上述目的，本发明采用了以下技术方案：

1)将发动机转速信号等间隔划分成长度为L的连续信号，然后对划分的每段长度为L的连续信号进行短时傅立叶变换，得到转速信号的时频分布；

2)根据转速信号的时频分布，采用基于代价函数的时频分布脊线提取方法计算瞬时频率；根据瞬时频率计算发动机的瞬时转速；

3)对瞬时转速波形按照不同气缸做功区间进行划分，在各自区间内分别计算转速波形的波动参数值、峰值和峭度，其中波动参数值定义为：

B＝|M_right-M_left|

其中M_right和M_left分别表示单缸各自区间转速波形右端点的幅值和左端点的幅值；

那么，由瞬时转速波形提取得到关于如下3N维特征向量F的特征数据：

F＝[B P K],B＝[B_i],P＝[P_i],K＝[K_i]；i＝1,2,…N

其中P为峰值，K为峭度，N为气缸数目，i按照气缸工作时序进行排列；

4)根据步骤3)获得的特征数据，采用主分量分析进行失火故障监测与诊断。

所述发动机转速信号的获取方法包括以下步骤：对采集的发动机原始转速信号通过低通滤波去除周期性脉冲型干扰信号，然后去除由发动机压缩上止点缺齿所造成的异常信号。

所述发动机原始转速信号为安装于发动机飞轮盘处的转速传感器所输出的电压信号。

所述低通滤波采用阶数为50～150的低通滤波器，该滤波器的截止频率为3000～5000Hz。

所述异常信号的去除方法具体包括以下步骤：对低通滤波后的转速信号进行峰值搜索，当相邻峰值之间的数据点数小于设定阈值时，去除该相邻峰值之间的所有数据点。

所述L的取值为4000～6000。

所述短时傅立叶变换按照下式进行：

$>{\mathrm{STFT}}_x(t,f)={\int }_{-\infty }^{\infty }x\left(t\right)h(t-\tau )e^{-j2\pi ft}dt$>

其中x(t)是长度为L的连续信号，h(t)是一个沿着时间轴滑动时间的窗函数，τ是时间平移因子，t表示信号时间，f表示信号频率。

所述瞬时频率的计算方法具体包括以下步骤：

2.1)构建如下式所示的代价函数：

CF_k＝-|S(k,a(k))|²+|a(k)-a(k-1)|²,k＝2,3,...,n

其中k为时间轴上的索引值，a(k)为k时刻的频率值，n为索引值的最大值，S为对应时刻短时傅立叶谱；

2.2)通过最小化代价函数，在时间轴上寻找谱峰模值较大并且变化光滑的局部最大值作为瞬时频率值。

所述瞬时转速按照下式进行计算：

S_n＝60C_ff_s/zN_ft

其中S_n为瞬时转速，C_f为瞬时频率，z为飞轮盘齿圈齿数，f_s为采样率，N_ft为傅里叶变换长度。

所述步骤4)具体包括以下步骤：

4.1)将特征数据定义为X_p∈R^m×d，R为特征数据集，m为样本的个数，d为监测变量的个数；将X_p分解为向量t_i和p_i的外积以及残余矩阵E_p之和：

$>X_p={\mathrm{\Sigma t}}_i{p}_i^T+E_p$>

其中t_i是包含样本之间信息的得分向量，p_i是包含变量之间信息的负荷向量；

4.2)根据特征数据分解结果，计算各个样本的监测指标T²或SPE的值以及监测指标对应的控制上限；

4.3)当监测指标T²或SPE的值超过控制上限时，对对应样本的各个气缸的特征数据对监测指标T²或SPE的贡献值进行计算，并根据贡献值大小，对发生失火故障的气缸进行定位。单气缸发生失火故障时，根据发生失火故障的气缸对应的特征数据计算得到的对监测指标的贡献值为最大。

本发明的有益效果体现在：

本发明方法首先基于短时傅立叶变换精确获得工作循环内反映各缸工作状态的瞬时转速波形，并在此基础上进行工作循环内波形的特征参数提取，提取反映各缸工作状态的高维特征向量。本发明以提取的高维特征向量为分析对象，利用主分量分析的监测指标对失火故障进行监测，并采用主分量分析的贡献图法对故障缸进行精确定位。

本发明方法的优点如下：

(1)基于短时傅立叶变换(STFT)的瞬时转速提取方法，获得内燃机瞬时转速精度高、方法抗干扰能力强。

(2)采用主分量分析(PCA)对失火故障进行监控，监控特征值和监测指标运算量小，对失火故障检测率高，利用贡献图法能准确定位故障缸。

(3)本发明方法直接采用正常工作数据进行建模，因此不需要事先获得故障数据，对发动机的失火故障诊断，具有较强的工程实用性。

附图说明

图1为转速传感器原始波形。

图2为滤波后转速传感器波形。

图3为转速传感器的时频分布。

图4为不同故障的整工作循环内的瞬时转速比。

图5为1缸失火故障监控图。

图6为1缸失火故障T²监测指标贡献图。

具体实施方式

下面结合附图和实施例对本发明作详细说明。

本发明以发动机台架的霍尔式转速传感器原始输出电压信号为失火故障诊断对象，由于发动机曲轴转速信号是一种典型的非平稳信号，而非平稳信号的瞬时频率因其具体的物理概念在内燃机瞬时转速计算方面也得到了广泛的应用。受此启发本发明首先采用基于短时傅立叶变换的瞬时转速提取方法获得准确的曲轴瞬时转速，在此基础上，对波形转速特征参数进行提取，采用基于主分量分析的多变量监测和诊断方法对发动机失火故障进行诊断。由于曲轴在一个工作循环内的转速波动信息能够反映各缸的工作状态，因此对各缸的高维特征向量进行深度挖掘，不但能准确的对内燃机故障状态进行识别，更能利用特征向量和各缸的对应关系进行故障缸的精确定位。

本发明采用型号为4G15S的汽油发动机台架进行实验，验证基于短时傅立叶变换和主分量分析的内燃机失火故障诊断方法的有效性。该发动机飞轮盘齿圈齿数为60。其中两个齿被去掉用于指示气缸压缩上止点位置。飞轮盘处(飞轮盘连接发动机曲轴，采集的就是发动机转速信号)安装有霍尔式转速传感器，设置采样率为200K，对转速信号进行了采集，采集的原始转速波形如图1所示，传感器输出值偏小处位置为气缸压缩上止点位置，此处由于齿圈缺齿，导致信号波形幅值明显偏小。由于受到飞轮盘齿圈齿形误差、分度误差等因素的影响，产生周期性脉冲型干扰信号，造成转速波形畸变。本发明采用阶数为100的FIR低通滤波器滤除脉冲干扰，滤波器窗函数设置为汉宁窗，截止频率设置为4000Hz。去除畸变干扰后转速波形如图2所示。为了去除气缸压缩上止点低电压输出对传感器整体输出电压的影响，对滤波后的转速波形进行峰值搜索，判别并去除压缩上止点缺齿所造成的低电压数据点。具体过程为：通过计算峰值之间的数据点数，当相邻峰值之间的点数小于设定阈值时(阈值设定依据为：非缺齿位置的转速信号相邻峰值之间的点数为w，则阈值设置为w/2)，通过删除相邻峰值之间数据点数小于阈值的该相邻峰值之间所有转速信号数据点，从而去除缺齿造成的低电压数据点。然后将去除低电压数据点的转速信号，首尾依次进行连接，实现转速信号重组。

在发动机一个工作循环过程中，其输出扭矩依次受当前做功气缸所影响，对于一个N缸内燃机，气缸i(i＝1,2…,N)在曲轴变化角范围内影响输出扭矩。TDC_i是气缸i(i＝1,2…,N)的曲轴起始角度。由于气缸i的输出扭矩可以影响到上述曲轴变化角内的瞬时转速，因此可以通过气缸i(i＝1,2…,N)在曲轴变化角范围内的瞬时转速，对气缸i(i＝1,2…,N)的工作状态做出评估。由于转速信号是典型的非平稳信号，因此采用时频分析算法STFT(短时傅立叶变换)，对重组的转速波形数据点进行时频分析，其具体过程如下：首先将被分析信号等间隔划分成长度L＝5000的连续信号，然后对每段信号x(t)进行如下短时傅立叶变换：

$>{\mathrm{STFT}}_x(t,f)={\int }_{-\infty }^{\infty }x\left(t\right)h(t-\tau )e^{-j2\pi ft}dt---\left(1\right)$>

其中h(t)是一个沿着时间轴滑动时间的窗函数，其设置为宽度l＝1000的汉宁窗，τ是时间平移因子，其设置为900，t表示信号时间，f表示信号频率，j表示复数单位。得到的转速信号时频分布如图3所示。从图3中可看出，转速信号的时频能量主要集中在1900Hz-2100Hz区间内。

由于干扰和其他误差因素的影响，信号在一个窗内的时频分布存在着多个峰值，直接采用谱峰最大值的方法计算瞬时频率会产生较大的误差。而由于转速波动在稳定工况条件下波动较小，可以认为相邻窗内的信号瞬时频率是连续的，因此本发明采用基于代价函数的时频分布脊线提取方法，精确计算瞬时频率。代价函数构建如下：

CF_k＝-|S(k,a(k))|²+|a(k)-a(k-1)|²,k＝2,3,...,n>

式中S为对应时刻短时傅立叶谱，其由公式1获得，k为时间轴上的索引值，对应公式1的时间t，a(k)为k时刻的频率值，n为索引值的最大值，其由信号长度决定。该代价函数由两部分构成，前者表示谱峰的模值，后者表示谱峰值的变化率，通过最小化代价函数，在时间轴上寻找谱峰模值较大并且变化光滑的局部最大值作为瞬时频率值，此时对应的瞬时转速为：

S_n＝60C_ff_s/zN_ft>

式中S_n为瞬时转速，单位为rpm(转/分)，C_f为瞬时频率，z为飞轮盘齿圈齿数。f_s为信号采样率(取值是200K)，N_ft为傅里叶变换长度。

在发动机正常工作以及气缸1、2、4分别出现失火故障时，对转速传感器信号进行了采集，采用以上方法提取瞬时转速，结果如图4所示。从图4中可看出，发动机正常工作时，随着活塞在气缸内的往复运动，瞬时转速也在有规律的波动，瞬时转速波形在气缸做功区间出现相应波峰，且波峰幅值近似相等。当发动机出现单缸失火故障时，其各缸工作状态的均匀性明显变化，故障缸对应波峰消失或幅值明显降低。在工作时序上，故障缸相邻气缸的工作状态也受到影响，其波峰幅值明显低于其他气缸谱峰，因此通过比较相邻气缸的工作状态，不但能对单缸失火故障进行有效检测，而且也能为故障缸的定位提供有用信息。

为了反映不同气缸对应的转速波形特征，在一个工作循环内，对瞬时转速波形按照不同气缸做功区间进行划分，根据瞬时转速的变化趋势，当瞬时转速由一个波谷上升到一个波峰，又回到一个波谷时，判断为一个做功区间，在各自区间内分别计算转速波形的波动参数值、峰值和峭度时域指标。其中波动参数值定义为：

B＝|M_right-M_left|>

式中M_right和M_left分别表示单缸各自做功区间波形的右端点和左端点的幅值。波动参数值和峰值指标分别反映不同气缸瞬时转速的波动和幅值大小分布情况，而峭度指标能够反映瞬时转速概率分布曲线的陡峭程度，对失火故障造成的瞬时转速波形趋势变化较为敏感。按照以上特征提取方法，对一个工作循环内的瞬时转速波形提取如下3N维特征向量：

F＝[B P K],B＝[B_i],P＝[P_i],K＝[K_i]；i＝1,2,…N>

其中P为峰值指标，K为峭度指标，N为气缸数目，i按照工作时序进行排列。

本发明实验所用发动机的点火顺序为1-3-4-2。发动机连接着测功机，通过切断单缸油路来模拟单缸失火故障，在不同的发动机转速和测功机功率模式下进行了实验。实验故障模式设置如表1所示。

表1.实验故障模式

根据式(5)定义的特征值，对转速2000rpm，正常和单缸分别失火(1缸，2缸和4缸)状态下的瞬时转速信号(瞬时转速波形)进行了特征提取，结果如表2所示，特征1～4分别对应做功区间1～4的瞬时转速波动参数，特征5～6分别对应做功区间1～2的瞬时转速峰值。表2中故障缸对应特征值用粗体加斜体表示，时序上相邻缸对应特征值用粗体表示。从表2中可以看出，故障缸对应波动参数值明显增大，而峰值指标明显降低，峭度指标增大。时序上相邻缸的相应指标值也发生相同趋势的变化，但变化幅度变小。

表2.瞬时转速信号特征提取结果(部分)

下面将采用基于多变量监测理论的主分量分析(PCA)方法，对发动机失火故障进行监测和诊断。对于包含d个变量的监测过程，其数据可定义为X_p∈R^m×d(R为特征数据集，m为特征值样本的个数，d为监测变量的个数(取值为3N维特征向量的特征个数)，一个工作循环包括两个样本，由一个样本求得一组特征值。PCA可将X_p分解为向量t_i和p_i的外积以及残余矩阵E_p之和：

$>X_p={\mathrm{\Sigma t}}_i{p}_i^T+E_p---\left(6\right)$>

式中t_i是包含样本之间信息的得分向量，p_i是包含变量之间信息的负荷向量，下标i表示特征编号，T为转置。进行过程监控时，首先基于正常历史数据进行PCA分析，建立正常数据的主分量模型，然后将监控过程中新得到的数据投影到两个正交的子空间上，并分别计算T²和SPE监测指标。当监测指标偏离正常值时，可以判断故障发生。T²监测指标定义为：

T²(k)＝t(k)^TΛ^-1t(k)>

式中Λ^-1是和保留的主分量相关的特征值构成的对角矩阵，t(k)为对应的特征向量矩阵。T²监测指标的控制上限UCL定义为：

$>UCL=\frac{r(m^2-r)}{m(m-r)}{F}_{\alpha }(r,m-r)---\left(8\right)$>

其中r是保留的主分量数，r按照保留分量数占原有分量数的70～95％(例如，85％)的标准进行设置(对特征值矩阵进行特征值求解，按照所求的特征值的大小进行排序，前几阶保留下来的较大的特征值所对应的特征值向量即为主分量)，m为样本数，F_α(r,m-r)是置信度为α，自由度为r和m-r的F分布。

SPE监测指标定义为：

$>SPE\left(k\right)=x{\left(k\right)}^T(I-P_r{P}_r^T)x\left(k\right)---\left(9\right)$>

式中P_r是主元模型中前r个负荷向量构成的矩阵，r为保留的主分量个数。SPE监测指标的控制上限Q_α定义为：

$>Q_{\alpha }={\theta }_1{[\frac{C_{\alpha }\sqrt{2{\theta }_2{h_0}^2)}}{{\theta }_1}+1+\frac{{\theta }_2{h}_0(h_0-1)}{{{\theta }_1}^2}]}^{1/h_0}---\left(10\right)$>

其中r为保留的主分量数，d为监测变量数，λ为特征向量(上标i为变量编号，下标j为特征值编号)。C_α为置信度为α的正态分布。

对照表1的实验故障模式，采用本发明方法对发动机单缸失火故障进行了诊断。按照式(5)提取了正常和故障状态下的瞬时转速特征值向量，按照式(7)和式(9)计算了T²和SPE监测指标的监控图，其中模式1的监控图如图5所示。图5中虚线为由正常状态所计算的在置信度95％条件下的监测指标阈值(即为控制上限)。图5中所示样本在第65个样本之前为正常状态，从第66个样本开始为一缸失火故障状态。从图5中可以看出，除个别样本外，T²和SPE监测指标值在故障发生时刻之前，均在监测指标阈值之下。而从故障发生时刻开始，监测指标值大大超过了监测指标阈值。通过计算，表1中各个模式下的故障检测率均为100％。为了验证本发明所提方法的误报率，采用正常状态下的样本进行了训练和测试，在置信度95％条件下，其T²和SPE监测指标的误报率分别为0.33％和0.79％，这表明，本发明方法能在低误报率条件下，有效对单缸失火故障进行监测。

在图5中，当T²和SPE监测指标值超过监测指标阈值时，对该时刻(具体指第66个样本)的对应的特征值对T²指标的贡献值进行了计算(采用SPE的贡献图实现故障缸定位，理论上也是可行的)，得到的贡献图如图6所示。从图6中可以看出1缸发生失火故障时，1缸所对应的特征值对故障造成的监测指标T²异常的贡献量最大，其次是时序上相邻的3缸的相关贡献量。因此，通过对故障发生时刻的监测指标贡献图进行分析，能准确的定位发生失火的故障缸。