作为一种有效的汽车排气处理手段，三元催化净化技术被广泛地应用于汽车排放系统，在实现汽车低污染排放方面发挥着重要作用^[1-3]。但在实际车辆行驶中，三元催化器易出现化学中毒、过热老化、机械损坏和结焦积碳等故障，致使排放水平急剧恶化，造成严重的环境污染。因此，针对三元催化转化器的故障诊断一直是汽车领域重点关注的问题^[4]。

近几十年来，众多学者对三元催化器做出了大量富有成效的研究^[5-9]，构建出三元催化器化学反应动力学模型^[5-6]、传热和流体数值分析模型^[7]、储放氧模型^[8-9]等诸多模型来监测分析三元催化器的性能。在模型构建过程中出现的模型复杂度和鲁棒性之间的矛盾长期困扰着研究者，即引入较多反应参数可以提高模型精确性，但会增加计算复杂度，并降低模型的鲁棒性。通过简化模型来降低复杂度提高模型的鲁棒性，却又无法详尽刻画潜在反应过程，使精确度降低，增加故障诊断的难度^[10]。当前，高速发展的现代传感技术已能实现对汽车尾气成分的准确检测^[3]，由此产生了海量的排放数据。这些海量的排放数据携带着三元催化器绝大部分物化反应信息，给基于尾气数据的三元催化器故障诊断奠定了基础。虽然当前利用SVM、神经网络、统计模式识别等智能方法对三元催化器的研究获得一定进展^[11-14]，但大多关注于对故障模型参数的辨识以及对三元催化器中部分废气的预测，未直接对尾气数据中蕴含的故障信息进行分析^[15-18]。

事实上，三元催化器的工作过程是一个涉及物化反应、传热、流体流动的复杂过程，加之车辆驾驶方式因人而异，车辆行驶环境复杂多变，所以尾气信号属于一种典型的非线性、随机性和耦合性的非平稳信号。考虑到尾气信号的复杂性且数据驱动的三元催化器故障诊断本质属于模式识别问题，选择合适的故障特征提取方法与分类(聚类)方法对提高三元催化器故障诊断准确性至关重要。分数阶傅里叶变换(FRFT, fractional fourier transform)是一种高精度的时频分析方法，可将原始数据空间中可分性差的信号映射到特定分数阶域来实现耦合信号的分离^[19]，这给尾气信号的特征提取提供一条有效途径。然而，在FRFT结果所构成的高维特征空间中，并不是任意阶次的FRFT都能最佳地刻画三元催化器故障特征。为此，引入粒子群优化算法(PSO)，借助其建模简单、收敛快速等优点，在构建合理的适应度函数基础上，定位所需的最优FRFT阶次。为进一步降低故障聚类的计算开销，采用分形计算^[20]和核熵成分分析^[21](KECA, kernel entropy component analysis)，在确保信息损失最小的前提下，将高维故障特征映射到较低维度^[22-23]。在诸多聚类分析方法中，模糊C均值聚类(FCM)^[24]在故障诊断领域应用最为广泛，但经典的FCM仅单一使用欧氏距离或角度作为测度来反映各特征向量之间的距离关系，使分类精度的提高受限。综合考虑这两种测度，通过构造含有欧式距离和角度距离测度来优化迭代目标函数，提高故障分类的精确性。

1 改进测度下的模糊C均值三元催化器故障诊断方法

基于改进测度的FCM三元催化器故障诊断方法流程如图 1所示。诊断过程主要包括原始信号获取，故障特征提取与故障聚类分析3个部分。在获得不同故障工况下的尾气信号之后，通过FRFT给出尾气信号从时域逐步变化到频域的所有特征，并采用PSO算法选择最优的FRFT阶次。这些最优阶次的FRFT特征经过分形算子转换成分形特征向量之后，通过KECA对候选故障特征实现维数约简。经过降维处理的故障特征向量被提交给改进测度的FCM进行故障聚类分析。

1.1 分数阶傅里叶变换下的故障特征提取与优化

对携带故障信息的尾气时域信号x(t)的p阶分数阶傅里叶变换(FRFT)可定义为

$ {X_p}\left( u \right) = \int_{ - \infty }^{ + \infty } {x\left( t \right){K_p}\left( {t, u} \right){\rm{d}}t}, $

(1)

式中K_p(t, u)为FRFT的核函数，可表示为

$ {K_p}\left( {t, u} \right) = \left\{ \begin{array}{l} \sqrt {\frac{{1 - j\cot \alpha }}{{2{\rm{ \mathit{ π} }}}}} {{\rm{e}}^{\frac{{j{\rm{ \mathit{ π} }}}}{2}\left[{\left( {{t^2} + {u^2}} \right)\cot \alpha-2tu\csc \alpha } \right]}}, \\ \alpha \ne n{\rm{ \mathit{ π} }}, \\ \delta \left( {t - u} \right), \alpha = 2n{\rm{ \mathit{ π} }}, \\ \delta \left( {t + u} \right), \alpha = \left( {2n - 1} \right){\rm{ \mathit{ π} }}, \end{array} \right. $

(2)

其中：p为变换的阶次，取值范围一般为0~1；p阶FRFT是在(t, ω)平面上旋转α角度所产生的坐标空间，α=πp/2为变换的旋转角度，一般取0~π/2的逆时针方向进行分析。由FRFT的对称性和周期性可知，其它角度的信号分析结果与0~π/2旋转角度相一致。随着p从0变到1，FRFT平滑地从原函数变化到傅里叶变换。因此，FRFT能给出信号从时域逐步变化到频域的所有特性。由于在不同FRFT变换域上的信号分布呈现不同的特性，当旋转到一定的分数阶，原本在时频平面上交错的信号将会分离。为实现故障尾气信号在FRFT时频平面上的最优分离，利用PSO算法对FRFT变换结果进行搜索。

PSO算法具体描述为：在p∈[0, 1]范围内，随机产生均匀分布、群体规模为N的粒子群，每个粒子代表一个FRFT变换结果。其中第i个粒子位置是v_i，飞行速度是v_i。根据粒子的适应度值，第k次迭代使粒子的自身历史最优位置为p_best，全局粒子最优位置为g_best，粒子通过这2个不断更新的极值来更新自己的速度和位置。迭代公式如下

$ v_i^{k + 1} = wv_i^k + {\eta _1}{\rm{rand}}\left( {} \right) \cdot \left( {{p_{{\rm{best}}}} - x_i^k} \right) + {\eta _2}{\rm{rand}}\left( {} \right) \cdot \left( {{g_{{\rm{best}}}} - x_i^k} \right), $

(3)

$ x_i^{k + 1} = x_i^k + v_i^{k + 1}, $

(4)

其中：k为当前进化代数；η₁、η₂为学习因子，一般取η₁=η₂=2，rand()为介于0到1的随机数；w为惯性权重。当达到最大迭代次数或适应度函数达到预设值时，停止迭代。此时的全局最优位置g_best即为所求的最优解。

采用散布矩阵来构建适应度函数J

$ \mathit{\boldsymbol{J}} = \frac{{{\rm{tr}}\left\{ {{\mathit{\boldsymbol{S}}_b}} \right\}}}{{{\rm{tr}}\left\{ {{\mathit{\boldsymbol{S}}_w}} \right\}}}, $

(5)

其中S_w和S_b分别为类内散布矩阵和类间散布矩阵。

$ {\mathit{\boldsymbol{S}}_w} = \sum\limits_{i = 1}^c {{P_i}E\left[{\left( {\mathit{\boldsymbol{x}}_j^i-{\mathit{\boldsymbol{\mu }}_i}} \right){{\left( {\mathit{\boldsymbol{x}}_j^i-{\mathit{\boldsymbol{\mu }}_i}} \right)}^{\rm{T}}}} \right]}, $

(6)

$ {\mathit{\boldsymbol{S}}_b} = \sum\limits_{i = 1}^c {{P_i}\left[{\left( {{\mathit{\boldsymbol{\mu }}_i}-{\mathit{\boldsymbol{\mu }}_0}} \right){{\left( {{\mathit{\boldsymbol{\mu }}_i}-{\mathit{\boldsymbol{\mu }}_0}} \right)}^{\rm{T}}}} \right]}, $

(7)

其中：c为故障类别数；P_i为第i类的先验概率；x_jⁱ表示第i类样本集中的第j个样本向量；E(*)为期望值函数；μ_i表示第i类样本集的均值向量；μ₀表示全部样本的均值向量。适应度函数J代表内离散度与间离散度的差别，J越大表明不同类别间的差别越大，相同类别的相似性越强。

1.2 核熵成分分析下的分形故障特征降维

由于FRFT给出的故障特征以长度较长的波形表示，不利于后续的故障聚类。因此，引入分形计算对FRFT表示的故障特征加以离散化，同时也降低特征向量的长度。分形盒维数的定义由下式给出

$ {d_V} = \mathop {\lim }\limits_{\varepsilon \to 0} \frac{{\ln {N_V}\left( \varepsilon \right)}}{{\ln \left( {1/\varepsilon } \right)}}, $

(8)

式中：V是Rⁿ的非空有界子集，由对应的最优FRFT特征幅值构成；N_V(ε)表示最大直径为ε且能覆盖V集合的最少个数。实际计算中取有限的ε(介于1和2之间的一个分数)，通过求取一系列ε和N_V(ε)，然后在双对数坐标中用最小二乘法拟合直线，所得直线的斜率即所求分形维数。

考虑到同一类型的三元催化器故障尾气数据包含CO、HC、NO_X、CO₂、H₂O、N₂等多种成分，这些气体成分数据所表征的故障特征向量够成一个高维故障特征空间。为降低故障聚类的计算开销，研究采用核熵成分分析(KECA)对用分形维表示的高维故障特征向量组进行维数约简。KECA以瑞利熵和Parzen窗密度估计为基础，其中二次瑞利熵定义为

$ H\left( p \right) = - \log \int {{p^2}\left( x \right){\rm{d}}x}, $

(9)

KECA的目标是利用最少的特征最大地呈现式(9)中输入数据的熵值。已知对数函数是单调函数，只需考虑V(p)=∫p²(x)dx值的大小。为了估计V(p)，引入Parzen窗密度估计，即

$ \hat p\left( x \right) = \frac{1}{N}\sum\limits_{{x_t} \in \mathit{\boldsymbol{D}}} {{k_\sigma }\left( {x, {x_t}} \right)}, $

(10)

进一步有

$ \mathit{\boldsymbol{\hat V}}\left( p \right) = \frac{1}{N}\sum\limits_{{x_i} \in \mathit{\boldsymbol{D}}} {\hat p\left( {{x_t}} \right)} = \frac{1}{N}\sum\limits_{{x_i} \in \mathit{\boldsymbol{D}}} {\frac{1}{N}} \sum\limits_{{x_i} \in \mathit{\boldsymbol{D}}} {{k_\sigma }\left( {{x_t}, {{x'}_t}} \right)} = \frac{1}{{{N^2}}}{\mathit{\boldsymbol{l}}^{\rm{T}}}\mathit{\boldsymbol{Kl}}, $

(11)

其中：K为由k(x_t, x′_t)组成的N×N核矩阵；l为元素均为1的N×1向量。将核矩阵特征分解，可得K=EDE^T，其中D为以特征值λ₁, …, λ_N为元素的对角矩阵，E为相应的特征向量e₁, …e_N构成的矩阵。将D和E矩阵代入式(11)，得到

$ \mathit{\boldsymbol{\hat V}}\left( p \right) = \frac{1}{N}\sum\limits_{{x_t} \in \mathit{\boldsymbol{D}}} {\hat p\left( {{x_t}} \right)} = \frac{1}{{{N_2}}}\sum\limits_{i = 1} {{{\left( {\sqrt {{\lambda _i}} \mathit{\boldsymbol{e}}_i^T\mathit{\boldsymbol{l}}} \right)}^2}} 。$

(12)

由此可知，上式中核矩阵的特征值、特征向量均对$ \mathit{\boldsymbol{\hat V}}\left( p \right)$有贡献。设U_k为高维故障特征数据通过非线性映射Φ得到的子空间，U_k由对$ \mathit{\boldsymbol{\hat V}}\left( p \right)$贡献最大的前k个特征值及特征向量组成，核熵成分分析的特征映射记为

$ {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_{{\rm{ECA}}}} = \mathit{\boldsymbol{U}}_k^{\rm{T}}\mathit{\boldsymbol{ \boldsymbol{\varPhi} }} = \sqrt {{\mathit{\boldsymbol{D}}_k}} \mathit{\boldsymbol{E}}_k^{\rm{T}}, $

(13)

对式(13)的求解即可实现对高维故障特征向量组的降维操作，该求解过程转化为求解最小化问题

$ {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_{{\rm{ECA}}}} = \sqrt {{\mathit{\boldsymbol{D}}_k}} \mathit{\boldsymbol{E}}_k^{\rm{T}}:\mathop {\min }\limits_{{\lambda _{\rm{l}}}, {e_{\rm{l}}}, \cdots {\lambda _N}, {e_N}} \mathit{\boldsymbol{\hat V}}\left( p \right) - {{\mathit{\boldsymbol{\hat V}}}_k}\left( p \right) = \mathop {\min }\limits_{{\lambda _{\rm{l}}}, {{\pmb e}_{\rm{l}}}, \cdots {\lambda _N}, {{\pmb e}_N}} \frac{1}{N}{\mathit{\boldsymbol{l}}^{\rm{T}}}\left( {K - {K_{{\rm{ECA}}}}} \right)\mathit{\boldsymbol{l}}, $

(14)

其中，K_ECA=Φ_ECA^TΦ_ECA=E_kD_kE_k^T。

1.3 改进相似测度下的模糊C均值故障特征聚类

假设维数约简之后的三元催化器故障特征样本集合为X={x₁, …, x_k, …, x_N}，相应的模糊分类矩阵U=[μ_ik]_c×N, 聚类中心向量C=[c₁, …, c_i, …, c_c]，其中c为聚类中心个数，μ_ik∈[0, 1]为任意故障特征点x_k相对于聚类中心c_i的隶属度，d_ik为故障特征点x_k到聚类中心c_i的距离测度。现有的FCM单一使用欧氏距离或余弦距离作为测度，在刻画各特征向量之间距离方面受到限制，进而影响聚类精确度。为提高聚类精确度，综合欧氏距离和余弦距离，提出一种新的测度G

$ {G_{ik}} = \sqrt {{D_{ik}}\left( {1 - {S_{ik}}} \right) \cdot \left[{\frac{{{D_{ik}}\left( {\max {D_{ik}}-\min {D_{ik}}} \right)}}{N}} \right]}, $

(15)

式中，S_ik为故障特征点x_k到聚类中心c_i的余弦距离测度，表示为

$ {S_{ik}} = \cos \angle \left( {{x_k}, {c_i}} \right)。$

(16)

D_ik为数据点x_k到聚类中心c_i的欧式距离，表示为

$ {D_{ik}} = \left\| {{x_k} - {c_i}} \right\|。$

(17)

综合考虑式(15)、(16)、(17)可得

$ J\left( {C, U} \right) = \sum\limits_{i = 1}^c {\sum\limits_{k = 1}^N {{{\left( {{\mathit{\boldsymbol{\mu }}_{ik}}} \right)}^m}{G_{ik}}} } = \sum\limits_{i = 1}^c {\sum\limits_{k = 1}^N {{{\left( {{\mathit{\boldsymbol{\mu }}_{ik}}} \right)}^m}\sqrt {{D_{ik}}\left( {1 - {S_{ik}}} \right) \cdot \left[{\frac{{{D_{ik}}\left( {\max {D_{ik}}-\min {D_{ik}}} \right)}}{N}} \right]} } }, $

(18)

其中，m为模糊加权参数、控制数据划分过程中的模糊程度，通常取2。对于给定的精度ε，改进的FCM聚类可按如下步骤实现

1) 计算故障特征的聚类中心

$ {c_i} = \frac{{\sum\limits_{k = 1}^N {{{\left( {\mathit{\boldsymbol{\mu }}_{ik}^{t - 1}} \right)}^m}{x_k}} }}{{\sum\limits_{k = 1}^N {{{\left( {\mathit{\boldsymbol{\mu }}_{ik}^{t - 1}} \right)}^m}} }}, 0 \le i \le c。$

(19)

2) 更新划分矩阵

$ \mathit{\boldsymbol{\mu }}_{ik}^{\left( t \right)} = \frac{1}{{\sum\limits_{j = 1}^c {{{\left( {\frac{{{G_{ik}}}}{{{G_{jk}}}}} \right)}^{\frac{2}{{m - 1}}}}} }}。$

(20)

3) 迭代计算

$ \left\| {{U^{\left( t \right)}} - {U^{\left( {t - 1} \right)}}} \right\| < \varepsilon 。$

(21)

2 仿真实验

数值实验按国家标准GB 18352.5：轻型汽车污染物排放限值及测量方法(中国第五阶段)进行尾气排放数据采集，采集过程分为4个低速运转循环和1个高速运转循环。其中，每个低速循环均采集195个数据样本，高速运转循环采集400个数据样本，每个数据样本长度为N=1 180。实际尾气数据采集中，考虑工作正常、化学中毒、过热老化、机械损坏、结焦积炭5种状态，每种状态下采集COH、THC、COL、NO_X、HC和CO₂ 6种气体数据。在相同工况及约束条件下，对上述5种状态分别进行20次独立实验，共得到100组排放数据。图 2给出了其中一组数据，该组数据包含了化学中毒、过热老化、机械损坏和结焦积碳4种故障状态下的三元催化器原始尾气信号。

为了清楚地说明三元催化器故障特征的提取与选择过程，以机械损坏故障数据集中的THC信号x_THC(t)为例。通过对信号x_THC(t)进行FRFT，得到该信号的时频域全局特征，如图 3所示。图 3中给出在机械损坏情况下，当分数阶次p从0变化到1时，x_THC(t)相应从时域到频域的一个连续平滑变化的过程。

为获得最优FRFT阶次，使用PSO算法对信号x_THC(t)的FRFT结果进行寻优搜索。在区间[0, 1]内随机产生均匀分布的30个阶次，以这30个阶次对应的FRFT变换结果初始化粒子群，由此，确定了30个粒子的初始位置。设置学习因子η=2，权重ω=0.75，迭代次数为100，分数阶次p的步长为0.002。寻优得到最佳FRFT的阶次p值为0.932，对应的FRFT特征由图 4(b)给出。图 4(a)、4(c)分别给出p=0(原始信号)和p=1(傅里叶变换)时的变换结果。从图 4中，可明显看出当p=0.932时，x_THC(t)信号能量最为集中，这种能量分布的可分性较强，特征更为突出。

对5种状态下的各气体信号分别进行FRFT，并采用PSO搜索最优阶次，构成5个故障特征数据集，每个故障特征数据集都采用三维数组表示，其维度为1180×6×20。此处，1180代表各类气体FRFT之后的长度，6代表气体种类，20代表独立实验次数。对上述5个故障特征数据集进行分形分析，得到5组6×20的分形故障特征，其均值如表 1所示。

从表 1可以看出5种模式下不同尾气排放数据的分形故障特征，在气体种类的横向(表 1的行方向)差异较为明显，在故障种类的纵向(表 1的列方向)除“机械损坏”与“正常”在“COH列”相近，以及“过热老化”与“结焦积碳”在“CO₂列”相近之外，其他各列差异明显。因此，采用分形来刻画5种状态的特征信息是可行的。

为降低后续故障聚类的计算复杂度，对以上6维候选故障特征运用KECA进行特征降维，核函数宽度取值为σ=5。得到5组3×20的降维故障特征。为评估KECA的有效性，将新的5组3×20故障特征的行向量进行任意抽取，构成一个3×100的数据矩阵。该数据矩阵的三维散点分布如图 5所示。从图 5中，可以明显看出分形故障特征经KECA提取之后，具有明显的聚类倾向。

对KECA降维后的5组3×20分形故障特征数据进行改进相似性测度下的FCM故障聚类。根据5种不同故障模式状态，聚类中心个数定为c=5，加权指数m=2，迭代终止容差ε=0.000 1。迭代计算至收敛得到聚类中心C={c₁, c₂, c₃, c₄, c₅}，分别对应正常、化学中毒、过热老化、机械损坏和结焦积炭5种状态

$ \mathit{\boldsymbol{C}} = \left[{\begin{array}{*{20}{c}} {0.5769}&{0.2977}&{0.4603}\\ {0.1807}&{0.2314}&{0.3348}\\ {0.1633}&{0.3876}&{0.2487}\\ {0.4951}&{0.4371}&{0.3285}\\ {0.6214}&{0.1307}&{0.1794} \end{array}} \right]。$

(22)

故障聚类结果如图 6所示。其中：图 6(a)为改进测度FCM故障聚类的3维空间分布；图 6(b)为图 6(a)的二维投影等高线分布。为验证研究方法的有效性，在相同条件下，分别采用欧式距离和余弦距离的FCM聚类，其故障聚类结果由图 6(c)和图 6(d)分别给出。从图 6可清楚看到，采用改进测度、欧式距离、余弦距离的FCM方法进行故障聚类，100组分形故障特征都围绕在5个聚类中心周围，正常模式和4种故障模式呈明显分离状态，且4种故障之间间距较大没有出现混叠现象。相比之下，选用改进测度与余弦距离测度的等高线比欧式距离测度更加符合特征点分布的走势。特别地，在故障特征点偏离聚类中心程度上，采用改进测度的FCM较之采用欧式距离和余弦距离的FCM要小。

为定量体现改进相似测度FCM聚类方法的优点，采用分区索引(SC)和Xie-Beni指标(XB)作为评价模糊聚类有效性指标标准。XB定义为

$ XB = \frac{{\sum\limits_{i = 1}^c {\sum\limits_{k = 1}^N {{{\left( {{u_{ik}}} \right)}^m}{d_{ik}}} } }}{{N \cdot \mathop {\min }\limits_{i \ne k} {d_{ik}}}}, $

(23)

SC定义为

$ SC = \sum\limits_{i = 1}^c {\frac{{\sum\limits_{k = 1}^N {{{\left( {{u_{ik}}} \right)}^m}{d_{ik}}} }}{{\sum\limits_{k = 1}^N {{u_{ik}}} \sum\limits_{j = 1}^c {{d_{ij}}} }}}, $

(24)

式中：XB指标以最小的类与类中心距离平方来衡量类间分离度，是类内紧密度与类间分离度的比值，其值越小表示聚类效果越好。SC表示聚类紧凑性与分离程度之和的占比，SC值越低表明划分效果越好。3种FCM故障聚类的定量评价指标由表 2给出。

为进一步验证方法的有效性，另采集5种模式下各60组数据作为测试样本，重复上述方法对测试样本进行特征提取，再利用最大贴近度原则进行故障识别。测试样本识别结果如表 3所示

从表 3可以得出，改进测度的FCM故障聚类准确率较之采用欧式距离和余弦距离的FCM聚类，其故障诊断准确率整体水平最高，平均故障诊断准确率分别为96%、89%和92.33%。

3 结论

研究针对三元催化器的故障诊断问题，提出一种利用汽车尾气排放数据对三元催化器进行故障诊断的方法。考虑到尾气信号的复杂性和数据驱动的三元催化器故障诊断本质属于模式识别问题，该方法主要包括原始信号获取、故障特征提取与故障聚类分析3个部分。在获得不同故障工况下的尾气信号之后，通过FRFT给出尾气信号从时域逐步变化到频域的所有特征，并采用PSO算法选择最优的FRFT阶次。这些最优阶次的FRFT特征经过分形算子转换成分形特征向量之后，通过KECA对候选故障特征实现维数约简。经过降维处理的故障特征向量被提交给改进测度的FCM进行故障聚类分析。数值实验结果表明，研究方法能有效对三元催化器的正常、化学中毒、过热老化、机械损坏和结焦积炭5种状态进行识别，较之采用欧式距离和余弦距离的FCM故障聚类方法，故障诊断准确率更高。