A Novel Method of Quality Abnormality Detection and Fault Quantitative Assessment for Industrial Processes
-
摘要: 质量异常检测(Quality abnormality detection, QAD)与故障量化评估(Fault quantitative assessment, FQA)作为工业过程监控的关键环节, 是故障诊断领域的研究热点. 本文提出了一种新的工业过程质量异常检测与故障量化评估方法. 首先, 采用弹性网络(Elastic net, EN)算法构建了质量相关的变量候选集, 借助典型相关分析(Canonical correlation analysis, CCA)构建了质量相关的特征向量, 并引入支持向量数据描述(Support vector data description, SVDD)实现质量异常检测. 其次, 从优化近邻点距离的角度提出了增强局部线性嵌入(Enhanced local linear embedding, ELLE)算法, 并提出了基于CCA-ELLE的质量异常故障量化评估方法. 最后, 通过田纳西−伊斯曼(Tennessee-Eastman, TE)过程进行仿真验证, 并与传统的方法进行对比分析, 实验结果验证了所提方法的优越性和有效性.Abstract: As the key link of industrial process monitoring, quality abnormality detection (QAD) and fault quantitative assessment (FQA) are the research hotspot in the field of fault diagnosis. In this paper, a novel framework of industrial process QAD and FQA is proposed. Firstly, the candidate set of quality related variables is constructed by elastic net (EN) algorithm. The quality related eigenvectors are constructed by canonical correlation analysis (CCA), and the support vector data description (SVDD) is introduced to detect quality anomaly. Secondly, an enhanced local linear embedding (ELLE) algorithm is proposed from the perspective of optimizing the distance between adjacent points, and a quantitative assessment method based on CCA- ELLE is proposed. Finally, the proposed framework is simulated and compared with some traditional methods by Tennessee-Eastman (TE) process simulation platform. The experimental results verify the superiority and effectiveness of the proposed method.
-
近年来, 多传感器信息融合技术得到广泛关注.经典Kalman滤波需要已知模型参数和噪声统计特性[1].而在实际应用中, 系统可能会出现模型参数或噪声统计特性未知情况.文献[2-4]研究了带未知噪声方差系统的自校正估计问题.文献[5]和[6]研究了模型参数和噪声方差未知系统的自校正融合估计问题.文献[7]针对带未知模型参数和噪声统计特性且噪声相关的多传感器系统, 分别应用RELS算法、Gevers-Wouters算法对未知模型参数和噪声统计特性进行辨识, 并提出了自校正融合估值器.文献[5-7]中在对未知模型参数进行融合处理时, 采用加权平均方法获得最终融合辨识器.该方法没有考虑不同传感器的局部参数辨识的差异, 不能保证融合后的参数估计精度都优于每个局部参数估计.上述文献所提出的辨识和估计算法都是基于完整的传感器观测数据, 而没有考虑数据的不完整现象.
在实际的网络化系统或传感器网络中, 由于传感器老化或故障、以及传感器的观测数据在通信传输过程中由于带宽有限, 可能会出现数据丢失、衰减、延迟等问题, 使得估值器所收到的传感器数据具有不完整现象[8].文献[9]研究了带多丢包系统的最优估计问题.文献[10]将文献[9]的单传感器系统推广到了多传感器系统, 对带有不同丢包率的多传感器系统, 提出了集中式和分布式融合估值器.文献[11-14]考虑数据在传输过程中存在的丢失和延迟现象, 设计了相应的最优估值器.其中文献[13]在文献[12]的基础上研究了带随机乘性噪声、多丢包和滞后系统的最优估计问题.文献[14]考虑了过程噪声和观测噪声具有一步自相关和互相关的情况.文献[15]和[16]对带有数据包丢失的多传感器系统, 应用协方差信息方法设计了分布式融合估值器.在上述文献[9-17]中数据丢失现象均由一组满足伯努利分布的随机变量所描述, 这种数据丢失可以看作是观测数据发生衰减现象的一种特殊情况.文献[18]和[19]考虑了带随机参数矩阵、相关噪声和衰减观测系统的最优和次优估计问题.文献[20]中研究了带衰减观测系统的Kalman滤波估计问题, 同时分析了误差协方差的有界性和稳态特性.文献[21-22]对带衰减观测的随机不确定系统, 分别设计了多传感器分布式和序贯融合估值器.上述文献都是在假设观测丢失率或衰减率已知的情况下得到的最优估计结果.而没有考虑观测丢失率或衰减率未知的自校正估计问题.目前, 有关同时考虑带有未知模型参数和未知衰减观测率多传感器系统的自校正融合估计问题的报导甚少.
基于以上文献分析, 本文将对带有未知模型参数和衰减观测率的多传感器随机系统, 应用相关函数和递推增广最小二乘算法, 分别在线辨识衰减观测的数学期望、方差和模型参数.应用线性无偏最小方差估计准则, 分别提出了分布式模型参数融合辨识器和自校正融合状态滤波器, 并分析算法的收敛性.
1. 问题阐述
考虑带衰减观测的多传感器随机系统
$$ \begin{equation} {\boldsymbol x(t+1)} = {\Phi}{\boldsymbol x(t)}+{\Gamma}{\boldsymbol w}(t) \end{equation} $$ (1) $$ \begin{equation} {{y}_{i}}(t) = {\mu_{i}(t)}{ h_{i}}{\boldsymbol x}(t)+{v_{i}(t)}, \ i = 1, 2, \cdots, L \end{equation} $$ (2) 其中状态$ \boldsymbol x(t)\in {\bf R}^n $, 观测$ {{y}_{i}}(t)\in {{\bf R}} $, 过程噪声$ \boldsymbol w(t)\in {\bf R}^r $, 观测噪声$ {{v}_{i}}(t)\in {{\bf R}} $, 下标$ i $表示第$ i $个传感器, $ L $表示传感器的个数. $ \left\{ {{\mu _i}(t)} \right\} $为一组在${[0, 1]}$区间取值的用来描述第$ i $个传感器衰减观测现象的标量随机变量.其中$ {\rm E}{[\mu_{i}(t)] = \alpha_{i}} $, $ {\rm Cov}[\mu_{i}(t)] = \sigma_{i}^2 $, $ {\rm E} $为数学期望符号, $ {\rm Cov} $为协方差符号. $ \left\{ {{\mu _i}(t)} \right\} $与其他随机变量不相关. $ \Phi $, $ \Gamma $, $ { h_{i}} $是适当维数矩阵.
假设1. $ {\boldsymbol w}(t) $和$ {v_{i}(t)} $为零均值、方差阵分别为$ Q_{{\boldsymbol w}} $和$ Q_{v_{i}} $的不相关白噪声.
假设2. 初值$ {\boldsymbol x(0)} $不相关于$ {\boldsymbol w}(t) $和$ {v_{i}(t)} $, 且$ {\rm{E}}\left\{ {\boldsymbol x(0)} \right\} = {\boldsymbol u_0}, {\rm{E}}\left\{ {\left[{\boldsymbol x(0) - {\boldsymbol u_0}} \right]{{\left[{\boldsymbol x(0) - {\boldsymbol u_0}} \right]}^{\rm{T}}}} \right\} = {P_0} $.其中$ {\rm T} $为转置号.
假设3. $ \Phi $为稳定矩阵, ($ \Phi, { h_{i}} $)是完全可观对, ($ \Phi, \Gamma $)是完全可控对.
假设4. $ \Phi $中部分参数未知, $ \left\{ {{\mu _i}(t)} \right\} $的数学期望$ {\alpha_i} $和方差$ \sigma_{i}^2 $未知.
问题是基于观测$ ({y}_{i}(1), \cdots, {y}_{i}(t)) $, $ i = 1, 2, $ $ \cdots, L $, 辨识$ \Phi $中未知参数、$ \left\{ {{\mu _i}(t)} \right\} $的数学期望$ {\alpha_i} $和方差$ \sigma_{i}^2 $, 并求$ \Phi $中未知参数的融合辨识器$ {{\rm{\hat \Phi }}_o}(t) $和状态$ {\boldsymbol x(t)} $的自校正融合滤波器$ {\hat{\boldsymbol {x}}_s(t|t)} $.
注1. 为了阅读方便, 这里对后文中经常遇到一些术语给予解释.局部滤波器, 即基于单个传感器的观测数据所获得的滤波器; 最优滤波器, 即系统的模型参数以及衰减观测期望和方差已知情况下, 获得的线性最小方差意义下的滤波器; 自校正滤波器, 即系统模型含有未知参数, 以及衰减观测期望和方差未知情况下, 通过辨识这些未知参数, 然后代入最优滤波算法中获得的滤波器.分布式融合滤波器, 即基于各个传感器的局部滤波器, 应用线性无偏最小方差矩阵加权融合估计算法[23]获得的融合滤波器.
2. 分布式最优融合滤波
当系统模型参数、$ \left\{ {{\mu _i}(t)} \right\} $的数学期望$ {\alpha_i} $和方差$ \sigma_{i}^2 $已知时, 应用线性无偏最小方差意义下的矩阵加权融合估计算法[23]可获得分布式最优融合滤波器.下面给出实现过程.
由式(2)可得
$$ \begin{equation} {{y}_{i}}(t) = {\alpha_i}{ h_{i}}{\boldsymbol x(t)}+{V_i(t)} \end{equation} $$ (3) 其中
$$ \begin{equation} {V_i(t)} = (\mu_{i}(t)-{\alpha_i}){h_{i}}{\boldsymbol x(t)}+v_i(t) \end{equation} $$ (4) 可计算其方差阵为
$$ \begin{equation} Q_{{V_i(t)}} = {\rm E}[V_i^2(t)] = {{{\sigma}}_{i}^2}{ h_{i}}X(t){ h_{i}^{\rm T}}+Q_{v_i} \end{equation} $$ (5) 状态二阶矩$ X(t) = {\rm E}[{\boldsymbol x}(t){\boldsymbol x}^{\rm T}(t)] $可递推计算如下:
$$ \begin{equation} X(t + 1) = {\rm{\Phi }}X(t){{\rm{\Phi }}^{\rm{T}}} + {\rm{\Gamma }}{Q_{\boldsymbol w}}{{\rm{\Gamma }}^{\rm{T}}} \end{equation} $$ (6) 初值为$ X(0) = {\boldsymbol \mu_0}{\boldsymbol \mu}_0^{\rm{T}}+ {P_0} $.由假设3可知$ X(t) $有界.
下面引理1给出了最优局部滤波器算法; 引理2给出了互协方差计算公式; 引理3给出了分布式最优加权融合滤波算法.
引理1[24]. 在假设1$ \sim $3下, 随机系统(1)和(3)基于每个传感器的观测有最优局部滤波器
$$ \begin{equation} {\hat {\boldsymbol x}_i}(t + 1|t + 1) = {\Psi _{fi}}(t + 1){\hat {\boldsymbol x}_i}(t|t)+{{K}_i}(t + 1){{ y}_i}(t + 1) \end{equation} $$ (7) $$ \begin{equation} {{ K}_i}(t + 1) = {\Sigma _i}(t + 1|t) F_i^{\rm{T}}Q_{{{{C}}_i}}^{ - 1}(t + 1) \end{equation} $$ (8) $$ \begin{equation} {\Sigma _i}(t + 1|t) = \Phi {P_i}(t|t){\Phi ^{\rm{T}}} + \Gamma {Q_{\boldsymbol w}}{\Gamma ^{\rm{T}}} \end{equation} $$ (9) $$ \begin{equation} {Q_{{{{C}}_i}}}(t + 1) = { F_i}{\Sigma _i}(t + 1|t)F_i^{\rm{T}} + {Q_{{V_i}}}(t + 1) \end{equation} $$ (10) $$ \begin{equation} {P_i}(t + 1|t + 1) = [{I_n} - {{ K}_i}(t + 1){ F_i}]{\Sigma _i}(t + 1|t) \end{equation} $$ (11) 其中$ { F_i} = {\alpha _i}{ h_i}, {\Psi _{fi}}(t + 1) = {\rm{[}}{I_n}-{{ K}_i}(t + 1){F_i}]\Phi $. $ {\hat {\boldsymbol x}_i}(t|t) $为第$ i $个传感器的局部滤波器, $ {K_i}(t + 1) $为相应的滤波增益, $ {P_i}(t|t) $为滤波误差方差阵.初值为$ {\hat {\boldsymbol x}_i}(0|0) = {\boldsymbol u_0}, {P_i}(0|0) = {P_0} $.
引理2[24]. 任意两个局部滤波误差之间的互协方差阵$ {P_{ij}}(t|t) = {\rm E}[{\tilde {\boldsymbol x}_i}(t|t)\tilde {\boldsymbol x}_j^{\rm T}(t|t)] $ (其中滤波误差$ {\tilde {\boldsymbol x}_i}(t|t) = {\boldsymbol x}(t)-{\hat {\boldsymbol x}_i}(t|t) $)可递推计算如下:
$$ \begin{align} {P_{ij}}&(t + 1|t + 1) = [{I_n} - { K_i}(t + 1){ F_i}]\times \\ &[\Phi{P_{ij}}(t|t){\Phi ^{\rm{T}}} + \Gamma {Q_{\boldsymbol w}}{\Gamma ^{\rm{T}}}]{[{I_n} - { K_j}(t + 1){ F_j}]^{\rm{T}}} \end{align} $$ (12) 初值为$ {P_{ij}}(0|0) = {P_0} $.
引理3[23]. 基于引理1的各局部滤波器和引理2的任意两个局部滤波误差之间的互协方差阵, 分布式最优矩阵加权融合滤波器可计算如下:
$$ \begin{equation} {\hat {\boldsymbol x}_o}(t|t) = \sum\limits_{i = 1}^L {W_i^{}(t){{\hat {\boldsymbol x}}_i}(t|t)} \end{equation} $$ (13) 加权矩阵计算为
$$ \begin{equation} [W_1^{}(t), \cdots , W_L^{}(t)]{\kern 1pt} {\kern 1pt} {\kern 1pt} = {(e_{}^{\rm{T}}P_{}^{ - 1}(t|t)e)^{ - 1}}e_{}^{\rm{T}}P_{}^{ - 1}(t|t) \end{equation} $$ (14) 其中$ e = [I_{n}, \cdots, I_{n}]^{\mathrm{T}}, P(t|t) = {\left[{{P_{ij}}(t|t)} \right]_{nL \times nL}} $是以$ {P_{ij}}(t|t) $为第$ (i, j) $元素的分块矩阵.融合滤波器的估计误差方差阵计算为
$$ \begin{equation} {P_o}(t|t) = {(e_{}^{\rm T}P_{}^{ - 1}(t|t)e)^{ - 1}} \end{equation} $$ (15) 且有$ {P_o}(t|t) \le {P_i}(t|t), i = 1, \cdots, L. $
3. 未知模型参数融合辨识器
上一节我们针对系统模型精确已知时给出了分布式最优融合估计算法.而在实际应用中, 系统模型可能含有未知参数.当$ \Phi $中含有未知参数时, 本节采用RELS算法辨识未知模型参数, 并对辨识得到的$ L $组参数估值进行加权融合, 获得模型参数的分布式融合辨识器.下面给出具体实现过程.
由式(1)可得
$$ \begin{equation} {\boldsymbol x}(t) = {({I_n} - {q^{ - 1}}\Phi )^{ - 1}}{q^{ - 1}}\Gamma {\boldsymbol w}(t) \end{equation} $$ (16) 式中$ q^{-1} $为单位滞后算子, 即$ {q^{ - 1}}{\boldsymbol x}(t) = {\boldsymbol x}(t - 1) $.将式(16)代入式(3)得
$$ \begin{equation} {y_i}(t) = {\alpha _i}{ h_i}{({I_n} - {q^{ - 1}}\Phi )^{ - 1}}{q^{ - 1}}\Gamma {\boldsymbol w}(t) + {V_i}(t) \end{equation} $$ (17) 将式(17)进一步化简得
$$ \begin{equation} A({q^{ - 1}}){ y_i}(t) = {\alpha _i}{B_i}({q^{-1}}){\boldsymbol w}(t)+ A({q^{ - 1}}){V_i}(t) \end{equation} $$ (18) 其中$ A({q^{-1}}) = \det ({I_n}-{q^{ - 1}}\Phi), {B_i}({q^{ - 1}}) = { h_i}{\rm{adj}}({I_n}-{q^{ - 1}}\Phi){q^{ - 1}}\Gamma $, 式中det和adj分别表示矩阵行列式和伴随矩阵. $ A(q^{-1}) $和$ B_i(q^{-1}) $具有如下多项式形式:
$$ \begin{equation} A({q^{ - 1}}) = 1 + {a_1}{q^{ - 1}} + \cdots + {a_{{n_a}}}{q^{ - {n_a}}} \end{equation} $$ (19) $$ \begin{equation} {B_i}({q^{ - 1}}) = {B_{i1}}{q^{ - 1}} + \cdots + {B_{i{n_{{B_i}}}}}{q^{ - {n_{{B_i}}}}} \end{equation} $$ (20) 其中$ {a_k}, k = {\rm{1}}, \cdots, {n_a} $和$ {B_{ik}}, k = {\rm{1}}, \cdots, {n_{{B_i}}} $是多项式系数. $ n_A, n_{B_i} $分别为$ A(q^{-1}) $和$ B_i(q^{-1}) $的阶次.式(18)等号右侧两个滑动平均过程可以等价为一个稳定的滑动平均过程$ {D_i}({q^{ - 1}}){\varepsilon _i}(t) $[24], 即
$$ \begin{equation} {D_i}({q^{ - 1}}){\varepsilon _i}(t) = {\alpha _i}{B_i}({q^{ - 1}}){\boldsymbol w}(t) + A({q^{ - 1}}){V_i}(t) \end{equation} $$ (21) 其中$ {\varepsilon _i}(t) $是零均值且带有未知噪声方差$ \sigma _{{\varepsilon _i}}^2 $的白噪声, $ {D_i}({q^{ - 1}}) $具有如下多项式形式:
$$ \begin{equation} {D_i}({q^{ - 1}}) = 1 + {d_{i1}}{q^{ - 1}} + \cdots + {d_{i{n_{{D_i}}}}}{q^{ - {n_{{D_i}}}}} \end{equation} $$ (22) 其中$ {d_{ik}}, k = 1, \cdots, {n_{{D_i}}} $是多项式$ {D_i}({q^{ - 1}}) $的系数, $ {n_{{D_i}}} $是$ {D_i}({q^{ - 1}}) $的阶次.
将式(18)重写为
$$ \begin{equation} A({q^{ - 1}}){y_i}(t) = {D_i}({q^{ - 1}}){\varepsilon _i}(t) \end{equation} $$ (23) 令$ \varphi _i^{\rm{T}}(t) = [-{y_i}(t- 1), \cdots, -{y_i}(t- {n_a}), {\hat \varepsilon _i}(t - 1), \cdots, {\hat \varepsilon _i}(t -{n_{{D_i}}})], {\vartheta _i} = [{a_1}, \cdots, {a_{{n_a}}}, {d_{i1}}, \cdots, {d_{i{n_{{D_i}}}}}{{\rm{]}}^{\rm{T}}} $, 则式(23)可表示为
$$ \begin{equation} {y_i}(t) = \varphi _i^{\rm{T}}(t){\vartheta _i} + {\varepsilon _i}(t) \end{equation} $$ (24) 参数$ {a_k}, k = 1, \cdots, {n_a}; {d_{ik}}, k = 1, \cdots, {n_{{D_i}}} $未知.
基于每个单传感器的观测数据, 应用RELS算法[24]可得到局部参数估计为
$$ \begin{equation} {\hat \vartheta _i}(t + 1) = {\hat \vartheta _i}(t) + {M_i}(t + 1){\hat \varepsilon _i}(t{\rm{ + }}1) \end{equation} $$ (25) $$ \begin{equation} {\hat \varepsilon _i}(t + 1) = {y_i}(t + 1) - \varphi _i^{\rm{T}}(t + 1){\hat \vartheta _i}(t) \end{equation} $$ (26) $$ \begin{equation} {M_i}(t + 1) = \frac{{{Z_i}(t){\varphi _i}(t + 1)}}{{1 + \varphi _i^{\rm{T}}(t + 1){Z_i}(t){\varphi _i}(t + 1)}} \end{equation} $$ (27) $$ \begin{equation} {Z_i}(t + 1) = [{I_{{n_a} + {n_{{D_i}}}}} - {M_i}(t + 1)\varphi _i^{\rm{T}}(t + 1)]{Z_i}(t) \end{equation} $$ (28) 代初值$ {\hat \vartheta _i}{\rm{(}}0) = 0, \; {Z_i}(0) = {\beta _i}I, {\beta _i} $为充分大的正数, 且规定$ {\hat \varepsilon _i}(j) = 0, \; {y_i}(j) = 0\; (j \le 0) $.
由文献[24]可知, RELS算法参数估计是一致的, 即$ {\hat \vartheta _i}(t) \to {\vartheta _i}, t \to \infty, w.p.1 $.符号"$ w.p.1. $"表示"以概率1".
由式(25)可得基于单传感器的局部参数估计误差$ {\tilde \vartheta _i}(t) = {\vartheta _i} - {\hat \vartheta _i}(t) $满足如下方程:
$$ \begin{align} {\tilde \vartheta _i}(t + 1) = \, & [{I_{{n_A} + {n_{{D_i}}}}} - {M_i}(t + 1)\varphi _i^{\rm{T}}(t + 1)]{\tilde \vartheta _i}(t) -\\ & {M_i}(t + 1){\varepsilon _i}(t + 1) \end{align} $$ (29) 于是, 任意两个局部参数估计之间的估计误差协方差阵$ {P_{{\vartheta _{ij}}}}(t) = {\rm{E}}[{\tilde \vartheta _i}(t)\tilde \vartheta _j^{\rm{T}}(t)] $可计算如下:
$$ \begin{align} {P_{{\vartheta _{ij}}}}(t + 1) = \, &[{I_{{n_A} + {n_{{D_i}}}}} - {M_i}(t + 1)\varphi _i^{\rm{T}}(t + 1)]{P_{{\vartheta _{ij}}}}(t) \\ &{[{I_{{n_A} + {n_{{D_j}}}}} - {M_j}(t + 1)\varphi _j^{\rm{T}}(t + 1)]^{\rm{T}}} +\\ & {M_i}(t + 1)\hat \sigma _{{\varepsilon _{ij}}}^2(t + 1)M_j^{\rm{T}}(t + 1) \end{align} $$ (30) 当$ i = j $时, $ {P_{{\vartheta _{ii}}}}(t) $即为局部参数估计误差方差阵$ {P_{{\vartheta _i}}}(t) $. $ {\varepsilon _i}(t) $与$ {\varepsilon _j}(t) $之间的互协方差$ \sigma _{{\varepsilon _{ij}}}^2 $可近似计算如下:
$$ \begin{equation} \hat \sigma _{{\varepsilon _{ij}}}^2(t) = \frac{1}{t}\sum\limits_{k = 1}^t {\hat \varepsilon _i^{}(k)\hat \varepsilon _j^{}(k)} \end{equation} $$ (31) 它可递推地计算为
$$ \begin{equation} \hat \sigma _{{\varepsilon _{ij}}}^2(t) = \hat \sigma _{{\varepsilon _{ij}}}^2(t - 1) + \frac{1}{t}[\hat \varepsilon _i^{}(t)\hat \varepsilon _j^{}(t) - \hat \sigma _{{\varepsilon _{ij}}}^2(t - 1)] \end{equation} $$ (32) 初值为$ \hat \sigma _{{\varepsilon _{ij}}}^2(0) = {y_i}(0){y_j}(0), i, j = 1, \cdots, L $.
令$ {\vartheta _A} = {[{a_1}, \cdots, {a_{{n_a}}}]^{\rm{T}}} $, 则有$ {\vartheta _{A_i}} = [{I_{{n_a}}}, 0]{\vartheta _i} $.于是, 我们有参数$ {\vartheta _A} $基于传感器$ i $的局部估计和估计误差协方差阵如下:
$$ \begin{equation} {\hat \vartheta _{Ai}}(t) = [{I_{{n_a}}}, 0]{\hat \vartheta _i}(t) \end{equation} $$ (33) $$ \begin{equation} {P_{{\vartheta _{Aij}}}}(t) = [{I_{{n_a}}}, 0]{P_{{\vartheta _{ij}}}}(t){[{I_{{n_a}}}, 0]^{\rm{T}}} \end{equation} $$ (34) 当$ i = j $时, $ {P_{{\vartheta _{Aii}}}}(t) $即为局部参数$ {\vartheta _A} $的估计误差方差阵$ {P_{{\vartheta _{Ai}}}}(t) $.
由式(18)可知, 参数$ {\vartheta _A} = {[{a_1}, \cdots, {a_{{n_a}}}]^{\rm{T}}} $是$ \Phi $中未知参数的函数.假设$ \Phi $中未知模型参数组成的列向量为$ {\Lambda ^{[{\rm{\Phi }}]}} \in {{\boldsymbol{\rm R}}^{{n_\Phi }}}, {n_\Phi } \le {n_a} $, 且可由$ {\vartheta _A} $唯一确定.设${\Lambda ^{[{\rm{\Phi }}]}}$与${\vartheta _A}$之间满足如下关系:
$$ \begin{equation} {\Lambda ^{[{\rm{\Phi }}]}} = f({\vartheta _A}) \end{equation} $$ (35) 其中$ f({\vartheta _A}) $为关于$ {\vartheta _A} $的线性或非线性函数.
1) 如果$ f({\vartheta _A}) $是线性函数, 我们将式(35)重写为
$$ \begin{equation} {\Lambda ^{{\rm{[\Phi }}]}} = S{\vartheta _A} + \gamma \end{equation} $$ (36) 式中$ S, \gamma $为适当维数的系数阵.
那么, 基于传感器$ i $的数据获得的$ \Phi $中未知模型参数在$ t $时刻的局部估计为
$$ \begin{equation} \hat \Lambda _i^{[{\rm{\Phi }}]}(t) = S{\hat \vartheta _{Ai}}(t) + \gamma \end{equation} $$ (37) 定义局部估计误差方差$ {P_{\Lambda _i^{[{\rm{\Phi }}]}}}(t) = {\rm{E}}[\tilde \Lambda _i^{[{\rm{\Phi }}]}(t){(\tilde \Lambda _i^{[{\rm{\Phi }}]}(t))^{\rm{T}}}] $, 任意两个局部参数估值器之间的估计误差互协方差为$ {P_{\Lambda _{ij}^{[{\rm{\Phi }}]}}}(t) = {\rm{E}}[\tilde \Lambda _i^{{\rm{[\Phi }}]}(t){(\tilde \Lambda _j^{[{\rm{\Phi]}}}(t))^{\rm{T}}}] $, 其中估计误差$ \tilde \Lambda _i^{[\Phi]}(t) = \Lambda _{}^{[\Phi]} - \hat \Lambda _i^{[\Phi]}(t) = S{\tilde \vartheta _{Ai}}(t) $, 则可获得
$$ \begin{equation} {P_{\Lambda _i^{[{\rm{\Phi }}]}}}(t) = S{P_{{\vartheta _{Ai}}}}(t){S^{\rm{T}}}, {P_{\Lambda _{ij}^{[{\rm{\Phi }}]}}}(t) = S{P_{{\vartheta _{Aij}}}}(t){S^{\rm{T}}} \end{equation} $$ (38) 2) 如果$ f({\vartheta _A}) $是非线性函数, 我们将$ f({\vartheta _A}) $在点$ {\hat \vartheta _{Ai}}(t - 1) $处进行线性化, 有
$$ \begin{equation} \Lambda _i^{[\Phi ]}(t) \approx S({\hat \vartheta _{Ai}}(t - 1)){\vartheta _A} + \gamma ({\hat \vartheta _{Ai}}(t - 1)) \end{equation} $$ (39) 那么, 基于传感器$ i $的数据获得的$ \Phi $中未知模型参数在$ t $时刻的局部估计为
$$ \begin{equation} \hat \Lambda _i^{[\Phi ]}(t) \approx S({\hat \vartheta _{Ai}}(t - 1)){\hat \vartheta _{Ai}}(t) + \gamma ({\hat \vartheta _{Ai}}(t - 1)) \end{equation} $$ (40) 类似地, 可以得到相应的估计误差方差$ {P_{\Lambda _i^{[{\rm{\Phi }}]}}}(t) $和互协方差$ P_{ij}^{{\Lambda ^{[{\rm{\Phi }}]}}}(t) $分别为
$$ \begin{align} &{P_{\Lambda _i^{[{\rm{\Phi }}]}}}(t) = S({\hat \vartheta _{Ai}}(t - 1)){{{P}}_{{\vartheta _{Ai}}}}(t){S^{\rm{T}}}({\hat \vartheta _{Ai}}(t - 1)) \\ &{P_{\Lambda _{ij}^{[{\rm{\Phi }}]}}}(t) = S({\hat \vartheta _{Ai}}(t - 1)){P_{{\vartheta _{Aij}}}}(t){S^{\rm{T}}}({\hat \vartheta _{Aj}}(t - 1)) \end{align} $$ (41) 其中$ i, j = 1, \cdots, L $.通过以上算法, 基于$ L $个传感器的数据可获得$ \Phi $中未知参数在时刻$ t $处的局部估值$ \hat \Lambda _i^{[\Phi]}(t), i, j = 1, \cdots, L $.由于对$ \Phi $中未知参数估计了$ L $次, 因此我们可应用线性无偏最小方差加权融合估计算法[23]将它们进行融合处理.下面定理1给出了未知参数分布式融合估计的结果.
定理1. 基于局部参数估计$ \hat \Lambda _i^{[{\rm{\Phi }}]}(t) $、局部参数估计误差方差阵$ {P_{\Lambda _i^{[{\rm{\Phi }}]}}}(t) $, 以及参数估计误差互协方差$ {P_{\Lambda _{ij}^{[{\rm{\Phi }}]}}}(t) $, 可得在线性无偏最小方差意义下的矩阵加权参数融合辨识器如下:
$$ \begin{equation} \hat \Lambda _o^{{\rm{[\Phi ]}}}(t) = \sum\limits_{i = 1}^L {{W_{\Lambda _i^{{\rm{[\Phi }}]}}}(t)\hat \Lambda _i^{[{\rm{\Phi }}]}(t)} \end{equation} $$ (42) 参数融合加权矩阵计算为
$$ \begin{align} [{W_{\Lambda _i^{[{\rm{\Phi }}]}}}(t), &\cdots , {W_{\Lambda _L^{[{\rm{\Phi }}]}}}(t)] = \\ &{(e_{{\Lambda ^{[\Phi ]}}}^{\rm{T}}P_{{\Lambda ^{[\Phi ]}}}^{ - 1}(t){e_{{\Lambda ^{[\Phi ]}}}})^{ - 1}}e_{{\Lambda ^{[\Phi ]}}}^{\rm{T}}P_{{\Lambda ^{[\Phi ]}}}^{ - 1}(t) \end{align} $$ (43) 其中$ {e_{{\Lambda ^{[\Phi]}}}} = {[{\begin{array}{*{20}{c}} {{I_{{n_\Phi }}}} & \cdots & {{I_{{n_\Phi }}}} \\ \end{array}}]^{\rm T}}_{{n_\Phi }L \times {n_\Phi }}, {P_{{\Lambda ^{[\Phi]}}}}(t) = {[{{P_{\Lambda _{ij}^{[{\rm{\Phi }}]}}}(t|t)}]_{{n_\Phi }L \times {n_\Phi }L}} $是以$ {P_{\Lambda _{ij}^{[{\rm{\Phi }}]}}}(t|t) $为第$ (i, j) $元素的分块矩阵.
参数融合辨识器的估计误差方差阵计算为
$$ \begin{equation} {P_{\Lambda _o^{[{\rm{\Phi }}]}}}(t) = {(e_{{\Lambda ^{[\Phi ]}}}^{\rm{T}}P_{{\Lambda ^{[\Phi ]}}}^{ - 1}(t){e_{{\Lambda ^{[\Phi ]}}}})^{ - 1}} \end{equation} $$ (44) 且$ {P_{\Lambda _o^{{\rm{[\Phi]}}}}}(t) \le {P_{\Lambda _i^{[{\rm{\Phi }}]}}}(t), i = 1, \cdots, L. $
基于以上算法可以获得未知模型参数融合辨识器$ \Lambda _o^{[{\rm{\Phi]}}}(t) $, 进而可获得$ \Phi $的融合估值$ {{\rm{\hat \Phi }}_o}(t) $.由前面的分析, 有$ {\hat \Phi _o}(t) \to \Phi, t \to \infty, w.p.1 $.
注2. 文献[5-7, 24]在进行模型参数融合辨识时, 将各传感器辨识得到的模型参数采用加权平均方法进行融合处理.该方法不能保证所获得的参数融合辨识器的估计精度不低于所有的局部参数估计.而本文采用线性无偏最小方差分布式矩阵加权融合算法[23]对各传感器辨识得到的参数进行融合处理.所获得的参数融合辨识器的估计精度不低于所有的局部参数估计.因此, 本文的矩阵加权融合的参数估计精度高于加权平均融合的参数估计精度.这在后面的仿真研究中也能看到.
4. $ \boldsymbol {\mu _i(t)} $的数学期望与方差辨识
当各传感器的衰减观测率未知时, 为了能应用第2节中的算法获得状态估计, 我们需要辨识描述衰减观测的随机变量$ \{ {\mu _i}(t)\} $的均值和方差.下面我们采用相关函数来辨识它们.
将模型参数融合辨识器$ {\hat \Phi _o}(t) $代入式(6)中, 得到
$$ \begin{equation} \hat X(t + 1) = {\hat \Phi _o}(t) \hat X(t)\hat \Phi _o^{\rm{T}}(t) + {\rm{\Gamma }}{Q_{\boldsymbol w}}{{\rm{\Gamma }}^{\rm{T}}} \end{equation} $$ (45) 由式(2)可计算零步相关函数$ R_i(t, 0) = {\rm E}[y_i^2(t)] $为
$$ \begin{equation} {R_i}(t, 0) = (\alpha _i^2 + \sigma _i^2){h_{i}}{{\hat X}}(t)h_{i}^{\rm{T}} + {{{Q}}_{{v_i}}} \end{equation} $$ (46) 一步相关函数$ {{{R}}_i}(t, {\rm{1}}) = {\rm{E}}[{y_i}(t){y_i}(t- {\rm{1}})] $为
$$ \begin{equation} {R_i}(t, 1) = \alpha _i^2{ h_{i}}{\hat \Phi _o}(t - 1)\hat X(t - 1) h_{i}^{\rm{T}} \end{equation} $$ (47) 零阶和一阶相关函数阵$ {{{R}}_i}(t, r), r = 0, {\rm{1}} $可通过如下采样相关函数$ {{{\hat R}}_i}(t, r) $逼近:
$$ \begin{align} {\hat R_i}(t, r) \approx &\frac{1}{t}\sum\limits_{s = 1}^t {{{ y}_i}(s){{y}_i}(s - r)} = {\hat R_i}(t - 1, r) + \\ & \frac{1}{t}[{y_i}(t){y_i}(t - r) - {\hat R_i}(t - 1, r)] \end{align} $$ (48) 初值为$ {{{\hat R}}_i}(0, 0) = {\rm E}[y_i^2{\rm{(0)]}}, \; {{{\hat R}}_i}(0, 1){\rm{ = 0}} $.
最后, 由式(47)可求出随机变量$ {\rm{\{ }}{\mu _i}{\rm{(}}t{\rm{)\} }} $的数学期望为
$$ \begin{equation} {\hat \alpha _i}(t) = \sqrt {\frac{{{{\hat R}_i}(t, 1)}}{{{{ h}_i}{{\hat \Phi }_o}(t - 1)\hat X(t - 1){ h}_i^{\rm{T}}}}} \end{equation} $$ (49) 将式(49)代入式(46)可得随机变量$ {\rm{\{ }}{\mu _i}{\rm{(}}t{\rm{)\} }} $的方差为
$$ \begin{equation} \hat \sigma _i^2(t) = \frac{{{{\hat R}_i}(t, 0) - {{{Q}}_{{v_i}}}}}{{ {h_{i}}{{\hat X}}(t)h_{i}^{\rm{T}}}} - \hat \alpha _i^2(t) \end{equation} $$ (50) 通过式(49)可获得$ {\alpha _i} $在$ t $时刻的估值.由第3节可知$ {\hat \Phi _o}(t) \to \Phi, t \to \infty $.由$ \Phi $的稳定性和式(45)可得$ \hat X(t) \to X(t), t \to \infty $.又由随机过程的遍历性可得$ {{{\hat R}}_i}(t, r) \to {{{R}}_i}(t, r), t \to \infty $.因而, 有估值$ {\hat \alpha _i}(t) $, $ \hat \sigma _i^2(t) $具有一致性, 即
$$ \begin{equation} {\hat \alpha _i}(t) \to {\alpha _i}, \hat \sigma _i^2(t) \to \sigma _i^2, t \to \infty , w.p.1 \end{equation} $$ (51) 进一步, 有$ {\hat Q_{{V_i}}}(t) = \hat \sigma _i^2(t){{{{ h}}}_i}\hat X(t){{ h}_i^{\rm{T}}} + {Q_{{v_i}}} \to {Q_{{V_i}}}(t) = \sigma _i^2{{{ h}}_i}X(t){{{ h}}_i^{\rm{T}}} + {Q_{{v_i}}}, t \to \infty, w.p.1 $.
将每时刻辨识的$ {\hat \Phi _o}(t) $、$ {\hat \alpha _i}(t) $和$ {\hat Q_{{V_i}}}(t) $替代第2节中的局部最优滤波器、互协方差阵和分布式融合滤波算法中的$ \Phi $、$ {\alpha _i} $和$ {Q_{{V_i}}}(t) $, 可得到相应的自校正状态估计算法.为了后文阐述方便, 记相应的自校正局部滤波器、预报误差方差、滤波误差协方差、增益、融合滤波器分别为$ {\hat {\boldsymbol x}_{si}}(t|t) $、$ {\hat \Sigma _{si}}(t|t - 1) $、$ {\hat P_{sij}}(t|t) $、$ {\hat K_{si}}(t) $和$ {\hat {\boldsymbol x}_s}(t|t) $.
5. 自校正滤波器的收敛性分析
下面的引理4给出了DESA方法, 应用DESA方法可证明自校正融合状态滤波器$ {\hat {\boldsymbol x}_s}(t|t) $的收敛性.
引理4[24]. 考虑动态误差系统
$$ \begin{equation} \delta (t) = T(t)\delta (t - 1) + u(t) \end{equation} $$ (52) 其中$ t \ge 0 $, 输出$ \delta (t) \in{{\bf R}}^n $, 输入$ u(t) \in {{\bf R}}^n $, 并且矩阵$ T(t) \in {{{\bf R}}^{n \times n}} $是一致渐近稳定的.若$ u(t) $是有界的, 则$ \delta {\rm{(}}t{\rm{)}} $是有界的.当$ t \to \infty $时, 若$ u(t) \to 0 $, 则$ \delta (t) \to 0 $.
记带已知参数$ \Phi $、$ {\alpha _i} $和$ \sigma _i^2 $的系统(1)和(3)为$ ({\rm{\Phi }}, {\rm{\Gamma}}, {F_i}, {Q_{\boldsymbol w}}, {Q_{{V_i}}}(t)) $, 相应的带时变参数$ {{\rm{\hat \Phi }}_o}(t) $、$ {\hat \alpha _i}(t) $和$ \hat \sigma_i^2(t) $的系统(1)和(3)为$ ({{\rm{\hat \Phi }}_o}(t), {\rm{\Gamma }}, {\hat F_i}(t), {Q_{\boldsymbol w}}, {\hat Q_{{V_i}}}(t)) $.
引理5[24]. 在假设1$ \, \sim\, $4下, 系统$ ({\rm{\Phi }}, {\rm{\Gamma }}, $ $ {F_i}, $ $ {Q_{\boldsymbol w}}, $ $ {Q_{{V_i}}}(t)) $和$ ({{\rm{\hat \Phi }}_o}(t), \Gamma, {\hat F_i}(t), {Q_{\boldsymbol w}}, {\hat Q_{{V_i}}}(t)) $是一致完全可观和一致完全可控的; 系统$ ({\rm{\Phi }}, \Gamma, {F_i}, {Q_{\boldsymbol w}}, {Q_{{V_i}}}(t)) $的最优局部滤波器的状态转移阵$ {\Psi _{fi}}(t) = [{I_n}-{K_i}(t){F_i}]{\rm{\Phi }} $是一致渐近稳定的, 且$ {K_i}(t) $是有界的; 系统$ ({{\rm{\hat \Phi }}_o}(t), \Gamma, {\hat F_i}(t), {Q_{\boldsymbol w}}, {\hat Q_{{V_i}}}(t)) $的自校正局部滤波器的状态转移阵$ {\hat \Psi _{sfi}}(t) = [{I_n}-{\hat K_{si}}(t){\hat F_i}{\rm{(}}t{\rm{)}}]{{\rm{\hat \Phi }}_o}(t) $是一致渐近稳定的, 且$ {\hat K_{si}}(t) $是有界的.而且$ \Delta {\hat K_i}(t) = {\hat K_{si}}(t)-{K_i}(t)\to 0 $和$ \Delta {\hat \Psi _{fi}}(t) = {\hat \Psi _{sfi}}(t)-{\Psi _{fi}}(t) \to 0 $, $ t \to\infty $.自校正局部滤波误差方差阵$ {\hat P_{si}}(t|t) $收敛于最优局部滤波误差方差阵$ {P_i}(t|t) $, 自校正滤波误差互协方差阵$ {\hat P_{sij}}(t|t) $收敛于最优滤波误差互协方差阵$ {P_{ij}}(t|t) $.
定理2. 在假设1$ \, \sim\, $4下, 自校正局部滤波器$ {\hat {\boldsymbol x}_{si}}(t|t) $收敛于最优局部滤波器$ {\hat {\boldsymbol x}_i}(t|t) $, 即
$$ \begin{equation} [{\hat {\boldsymbol x}_{si}}(t|t) - {\hat {\boldsymbol x}_i}(t|t)] \to 0, t \to \infty \end{equation} $$ (53) 证明. 由式(7)可得自校正滤波器为
$$ \begin{equation} {\hat {\boldsymbol x}_{si}}(t|t) = {\hat \Psi _{sfi}}(t){\hat {\boldsymbol x}_{si}}(t - 1|t - 1)+{\hat K_{si}}(t){y_i}(t) \end{equation} $$ (54) 注意到$ {\hat K_{si}}(t) $和$ {y_i}(t) $是有界的, $ {\hat \Psi _{sfi}}(t) $是一致渐近稳定的, 应用引理4得到$ {\hat {\boldsymbol x}_{si}}(t|t) $是有界的.令$ {\delta _i}(t) = {\hat {\boldsymbol x}_{si}}(t|t) - {\hat {\boldsymbol x}_i}(t|t) $, 式(54)减式(7)有动态误差方程
$$ \begin{equation} {\delta _i}(t) = {\Psi _{fi}}(t){\delta _i}(t - 1) + {u_i}(t) \end{equation} $$ (55) 其中$ {u_i}(t) = \Delta {\hat \Psi _{fi}}(t){\hat {\boldsymbol x}_{si}}(t - 1|t - 1) + \Delta {\hat K_i}(t){y_i}(t) $.根据$ {\hat {\boldsymbol x}_{si}}(t|t) $和$ { y_i}(t) $的有界性, 由$ \Delta {\hat K_i}(t)\to 0 $和$ \Delta{\hat \Psi _{fi}}(t){\to0} $, 有$ {u_i}(t)\to0 $.对式(55)应用引理4, 当$ t \to \infty $时, 有$ {\delta _i}(t) \to 0 $, 即式(53)成立.
定理3. 在假设1$ \, \sim\, $4下, 自校正加权融合滤波器$ {\hat {\boldsymbol x}_s}(t|t) $收敛于最优加权融合滤波器$ {\hat {\boldsymbol x}_o}(t|t) $, 即
$$ \begin{equation} [{\hat {\boldsymbol x}_s}(t|t) - {\hat {\boldsymbol x}_o}(t|t)] \to 0, t \to \infty \end{equation} $$ (56) 证明. 由引理3和引理5有$ {W_i}(t) $有界且$ \Delta {\hat W_i}(t) = [{\hat W_{si}}(t) -{W_i}(t)] \to 0 $.由式(53)以及$ {\hat {\boldsymbol x}_{si}}(t|t) $的有界性可得
$$ \begin{align} {\hat {\boldsymbol x}_s}(t|t) & - {\hat {\boldsymbol x}_o}(t|t) = \sum\limits_{i = 1}^L {{W_i}(t)[{{\hat {\boldsymbol x}}_{si}}(t|t) - {{\hat {\boldsymbol x}}_i}(t|t)} ]+ \\ &\sum\limits_{i = 1}^L {\Delta {{\hat W}_i}(t){{\hat {\boldsymbol x}}_{si}}(t|t)} \to 0 \end{align} $$ (57) 即式(56)成立.
6. 仿真例子
考虑带三传感器系统(1)和(2), 其中系数阵为$ \Phi = \left[{\begin{array}{*{20}{c}} {{a_{11}}} & {{a_{12}}} \\ {0.4} & { - 0.8} \\ \end{array}} \right], \Gamma = \left[{\begin{array}{*{20}{c}} {0.5} \\ {0.6} \\ \end{array}} \right], {{h}_1} = \left[{\begin{array}{*{20}{c}} {0.5} & {1.2} \\ \end{array}} \right] $, $ {{ h}_2} = \left[{\begin{array}{*{20}{c}} {0.6} & {1.9} \\ \end{array}}\right], {{ h}_3} = \left[{\begin{array}{*{20}{c}} {1.4} & 2 \\ \end{array}} \right], $噪声方差为$ {Q_{\boldsymbol w}} = 3, {Q_{{v_1}}} = 2, {Q_{{v_2}}} = 0.4, {Q_{{v_3}}} = 1, \left\{ {{\mu _i}(t)} \right\} $, $ i = 1, 2, 3 $为在$ [0, 1] $区间取值的标量随机变量.取初值$ {\hat{\boldsymbol x}_i}(0|0) = 0, {P_{ij}}(0|0) = 0.1{I_2}, i, j = 1, 2, 3 $.假设$ {a_{11}}, {a_{12}}, \{ {\mu _i}(t)\} $的数学期望$ {\alpha _i} $与方差$ \sigma _i^2, i = 1, 2, 3 $, 未知.目的是辨识未知模型参数$ {a_{11}} $和$ {a_{12}} $、期望$ {\alpha _i} $和方差$ \sigma _i^2, i = 1, 2, 3 $并求自校正融合状态滤波器.
在仿真中假设未知模型参数$ {a_{11}} = 0.6 $, $ {a_{12}} = -0.2 $, $ {\mu _i}(t), i = 1, 2, 3 $, 的概率分布分别为$ P{\rm{\{ }}{\mu _1}(t){\rm{ = 0.3\} = 0.3}} $, $ P{\rm{\{ }}{\mu _1}(t){\rm{ = 0.5\} = 0.2}} $, $ P{\rm{\{ }}{\mu _1}(t){\rm{ = 1\} = 0.5}} $, $ P{\rm{\{ }}{\mu_2}(t){\rm{ = 0.4\} = 0.4}} $, $ P{\rm{\{ }}{\mu_2}(t) {\rm{ = 0.7\} = 0.3}} $, $ P{\rm{\{ }}{\mu_2}(t){\rm{ = 0.9\} = 0.3}} $, $ P{\rm{\{ }}{\mu_3}(t){\rm{ = 0.1\}}} $ $ {{ = 0.2}} $, $ P{\rm{\{ }}{\mu_3}(t){\rm{ = 0.6\} = 0.6}} $, $ P{\rm{\{ }}{\mu_3}(t){\rm{ = 0.9\} = 0.2}} $.我们可以计算$ {\mu _i}(t) $, $ i = 1, 2, 3 $的数学期望和方差分别为$ {\alpha _1} = 0.69 $, $ {\alpha _2}{\rm{ = }}0.64 $, $ {\alpha_3} = 0.56 $, $ \sigma _1^2 = 0.1009 $, $ \sigma_2^2 = 0.0444 $, $ \sigma_3^2 = 0.0664 $.
为了与文献[5-7, 24]中的参数加权平均融合算法相比较.如下给出了模型参数$ {a_{11}} $和$ {a_{12}} $的局部估计误差方差、加权平均估计误差方差以及分布式加权融合估计误差方差算法.
a) 局部参数估计和分布式加权融合估计:
根据前面的第3节, 可知未知模型参数与$ {\vartheta _A} = {[{a_1}, {a_2}]^{\rm{T}}} $有如下线性关系
$$ \begin{equation} {\Lambda ^{[\Phi ]}} = S{\vartheta _A} + \gamma \end{equation} $$ (58) 其中$ {\Lambda ^{[\Phi]}} = \left[{\begin{array}{*{20}{c}} {{a_{11}}} \\ {{a_{12}}} \\ \end{array}} \right] $, $ S = \left[{\begin{array}{*{20}{c}} { - 1} & 0 \\ {{ {{ - {a_{22}}} \over {{a_{21}}}}}} & {{ {{ - 1} \over {{a_{21}}}}}} \\ \end{array}} \right] $, $ \gamma = \left[{\begin{array}{*{20}{c}} {{a_{22}}} \\ {{ {{a_{_{22}}^2} \over {{a_{21}}}}}} \\ \end{array}} \right]. $
${\Lambda ^{[{\rm{\Phi }}]}}$的局部估值为
$$ \begin{equation} \hat \Lambda _i^{[\Phi ]}(t) = S{\hat \vartheta _{Ai}}(t) + \gamma , i = 1, 2, 3 \end{equation} $$ (59) $ {\hat \vartheta _{Ai}}(t) $可由最小二乘算法辨识, 将$ S, {\hat \vartheta _{Ai}}(t), \gamma $代入式(38)中, 可分别获得局部估值$ \hat \Lambda _i^{[\Phi]}(t) $和局部估计误差方差$ {P_{\Lambda _i^{[\Phi]}}}(t) $以及误差互协方差$ {P_{\Lambda _{ij}^{[\Phi]}}}(t), i, j = 1, 2, 3 $.再根据定理1中参数融合估计算法, 可得参数融合估计误差方差$ {P_{\Lambda _0^{[{\rm{\Phi }}]}}}(t) $.
b) 参数的加权平均融合估计:
${\Lambda ^{[{\rm{\Phi }}]}}$的加权平均估值为
$$ \begin{equation} \bar \Lambda _{}^{[\Phi ]}(t) = { {1 \over 3}}\sum\limits_{i = 1}^3 {\hat \Lambda _i^{[\Phi ]}(t)} \end{equation} $$ (60) 定义加权平均估计误差方差为$ {P_{\bar \Lambda _{}^{[\Phi]}}}(t) = {\rm{E}}[\tilde {\bar \Lambda}_{}^{[\Phi]}(t)(\tilde {\bar \Lambda} _{}^{[\Phi]}(t)){^{\rm{T}}}] $, 其中$ \tilde {\bar \Lambda} _{}^{[\Phi]}(t) = {\Lambda ^{[\Phi]}} -\bar \Lambda _{}^{[\Phi]}(t) = { {1 \over 3}}\sum_{i = 1}^3 {(S{{\tilde \vartheta }_{Ai}}(t))} $.于是, 加权平均估计误差方差可计算为
$$ \begin{equation} {P_{\bar \Lambda _{}^{[\Phi ]}}}(t) = { {1 \over 9}}S\left(\sum\limits_{i = 1}^3 {\sum\limits_{j = 1}^3 ({P_{{\vartheta _{Aij}}}}}(t))\right){S^{\rm{T}}} \end{equation} $$ (61) 图 1给出了应用部分3参数融合辨识算法获得的未知模型参数融合辨识结果.由图可知随着时间的增长, 辨识结果收敛于真值.
图 2和图 3分别给出了未知模型参数$ {a_{11}} $和$ {a_{12}} $的局部估计误差方差、加权平均估计误差方差和分布式融合估计误差方差的比较结果.由图可知分布式加权融合辨识误差方差小于各局部辨识误差方差和加权平均辨识误差方差.图中$ {S_i}, i = 1, 2, 3 $表示第$ i $个传感器的局部辨识的误差方差, DWF表示分布式加权融合辨识的误差方差, WAEV表示加权平均融合辨识的误差方差.
图 4和图 5给出了应用部分4中辨识算法分别对不同传感器的随机变量$ \{ {\mu _i}(t)\}, i = 1, 2, 3 $的数学期望和方差进行辨识的结果.曲线表示辨识结果, 直线表示相应的真值.由图可知随着时间的增长, 辨识结果收敛于真值.图 6和图 7给出了自校正融合状态滤波器, 可见自校正融合估计具有有效性.
图 4 $\mu_{i}(t)$的数学期望辨识Fig. 4 Identification of Mathematical expectation of $\mu_{i}(t)$图 8和图 9给出了局部和融合的最优与自校正状态估计误差方差图.由图可见, 各局部自校正误差方差收敛于局部最优误差方差, 自校正融合误差方差收敛于最优融合误差方差, 即自校正滤波器具有渐近最优性.而且自校正融合滤波器比各局部自校正滤波器具有更高精度.图中$ {S_i}, i = 1, 2, 3 $表示第$ i $个传感器的局部自校正估计误差方差, SF表示自校正融合估计误差方差, 直线表示相应的最优方差.
图 8 局部、融合最优与自校正状态分量1的滤波误差方差Fig. 8 Variance of the first state component of local, fusion optimal and self-tuning filters
图 9 局部、融合最优与自校正状态分量2的滤波误差方差Fig. 9 Variance of the second state component of local, fusion optimal and self-tuning filters目前参考文献[5-7, 24]中的自校正滤波算法大都没有考虑传感器的衰减观测现象.图 10给出了传感器存在衰减观测而没有给予考虑的自校正融合滤波器与本文考虑衰减观测的自校正融合滤波器在30次蒙特卡洛实验下均方误差迹的比较.可见, 在传感器存在衰减观测时, 本文考虑衰减观测的自校正融合滤波器具有更高的精度.
图 10 考虑衰减观测与没有考虑衰减观测自校正融合滤波器的均方误差的迹Fig. 10 Trace of mean square error of the self-tuning fusion filters with/without considering fading measurements7. 结论
对带未知模型参数和衰减观测率的多传感器随机系统, 应用RELS算法和相关函数分别对未知模型参数、描述衰减观测现象的随机变量的数学期望和方差进行在线实时辨识, 提出了线性无偏最小方差矩阵加权融合模型参数辨识器.与已有文献的加权平均融合模型参数辨识算法相比, 本文所提出的线性无偏最小方差矩阵加权融合参数辨识算法具有更高的估计精度.将实时辨识的模型参数、数学期望和方差代入到最优局部和融合状态估计算法中获得了相应的自校正状态滤波算法.利用DESA方法证明了自校正状态滤波器收敛于最优状态滤波器.与现有文献的带未知模型参数的自校正估计算法相比, 本文还考虑了传感器的衰减观测现象, 并给出了采用相关函数辨识衰减观测的数学期望和方差的算法.
-
表 1 TE过程变量
Table 1 Process variables in the TE process
变量描述 变量描述 1) 物料 A 流量 12) 气/液分离器液位 2) 物料 D 流量 13) 气/液分离器压力 3) 物料 E 流量 14) 气/液分离器出口流量 4) 物料 C 流量 15) 汽提塔液位 5) 压缩机返回物料流量 16) 汽提塔压力 6) 反应器给料流量 17) 汽提塔塔底流量 7) 反应器压力 18) 汽提塔温度 8) 反应器液位 19) 汽提塔蒸汽流量 9) 反应器温度 20) 压缩机功率 10) 排空物料流量 21) 反应器冷却水出口温度 11) 气/液分离器温度 22) 冷凝器冷却水出口温度 
下载: 导出CSV
表 2 验证数据集
Table 2 Data sets used for validation
故障程度 样本数量 正常状态样本 (TE 标准数据) 500 (训练集) 故障数据 (TE 标准数据) 960 (测试集) 正常状态样本 (生成数据 FS = 0) 500 + 160 (训练集 + 测试集) 故障程度 1 (生成数据 FS = 0.2) 160 (测试集) 故障程度 2 (生成数据 FS = 0.4) 160 (测试集) 故障程度 3 (生成数据 FS = 0.6) 160 (测试集) 故障程度 4 (生成数据 FS = 0.8) 160 (测试集) 故障程度 5 (生成数据 FS = 1.0) 160 (测试集)
下载: 导出CSV
表 3 两种方法的性能比较
Table 3 Comparison of the two methods
方法 误报率 (%) 检测率 (%) 建模运行时间 (s) KPLS 3.75 98.375 0.307 CCA-SVDD 0 97.75 0.033
下载: 导出CSV
-
[1] 彭开香, 马亮, 张凯. 复杂工业过程质量相关的故障检测与诊断技术综述. 自动化学报, 2017, 43(3): 349-365Peng Kai-Xiang, Ma Liang, Zhang Kai. Review of quality-related fault detection and diagnosis techniques in complex industrial processes. Acta Automatica Sinica, 2017, 43(3): 349-365 [2] Krzanowski W J. A user's guide to principal components. Journal of the Royal Statistical Society: Series A, 1992, 155(1): 173-174 doi: 10.2307/2982678 [3] Qin S J. Survey on data-driven industrial process monitoring and diagnosis. Annual Reviews in Control, 2012, 36(2): 220-234 doi: 10.1016/j.arcontrol.2012.09.004 [4] 马洁, 李钢, 陈默. 基于非线性故障重构的旋转机械故障预测方法. 自动化学报, 2014, 40(9): 2045-2049Ma Jie, Li Gang, Chen Mo. Fault prediction method of rotating machinery based on nonlinear fault reconstruction. Acta Automatica Sinica, 2014, 40(9): 2045-2049 [5] 王静, 胡益, 侍洪波. 基于GMM的间歇过程故障检测. 自动化学报, 2015, 41(5): 899-905Wang Jing, Hu Yi, Shi Hong-Bo. Intermittent process fault detection based on GMM. Acta Automatica Sinica, 2015, 41(5): 899-905 [6] Zhou D H, Li G, Qin S J. Total projection to latent structures for process monitoring. AIChE Journal, 2010, 56(1): 168-178 [7] Ding S X, Yin S, Peng K X, Hao H Y, Shen B. A novel scheme for key performance indicator prediction and diagnosis with application to an industrial hot strip mill. IEEE Transactions on Industrial Informatics, 2013, 9(4): 2239-2247 doi: 10.1109/TII.2012.2214394 [8] Huang J P, Yan X F. Quality-driven principal component analysis combined with kernel least squares for multivariate statistical process monitoring. IEEE Transactions on Control System Technology, 2019, 27(6): 2688-2695 doi: 10.1109/TCST.2018.2865130 [9] Wang G, Yin S. Quality-related fault detection approach based on orthogonal signal correction and modified PLS. IEEE Transactions on Industrial Informatics, 2015, 11(2): 398-405 [10] Jiang Y C, Yin S. Recent advances in key-performance-indicator oriented prognosis and diagnosis with a matlab toolbox: DB-KIT. IEEE Transactions on Industrial Informatics, 2019, 15(5): 2849-2858 doi: 10.1109/TII.2018.2875067 [11] 尚林源, 田学民, 曹玉苹, 蔡连芳. 基于PLS交叉积矩阵非相似度分析的MPC性能监控与诊断. 自动化学报, 2017, 43(2): 271-279Shang Lin-Yuan, Tian Xue-Min, Cao Yu-Ping, Cai Lian-Fang. Performance monitoring and diagnosis of MPC based on PLS cross product matrix non-similarity analysis. Acta Automatica Sinica, 2017, 43(2): 271-279 [12] Chen Z, Ding S X, Zhang K. Canonical correlation analysis-based fault detection methods with application to alumina evaporation process. Control Engineering Practice, 2016, 46: 51-58 doi: 10.1016/j.conengprac.2015.10.006 [13] Zhang K, Peng K X, Ding S X, Chen Z W, Yang X. A correlation-based distributed fault detection method and its application to a hot tandem rolling mill process. IEEE Transactions on Industrial Electronics, 2013, 9(4): 2226-2238 doi: 10.1109/TII.2013.2243743 [14] David M J T, Robert P W D. Support vector domain description. Pattern Recognition Letters, 1999, 20(11): 1191-1199 [15] Khediri I B, Weihs C, Limam M. Kernel k-means clustering based local support vector domain description fault detection of multimodal processes. Expert Systems with Applications, 2012, 39(2): 2166-2171 doi: 10.1016/j.eswa.2011.07.045 [16] Zhang C F, Peng K X, Dong J. A novel plant-wide process monitoring framework based on distributed Gap-SVDD with adaptive radius. Neurocomputing, 2019, 350: 1-12 doi: 10.1016/j.neucom.2019.04.026 [17] Zhu K, Mei F, Zheng J. Adaptive fault diagnosis of HVCBs based on P-SVDD and P-KFCM. Neurocomputing, 2017, 240: 127-136 doi: 10.1016/j.neucom.2017.02.042 [18] Zhu Q, Qin S J. Supervised diagnosis of quality and process faults with canonical correlation analysis. Industrial and Engineering Chemistry Research, 2019, 58(26): 11213-11223 doi: 10.1021/acs.iecr.9b00320 [19] Ma L, Dong J, Peng K X. A novel key performance indicator oriented hierarchical monitoring and propagation path identification framework for complex industrial processes. ISA Transactions, 2020, 96(1): 1-13 [20] Luo L J, Lovelett R J, Ogunnaike B A. Hierarchical monitoring of industrial processes for fault detection, fault grade evaluation, and fault diagnosis. AIChE Journal, 2017, 63(7): 2781-2795 doi: 10.1002/aic.15662 [21] Song B, Zhou X G, Shi H B, Yang T. Performance indicator oriented concurrent subspace process monitoring method. IEEE Transactions on Industrial. Electronics, 2018, 65(2): 1508-1571 doi: 10.1109/TIE.2017.2733443 [22] Yang T, Shi H B, Song B, Tan S. Operating performance assessment and non-optimal cause identification for chemical process. Canadian. Journal of Chemical. Engineering, 2017, 97(1): 1475-1487 [23] Guo S, Sun Y, Wu F, Li Y. Integrating laplacian eigenmaps feature space conversion into deep neural network for equipment condition assessment. Automatic Control and Computer Sciences, 2018, 52(6): 465-475 doi: 10.3103/S0146411618060056 [24] Yan H, Liu K, Zhang X. Multiple sensor data fusion for degradation modeling and prognostics under multiple operational conditions. IEEE Transactions on Reliability, 2016, 65(3): 1416-1426 doi: 10.1109/TR.2016.2575449 [25] Sun C, Wang P, Yan R Q, Gao R X, Chen X F. Machine health monitoring based on locally linear embedding with kernel sparse representation for neighborhood optimization bearing fault diagnosis. Mechanical System and Signal Processing, 2019, 114: 25-34 doi: 10.1016/j.ymssp.2018.04.044 [26] Atamuradov V, Medjaher K, Camci F, Dersin P, Zerhouni N. Railway point machine prognostics based on feature fusion and health state assessment. IEEE Transactions on Instrumentation and Measurement, 2019, 68(8): 2691-2704 doi: 10.1109/TIM.2018.2869193 [27] Zhang F, Sun K, Wu X L. A novel variable selection algorithm for multi-layer perceptron with elastic net. Neurocomputing, 2019, 361: 110-118 doi: 10.1016/j.neucom.2019.04.091 [28] Downs J, Vogel E F. A plant-wide industrial process control problem. Computers and Chemical Engineering, 1993, 17(3): 245-255 doi: 10.1016/0098-1354(93)80018-I -
下载:
计量
- 文章访问数: 1230
- HTML全文浏览量: 173
- PDF下载量: 417
- 被引次数: 0
