自动化学报  2017, Vol. 43 Issue (8): 1329-1338   PDF    
基于异步IMM融合滤波的网络化系统故障诊断
胡艳艳1, 金增旺1, 薛晓玲1, 孙长银2     
1. 北京科技大学自动化学院北京市工业波谱成像工程技术研究中心 北京 100083;
2. 东南大学自动化学院 南京 210096
摘要: 针对一类带随机丢包的异步多传感器网络化系统,提出了基于网络化异步交互式多模型(Interacting multiple model,IMM)融合滤波的故障诊断方法.考虑不同传感器通道具有不同丢包概率的情况,将未知的故障幅值看作扩维的系统状态,利用提出的网络化异步IMM融合滤波算法对由系统正常模型和各种可能的故障模型构成的模型集进行滤波,根据模型概率进行故障检测和定位,同时得到故障幅值和系统状态的联合估计.提出的方法避免了传统IMM故障诊断方法模型集设计中故障大小难以确定的问题,适用于具有任意采样速率和任意初始采样时刻的异步多传感器网络化系统,并且通过融合多个传感器的信息提高了故障诊断的准确性.仿真实例验证了所提出方法的可行性和有效性.
关键词: 故障诊断     网络化系统     异步融合     交互式多模型滤波    
Fault Diagnosis for Networked Systems By Asynchronous IMM Fusion Filtering
HU Yan-Yan1, JIN Zeng-Wang1, XUE Xiao-Ling1, SUN Chang-Yin2     
1. Beijing Engineering Research Center of Industrial Spectrum Imaging, School of Automation and Electrical Engineering, University of Science and Technology Beijing, Beijing 100083;
2. School of Automation, Southeast University, Nanjing 210096
Manuscript received : November 14, 2016, accepted: February 3, 2017.
Foundation Item: Supported by National Natural Science Foundation of China (61304105, 61520106009, 61533008)
Author brief: HU Yan-Yan    Lecturer at the School of Automation and Electrical Engineering, University of Science and Technology Beijing. She received her Ph. D. degree from the Department of Automation, Tsinghua University in 2011. Her research interest covers fault diagnosis, fault prediction, and information fusion.E-mail:huyanyan@ustb.edu.cn;
JIN Zeng-Wang    Ph. D. candidate at the School of Automation and Electrical Engineering, University of Science and Technology Beijing. He received his bachelor degree from the School of Automation and Electrical Engineering, University of Science and Technology Beijing in 2013. His research interest covers estimation fusion, fault diagnosis and prediction, and event-triggered systems.E-mail:b20130374@xs.ustb.edu.cn;
XUE Xiao-Ling     Master student at the School of Automation and Electrical Engineering, University of Science and Technology Beijing. She received her bachelor degree from Tianjin Normal University in 2015. Her research interest covers information fusion and fault diagnosis.E-mail:xuexiaoling@xs.ustb.edu.cn
Corresponding author. SUN Chang-Yin    Professor at the School of Automation and Electrical Engineering, University of Science and Technology Beijing. He received his master and Ph. D. degrees from the School of Automation, Southeast University in 2001 and 2003, respecttively. His research interest covers intelligent control, flight control, pattern recognition, and optimal theory. Corresponding author of this paper.E-mail:cysun@seu.edu.cn
Recommended by Associate Editor WEN Cheng-Lin
Abstract: A fault diagnosis method is proposed based on networked asynchronous interacting multiple model (IMM) fusion filtering for a kind of networked systems with multiple asynchronous sensors and stochastic packet dropouts. The proposed networked asynchronous IMM fusion filtering algorithm is used to perform fusion filtering for the model set consisting of the normal model and all kinds of possible fault models of the system, where different arriving probabilities are considered for different sensor communication channels and the unknown fault amplitude is taken as the augmented system state. Fault detection and location are achieved based on the model possibilities, and the estimates of system state and fault amplitude can be obtained simultaneously. The proposed method avoids the problem of determining fault amplitude in the model set design of traditional IMM fault diagnosis approaches, improves the accuracy of fault diagnosis by fusing information from multiple sensors, and can be used to asynchronous multi-sensor networked systems with arbitrary sampling rates and arbitrary initial sampling time instants. The feasibility and effectiveness of the proposed algorithm are illustrated by simulation examples.
Key words: Fault diagnosis     networked systems     asynchronous fusion     interacting multiple model (IMM) filtering    

随着对系统可靠性和安全性要求的不断提高, 复杂系统故障诊断问题受到了越来越多的关注, 特别是在航天系统和石油化工等安全至关重要的领域[1-2].发生故障的系统往往既含有连续时间的状态演化, 又有离散时间的参数或者结构的变化, 是一个典型的混杂系统.交互式多模型算法作为一种高效的多模型滤波估计方法, 被广泛用于各类系统的故障诊断[3].

Zhang等在文献[4]中给出了利用交互式多模型(Interacting multiple model, IMM)滤波方法进行动态系统故障诊断的框架.基本思想是建立故障候选模型集, 并利用观测数据计算更新模型集中各个模型的模型概率, 并依此进行故障诊断.与其他传统的基于模型的故障诊断方法相比, IMM故障诊断方法不但可以实现故障的检测和定位, 而且可以同时提供关于故障程度大小的信息.文献[3]主要考虑乘性故障, 且在故障模型集设计中, 需要将故障大小作为给定的模型参数.但实际中, 即使同一故障也会有不同的故障发生程度.因此, 为了能够包含尽可能多的故障情况, 候选模型集中故障模型的个数就会较多, 从而带来很大的计算负担.为了解决上述问题, Ru等[5]针对乘性故障提出了分层IMM故障诊断方法和将IMM与最大似然估计相结合的故障诊断方法.分层IMM故障诊断方法先利用第一层IMM检测并定位故障, 再用第二层IMM确定故障的大小.而在IMM极大似然方法中, 第二层IMM滤波被一个极大似然估计器代替, 并且基于极大似然估计结果得到的新的故障模型也被添加到候选模型集中, 因此该方法是一种变结构的IMM方法. Ru等在文献[6]中对后一种方法进行了更详细的讨论.为了提高故障诊断的快速性和准确性, Kim等[7]利用模糊逻辑对多模型的转移概率进行调整, 并将上述方法用于飞行器执行器故障的诊断.文献[8]同样针对IMM故障诊断方法中模型转移概率的确定问题, 提出了利用在线测量对模型转移概率进行修正的一种改进的IMM故障诊断方法.此外, 基于IMM的故障诊断方法还被用于主动容错控制[9]、卫星姿态控制[10]以及线控车辆的故障诊断[11].

上述工作考虑的都是单传感器系统, 而实际系统中为了减少不确定性, 提高测量精度, 经常使用多个传感器甚至是多个异类的传感器对系统进行观测[12-13].针对具有广播式自动相关监视传感器和多点定位传感器的空中交通管制系统, Liu等首先根据IMM滤波残差的统计特性对两个传感器的故障进行诊断, 进而通过融合两个传感器的信息实现了对飞行器状态的融合估计[14].文献[14]考虑的是同步多传感器系统, 然而实际系统中多传感器, 特别是多个异质的传感器, 由于采样速率不同或初始采样时间不同等原因往往会导致其测量在时间上并不是同步的, 而是异步的.文献[15]针对具有多个异步传感器的动态系统发生执行器故障的情况, 提出了首先利用异步IMM算法进行故障检测和定位, 进而利用扩维的卡尔曼滤波进行故障幅值估计的方法.

综上所述, 目前已有基于IMM的故障诊断方法主要考虑的是单传感器系统, 针对多传感器特别是异步多传感器系统的IMM故障诊断方法的研究还相对较少.此外, 随着传感器网络的广泛应用, 网络化系统的故障诊断问题在近年来也引起了大量关注[16-17].网络的引入带来了各种不确定性, 例如随机延迟和丢包现象等, 因此网络化系统的IMM故障诊断方法, 特别是异步多传感器网络化系统的IMM故障诊断方法面临新的挑战, 有待进一步研究.

本文研究了一类带随机数据包丢失的异步多传感器网络化系统的IMM故障诊断问题.考虑传感器的采样速率和初始采样时间都是任意的情况, 利用具有不同分布的伯努利分布描述不同传感器通道上的随机数据包丢失.主要完成了: 1) 将故障赋值作为未知量与系统状态一起构成扩维的状态向量, 并得到扩维的故障模型. 2) 对由各种可能故障情况下的故障模型与正常模型一起构成的多模型集, 利用网络化异步IMM融合滤波算法进行融合滤波. 3) 依据滤波得到的模型概率更新进行故障的检测和定位, 同时得到系统状态和幅值的融合估计.本文的主要创新点可归纳为:1) 提出了网络化异步IMM融合滤波方法, 考虑了传感器到融合中心之间的随机丢包现象, 且认为不同传感器到融合中心之间具有不同的丢包率;2) 与传统IMM故障诊断方法将故障大小作为给定模型参数不同, 本文方法在模型集设计中将故障幅值作为未知量与系统状态一起构成扩维状态向量同时进行滤波估计, 因此对一种故障只用建立一个故障候选模型, 避免了传统方法中故障大小取值难以确定, 或为了包含尽可能多的故障情况而造成的模型集中候选模型数量过多、计算量过大等问题.

1 问题描述

考虑如下一类连续时间线性随机动态系统

$\begin{align} & \overset{.}{\mathop{\mathit{\boldsymbol{x}}}}\,(t)=\ A(t)\mathit{\boldsymbol{x}}(t)+B(t)\mathit{\boldsymbol{u}}(t)+ \\ & \quad \ \quad \ G(t)\mathit{\boldsymbol{w}}(t)+F(t)\mathit{\boldsymbol{f}}(t) \\ \end{align}$ (1)

其中, $\pmb{x}(t)\in {\bf R}^{d_x}$表示$d_x$维的系统状态, $\pmb{u}(t)\in {\bf R}^{d_u}$表示$d_u$维的执行器输入, $\pmb{f}(t)\in {\bf R}^{d_f}$是故障信号, $\pmb{w}(t)$ $\in$ ${\bf R}^{d_w}$$d_w$维的系统过程噪声, $A(t)$, $B(t)$, $G(t)$$F(t)$为相应维数的矩阵.

设有N个异步传感器同时对系统(1)进行观测, 这些传感器具有不同的采样速率和初始采样时间, 并且可以是非均匀采样的.在时刻$t_{(n)}^{l}$, 第n个传感器的第l个测量表示为

$\underline{\mathit{\boldsymbol{z}}}_{(n)}^{l}=\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{H}_{(n)}^{l}\mathit{\boldsymbol{x}}(t_{(n)}^{l})+\underline{\mathit{\boldsymbol{v}}}(t_{(n)}^{l})$ (2)

其中, n表示传感器编号, $\underline{\pmb{z}}_{(n)}^{l} \in {\bf R}^{d_{\underline{z}}^{(n)}}$, ${d_{\underline{z}}^{(n)}}$是测量的维数, $\underline{H}_{(n)}^{l}$表示测量矩阵, $\underline{\pmb{v}}(t_{(n)}^{l})$是测量噪声.

传感器测量$\underline{\pmb{z}}_{(n)}^{l}$通过网络被传输到融合中心, 并且网络存在随机丢包现象.融合中心接收到的测量为

$\underline{\mathit{\boldsymbol{y}}}_{(n)}^{l}=\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{\theta }_{(n)}^{l}\underline{\mathit{\boldsymbol{z}}}_{(n)}^{l}+\left( 1-\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{\theta }_{(n)}^{l} \right)\underline{\mathit{\boldsymbol{y}}}_{(n)}^{l-1}$ (3)

其中, 随机变量$\underline{\theta}_{(n)}^{l}$取值为0或1, 表示融合中心是否接收到测量值$\underline{\pmb{z}}_{(n)}^{l}$. $\underline{\theta}_{(n)}^{l}=0$表示发生了数据丢包现象, 此时融合中心使用从传感器j接收到的离当前时刻最近的测量值, 即$\underline{\pmb{y}}_{(n)}^{l}=\underline{\pmb{y}}_{(n)}^{l-1}$.假设$\underline{\theta}_{(n)}^{l}$服从贝努利分布, 即

$\begin{align} &Pr\left\{ \underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{\theta }_{(n)}^{l}=1 \right\}={{\beta }_{(n)}}P \\ &Pr\left\{ \underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{\theta }_{(n)}^{l}=0 \right\}=1-{{\beta }_{(n)}} \\ \end{align}$

其中, $\beta_{(n)}\in[0,1]$是一个已知标量, 表示第n个传感器通信通道数据包的到达率, 这里考虑不同的传感器通道具有不同丢包概率的情况.

假设1.  系统过程噪声$\pmb{w}(t)$和传感器测量噪声$\underline{\pmb{v}}(t_{(n)}^{l})$是零均值的高斯白噪声, 且传感器的测量噪声之间互相独立, 其协方差矩阵分别为$E\{\pmb{w}(t)\pmb{w}^{\rm T}(t')\}$ $=Q(t)\delta_{tt'}$, $E\{\underline{\pmb{v}}(t_{(n)}^{l})\underline{\pmb{v}}^{\rm T}(t_{(n')}^{l'})\}=R_{(n)}^{l}\delta_{nn'}\delta_{ll'}$, 其中$\delta$是克罗内克函数, $\pmb{w}(t)$, $\underline{\pmb{v}}(t_{(n)}^{l})$和初始状态$\pmb{x}(0) $是相互独立的.

假设2.  不同传感器通信通道的随机丢包是互相独立的, 且相同传感器通信通道在不同时刻的随机丢包也是互相独立的, 即随机变量$\underline{\theta}_{(n)}^{l}$$\underline{\theta}_{(n')}^{l'}$n $\neq n'$$l \neq l'$情况下是相互独立的, 也就是说$E\{\underline{\theta}_{(n)}^{l}(\underline{\theta}_{(n')}^{l'})^{\rm T}\}$ $=\beta_{(n)}\delta_{nn'}\delta_{ll'}$.同时, 随机丢包$\underline{\theta}_{(n)}^{l}$与过程噪声$\pmb{w}(t)$、传感器测量噪声$\underline{\pmb{v}}(t_{(n)}^{l})$和系统初始状态$\pmb{x}(0) $相互独立.

$t_{k-1}$表示融合中心上一个融合时刻, $t_k$表示当前融合时刻, 则当$k=1$时, $t_{k-1}=0$表示初始时刻.因为传感器是异步工作的, 所以在一个融合周期内, 单个传感器可能有不止一个的传感器测量.令$N_k^{(n)}$表示传感器n在融合区间$(t_{k-1}, t_k]$内的测量数, 则$N_k^{(n)}=0$表示传感器n在当前融合区间内没有测量.另外, 记$N_k$为融合区间$(t_{k-1}, t_k]$N个异步传感器的测量总个数, 则有$N_k=\sum_{n=1}^{N}N_k^{(n)}$.

对第k个融合区间内$N_k$个异步传感器测量按时间顺序排列后得到测量序列$\left\{\pmb{y}_k^i\right\}_{i=1}^{N_k}$.其中, 若排序后的第i个测量$\pmb{y}_k^i$是第n个传感器的第l个测量, 则可以表示为$\pmb{y}_k^i=\underline{\pmb{y}}_{(n)}^{l}$.用$t_k^i$表示测量$\pmb{y}_k^i$的采样时间, 则$t_k^i=t_{(n)}^{l}$, 并且有$t_{k-1} \leq t_k^1\leq t_k^2\leq$ $\cdots \leq t_k^{N_k} \leq t_k$.当两个传感器测量是同一时刻的测量时, 上式中的等号成立. 图 1刻画了$(k, i)$$(n, l)$之间的对应关系, 其中$i \in\{1, 2, \cdots, N_k\}$, $ n$ $\in$ $\{1$, $2$, $\cdots$, $N\}$, $l\in {\bf N}$.由式(2)、式(3) 和上述对应关系得到

图 1 融合区间$(t_{k-1}, t_k]$内异步多传感器网络化测量 Figure 1 Networked measurements from asynchronous multi-sensors during fusion interval $(t_{k-1}, t_k]$
$\begin{align} & \mathit{\boldsymbol{y}}_{k}^{i}=\theta _{k}^{i}\mathit{\boldsymbol{z}}_{k}^{i}+(1-\theta _{k}^{i})\mathit{\boldsymbol{y}}_{k}^{{{i}^{-}}} \\ \end{align}$ (4)
$\begin{align} \mathit{\boldsymbol{z}}_{k}^{i}=H_{k}^{i}\mathit{\boldsymbol{x}}(t_{k}^{i})+\mathit{\boldsymbol{v}}(t_{k}^{i}) \\ \end{align}$ (5)

其中, $ \pmb{z}_k^i=\underline{\pmb{z}}_{(n)}^{l}$, $\pmb{y}_k^{i^-}=\underline{\pmb{y}}_{(n)}^{l-1}$, $H_k^i=\underline{H}_{(n)}^{l}$, $\pmb{v}_k^i=\underline{\pmb{v}}_{(n)}^{l}$, $\theta_k^i$ = $\underline{\theta}_{(n)}^{l}$. $d_z^i$是测量$\pmb{z}_k^i$的维数, 且有$d_z^i=d_{\underline{z}}^{(n)}$.从假设1和假设2可知, $\pmb{v}_k^i$是零均值高斯白噪声, 且有其协方差矩阵为

$E\{\pmb{v}_k^i(\pmb{v}_k^{j})^{\rm T}\}=\underline{R}_{(n)}^{l}\delta_{ij}=R_k^i\delta_{ij}$ (6)

$\theta_k^i$$\theta_k^{j}$相互独立, 且

$E\{\theta_k^i(\theta_k^{j})^{\rm T}\}=\underline{\beta}_{(n)}^{l}\delta_{ij} =\beta_k^i\delta_{ij}$ (7)

其中, $i, j \in {1, 2, \cdots, N_k}$.

本文提出的故障诊断方法的核心思想是:利用融合中心收到的异步多传感器测量对各种可能故障情况下的故障模型和正常情况下的系统模型一起构成的模型集进行网络化异步IMM融合滤波, 并基于滤波得到的模型概率和状态估计同时进行故障的联合检测、定位和估计, 即判断故障是否发生、确定故障发生的位置以及估计故障幅值的大小.为此, 首先给出网络化IMM融合滤波的模型集设计.

2 模型集设计

考虑系统(1) 某一时刻仅发生单个故障的情况, 且假设故障幅值保持不变或者变化相对缓慢.则当系统第m维发生故障时, $m=1, \cdots, d_f$, 故障$\pmb{f}(t)$可表示为$\pmb{f}(t)=\pmb{e}_mf$, 其中标量f表示故障幅值, $\pmb{e}_m$是单位矩阵的第m列, 表示故障的方向, 则发生第m种故障情况下对应的系统模型为

$\begin{align} & \overset{.}{\mathop{\mathit{\boldsymbol{x}}}}\,(t)=A(t)\mathit{\boldsymbol{x}}(t)+B(t)\mathit{\boldsymbol{u}}(t)+ \\ & \quad \quad G(t)\mathit{\boldsymbol{w}}(t)+F(t){{\mathit{\boldsymbol{e}}}_{m}}f \\ \end{align}$ (8)

对应的离散时间状态转移方程为

$\begin{align} & \mathit{\boldsymbol{x}}({{t}_{2}})=\Phi ({{t}_{2}},{{t}_{1}})\mathit{\boldsymbol{x}}({{t}_{1}})+\mathit{\boldsymbol{u}}({{t}_{2}},{{t}_{1}})+ \\ & \quad \quad \quad {{\mathit{\boldsymbol{ }}\!\!\xi\!\!\rm{ }}^{m}}({{t}_{2}},{{t}_{1}})f+\mathit{\boldsymbol{w}}({{t}_{2}},{{t}_{1}}) \\ \end{align}$ (9)

其中, $\Phi(t_2, t_1) $为系统(1) 从$t_1$$t_2$的状态转移矩阵,

$\mathit{\boldsymbol{u}}({{t}_{2}},{{t}_{1}})=\int_{{{t}_{1}}}^{{{t}_{2}}}{\Phi }({{t}_{2}},\tau )B(\tau )\mathit{\boldsymbol{u}}(\tau )\rm{d}\tau $ (10)
${{\mathit{\boldsymbol{ }}\!\!\xi\!\!\rm{ }}^{m}}({{t}_{2}},{{t}_{1}})=\int_{{{t}_{1}}}^{{{t}_{2}}}{\Phi }({{t}_{2}},\tau )F(\tau ){{\mathit{\boldsymbol{e}}}_{m}}\rm{d}\tau $ (11)
$\mathit{\boldsymbol{w}}({{t}_{2}},{{t}_{1}})=\int_{{{t}_{1}}}^{{{t}_{2}}}{\Phi }({{t}_{2}},\tau )G(\tau )w(\tau )\rm{d}\tau $ (12)

且有

$\begin{align} & \Omega ({{t}_{2}},{{t}_{1}})=E\left\{ \mathit{\boldsymbol{w}}({{t}_{2}},{{t}_{1}}){{\mathit{\boldsymbol{w}}}^{\rm{T}}}({{t}_{2}},{{t}_{1}}) \right\}= \\ & \quad \quad \quad \quad \quad \int_{{{t}_{1}}}^{{{t}_{2}}}{\Phi }({{t}_{2}},\tau )G(\tau )Q(\tau ){{G}^{\rm{T}}}(\tau ){{\Phi }^{\rm{T}}}({{t}_{2}},\tau )\rm{d}\tau \\ \end{align}$ (13)

则相应地从上一个融合时刻$t_{k-1}$到当前融合时刻$t_k$

$\begin{array}{l} \mathit{\boldsymbol{x}}({t_k}) = \Phi ({t_k},{t_{k - 1}})\mathit{\boldsymbol{x}}({t_{k - 1}}) + \mathit{\boldsymbol{u}}({t_k},{t_{k - 1}}) + \\ \quad \quad \quad {\mathit{\boldsymbol{\xi }}^m}({t_k},{t_{k - 1}})f + \mathit{\boldsymbol{w}}({t_k},{t_{k - 1}}) \end{array}$ (14)

定义

${\mathit{\boldsymbol{\overline x}} _k} = \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{x}}({t_k})}\\ f \end{array}} \right]$ (15)
$\bar \Phi _{k - 1}^m = \left[ {\begin{array}{*{20}{c}} {\Phi ({t_k},{t_{k - 1}})} & {{\mathit{\boldsymbol{\xi }}^m}({t_k},{t_{k - 1}})}\\ 0 & 1 \end{array}} \right]$ (16)
$\mathit{\boldsymbol{\overline u}} _{k - 1}^x = \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{u}}({t_k},{t_{k - 1}})}\\ 0 \end{array}} \right]$ (17)
${\mathit{\boldsymbol{\overline w}} _{k - 1}} = \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{w}}({t_k},{t_{k - 1}})}\\ 0 \end{array}} \right]$ (18)

则根据式(14) 和常值故障或缓变故障的假设, 经状态扩维有

${\mathit{\boldsymbol{\overline x}} _k} = \bar \Phi _{k - 1}^m{\mathit{\boldsymbol{\overline x}} _{k - 1}} + \mathit{\boldsymbol{\overline u}} _{k - 1}^x + {\mathit{\boldsymbol{\overline w}} _{k - 1}}$ (19)

且有

${{\bar \Omega }_{k - 1}} = E\left\{ {{{\mathit{\boldsymbol{\overline w}} }_{k - 1}}\mathit{\boldsymbol{\overline w}} _{k - 1}^{\rm{T}}} \right\} = {\rm{diag}}\left\{ {\Omega ({t_k},{t_{k - 1}}),0} \right\}$ (20)

式(19) 即为第m种故障对应的故障模型, 令$m= 0$表示系统正常运行时的状态模型, 此时$\pmb{\xi}^0=$ $\pmb{0}$, f没有实际意义, 只是为了使系统正常情况下的模型和故障模型在形式上能够统一.则当$m=0, 1$, $\cdots$, $d_f$时, 对应的$d_f+1$个模型一起构成后续IMM融合滤波的模型集.

注1.从上述模型集的构造过程中可以看出, 本文是将系统故障分为故障幅值f和故障方向$\pmb{e}_m$两部分, 利用故障方向来区分不同的故障, 而将故障幅值f作为未知量, 与系统状态$\pmb{x}$一起构成扩维的状态向量$\bar{ \pmb{x}}$进行滤波估计.因此, 与传统IMM故障诊断方法不同, 本文提出的方法对一种故障只需建立一个候选故障模型.

注2.虽然本文为了论述简单, 假设系统在同一时刻仅发生单个故障, 但本文提出的方法可以直接推广到系统同时发生多个故障的情况.

注3.本文假设故障是定常的或者慢时变的, 对于时变故障, 可采用基于强跟踪滤波器[18]的IMM滤波算法, 利用强跟踪滤波器对模型不确定性极强的鲁棒性和对状态突变具有的强跟踪能力实现对时变故障的跟踪估计.

3 网络化异步IMM融合滤波算法

众所周知, 多模型滤波的基本思想是同时运行一组基本滤波器, 其中每个基本滤波器基于系统的一个候选模型.与其他多模型算法相比, IMM能够用较小的计算量获得较好的滤波估计效果, 关键在于IMM对每个基本滤波器在每个滤波周期的初始加入了一个有效的输入混合步骤, 即滤波器的重启过程.与标准IMM算法类似, 网络化异步IMM融合滤波算法同样包括滤波器重启、故障情况下滤波、故障概率更新和输出融合四个步骤.

3.1 滤波器重启

$\pmb{Y}_k:=\left\{\pmb{y}_k^i\right\}_{i=1}^{N_k}$, $\pmb{Y}^k:=\left\{\pmb{Y}_l\right\}_{l=1}^{k}$, 则$\pmb{Y}^k$表示到当前融合时刻$t_k$为止融合中心收到的所有异步传感器累积测量集.

定义

$\mathit{\boldsymbol{\widehat x}}_k^m = E\left\{ {{{\mathit{\boldsymbol{\overline x}} }_k}|\rlap{--} m({t_k}) = m,{\rm{ }}{\mathit{\boldsymbol{Y}}^k}} \right\}$ (21)
$P_k^m = Cov\left\{ {{{\mathit{\boldsymbol{\overline x}} }_k}|\rlap{--} m({t_k}) = m,{\mathit{\boldsymbol{Y}}^k}} \right\}$ (22)
$\mu _k^m = Pr\left\{ {\rlap{--} m({t_k}) = m|{\mathit{\boldsymbol{Y}}^k}} \right\}$ (23)

其中, $\rlap{---}m(t_k)=m$表示$t_k$时刻第m个模型有效, 则$\mu_k^m$表示$t_k$时刻模型m为有效模型的后验概率.则重启后第m个基本滤波器的滤波初始值为[3]

$\begin{align} \check{\pmb{{x}}}_{k-1}^m&=\sum\limits_{i=0}^{d_f}\hat{\pmb{{x}}}_{k-1}^i\check{\mu}_{k-1}^{i|m}\end{align}$ (24)
$\begin{align}\check{P}_{k-1}^m&=\sum\limits_{i=0}^{d_f}[P_{k-1}^i+(\check{\pmb{{x}}}_{k-1}^m-\hat{\pmb{{x}}}_{k-1}^i) (\cdot)^{\rm T}]\check{\mu}_{k-1}^{i|m}\end{align}$ (25)
$\begin{align}\check{\mu}_{k-1}^{i|m}&=\frac{\lambda^{i, m}\mu^i_{k-1}}{\sum\limits_{i=0}^{d_f} \lambda^{i, m}\mu^i_{k-1}}\end{align}$ (26)

其中, $(\cdot)$表示与之前括弧中相同的内容, $m=0, 1$, $\cdots$, $d_f$, $\lambda^{i, m}=Pr\{\rlap{---}m(t_k)=m|\rlap{---}m(t_{k-1})=i\}$是系统从$t_{k-1}$时刻模型i转移到$t_k$时刻模型m的转移概率, 且$\Lambda=(\lambda^{i, m})$为相应的一步转移概率矩阵, 并有$0\leq \lambda^{i, m} \leq 1$, $\Sigma_m \lambda^{i, m} =1 $.

3.2 故障条件下的滤波

引理1.定义

$\begin{align} \bar{\pmb{y}}_k=\left[(\pmb{y}_k^1) ^{\rm T}, (\pmb{y}_k^2) ^{\rm T}, \cdots, (\pmb{y}_k^{N_k})^{\rm T}\right]^{\rm T} \label{eq:y_kstacked} \end{align}$ (27)

则在给定第m种故障情况下, $\pmb{y}_k$可以看作融合时刻$t_k$关于扩维的状态$\bar{\pmb{x}}_k$的等价测量, 其等价测量方程为

$\begin{align}\label{eq:bary} \bar{\pmb{y}}_k=\Theta_k\bar{H}^m_k\bar{\pmb{x}}_k+\Theta_k\bar{\pmb{u}}_k+\Theta_k\bar{\pmb{\eta}}_k+(I-\Theta_k)\bar{\pmb{y}}_{k}^{-} \end{align}$ (28)

其中,

${{\Theta }_{k}}=\text{diag}\{\theta _{k}^{i}{{I}_{d_{z}^{i}}}\}_{i=1}^{{{N}_{k}}}$ (29)
$\bar{H}_{k}^{m}={{\left[ {{\left( \bar{H}_{k}^{m,1} \right)}^{\text{T}}},{{\left( \bar{H}_{k}^{m,2} \right)}^{\text{T}}},\cdots ,{{\left( \bar{H}_{k}^{m,{{N}_{k}}} \right)}^{\text{T}}} \right]}^{\text{T}}}$ (30)
${{\overline{\mathit{\boldsymbol{u}}}}_{k}}={{\left[ {{\left( \mathit{\boldsymbol{u}}_{k}^{1} \right)}^{\rm{T}}},{{\left( \mathit{\boldsymbol{u}}_{k}^{2} \right)}^{\rm{T}}},\cdots ,{{\left( \mathit{\boldsymbol{u}}_{k}^{{{N}_{k}}} \right)}^{\rm{T}}} \right]}^{\rm{T}}}$ (31)
${\mathit{\boldsymbol{\overline \eta }} _k} = {\left[ {{{\left( {\mathit{\boldsymbol{\eta }}_k^1} \right)}^{\rm{T}}},{{\left( {\mathit{\boldsymbol{\eta }}_k^2} \right)}^{\rm{T}}}, \cdots ,{{\left( {\mathit{\boldsymbol{\eta }}_k^{{N_k}}} \right)}^{\rm{T}}}} \right]^{\rm{T}}}$ (32)
$\mathit{\boldsymbol{\overline y}} _k^ - = {\left[ {{{\left( {\mathit{\boldsymbol{y}}_k^{{1^ - }}} \right)}^{\rm{T}}},{{\left( {\mathit{\boldsymbol{y}}_k^{{2^ - }}} \right)}^{\rm{T}}}, \cdots ,{{\left( {\mathit{\boldsymbol{y}}_k^{N_k^ - }} \right)}^{\rm{T}}}} \right]^{\rm{T}}}$ (33)

并且

$\bar H_k^{m,i} = \Pi _k^i\left[ {\begin{array}{*{20}{c}} I&{ - {\boldsymbol{\xi }^m}({t_k},t_k^i)} \end{array}} \right]$ (34)
$\mathit{\boldsymbol{\eta }}_k^i = \mathit{\boldsymbol{v}}_k^i - \Pi _k^i\mathit{\boldsymbol{w}}({t_k},t_k^i)$ (35)
$\mathit{\boldsymbol{u}}_k^i = - \Pi _k^i\mathit{\boldsymbol{u}}({t_k},t_k^i)$ (36)
$\Pi _k^i = H_k^i{\Phi ^{ - 1}}({t_k},t_k^i)$ (37)

证明.由式(9) 可知

$\begin{align}\label{eq:xt_ktki} \pmb{x}(t_k)=&\ \Phi(t_k, t_k^{i})\pmb{x}(t_k^{i})+\pmb{u}(t_k, t_k^{i}) + \nonumber\\ &\ \pmb{\xi}^m(t_k, t_k^{i})f+\pmb{w}(t_k, t_k^{i}) \end{align}$ (38)

将式(38) 代入测量方程(5), 可得

$\begin{align}\label{eq:z_k^itrans} \pmb{z}_k^{i}=&\ H_k^{i}\Phi^{-1}(t_k, t_k^{i})\left[\pmb{x}(t_k)-\pmb{u}(t_k, t_k^{i}) -\right.\notag\\ &\ \left.\pmb{\xi}^m(t_k, t_k^{i})f-\pmb{w}(t_k, t_k^{i})\right]+\pmb{v}_k^{i}=\notag\\ &\ \bar{H}_k^{m, i}\bar{\pmb{x}}^m_k+\pmb{u}_k^{i}+\pmb{\eta}_k^{i} \end{align}$ (39)

将式(39) 代入式(4), 可得

$\begin{align} \pmb{y}_k^i=\theta_k^i \bar{H}_k^{m, i}\bar{\pmb{x}}_k+\theta_k^i\pmb{u}_k^{i}+\theta_k^i\pmb{\eta}_k^{i}+(1-\theta_k^i)\pmb{y}_k^{i^-} \end{align}$ (40)

则式(27) 可由式(40) 及定义(29) ~ (33) 通过扩维直接得到.

引理2.等价测量噪声$\bar{\pmb{\eta}}_k$是均值为零的高斯白噪声, 其协方差矩阵$\bar{R}_k=E\{\bar{\pmb{\eta}}_{k}\bar{\pmb{\eta}}_{k}^{\rm T}\}=(\bar{R}_k^{i, j})$, 其中$\bar{R}^{i, j}_k$表示$\bar{R}_k$中第ij列的子矩阵块, i, $j\in\{1$, $2$, $\cdots$, $N_k\}$, 且有

$\begin{align}\label{eq:eta_k^ieta_k^j} \bar{R}_k^{i, j}=&\ E\{\pmb{\eta}_k^i(\pmb{\eta}_k^j)^{\rm T}\}=\notag\\ &\ \begin{cases} R_k^i+\Pi_k^i\Omega(t_k, t_k^i)(\Pi_k^i)^{\rm T},&i=j\\ \Pi_k^i \Omega(t_k, t_k^{max\{i, j\}})(\Pi_k^j)^{\rm T},&i\neq j \end{cases} \end{align}$ (41)

同时, 由于共同的过程噪声的影响, 等价测量噪声$\bar{\pmb{\eta}}_k$和过程噪声$\pmb{w}(t_k, t_{k-1})$也是相关的, 且有

$\begin{align}\label{eq:omegaandeta} \Psi_k=&\ E\left\{\pmb{w}(t_k, t_{k-1})(\bar{\pmb{\eta}}_k)^{\rm T}\right\}=\notag\\&-\left[\Omega(t_k, t_k^1) (\Pi_k^1) ^{\rm T}, \Omega(t_k, t_k^2) (\Pi_k^2) ^{\rm T}, \cdots, \right.\notag\\ &\left. \Omega(t_k, t_k^{N_k})(\Pi_k^{N_k})^{\rm T}\right] \end{align}$ (42)

证明.由假设1、式(32)、式(35) 和式(12) 可知, $\bar{\pmb{\eta}}_k$是均值为零的高斯白噪声, 且

$\begin{align} E\left\{\pmb{\eta}_k^i(\pmb{\eta}_k^i)^{\rm T}\right\}= R_k^i+\Pi_k^i\Omega(t_k, t_k^i)(\Pi_k^j)^{\rm T} \end{align}$ (43)

同时, 当$i\neq j$时, 有

$\begin{align}\label{eq:etaeta} &E\left\{\pmb{\eta}_k^i(\pmb{\eta}_k^j)^{\rm T}\right\}=\notag\\ &\qquad E\left\{[\pmb{v}_k^{i}-\Pi_k^i\pmb{w}(t_k, t_k^{i})][\pmb{v}_k^{j}-\Pi_k^j\pmb{w}(t_k, t_k^{j})]^{\rm T}\right\}=\notag\\ &\qquad\Pi_k^iE\left\{\pmb{w}(t_k, t_k^i)\pmb{w}^{\rm T}(t_k, t_k^j)\right\}(\Pi_k^j)^{\rm T}=\notag\\ &\qquad\Pi_k^i \Omega\left(t_k, t_k^{\max\{i, j\}}\right)(\Pi_k^j)^{\rm T} \end{align}$ (44)

类似地, 由式(35) 和式(12) 可得

$\begin{array}{l} E\left\{ {\mathit{\boldsymbol{w}}({t_k},{t_{k - 1}}){{(\mathit{\boldsymbol{\eta }}_k^i)}^{\rm{T}}}} \right\} = \\ \qquad E\left\{ {\mathit{\boldsymbol{w}}({t_k},{t_{k - 1}}){{[\mathit{\boldsymbol{v}}_k^i - \Pi _k^i\mathit{\boldsymbol{w}}({t_k},t_k^i)]}^{\rm{T}}}} \right\} = \\ \qquad - \Omega ({t_k},t_k^i){(\Pi _k^i)^{\rm{T}}} \end{array}$ (45)

定理1.根据引理1和引理2, 对由扩维状态方程(19) 和等价测量方程(28) 组成的系统, $t_k$时刻扩维状态$\bar{\pmb{x}}_k$在第m种故障条件下, 基于$(t_{k-1}, t_k]$内异步网络化测量$\pmb{Y}_k$的融合估计$\hat{{\pmb{x}}}_k^m$及其误差协方差$P_k^m$的更新方程为

$\begin{align} \hat{{\pmb{x}}}_{k|k-1}^m=E\left\{\bar{\pmb{x}}_k|\rlap{---}m(t_k)=m, \pmb{Y}^{k-1}\right\}=\\ \qquad\qquad{\Phi}_{k-1}^m\check{\pmb{x}}^m_{k-1}+\bar{\pmb{u}}^x_{k-1} \end{align}$ (46)
$\begin{align} P_{k|k-1}^m=Cov\left\{\bar{\pmb{x}}_k|\rlap{---}m(t_k)=m, \pmb{Y}^{k-1}\right\}=\\ \qquad\qquad\bar{\Phi}_{k-1}^m\check{\pmb{P}}^m_{k-1}(\bar{\Phi}_{k-1}^m)^{\rm T}+\bar{\Omega}_{k-1} \end{align}$ (47)
$\begin{align} \hat{{\pmb{x}}}_{k}^m=\hat{{\pmb{x}}}_{k|k-1}^m+\Gamma_k^m \left (C_k^m\right)^{-1}\tilde{{\pmb{y}}}_{k|k-1}^m \end{align}$ (48)
$\begin{align} P_{k}^m=P_{k|k-1}^m-\Gamma_k^m \left (C_k^m\right)^{-1}(\Gamma_k^m)^{\rm T} \end{align}$ (49)

其中,

$\begin{array}{l} \mathit{\boldsymbol{\widetilde y}}_{k|k - 1}^m = {\mathit{\boldsymbol{\overline y}} _k} - E\left\{ {{{\mathit{\boldsymbol{\overline y}} }_k}|\rlap{--} m({t_k}) = m,{Y^{k - 1}}} \right\} = \\ \quad \quad \quad {\mathit{\boldsymbol{\overline y}} _k} - {\Xi _k}\bar H_k^m\mathit{\boldsymbol{\widehat x}}_{k|k - 1}^m - {\Xi _k}{\mathit{\boldsymbol{\overline u}} _k} - (I - {\Xi _k})\mathit{\boldsymbol{\overline y}} _k^ - \end{array}$ (50)
$\Gamma _k^m = P_{k|k - 1}^m{(\bar H_k^m)^{\rm{T}}}{\Xi _k} + {\bar \Psi _k}{\Xi _k}$ (51)
${\Xi _k} = E\{ {\Theta _k}\} = {\rm{diag}}\{ \beta _k^i{I_{d_z^i}}\} _{i = 1}^{{N_k}}$ (52)
${\bar \Psi _k} = {\left[ {\Psi _k^{\rm{T}},{\rm{ }}{{\bf{0}}_{\sum\limits_{i = 1}^{{N_k}} {d_z^i \times 1} }}} \right]^{\rm{T}}}$ (53)
$C_k^m = C_{k,1}^m + C_{k,2}^m + C_{k,3}^m + C_{k,4}^m + {(C_{k,4}^m)^{\rm{T}}}$ (54)
$\begin{array}{l} C_{k,1}^m = {\Xi _k}\bar H_k^mP_{k|k - 1}^m{(\bar H_k^m)^{\rm{T}}}{\Xi _k} + ({\Xi _k} - \Xi _k^2) \times \\ \quad \quad \quad {\rm{diag}}\{ \bar H_k^{m,i}P_{k|k - 1}^m{(\bar H_k^{m,i})^{\rm{T}}}\} _{i = 1}^{{N_k}} \end{array}$ (55)
$\begin{array}{l} C_{k,2}^m = ({\Xi _k} - \Xi _k^2){\rm{diag}}\{ (\bar H_k^{m,i}\mathit{\boldsymbol{\hat x}}_{k|k - 1}^m + \mathit{\boldsymbol{\overline u}} _k^i - \mathit{\boldsymbol{\overline y}} _k^{i - }) \times \\ \quad \quad \quad {(\bar H_k^{m,i}\mathit{\boldsymbol{\hat x}}_{k|k - 1}^m + \mathit{\boldsymbol{\overline u}} _k^i - \mathit{\boldsymbol{\overline y}} _k^{i - })^{\rm{T}}}\} _{i = 1}^{{N_k}} \end{array}$ (56)
$\begin{array}{l} C_{k,3}^m = {\Xi _k}{{\bar R}_k}{\Xi _k} + ({\Xi _k} - \Xi _k^2) \times \\ \quad \quad {\rm{diag}}\{ R_k^i + \Pi _k^i\Omega ({t_k},t_k^i){(\Pi _k^i)^{\rm{T}}}\} _{k = i}^{{N_k}} \end{array}$ (57)
$\begin{array}{l} C_{k,4}^m = Cov\{ {\Theta _k}\bar H_k^m{\mathit{\boldsymbol{\overline w}} _{k - 1}},{\Theta _k}{\mathit{\boldsymbol{\overline \eta }} _k}\} = \\ \quad \quad \quad {\Xi _k}\bar H_k^m{{\bar \Psi }_k}{\Xi _k} + (\Xi _k^2 - {\Xi _k}) \times \\ \quad \quad \quad {\rm{diag}}\{ \bar H_k^{m,i}{[\Pi _k^i\Omega ({t_k},t_k^i),{{\bf{0}}_{d_z^i \times 1}}]^{\rm{T}}}\} _{k = i}^{{N_k}} \end{array}$ (58)

证明.对离散系统方程(19) 和等价测量方程(28) 利用递归形式的滤波估计器可得[19]

$\begin{array}{l} \mathit{\boldsymbol{\widehat x}}_k^m = \;\mathit{\boldsymbol{\widehat x}}_{k|k - 1}^m + \;\\ \quad \quad \quad Cov\{ \mathit{\boldsymbol{\widetilde x}}_{k|k - 1}^m,\mathit{\boldsymbol{\widetilde y}}_{k|k - 1}^m|\rlap{--} m({t_k}) = m,{Y^{k - 1}}\} \times \;\\ \quad \quad \quad {\left( {Cov\{ \mathit{\boldsymbol{\widetilde y}}_{k|k - 1}^m|\rlap{--} m({t_k}) = m,{Y^{k - 1}}\} } \right)^{ - 1}}\mathit{\boldsymbol{\widetilde y}}_{k|k - 1}^m \end{array}$ (59)
$\begin{array}{l} P_k^m = P_{k|k - 1}^m - \\ \quad \quad \quad Cov\left\{ {\mathit{\boldsymbol{\widetilde x}}_{k|k - 1}^m,\mathit{\boldsymbol{\widetilde y}}_{k|k - 1}^m|\rlap{--} m({t_k}) = m,{Y^{k - 1}}} \right\} \times \\ \quad \quad \quad \;{\left( {Cov\{ \mathit{\boldsymbol{\widetilde y}}_{k|k - 1}^m|\rlap{--} m({t_k}) = m,{Y^{k - 1}}\} } \right)^{ - 1}} \times \\ \quad \quad \quad \left. {Cov\{ \mathit{\boldsymbol{\widetilde y}}_{k|k - 1}^m,\mathit{\boldsymbol{\widetilde x}}_{k|k - 1}^m|\rlap{--} m({t_k}) = m,{\mathit{\boldsymbol{Y}}^{k - 1}}} \right\} \end{array}$ (60)

其中,

$\begin{align}\label{eq:tildex} \tilde{\pmb{{x}}}^m_{k|k-1}=&\ \bar{{\pmb{x}}}_k-\hat{{\pmb{x}}}_{k|k-1}^m=\notag\\ &\ \bar{\Phi}_{k-1}^m(\bar{\pmb{x}}_{k-1}-\check{\pmb{{x}}}^m_{k-1})+\bar{\pmb{w}}_{k-1}=\notag\\ &\ \bar{\Phi}_{k-1}^m\tilde{\pmb{x}}^m_{k-1}+\bar{\pmb{w}}_{k-1} \end{align}$ (61)

以下本定理的证明中, 在不引起歧义的前提下, 为了描述简便, 将条件方差$Cov\{\cdot|\rlap{---}m(t_k)=m$, ${\pmb Y}^{k-1}\}$统一简写为无条件方差$Cov\{\cdot\}$.

根据引理1, 将式(28) 代入式(50), 可得

$\begin{align}\label{eq:tildebmyk_k-1} \tilde{{\pmb{y}}}_{k|k-1}^m=&\ \Theta_k\bar{H}^m_k\tilde{\pmb{{x}}}^m_{k|k-1}+(\Theta_k-\Xi_k) \bar{H}^m_k\hat{\pmb{{x}}}^m_{k|k-1} +\notag\\ &\ (\Theta_k-\Xi_k)(\bar{\pmb{u}}_k-\bar{\pmb{y}}_{k}^{-})+\Theta_k\bar{\pmb{\eta}}_k \end{align}$ (62)

由式(62)、式(61) 和假设2可知

$\begin{array}{l} \Gamma _k^m = \;Cov\{ \mathit{\boldsymbol{\widetilde x}}_{k|k - 1}^m,\mathit{\boldsymbol{\widetilde y}}_{k|k - 1}^m\} = \\ \quad \quad \;Cov\{ \mathit{\boldsymbol{\widetilde x}}_{k|k - 1}^m,{\Theta _k}\bar H_k^m\mathit{\boldsymbol{\widetilde x}}_{k|k - 1}^m + {\Theta _k}{\mathit{\boldsymbol{\overline \eta }} _k}\} = \\ \quad \quad \;P_{k|k - 1}^m{(\bar H_k^m)^{\rm{T}}}{\Xi _k} + Cov\{ \mathit{\boldsymbol{\widetilde x}}_{k|k - 1}^m,{\Theta _k}{\mathit{\boldsymbol{\overline \eta }} _k}\} \end{array}$ (63)

并且由引理2和式(18) 可得

$\begin{align}\label{eq:covxeta} Cov&\{\tilde{{\pmb{x}}}_{k|k-1}^m, \Theta_k\bar{\pmb{\eta}}_k\}=\notag\\ &Cov\{\bar{\Phi}_{k-1}^m\tilde{\pmb{x}}^m_{k-1}+\bar{\pmb{w}}_{k-1}, \Theta_k\bar{\pmb{\eta}}_k\}=\notag\\ &Cov\{\bar{\pmb{w}}_{k-1}, \Theta_k\bar{\pmb{\eta}}_k\}=\notag\\ &\bar{\Psi}_k\Xi_k \end{align}$ (64)

其中, $\bar{\Psi}_k=[{\Psi}_k^{\rm T} , \pmb{0}_{\sum_{i=1}^{N_k}{d_z^i\times 1}}]^{\rm T}$.因此, 将式(64) 代入式(63), 即可得式(51).

另外, 由式(62) 可知

$\begin{array}{l} C_k^m = \;Cov\{ \mathit{\boldsymbol{\widetilde y}}_{k|k - 1}^m\} = \\ \quad \quad \quad Cov\{ {\Theta _k}\bar H_k^m\mathit{\boldsymbol{\widetilde x}}_{k|k - 1}^m + {\Theta _k}{\mathit{\boldsymbol{\overline \eta }} _k} + ({\Theta _k} - {\Xi _k}) \times \\ \quad \quad \quad \;(\bar H_k^m\mathit{\boldsymbol{\widehat x}}_{k|k - 1}^m + {\mathit{\boldsymbol{\overline u}} _k} - \mathit{\boldsymbol{\overline y}} _k^ - )\} \end{array}$ (65)

其中, 由假设1、假设2及正交性原理可知

$Cov\{ {\Theta _k}\bar H_k^m\mathit{\boldsymbol{\widetilde x}}_{k|k - 1}^m,({\Theta _k} - {\Xi _k})\bar H_k^m\mathit{\boldsymbol{\widehat x}}_{k|k - 1}^m\} = 0$ (66)
$Cov\{ {\Theta _k}\bar H_k^m\mathit{\boldsymbol{\widetilde c}}_{k|k - 1}^m,({\Theta _k} - {\Xi _k})({\mathit{\boldsymbol{\overline u}} _k} - \mathit{\boldsymbol{\overline y}} _k^ - )\} = 0$ (67)
$Cov\{ {\Theta _k}{\mathit{\boldsymbol{\overline \eta }} _k},({\Theta _k} - {\Xi _k})\bar H_k^m\mathit{\boldsymbol{\widehat x}}_{k|k - 1}^m\} = 0$ (68)
$Cov\{ {\Theta _k}{\mathit{\boldsymbol{\overline \eta }} _k},({\Theta _k} - {\Xi _k})({\mathit{\boldsymbol{\overline u}} _k} - \mathit{\boldsymbol{\overline y}} _k^ - )\} = 0$ (69)

因此有式(54) 成立, 并且

$\begin{array}{l} C_{k,1}^m = Cov\{ {\Theta _k}\bar H_k^m\mathit{\boldsymbol{\widetilde x}}_{k|k - 1}^m\} = \\ \quad \quad {\Xi _k}\bar H_k^mP_{k|k - 1}^m{(\bar H_k^m)^{\rm{T}}}{\Xi _k} + \\ \quad \quad ({\Xi _k} - \Xi _k^2){\rm{diag}}\{ \bar H_k^{m,i}P_{k|k - 1}^m{(\bar H_k^{m,i})^{\rm{T}}}\} _{i = 1}^{{N_k}} \end{array}$ (70)
$\begin{array}{l} C_{k,2}^m = Cov\{ ({\Theta _k} - {\Xi _k})(\bar H_k^m\mathit{\boldsymbol{\widehat x}}_{k|k - 1}^m + {\mathit{\boldsymbol{\overline u}} _k} - \mathit{\boldsymbol{\overline y}} _k^ - )\} = \\ \;\quad \quad ({\Xi _k} - \Xi _k^2){\rm{diag}}\{ (\bar H_k^{m,i}\mathit{\boldsymbol{\hat x}}_{k|k - 1}^m + \mathit{\boldsymbol{\overline u}} _k^i - \mathit{\boldsymbol{\overline y}} _k^{i - }) \times \\ \;\quad \quad \;{(\bar H_k^{m,i}\mathit{\boldsymbol{\hat x}}_{k|k - 1}^m + \mathit{\boldsymbol{\overline u}} _k^i - \mathit{\boldsymbol{\overline y}} _k^{i - })^{\rm{T}}}\} _{i = 1}^{{N_k}} \end{array}$ (71)
$C_{k,3}^m = Cov\{ {\Theta _k}{\mathit{\boldsymbol{\overline \eta }} _k}\} $ (72)
$C_{k,4}^m = Cov\{ {\Theta _k}\bar H_k^m\mathit{\boldsymbol{\widetilde x}}_{k|k - 1}^m,{\Theta _k}{\mathit{\boldsymbol{\overline \eta }} _k}\} $ (73)

进一步, 由引理2可知

$\begin{array}{l} C_{k,3}^m = \;Cov\{ {\Theta _k}{\mathit{\boldsymbol{\overline \eta }} _k}\} = \\ \;\quad \quad {\Xi _k}{{\bar R}_k}{\Xi _k} + ({\Xi _k} - \Xi _k^2) \times \\ \quad \quad {\rm{diag}}\{ R_k^i + \Pi _k^i\Omega ({t_k},t_k^i){(\Pi _k^i)^{\rm{T}}}\} _{k = i}^{{N_k}} \end{array}$ (74)

由式(64) 可知

$\begin{array}{l} C_{k,4}^m = Cov\{ {\Theta _k}\bar H_k^m\mathit{\boldsymbol{\widetilde x}}_{k|k - 1}^m,{\Theta _k}{\mathit{\boldsymbol{\overline \eta }} _k}\} = \\ \quad \quad Cov\{ {\Theta _k}\bar H_k^m{\mathit{\boldsymbol{\overline w}} _{k - 1}},{\Theta _k}{\mathit{\boldsymbol{\overline \eta }} _k}\} = \\ \quad \quad {\Xi _k}\bar H_k^m{{\bar \Psi }_k}{\Xi _k} + (\Xi _k^2 - {\Xi _k}) \times \\ \quad \quad {\rm{diag}}\{ \bar H_k^{m,i}{[\Pi _k^i\Omega ({t_k},t_k^i),{{\bf{0}}_{d_z^i \times 1}}]^{\rm{T}}}\} _{k = i}^{{N_k}} \end{array}$ (75)
3.3 故障概率更新

类似于标准IMM算法, 在$t_k$时刻, 第m种故障模型有效的模型概率为

$\begin{array}{l} \mu _k^m = Pr\left\{ {\rlap{--} m({t_k}) = m|{\mathit{\boldsymbol{Y}}^k}} \right\} = \;\\ \quad \quad \frac{{Pr\left\{ {\rlap{--} m({t_k}) = m|{\mathit{\boldsymbol{Y}}^{k - 1}}} \right\}p\left( {{{\mathit{\boldsymbol{\overline y}} }_k}|\rlap{--} m({t_k}) = m,{\mathit{\boldsymbol{Y}}^{k - 1}}} \right)}}{{\sum\limits_{i = 0}^{{d_f}} P r\left\{ {\rlap{--} m({t_k}) = i|{\mathit{\boldsymbol{Y}}^{k - 1}}} \right\}p\left( {{{\mathit{\boldsymbol{\overline y}} }_k}|\rlap{--} m({t_k}) = i,{\mathit{\boldsymbol{Y}}^{k - 1}}} \right)}} \end{array}$ (76)

其中,

$Pr\left\{ {\rlap{--} m({t_k}) = m|{\mathit{\boldsymbol{Y}}^{k - 1}}} \right\} = \sum\limits_{i = 0}^{{d_f}} {\mu _{k - 1}^i} {\lambda ^{i,m}}$ (77)
$\begin{array}{l} p\left( {{{\mathit{\boldsymbol{\overline y}} }_k}|\rlap{--} m({t_k}) = m,{\mathit{\boldsymbol{Y}}^{k - 1}}} \right) = \\ \quad \quad {(2\pi )^{ - \frac{1}{2}\sum\limits_{i = 1}^{{N_k}} {d_z^i} }}{\left| {C_k^m} \right|^{ - \frac{1}{2}}} \times \\ \quad \quad {\rm{exp}}\left\{ { - \frac{1}{2}{{(\mathit{\boldsymbol{\widetilde y}}_{k|k - 1}^m)}^{\rm{T}}}{{(C_k^m)}^{ - 1}}\mathit{\boldsymbol{\widetilde y}}_{k|k - 1}^m} \right\} \end{array}$ (78)
3.4 输出融合

由全概率公式可知

${\mathit{\boldsymbol{\widehat x}}_k} = E\left\{ {{{\mathit{\boldsymbol{\overline x}} }_k}|{\mathit{\boldsymbol{Y}}^k}} \right\} = \sum\limits_{m = 0}^{{d_f}} {\mathit{\boldsymbol{\widehat x}}_k^m} \mu _k^m$ (79)
$\begin{array}{l} {P_k} = E\left\{ {({{\mathit{\boldsymbol{\widehat x}}}_k} - {{\mathit{\boldsymbol{\overline x}} }_k}){{({{\mathit{\boldsymbol{\widehat x}}}_k} - {{\mathit{\boldsymbol{\overline x}} }_k})}^{\rm{T}}}} \right\} = \\ \quad \quad \sum\limits_{m = 0}^{{d_f}} {\left[ {P_k^m + \left( {{{\mathit{\boldsymbol{\widehat x}}}_k} - \mathit{\boldsymbol{\widehat x}}_k^m} \right){{\left( {{{\mathit{\boldsymbol{\widehat x}}}_k} - \mathit{\boldsymbol{\widehat x}}_k^m} \right)}^{\rm{T}}}} \right]} \mu _k^m \end{array}$ (80)
4 异步多传感器网络化系统故障融合诊断

故障模型概率反映了当前模型是有效模型的概率, 也即系统发生相应故障的概率, 因此可基于故障模型概率实现系统的故障诊断.具体步骤如下:

步骤1. 对每种故障情况, 即$m=0, 1, \cdots, d_f$分别建立如式(19)所示的系统模型, 并利用建立的$d_f+1$个模型构造模型集$\mathbb{M}$, 其中$m=0$对应系统无故障发生时的正常情况.

步骤2. 在当前融合时刻k, 对构建的模型集$\mathbb{M}$利用定理1中提出的网络化异步IMM融合滤波算法进行多模型滤波, 得到系统状态和故障幅值的联合估计${\hat {\pmb{x}}_k}$以及每个模型的模型概率$\mu_k^m$, $m=0, 1$, $\cdots$, $d_f$.

步骤3. 如果有

$\mu _k^i = \max \{ \mu _k^m\} _{m = 1}^{{d_f}} > {\mu _T}$ (81)

则可以判定有故障发生, 且发生的是第i种故障.相应模型条件下的联合估计${\hat {\pmb{x}}_k}^i$同时给出了故障幅值的融合估计, 从而实现了故障的联合检测、定位和估计.其中, ${\mu _T}$是一个预先给定的概率门限值, 并且$ {\mu _T}$ $\in$ $(0, 1]$.实际中为了减少故障虚报率, 可设定时间窗内L步的模型概率连续超过阈值$\mu _T$才判定故障发生.

步骤4. 令$k = k + 1$, 并转至步骤2.

5 仿真实验

考虑式(1) 描述的动态系统, 相关参数设置如下

$\begin{array}{l} A(t) = \left[ {\begin{array}{*{20}{c}} 0&1&0&0\\ 0&0&0&{ - \omega }\\ 0&0&0&1\\ 0&\omega &0&0 \end{array}} \right]\\ Q(t) = \left[ {\begin{array}{*{20}{c}} {0.01}&0\\ 0&{0.01} \end{array}} \right]\\ B(t) = G(t) = F(t) = {\left[ {\begin{array}{*{20}{c}} 0&{ - 1}&0&0\\ 0&0&0&1 \end{array}} \right]^{\rm{T}}} \end{array}$

其中, 角速度$\omega=0.2$, 控制输入$u(t)=[0.5~~~0.5]^{\rm T}$.正常情况下, $\pmb{f}(t)= 0$, 目标在二维平面中以固定角速度做匀速圆周运动.

运动目标由三个异步传感器进行观测, 其初始采样时间为$t_{(1) }=0.2$ s, $t_{(2) }=0.3$ s和$t_{(3) }=0.4$ s, 对应的采样周期为$T_{(1) }=0.6$ s, $T_{(2) }=0.6$ s, $T_{(3) }$在集合{0.8 s, 0.9 s, 1 s}中随机取值, 可见三个传感器是异步工作的, 且传感器3是非均匀采样的.传感器的测量矩阵和测量噪声的协方差矩阵分别为

$\begin{array}{l} {H_{(1)}} = {H_{(2)}} = {H_{(3)}} = \left[ {\begin{array}{*{20}{c}} 1&0&0&0\\ 0&0&1&0 \end{array}} \right]\\ {R_{(1)}} = {R_{(2)}} = {R_{(3)}} = \left[ {\begin{array}{*{20}{c}} {25}&0\\ 0&{25} \end{array}} \right] \end{array}$

传感器测量通过无线通信网络传输至融合中心, 但通信网络存在随机丢包现象, 丢包概率分别为$0.1$, $0.2$$0.2$, 即$\beta_{(1) }=0.9$, $\beta_{(2) }=0.8$, $\beta_{(3) }=0.8$.融合中心的融合周期$T=1$ s.

本例仅考虑单一执行器发生常值故障的情况, 系统中存在两个执行器, 假定第2个执行器在$\tau_f=$ $16.5$ s时刻发生故障, 其故障幅值$f=2$.采用本文提出的基于网络化异步IMM融合滤波的方法对上述系统进行故障诊断.分别对两个执行器发生故障的情况和系统正常情况下的故障情况建立模型, 构成包含三个模型的多模型集进行多模型滤波.其中模型的一步转移概率取值为

$\Lambda = \left[ {\begin{array}{*{20}{c}} {0.9}&{0.05}&{0.05}\\ {0.05}&{0.9}&{0.05}\\ {0.05}&{0.05}&{0.9} \end{array}} \right]$

图 2给出了三个模型的模型概率演化曲线, 可以看到在发生故障之前, 正常模型的概率明显占优.随着故障的发生, 正常模型的概率急剧下降, 执行器2故障模型的概率急剧上升, 并迅速超过了故障阈值$0.5$.因此, 可以判定系统发生了故障, 且是执行器2发生了故障, 同时可以得到故障幅值的估计, 如图 3所示. 图 3同时给出了利用单个传感器数据和所提出的网络化异步融合算法融合三个传感器数据后所得到的故障幅值估计曲线.为了更直观地反映融合算法的效果, 图 4图 5分别给出了500次蒙特卡洛仿真得到的系统状态估计和故障幅值估计的均方根误差(Root mean square error, RMSE)[2]曲线.由图 4图 5可知, 与基于单个传感器数据得到的结果相比较, 融合算法得到的状态和故障幅值估计具有更小的RMSE误差, 并且对故障具有最好的鲁棒性.上述结果验证了本文提出的基于异步IMM融合滤波的网络化系统故障诊断方法的可行性和有效性.此外, 从图 4图 5中还可以看出, 传感器1比传感器2具有更好的估计效果.这是因为传感器1和传感器2虽然具有相同的采样速率, 但是前者有更小的丢包率.传感器3的丢包率与传感器2相同但采样率更低, 因此效果最差.这反映了融合估计效果与丢包率和采样速率之间的关系, 也与人们的直观结论相符.

图 2 模型后验概率曲线 Figure 2 The posterior probability curves of models
图 3 故障幅值的估计曲线 Figure 3 Estimation curves of fault amplitude
图 4 状态估计的均方根误差曲线 Figure 4 RMSE curves of state estimation
图 5 故障幅值估计的均方根误差曲线 Figure 5 RMSE curves of fault amplitude
6 结论

本文针对具有多个异步传感器的网络化系统, 提出了基于异步IMM融合滤波的网络化系统故障诊断方法.提出的方法考虑了网络传输中存在的随机丢包现象, 实现了故障的联合检测、分离和估计.由于提出的故障诊断方法是依据IMM融合滤波得到的故障模型概率进行故障诊断的, 与其他基于模型的故障诊断方法相比, 其检测阈值的含义更加明确, 取值更容易确定.此外, 本文提出的方法避免了传统基于IMM滤波的故障诊断方法模型集设计中作为模型参数的故障大小取值难以确定的问题, 并且实现了系统状态和故障幅值基于异步传感器测量的联合融合估计.本文提出的方法适用于具有任意采样速率和任意初始采样时间的异步多传感器系统, 并且不要求传感器是均匀采样的.但是, 本文仅针对线性系统考虑了网络中的丢包现象, 对于非线性系统或网络同时存在丢包和延迟等现象情况下的IMM故障诊断方法还需进一步的研究[20].

参考文献
1
Yan Rong-Yi, He Xiao, Zhou Dong-Hua. Robust diagnosis of intermittent faults for linear stochastic systems subject to time-varying perturbations. Acta Automatica Sinica, 2016, 42(7): 1004-1013.
( 鄢镕易, 何潇, 周东华. 一类存在参数摄动的线性随机系统的鲁棒间歇故障诊断方法. 自动化学报, 2016, 42(7): 1004-1013.)
2
Xu Xiao-Bin, Zhang Zhen, Li Shi-Bao, Wen Cheng-Lin. Fault diagnosis based on fusion and updating of diagnosis evidence. Acta Automatica Sinica, 2016, 42(1): 107-121.
( 徐晓滨, 张镇, 李世宝, 文成林. 基于诊断证据静态融合与动态更新的故障诊断方法. 自动化学报, 2016, 42(1): 107-121.)
3
Li X R, Jilkov V P. Survey of maneuvering target tracking. Part V. Multiple-model methods. IEEE Transactions on Aerospace and Electronic Systems, 2005, 41(4): 1255-1321. DOI:10.1109/TAES.2005.1561886
4
Zhang Y M, Li X R. Detection and diagnosis of sensor and actuator failures using IMM estimator. IEEE Transactions on Aerospace and Electronic Systems, 1998, 34(4): 1293-1313. DOI:10.1109/7.722715
5
Ru J F, Li X R. Interacting multiple model algorithm with maximum likelihood estimation for FDI. In:Proceedings of the 2003 IEEE International Symposium on Intelligent Control. Houston, TX, USA:IEEE, 2003. 661-666
6
Ru J F, Li X R. Variable-structure multiple-model approach to fault detection, identification, and estimation. IEEE Transactions on Control Systems Technology, 2008, 16(5): 1029-1038. DOI:10.1109/TCST.2007.916318
7
Kim S, Choi J, Kim Y. Fault detection and diagnosis of aircraft actuators using fuzzy-tuning IMM filter. IEEE Transactions on Aerospace and Electronic Systems, 2008, 44(3): 940-952. DOI:10.1109/TAES.2008.4655354
8
Zhao S Y, Huang B, Liu F. Fault detection and diagnosis of multiple-model systems with Mismodeled transition probabilities. IEEE Transactions on Industrial Electronics, 2015, 62(8): 5063-5071. DOI:10.1109/TIE.2015.2402112
9
Zhang Y M, Jiang J. Integrated active fault-tolerant control using IMM approach. IEEE Transactions on Aerospace and Electronic Systems, 2001, 37(4): 1221-1235. DOI:10.1109/7.976961
10
Tudoroiu N, Khorasani K. Fault detection and diagnosis for satellite's attitude control system (ACS) using an interactive multiple model (IMM) approach. In:Proceedings of the 2005 IEEE Conference on Control Applications. Toronto, Canada:IEEE, 2005. 1287-1292
11
Judalet V, Glaser S, Gruyer D, Mammar S. IMM-based sensor fault detection and identification for a drive-by-wire vehicle. In:Proceedings of the 9th IFAC Symposium on Fault Detection, Supervision and Safety for Technical Processes. Paris, France:IFAC, 2015. 1158-1164
12
Han Chong-Zhao, Zhu Hong-Yan, Duan Zhan-Sheng. Multi-Source Information Fusion. 2nd edition. Beijing: Tsinghua University Press, 2010, 455-496.
( 韩崇昭, 朱洪艳, 段战胜. 多源信息融合. 第2版. 北京: 清华大学出版社, 2010, 455-496.)
13
Pan Quan. Multi-Source Information Fusion Theory and Its Applications. Beijing: Tsinghua University Press, 2013, 1-21.
( 潘泉. 多源信息融合理论及应用. 北京: 清华大学出版社, 2013, 1-21.)
14
Liu W Y, Wei J, Liang M C, Cao Y, Hwang I. Multi-sensor fusion and fault detection using hybrid estimation for air traffic surveillance. IEEE Transactions on Aerospace and Electronic Systems, 2013, 49(4): 2323-2339. DOI:10.1109/TAES.2013.6621819
15
Hu Y Y, Zhou D H. Actuator fault diagnosis for dynamic systems with multiple asynchronous sensors using the IMM approach. In:Proceedings of the 3rd International Conference on Measuring Technology and Mechatronics Automation. Shanghai, China:IEEE, 2011. 318-321
16
He X, Wang Z D, Liu Y, Zhou D H. Least-squares fault detection and diagnosis for networked sensing systems using a direct state estimation approach. IEEE Transactions on Industrial Informatics, 2013, 9(3): 1670-1679. DOI:10.1109/TII.2013.2251891
17
Dong H L, Wang Z D, Ding S X, Gao H J. On H estimation of randomly occurring faults for a class of nonlinear time-varying systems with fading channels. IEEE Transactions on Automatic Control, 2016, 61(2): 479-484. DOI:10.1109/TAC.2015.2437526
18
He X, Wang Z D, Wang X F, Zhou D H. Networked strong tracking filtering with multiple packet dropouts:algorithms and applications. IEEE Transactions on Industrial Electronics, 2014, 61(3): 1454-1463. DOI:10.1109/TIE.2013.2261038
19
Li X R. Optimal linear estimation fusion-part Ⅶ:dynamic systems. In:Proceedings of the 6th International Conference of Information Fusion. Cairns, Australia:IEEE, 2003. 455-462
20
Peng Kai-Xiang, Zhang Kai, Li Gang. Online contribution rate based fault diagnosis for nonlinear industrial processes. Acta Automatica Sinica, 2014, 40(3): 423-430.
( 彭开香, 张凯, 李钢. 基于贡献率法的非线性工业过程在线故障诊断. 自动化学报, 2014, 40(3): 423-430.)