现代工业过程规模越来越庞大，系统内部各部件联系非常紧密，随着对现代复杂设备故障诊断精度的要求提高，设备的可靠性及安全性备受学术界的广泛关注^[1-2]。基于模型的故障诊断方法是一种能够有效反映物理系统故障本质特征的方法，可充分利用系统内部的深层信息，有利于系统的故障隔离和辨识，获得了广泛研究。但因现代工业系统的复杂性及其部件工作特性的耦合影响，建立精确的系统模型极为困难，模型偏移误差与参数偏差是必然存在，同时受环境与人为干扰因素的影响，大大降低了故障诊断的准确性。

粒子滤波(Particle Filter, PF)是基于蒙特卡罗采样思想来实现状态的近似估计，理论上可以适用于任何非线性非高斯系统，因此在故障诊断领域获得了广泛的研究和应用^[3-4]。但是目前存在的PF算法及其改进算法，在应用中假设噪声统计特性均为已知的高斯分布^[5-6]。事实上，噪声统计特性在实际应用中很可能是未知的，且可能随系统本身或者量测环境的变化而变化，现有的PF算法难以满足诊断系统精度的要求。代价评估粒子滤波(Cost Reference Particle Filter, CRPF)不需要已知系统噪声的统计特性，便可实现状态估计^[7]，因此在一定程度上解决了复杂噪声背景下的状态估计问题，该方法在故障诊断和目标跟踪等领域得到了良好的应用。胡振涛和潘泉等^[8]提出了一种基于CRPF的残差似然比检验故障诊断算法，克服了外界随机扰动对于滤波精度的不利影响。卢锦等^[9]针对天波雷达背景噪声强度大和统计特性未知等特点，采用CRPF估计天波超视距雷达目标状态。Lim^[10]采用CRPF实现未知噪声统计特性的非线性系统动态状态估计。Yu^[11]提出了一种H_∞和CRPF结合的滤波算法，提高了未知非高斯噪声的非线性动力学系统的状态估计精度。

CRPF是以粒子的代价或风险为参考来实现重采样及状态估计，代价和风险的计算不需要已知噪声统计特性^[9]，因此，在解决非线性、噪声统计特性未知问题中得到了良好应用。但在实际应用中仍然存在不足，CRPF以当前时刻的状态预测值为均值，在一定的高斯分布范围内随机选取粒子完成状态更新^{[10, 12]}，但随着时间的不断推移，状态转移密度协方差逐渐趋于恒定值，对强噪声扰动或时变噪声的修正能力将减弱甚至消失，大大降低了算法对噪声环境的适应能力，将导致较大的误诊和漏诊。对复杂噪声环境下的故障诊断方法的研究是现代复杂系统故障诊断与故障预测研究中必须解决的基础科学问题，因此，本文研究对提高诊断系统的精确性，进一步提高现代复杂设备故障诊断的可靠性具有重要意义。

本文根据噪声的变化情况自适应调整状态转移密度方差，确定状态更新的合理区域，增强算法对时变噪声的适应能力；同时，利用统计方法设计残差判别函数的自适应阈值，提高故障检测的准确性。在此基础上，对强噪声背景下的非线性非高斯系统进行了故障诊断仿真分析。

1 非线性系统状态空间模型

考虑非线性非高斯随机状态空间模型如下：

式中：k为采样时刻；x_k为k时刻系统状态向量；u_k为k时刻系统输入；θ为系统参数变量；y_k为k时刻系统状态的量测向量；g(·)和h(·)分别为系统状态转移函数和量测函数，它们都可以是非线性函数；v_k和w_k分别为系统状态噪声和量测噪声，其统计特性未知。

假设系统可能存在S种故障模式，则共有S+1个系统模型，其中1个模型为无故障模式。采用多模型方法表示的非线性系统状态空间模型如下：

系统量测的预测输出为

式中：x_k|k－1ⁱ为系统状态的一步预测值；ŷ_kⁱ为系统量测的预测输出值。

系统残差可表示为

残差是系统实际输出和依据模型预测输出的差值，反映了实际系统和数学模型之间的不一致程度，主要由系统噪声和故障确定，当系统无故障时，残差的变化反映了系统中的噪声和干扰因素对系统输出的影响。当系统处于强噪声干扰环境，噪声会引起残差较明显的变化，可能引起故障误报或漏报。可见，获得准确的残差是提高故障诊断准确性的关键。因此，利用CRPF对未知噪声系统的状态估计优势，克服强噪声的影响而得到准确的系统状态的估计值，当估计值越接近真实值，则会减小噪声对残差信号的影响，提高故障诊断的准确性。

2 代价评估粒子滤波

在CRPF算法中用代价函数和风险函数表示粒子性能质量，同时引入遗忘因子，并基于代价最小化原则给出一种噪声统计特性未知情况下的粒子权重评价方法。

代价函数定义为

将式(5)简写为：c_k^i_p=λc_k－1^i_p+Δc_k^i_p。式中：i_p为粒子索引；λ(0≤λ≤1)为遗忘因子；Δc_k^i_p=‖y_k－h(x_k^i_p)‖^q(q≥1)为代价增量。

风险函数定义为

概率质量函数(PMF)也称类权值，计算如下：

式中：δ, β>0，δ是为了保证分母不为0；N为粒子数。

按照以上参数定义，CRPF算法通过风险估计、选择、粒子传递和代价更新步骤，递推计算获得状态估计。算法步骤如下：

步骤1    在k=0时刻，从先验分布中获取N个样本，x₀^i_p~P₀(x₀)，设定初始时刻粒子的代价c₀^i_p和状态转移密度方差σ₀^{2, i_p}，则k=0时刻样本和代价集合为{x₀^i_p, c₀^i_p}^N_ip=1，i_p=1, 2, …, N。

步骤2    递归更新。

1) 分别按照式(6)和式(7)计算风险函数R_k^i_p和概率质量函数。

2) 重采样。根据的大小，随机选择N个粒子，生成粒子代价集。

3) 粒子更新。

式中：I为维数与x相同的单位矩阵。

状态转移密度的方差为

式中：dim[x]表示取x的维数。

可见，在CRPF算法中，σ_k^{2, i_p}是在线调整的。

4) 按照式(5)计算代价函数c_k^i_p，按式(10)计算π_k^i_p，并归一化π_k^i_p。

步骤3    状态估计。

从以上算法步骤可见，在CRPF算法的递推估计中不涉及对过程噪声与量测噪声的计算，因此，在计算中不需要已知噪声统计特性，改善了PF中由于外界随机干扰对依据量测似然度评价的不利影响。从式(8)和式(9)的粒子更新过程可知，粒子的时间更新是在确定的概率区域中随机选取新粒子，此区域的方差σ_k^{2, i_p}是根据当前时刻的状态估计值与前一时刻的方差递推更新，当k→∞时，，，这时σ_k^{2, i_p}→σ_k－1^{2, i_p}，可见随着时间的推移，系统状态的当前值失去对方差的调节作用，其方差趋于恒定，对于时变噪声或噪声波动较大的情况下，会大大降低算法的准确度。

3 自适应CRPF故障诊断方法 3.1 状态转移密度的方差自适应更新

由对CRPF的递推估计过程的分析可知，状态转移密度方差直接影响着状态预测的准确性，但在CRPF算法中，方差的递归更新随着时间的推移最终趋于恒定值，对于时变噪声出现随机波动的情况下，在确定的小范围中很有可能无法获得准确的状态预测。因此，本文考虑当前观测值与先验状态之间的互相关系数的大小，自适应调整状态转移密度方差，增强算法对时变噪声的适应能力。

测量值和先验状态之间的相关系数表示为

状态转移密度的方差按照式(12)进行修正：

式中：ω_{1, k}和ω_{2, k}为方差的自适应调节系数，ω_{1, k}=γ_k，ω_{2, k}=1－γ_k。相关系数大，说明状态噪声小，则ω_{1, k}>ω_{2, k}，方差主要由前一时刻的方差σ_k－1^{2, i_p}决定，这样可适当增强粒子的扰动；相关系数小，说明噪声扰动较大，这时ω_{1, k} < ω_{2, k}，方差主要由确定，利用当前时刻的噪声扰动对方差进行修正，提高算法对时变噪声背景下状态估计的准确性。

3.2 自适应阈值设计

残差是判断是否发生故障的主要依据，若k时刻的残差r_k>r_th，r_th为所设定的阈值，则说明系统发生了故障。最简单的阈值选择方法是将阈值设定为一个固定的常数，但是实际系统不可避免地存在建模误差、复杂噪声及干扰等不确定性因素，在这种情况下若仍采用传统的固定阈值来判定故障是否发生，很容易出现故障的误报和漏报。因此，依据系统不确定性的变化，采用自适应阈值可以避免上述情况的发生。文献[13-14]将统计学中参数置信区间的思想应用于自适应阈值的设计，假设残差的统计特性符合正态分布，残差的均值和方差求取方法分别为

式中：η(u_j, t_k)为输入为u_j时t_k时刻残差的均值；σ_v²(u_j, t_k)为输入为u_j时t_k时刻残差的方差；n为系统重复运行并测量系统残差的次数。

置信度为1－α的残差的置信区间表示为

式中：α为置信水平，通常设定置信度为95%~98%。若设α=0.03，则置信度为97%，z=2.17，由式(16)求得自适应阈值为

但采用这种方法计算自适应阈值，为了得到当前时刻残差的均值和方差，在算法迭代的每一步都需要系统重复运行n次，这样计算量近似成指数增长，很难满足实际系统的实时性要求。本文受Hashemi和Pisu^[15]的启发，引入滑动窗通过求区间均值来解决上述问题。将残差的标准差表示为

式中：μ_k为残差的均值；M为滑动窗宽度。

由式(17)和式(18)可见，在每一个采样周期，只需要在前一时刻的基础上通过当前时刻的残差对均值与方差进行更新，无需重复多次运行系统，大大减少运行时间。同时，对于随机噪声来说，在一定的时间内对某一个采样值重复执行多次，与对不同时间点的多个采样值重复执行多次，在该时间范围内噪声的分布趋于一致，因此，令σ_v²(u_j, t_k)=(STD_k)²，η(u_j, t_k)=μ_k，采用式(17)和式(18)分别计算残差的方差与均值，再由式(16)求得每一时刻的残差自适应阈值。

3.3 算法步骤

步骤1    初始化x₀^i_p~P₀(x₀)、c₀^i_p和σ₀^{2, i_p}，生成初始样本代价集合{x₀^i_p, c₀^i_p}^N_ip=1。

步骤2    递归更新：k=1, 2, …, K(总时间步数)。

1) for i_p=1, 2, …, N，计算R_k^i_p和。

2) 重采样。根据的大小，随机选择N个粒子，生成粒子代价集。

3) 根据式(8)更新状态，由式(11)计算互相关系数，然后根据式(12)修正状态转移密度协方差。

4) 计算代价函数：，按式(10)计算并归一化π_k^i_p。

步骤3    估计状态。

步骤4    故障检测和隔离。将获得的状态估计代入量测方程，计算得到ŷ_k，由式(4)计算残差，由式(16)~式(18)计算残差自适应阈值，若r_k>r_{th, k}，则有故障发生; 否则，无故障发生。然后，进一步进行故障隔离。

4 实例分析与讨论

采用一个160 MW燃油机组的动态模型^[16]为对象，以燃油调节阀开度、汽轮机调节阀开度和给水调节阀开度为输入，选取汽包压力、汽包液体密度为状态变量，汽包水位作为观测变量，分别在伽马噪声和高斯混合噪声影响下对汽包水位传感器的3种故障进行实例分析，验证本文方法的性能。该燃油机组模型的离散方程为

式中：

其中：x₁为汽包压力；x₂为汽包液体密度；u₁为燃油调节阀开度；u₂为汽轮机调节阀开度；u₃为给水调节阀开度；y为汽包水位；v_k为状态噪声；w_k为量测噪声，它们都是统计特性未知的非高斯噪声。

针对汽包水位传感器故障进行实验，假设传感器有3种故障模态：

故障模态1    恒偏差故障，y_k=0.05B+0.5+w_k。

故障模态2    恒增益故障，y_k=0.08B+w_k。

故障模态3    卡死故障，y_k=w_k。

设置初始状态x₀=[108  428]^T，离散步长Δt=0.1 s，粒子数N=500，传感器采样频率为1 Hz，其他参数：a₁₁=0.001 8，a₁₂=0.9，a₁₃=0.15，a₂₁=141，a₂₂=1.1，a₂₃=0.19，b₁=0.131，b₂=0.068，b₃=0.001 54，b₄=0.8，b₅=25.6，b₆=1.039 4，b₇=0.001 23，b₈=0.854，b₉=0.147，b₁₀=45.59，b₁₁=2.514，b₁₂=2.096，u_t=[0.3  0.4  0.5]。状态初始先验分布x₀^i_p~N(x₀, Σ₀), Σ₀=diag(0.01, 0.01)，c₀^i_p=0，σ₀^{2, i_p}=[0.01  0.01]^T，滑动窗宽度M=20。

本文方法的评价指标选择平均绝对误差、故障误报率和漏报率，评价指标定义如下:

平均绝对误差MAE为

式中：x_sk和分别为第s次仿真第k步状态的实际值和估计值；总的仿真次数为N_s=50；一次仿真中的时间步数T=3 000。

故障漏报率p_m和误报率p_f为

式中：F为系统总共的运行次数；E为发生故障时残差小于阈值的采样点总数；G为系统一次运行的总采样点数；C为系统无故障发生且残差值大于阈值的采样点总数；D为系统无故障时的采样点总数。

4.1 实验1

本文所设定的偏差故障为0.5的小偏差，增益故障的增益系数变化仅为0.03，噪声对这种小故障的检测和诊断结果影响更为严重，采用伽马噪声和高斯混合噪声来模拟工程实际中的复杂噪声，采用CRPF算法和自适应代价评估粒子滤波(ACRPF)算法对汽包压力x₁和汽包液体密度x₂进行估计，验证本文方法在强噪声环境下的状态估计精度。

采样周期为0.1 s，步长为3 000，总仿真时长为300 s，粒子数为500。设定0≤t≤100 s时，系统处于正常模态；100 s < t≤200 s时，系统处于故障模态1；200 s < t≤300 s时，系统处于故障模态2。

实验中所采用的噪声如下:

伽马噪声：

高斯混合噪声：

图 1和图 2分别为伽马噪声和高斯混合噪声背景下，CRPF算法和改进的ACRPF算法对汽包压力x₁和汽包液体密度x₂的估计误差，0~100 s系统处于正常工作状态，2种算法的估计误差都比较小。可见，由于ACRPF算法能够根据状态估计误差自适应调整状态转移密度方差，增强了对强噪声扰动的自适应修正能力，因此，对x₁和x₂的估计精度均高于改进前的CRPF算法。尤其对于x₂，在100 s系统出现故障，采用原有CRPF算法估计误差急剧增加，而本文改进算法的估计精度得到了大幅提高。表 1为图 1和图 2中x₁和x₂的平均绝对误差，数据同样说明了本文算法的状态估计准确性均得到了明显提高。

4.2 实验2

分别在伽马噪声和高斯混合噪声影响下对上述系统进行故障诊断，对本文ACRPF故障诊断方法的性能进行验证。设残差r₀=|y-ŷ_正常|，若r₀>r_th，说明有故障发生；然后分别计算r₁=|y-ŷ_故障1|，r₂=|y-ŷ_故障2|，r₃=|y-ŷ_故障3|，若r_i→0(i=1, 2, 3)，则说明发生了模态为i的故障。系统工作过程中的故障状态设定和噪声形式同实验1。

图 3为文献[13-14]的自适应阈值与本文自适应阈值分别在伽马噪声和高斯混合噪声下的对比。对比2种自适应阈值曲线，其变化趋势和范围均非常接近；表 2中，第2列和第3列、第4列和第5列相应的漏报率非常接近。可见，2种自适应阈值的漏报率基本相同，而误报率取近似值均为0。因此，采用本文改进的自适应阈值完全可等价文献[13-14]中的阈值，而本文自适应阈值的时间计算复杂度要比文献[13-14]小一个数量级，大大节约了计算时间。

图 4和图 5分别为伽马噪声和高斯混合噪声背景下对160 MW燃油机组的2种故障进行检测和隔离的结果。由残差r₀可判断，系统在100s以前处于正常工作状态，100 s以后均处于故障状态；为了充分说明故障检测的有效性，在图 4和图 5中0~100 s时间段内，选择残差和自适应阈值的极小值点分别标注了残差r₀的坐标值以及对应的自适应阈值的坐标值。图 4中，在时刻17.3 s对应的残差r₀的值为0.002 63，自适应阈值为0.008 872；在时刻98.1 s对应的r₀为0.002 21，自适应阈值为0.01。图 5中，在时刻3.4 s对应的r₀为0.002 876，自适应阈值为0.006 455。可见，r₀的值均小于对应时刻的自适应阈值，在图 5中的25 s和84.3 s满足同样的结论，说明系统处于正常工作状态。在时间段100~200 s，比较r₁~r₃可见，r₁更接近于0，说明系统发生了模态1故障，即恒偏差故障。为了更清楚比较残差大小，在图 4的199.3 s和图 5的180、199.6 s处分别标注了图中最接近的残差r₁和r₂的值，均满足r₁ < r₂；在时间段200~300 s，r₂最接近0，说明系统发生了模态2故障，即恒增益故障。虽然在图 4的250.1 s和图 5的251.5 s处r₀的值较小，但依然大于自适应阈值。可见，ACRPF故障诊断方法在不同噪声下均可实现故障的准确检测和隔离。

表 2对比了CRPF故障诊断方法和ACRPF故障诊断方法的漏报率，实验中系统运行50次，每次运行的采样点数为3 000，通过对所有采样点在50次运行中出现漏报的采样点数进行统计，并利用漏报率公式计算得到漏报率。总体来看，ACRPF故障诊断方法在不同噪声下的漏报率均小于CRPF故障诊断方法的漏报率，分别对比表 2的前2列和后2列的数据，可见本文的改进自适应阈值和文献[13-14]所得到的故障漏报率基本相同，但本文改进自适应阈值的时间计算复杂度低。结果表明，ACRPF故障诊断方法在强噪声背景下，可提高故障诊断准确性，减少系统的运行时间。

5 结论

1) CRPF算法可在未知噪声统计特性的条件下实现对非线性非高斯系统状态的递推估计，为复杂非线性系统的故障诊断提供了一种可借鉴的方法。

2) 本文通过分析CRPF算法对强噪声干扰下状态估计的性能，设计相关性判别函数实现对CRPF的状态转移密度方差的自适应调整，克服了CRPF随着时间的递推，对强噪声及时变噪声的修正能力逐渐消失而导致的误差变大甚至发散的问题。

3) 引入滑动窗利用求区间均值的方法对残差的自适应阈值进行了改进，解决了在计算残差均值和方差时系统重复运行而导致的耗时问题。

4) 应用本文方法，在不同的强噪声背景下对非线性非高斯系统进行了故障检测和隔离，对漏报率和误报率进行了分析，提高了故障诊断准确性并可减少运行时间，具有重要的理论意义和工程应用价值。