1
基于CESCKF的状态估计
本节主要针对一类带有白噪声的非线性系统采用CESCKF进行状态估计.
考虑如下非线性系统:
1 \begin{document}$ \begin{align} \begin{cases} x_k=f(x_{x-1}, u_{k-1})+\nu_{k-1}, \\ y_k=h(x_k, u_k)+g_k+\omega_k. \end{cases} \end{align} $\end{document}
其中: \begin{document}$ x_k\in {\cal R}^{n_x} $\end{document} 为第\begin{document}$ k $\end{document} 次采样的系统状态向量, \begin{document}$ f:{\cal R}^{n_x}\times {\cal R}^{n_u}\longrightarrow {\cal R}^{n_x} $\end{document} 和\begin{document}$ h:{\cal R}^{n_x}\times {\cal R}^{n_u}\longrightarrow {\cal R}^{n_y} $\end{document} 为已知非线性函数, \begin{document}$ u_k\in {\cal R}^{n_u} $\end{document} 为控制输入, \begin{document}$ y_k\in {\cal R}^{n_y} $\end{document} 为测量输出, 过程高斯噪声序列\begin{document}$ \{\nu_k\} $\end{document} 和测量高斯噪声序列\begin{document}$ \{\omega_k\} $\end{document} 是分别具有零均值与协方差\begin{document}$ \varXi_{\nu, k} $\end{document} 和\begin{document}$ \varXi_{\omega, k} $\end{document} 的独立序列, \begin{document}$ g_k $\end{document} 为系统的输出传感器故障.
CESCKF算法[6 , 8 ] 如下:
step 1:初始化
2 \begin{document}$ \begin{align} x_0\thicksim N(\hat{x}_0, P_0), \; \hat{x}_0={E}[x_0], \; S_0={\rm chol}(P_0). \end{align} $\end{document}
其中: \begin{document}$ {x_0} $\end{document} 服从均值为\begin{document}$ {{{\hat x}_0}} $\end{document} 、方差为\begin{document}$ {{P_0}} $\end{document} 的正态分布, \begin{document}$ {E} $\end{document} 为数学期望, chol\begin{document}$ (\star) $\end{document} 表示对矩阵\begin{document}$ \star $\end{document} 进行Cholesky分解得到的上三角矩阵.
step 2:先验状态估计. 估计第\begin{document}$ k $\end{document} 次先验状态估计\begin{document}$ \hat{x}_{k|k-1} $\end{document} , 即
3 \begin{document}$ \begin{align} \hat{x}_{k| k-1}=\sum^{2^{n+1}+1}_{i=1}{w_i x^{\ast}_{i, k|k-1}}, \; i=1, 2, \ldots, 2^{n+1}+1. \end{align} $\end{document}
其中: \begin{document}$ {w_i} $\end{document} 为权值, 取值为
\begin{document}$ \begin{align*} {w_i} = \begin{cases} \dfrac{{(2\delta _1^2 - 1)(2\delta _1^2 + 3)}}{{8\delta _1^4}}, \; i = 1;\\[9pt] \dfrac{{{{(2\delta _1^2 - 1)}^2}}}{{{2^{n + 3}}\delta _1^4}}, \; i = 2, 3, \ldots , {2^n} + 1;\\[9pt] \dfrac{1}{{{2^{n + 2}}\delta _1^4}}, \; i = {2^n} + 2, \ldots , {2^{n + 1}} + 1. \end{cases} \end{align*} $\end{document}
\begin{document}$ {\delta _1} $\end{document} 为可调参数, 且满足\begin{document}$ \delta _1^2 > 1/2 $\end{document} , \begin{document}$ {\delta _2} = \delta _1\sqrt {4\delta _1^2 -2}/ (2\delta _1^2 - 1 ) $\end{document} . \begin{document}$ x^{\ast}_{i, k|k-1} $\end{document} 为传播嵌入式容积点, 由下式得出:
4 \begin{document}$ \begin{align} & x_{i, k-1|k-1}=S_{k-1|k-1} \varepsilon_i + \hat{x}_{k-1|k-1}, \end{align} $\end{document}
5 \begin{document}$ \begin{align} & x^{\ast}_{i, k|k-1}=f(x_{i, k-1|k-1}, u_{k-1}). \end{align} $\end{document}
\begin{document}$ x_{i, k-1|k-1} $\end{document} 为嵌入式容积点, \begin{document}$ S_{k-1|k-1} $\end{document} 和\begin{document}$ \hat{x}_{k-1|k-1} $\end{document} 分别为\begin{document}$ k-1 $\end{document} 时刻误差协方差均方根和后验状态估计值. \begin{document}$ {\varepsilon _i} $\end{document} 为容积点, 取值为
\begin{document}$ \begin{align*} \varepsilon _i = \begin{cases} [0]_i^{\rm C}, \; i = 1;\\ \sqrt 2 [{\delta _2}]_i^{\rm CE}, \; i = 2, 3, \ldots , {2^n} + 1;\\ \sqrt 2 [{\delta _1}]_i^{\rm CE}, \; i = {2^n} + 2, \ldots , {2^{n + 1}} + 1. \end{cases} \end{align*} $\end{document}
上角标C和CE分别表示CKF算法与CECKF算法, \begin{document}$ {[0]^{\rm C}} $\end{document} 表示CKF的零点, \begin{document}$ {[{\delta _1}]^{\rm CE}} $\end{document} 和\begin{document}$ {[{\delta _2}]^{\rm CE}} $\end{document} 为CECKF的两种采样点.
估计第\begin{document}$ k $\end{document} 次先验测量估计值\begin{document}$ \hat{y}_{k|k-1} $\end{document} 及残差\begin{document}$ r_k $\end{document} 为
6 \begin{document}$ \begin{align} &\hat{y}_{k|k-1}=\sum^{2^{n+1}+n}_{i=1}{w_iy^{\ast}_{i, k|k-1}}, \; i=1, 2, \ldots, 2^{n+1}+1; \end{align} $\end{document}
7 \begin{document}$ \begin{align} & r_k=y_k-\hat{y}_{k|k-1}. \end{align} $\end{document}
其中: \begin{document}$ y_k $\end{document} 为传感器\begin{document}$ k $\end{document} 时刻的测量实际值; \begin{document}$ y^{\ast}_{i, k|k-1} $\end{document} 为传播嵌入式容积点, 可表示为
8 \begin{document}$ \begin{align} &x_{i, k|k-1}=S_{xx, k|k-1}\varepsilon_i+\hat{x}_{k|k-1}, \end{align} $\end{document}
9 \begin{document}$ \begin{align} &y^{\ast}_{i, k|k-1}=h(x_{i, k|k-1}, u_k). \end{align} $\end{document}
这里先验误差协方差均方根\begin{document}$ S_{xx, k|k-1} $\end{document} 为
10 \begin{document}$ \begin{align} S_{xx, k|k-1}={\rm qr}([M_{k|k-1}\quad S_{\nu, k-1}]). \end{align} $\end{document}
qr\begin{document}$ (\star) $\end{document} 表示对矩阵\begin{document}$ \star $\end{document} 进行三角分解得到的下三角矩阵; \begin{document}$ S_{v, k-1} $\end{document} 为过程高斯噪声协方差\begin{document}$ \varXi_{\nu, k-1} $\end{document} 的平方根矩阵, 即\begin{document}$ \varXi_{\nu, k-1}=S_{\nu, k-1}S_{\nu, k-1}^{\rm T} $\end{document} ; 矩阵\begin{document}$ M_{k|k-1} $\end{document} 为
11 \begin{document}$ \begin{align} M_{k|k-1}= &[\sqrt{w_1}(x^{\ast}_{1, k|k-1}-\hat{x}_{k|k-1}), \ldots, \\ &\sqrt{w_{2^{n+1}+1}}(x^{\ast}_{2^{n+1}+1, k|k-1}-\hat{x}_{k|k-1})]. \end{align} $\end{document}
step 3:后验测量更新. 估计第\begin{document}$ k $\end{document} 次后验状态估计\begin{document}$ \hat{x}_{k|k} $\end{document} , 有
12 \begin{document}$ \begin{align} {\hat x_{k| k }} = {\hat x_{k| {k - 1}}} + {K_{c, k}}( {{y_k} - {{\hat y}_{k| {k - 1} }}} ), \end{align} $\end{document}
其中\begin{document}$ {K_{c, k}} $\end{document} 为卡尔曼增益, 由下式计算:
13 \begin{document}$ \begin{align} {K_{c, k}} = (P_{xy, k|k - 1}/S^{\rm T}_{yy, k| k - 1})/S_{yy, k|k - 1}. \end{align} $\end{document}
这里: 交叉协方差矩阵\begin{document}$ {P_{xy, }}_{k| {k - 1} } $\end{document} 和下三角矩阵\begin{document}$ {S_{yy, k| {k - 1}}} $\end{document} 分别为
14 \begin{document}$ \begin{align} &{P_{xy, }}_{k| {k - 1} } = {R_{k| {k - 1} }}{N^{\rm T}}_{k| {k - 1} }, \end{align} $\end{document}
15 \begin{document}$ \begin{align} &{S_{yy, k| {k - 1} }} = {\rm qr}( [N_{k| k - 1}\; \; S_{\omega , k}] ). \end{align} $\end{document}
\begin{document}$ {S_{\omega , k}} $\end{document} 为测量高斯噪声协方差\begin{document}$ {\varXi _{\omega , k}} $\end{document} 的平方根矩阵, 即\begin{document}$ {\varXi _{\omega , k}} = {S_{\omega , k}}S_{\omega , k}^{\rm T} $\end{document} . 矩阵\begin{document}$ {R_{k| {k - 1}}} $\end{document} 和\begin{document}$ {N_{k| {k - 1} }} $\end{document} 分别为
16 \begin{document}$ \begin{align} &{R_{k| {k - 1} }} = [ {\sqrt {{w_1}} ( {{x_{1, k| {k - 1}}} - {{\hat x}_{k| {k - 1} }}} ), \ldots , } \\ &\quad \quad \quad \quad \; {\sqrt {{w_{{2^{n + 1}} + 1}}} ( {{x_{( {{2^{n + 1}} + 1} ), k| {k - 1}}} - {{\hat x}_{k| {k - 1} }}} )} ], \end{align} $\end{document}
17 \begin{document}$ \begin{align} &{N_{k| {k - 1}}} = [ {\sqrt {{w_1}} ( {{y^ * }_{1, k| {k - 1}} - {{\hat y}_{k| {k - 1}}}} ), \ldots , }\\ &\quad \quad \quad \quad \; \sqrt {w_{2^{n + 1} + 1}} ( y^*_{( {2^{n + 1}} + 1), k| {k - 1}}-{\hat y}_{k| k - 1} ) ]. \end{align} $\end{document}
误差协方差的均方根为
18 \begin{document}$ \begin{align} {S_{k| k}} = {\rm qr}( {[{{R_{k| {k - 1}}}- {K_{c, k}}{N_{k| {k - 1} }}}\; \; {{K_{c, k}}{S_{\omega , k}}}]} ). \end{align} $\end{document}
跳转到step 2, 开始下一次滤波, 直到滤波周期结束.
2
SS-DEWMA图的多目标优化(MO-SS-DEWMA图)及故障检测
本节首先介绍SS-DEWMA图; 然后, 构造残差评价指标MDR和FAR与SS-DEWMA图控制宽度和平滑参数的函数, 通过分析两个参数对MDR和FAR的影响, 得出SS-DEWMA图难以选取合适的参数使得MDR和FAR同时达到最小. 为了解决上述问题, 在SS-DEWMA图中引入MO-PSO算法, 以MDR和FAR为两个目标函数, 得到离线优化的参数, 从而得到MO-SS-DEWMA图, 进而对故障进行在线检测. 其中, 采用小波分析算法对SS-DEWMA图进行多尺度分解、阈值去噪和重构, 减弱噪声对系统的影响.
2.1
SS-DEWMA图
SS-DEWMA图采用平方统计量之和, 可同时检测数据的过程均值与方差的变化. 由式(7)得到的第\begin{document}$ k $\end{document} 个残差向量\begin{document}$ {r_k} = {[ {{r_{k, 1}}, {r_{k, 2}}, \ldots , {r_{k, {n_y}}}} ]^{\rm T}} $\end{document} , \begin{document}$ {n_y} $\end{document} 表示样本向量维数, \begin{document}$ {r_k} $\end{document} 向量服从正态分布, 即\begin{document}$ {r_k} \sim N( \mu+ $\end{document} \begin{document}$ a\sigma , {b^2}{\sigma ^2}) $\end{document} . 其中: \begin{document}$ \mu $\end{document} 和\begin{document}$ \sigma $\end{document} 分别为无故障数据样本的均值和标准差; \begin{document}$ a $\end{document} 和\begin{document}$ b $\end{document} 分别表征数据样本的均值故障和方差故障参数, \begin{document}$ a = 0 $\end{document} , \begin{document}$ b = 1 $\end{document} 即为无故障, \begin{document}$ a \ne 0 $\end{document} 表示存在均值故障(如位移偏移、漂移故障等), \begin{document}$ b > {\rm{1}} $\end{document} 表示存在方差故障; \begin{document}$ k=1, 2, \ldots, N_r $\end{document} , \begin{document}$ {N_r} $\end{document} 为样本的大小. 于是, \begin{document}$ {n_y} $\end{document} 维向量\begin{document}$ [r_{k, 1}, r_{k, 2}, \ldots, r_{k, n_y}]^{\rm T} $\end{document} 组成样本向量为
\begin{document}$ \begin{align*} &\{[r_{1, 1}, r_{1, 2}, \ldots, r_{1, n_y}]^{\rm T}, \; [r_{2, 1}, r_{2, 2}, \ldots, r_{2, n_y}]^{\rm T}, \ldots, \\ &[r_{N_r, 1}, r_{N_r, 2}, \ldots, r_{N_r, n_y}]^{\rm T}\}. \end{align*} $\end{document}
SS-DEWMA图算法[13 ] 如下:
step 1:根据向量\begin{document}$ [r_{k, 1}, r_{k, 2}, \ldots, r_{k, n_y}]^{\rm T} $\end{document} 定义如下两个独立统计量\begin{document}$ {U_k} $\end{document} 和\begin{document}$ {V_k} $\end{document} :
19 \begin{document}$ \begin{align} & {U_k} = \dfrac{{{{\bar r}_k} - \mu }}{{{\sigma / {\sqrt {{n_y}} }}}}, \end{align} $\end{document}
20 \begin{document}$ \begin{align} & {V_k} = {\varPhi ^{ - 1}}\Big\{ {H\Big( {\dfrac{{({n_y} - 1)S_k^2}}{{{\sigma ^2}}};{n_y} - 1}\Big)}\Big\}. \end{align} $\end{document}
其中: \begin{document}$ {\bar r_k} = ( {{r_{k, 1}} + \ldots + {r_{k, {n_y}}}} )/{n_y} $\end{document} 表示第\begin{document}$ k $\end{document} 个样本向量均值; \begin{document}$ S_k^2 = \sum\limits_{i = 1}^{{n_y}} {{{{{( {{r_{k, j}} - {{\bar r}_k}} )}^2}} / {( {{n_y} - 1} )}}} $\end{document} , \begin{document}$ j = 1, 2, \ldots, $\end{document} \begin{document}$ {n_y} $\end{document} , 表示第\begin{document}$ k $\end{document} 个样本向量方差; \begin{document}$ {\varPhi ^{ - 1}} $\end{document} 表示标准正态分布函数的逆函数; \begin{document}$ \Big( {\dfrac{{({n_y} - 1)S_k^2}}{{{\sigma ^2}}};{n_y} - 1} \Big) $\end{document} 表示具有\begin{document}$ ( {{n_y} - 1} ) $\end{document} 自由度的卡方分布函数, \begin{document}$ \dfrac{{({n_y} - 1)S_k^2}}{{{\sigma ^2}}} $\end{document} 是独立卡方随机变量.
step 2:基于\begin{document}$ {U_k} $\end{document} 和\begin{document}$ {V_k} $\end{document} 定义两个EWMA统计量如下:
21 \begin{document}$ \begin{align} & {Y_k} = (1 - \lambda ){Y_{k - 1}} + \lambda {U_k}, \end{align} $\end{document}
22 \begin{document}$ \begin{align} & {Z_k} = (1 - \lambda ){Z_{k - 1}} + \lambda {V_k}. \end{align} $\end{document}
其中: \begin{document}$ \lambda $\end{document} 为平滑参数, \begin{document}$ 0 < \lambda < 1 $\end{document} . \begin{document}$ {Y_k} $\end{document} 和\begin{document}$ {Z_k} $\end{document} 的初值设置为\begin{document}$ {Y_0} = {Z_0} = 0 $\end{document} .
step 3:基于\begin{document}$ {Y_k} $\end{document} 和\begin{document}$ {Z_k} $\end{document} 定义两个DEWMA统计量如下:
23 \begin{document}$ \begin{align} & {W_k} = (1 - \lambda ){W_{k - 1}} + \lambda {Y_k}, \end{align} $\end{document}
24 \begin{document}$ \begin{align} & {Q_k} = (1 - \lambda ){Q_{k - 1}} + \lambda {Z_k}, \end{align} $\end{document}
其中\begin{document}$ {W_k} $\end{document} 和\begin{document}$ {Q_k} $\end{document} 的初值设置为\begin{document}$ {W_0} = {Q_0} = 0 $\end{document} .
step4: SS-DEWMA图定义如下:
25 \begin{document}$ \begin{align} {L_k} = W_k^2 + Q_k^2. \end{align} $\end{document}
step5: SS-DEWMA图的上界UCL定义为
26 \begin{document}$ \begin{align} {\rm UCL} = E( {\{ {{L_k}} \}} ) + {K_{\rm SD}}\sqrt {\sigma _{{L_k}}^2} . \end{align} $\end{document}
其中: \begin{document}$ {K_{\rm SD}} $\end{document} 为UCL的控制宽度, \begin{document}$ {K_{{\rm{SD}}}} > 0 $\end{document} ; \begin{document}$ E( {\{ {{L_k}} \}} ) $\end{document} 和\begin{document}$ \sigma _{{L_k}}^2 $\end{document} 分别为序列\begin{document}$ {\{ {{L_k}} \}} $\end{document} 的均值和方差.
为减少噪声对实际数据的影响, 利用小波分析方法[14 ] 对SS-DEWMA输出序列\begin{document}$ {\{ {{L_k}} \}} $\end{document} 进行多尺度分解、去噪和重构, 得到新的序列\begin{document}$ {\{ {{l_k}} \}} $\end{document} (参见3.2节). 相应地, 式(26)重写为
27 \begin{document}$ \begin{align} {\rm UCL} = E( {\{ {{l_k}} \}} ) + {K_{\rm SD}}\sqrt {\sigma _{{l_k}}^2}. \end{align} $\end{document}
在上述SS-DEWMA图算法中, \begin{document}$ {U_k} $\end{document} 和\begin{document}$ {V_k} $\end{document} 是两个分别表征过程均值和方差的独立统计量, 当\begin{document}$ a = 0 $\end{document} , \begin{document}$ b = 1 $\end{document} 时, \begin{document}$ \mu $\end{document} 和\begin{document}$ \sigma^2 $\end{document} 分别为样本的均值和方差. 由式(19)和(20)可知, \begin{document}$ {U_k} $\end{document} 和\begin{document}$ {V_k} $\end{document} 都服从标准正态分布, 即\begin{document}$ {U_k} \sim N(0, 1) $\end{document} , \begin{document}$ {V_k} \sim N(0, 1) $\end{document} , 并且, 其分布都不依赖于样本序列维数\begin{document}$ {n_y} $\end{document} , 此时, 该过程在统计上处于可控状态; 而\begin{document}$ a \ne 0 $\end{document} 或\begin{document}$ b \ne 1 $\end{document} 的情况分别表示样本的均值或方差发生改变. 这样, \begin{document}$ {U_k} $\end{document} 和\begin{document}$ {V_k} $\end{document} 构造成一个单一统计图来监控过程均值和方差的变化.
对式(21)进行递归, 并考虑到\begin{document}$ {Y_0}=0 $\end{document} , 可得
28 \begin{document}$ \begin{align} {Y_k} =\lambda \sum\limits_{j = 1}^k {{{(1 - \lambda )}^{k - j}}} {U_j}. \end{align} $\end{document}
由于\begin{document}$ {U_k} \sim N(0, 1) $\end{document} , 有\begin{document}$ {Y_k} \sim N( {0, \sigma _{{Y_k}}^2} ) $\end{document} , 其中, 方差为
29 \begin{document}$ \begin{align} \sigma _{{Y_k}}^2 = {{\lambda ( {1 - {{( {1 - \lambda } )}^{2k}}} )} / {( {2 - \lambda } )}}. \end{align} $\end{document}
同理可得\begin{document}$ {Z_k} \sim N( {0, \sigma _{{Z_k}}^2} ) $\end{document} , \begin{document}$ \sigma _{{Z_k}}^2 = \sigma _{{Y_k}}^2 $\end{document} . 同样, 对式(23)进行递归, 并考虑到\begin{document}$ {W_0} = 0 $\end{document} , 可得
30 \begin{document}$ \begin{align} {W_k} = {\lambda ^2}\sum\limits_{j = 1}^k {( {k - j + 1} ){{( {1 - \lambda } )}^{k - j}}{U_j}} . \end{align} $\end{document}
因此有\begin{document}$ {W_{{k}}} \sim N(0, \sigma _{{W_k}}^2) $\end{document} , 其中, 方差为
31 \begin{document}$ \begin{align} \sigma _{{W_k}}^2 = {\lambda ^4}\sum\limits_{j = 1}^k {{{( {k - j + 1} )}^2}{{( {1 - \lambda } )}^{2( {k - j} )}}} . \end{align} $\end{document}
同理可得\begin{document}$ {Q_k} \sim N(0, \sigma _{{Q_k}}^2) $\end{document} , \begin{document}$ \sigma _{{Q_k}}^2 = \sigma _{{W_k}}^2 $\end{document} .
2.2
SS-DEWMA图的多目标优化
2.2.1
构造FAR和MDR函数并分析其影响规律
分别建立MDR和FAR与SS-DEWMA图的两个参数的计算公式, 如下:
32 \begin{document}$ \begin{align} &{\rm FAR} = {\rm prob}\{ {{l_k} > {\rm UCL}|{g_k} = 0} \}, \end{align} $\end{document}
33 \begin{document}$ \begin{align} &{\rm MDR} = 1 - {\rm FDR}=1 - {\rm prob}\{ {{l_k} > {\rm UCL}|{g_k} \ne 0} \}, \end{align} $\end{document}
34 \begin{document}$ \begin{align} &{\rm FDR} = {\rm prob}\{ {{l_k} > {\rm UCL}|{g_k} \ne 0} \}. \end{align} $\end{document}
其中: FDR是故障检测率; \begin{document}$ {g_k} = 0 $\end{document} 表示传感器无故障, \begin{document}$ {g_k} \ne 0 $\end{document} 表示传感器出现故障; prob\begin{document}$ \{ * \} $\end{document} 表示序列\begin{document}$ \{ {{l_k}} \} $\end{document} 中符合条件\begin{document}$ * $\end{document} 的样本个数占序列\begin{document}$ \{ {{l_k}} \} $\end{document} 总样本个数的百分比. 由式(32) \begin{document}$ \sim $\end{document} (34)可知, MDR和FAR值的大小由序列\begin{document}$ \{ {{l_k}}\} $\end{document} 和UCL共同决定.
由式(26)可知, UCL的值与序列\begin{document}$ \{ {{L_k}} \} $\end{document} 的均值\begin{document}$ E( {\{ {{L_k}} \}} ) $\end{document} 、方差\begin{document}$ \sigma _{{L_k}}^2 $\end{document} 以及控制宽度\begin{document}$ {K_{\rm SD}} $\end{document} 有关. 根据式(25), 并考虑\begin{document}$ \sigma _{{Q_k}}^2 = \sigma _{{W_k}}^2 $\end{document} , 有\begin{document}$ \dfrac{{{L_k}}}{{\sigma _{{W_k}}^2}} = \dfrac{{W_k^2}}{{\sigma _{{W_k}}^2}} + \dfrac{{Q_k^2}}{{\sigma _{{Q_k}}^2}} $\end{document} , 则\begin{document}$ \dfrac{{{L_k}}}{{\sigma _{{W_k}}^2}} $\end{document} 服从自由度为2的卡方分布, 即\begin{document}$ \dfrac{{{L_k}}}{{\sigma _{{W_k}}^2}} \sim \chi _2^2 $\end{document} . 进而有\begin{document}$ E( {\{ {{L_k}} \}} ) = 2\sigma _{{W_k}}^2 $\end{document} , \begin{document}$ \sigma _{{L_k}}^2 = 4\sigma _{{W_k}}^4 $\end{document} . 因此, 式(26)可写为
35 \begin{document}$ \begin{align} {\rm UCL} = 2\sigma _{{W_k}}^2( {1 + {K_{\rm SD}}} ). \end{align} $\end{document}
当\begin{document}$ \lambda $\end{document} 增大时, 由式(21) \begin{document}$ \sim $\end{document} (25)可知序列\begin{document}$ \{ {{L_k}} \} $\end{document} 实时性变强, 但同时受随机干扰程度也会增大, 平稳性下降. 由式(31)可知, 当\begin{document}$ \sigma _{{W_k}}^2 $\end{document} 增大时, 在控制宽度\begin{document}$ {K_{\rm SD}} $\end{document} 不变的情况下, 由式(35)可知UCL将会增大; 相反, 当\begin{document}$ \lambda $\end{document} 减小时, 不能及时反映当前数据信息, 但平稳性增加, \begin{document}$ \sigma _{{W_k}}^2 $\end{document} 减小, 上界UCL将会减小, 但并不能直接分析出对故障检测误报率FAR和漏报率MDR的影响. 当平滑参数\begin{document}$ \lambda $\end{document} 不变时, 若增大控制宽度\begin{document}$ {K_{\rm SD}} $\end{document} , 则上界UCL将会增大, 对故障检测的FAR将会降低, 但可能增大MDR; 相反, 若减小控制宽度\begin{document}$ {K_{\rm SD}} $\end{document} , 则会导致上界UCL减小, 将会减小MDR, 但FAR可能增大.
为更直观地说明SS-DEWMA图的两个参数对MDR和FAR的影响, 将\begin{document}$ {K_{{\rm{SD}}}} \in({0, 10}) $\end{document} 和\begin{document}$ \lambda \in ( {0, 1} ) $\end{document} 都分为20等分进行仿真. 在每组\begin{document}$ ( {\lambda , {K_{\rm SD}}} ) $\end{document} 参数下, 采用SS-DEWMA图方法对存在输出传感器偏移故障的数据样本进行故障检测, 并逐一计算出各组参数取值情况下的MDR和FAR, 绘制成三维曲面图, 如图 1 所示.
1
不同参数下的MDR和FAR三维图
由图 1 可见, \begin{document}$ \lambda $\end{document} 和\begin{document}$ {K_{\rm SD}} $\end{document} 的选取影响MDR和FAR的大小, 并且其影响规律较为复杂, 很难确定选取参数\begin{document}$ \lambda $\end{document} 和\begin{document}$ {K_{\rm SD}} $\end{document} 的值可使MDR和FAR同时达到最小.
2.2.2
SS-DEWMA图参数的离线优化
本节利用多目标算法可以得到一组帕累托(Pareto)最优解使得多个目标函数同时满足要求的特性, 在SS-DEWMA图中引入MO-PSO[15 ] 算法, 并以MDR和FAR为两个目标函数, 用来解决MDR和FAR难以同时达到最小的问题.
SS-DEWMA图通过MO-PSO算法在\begin{document}$ ( {\lambda, {K_{\rm SD}}} ) $\end{document} 的取值范围内随机生成的粒子位置信息表示参数的不同取值. 再使用帕累托支配原则对当前粒子的MDR和FAR进行评价, 将粒子中的非劣解拷贝到Archive集中. 然后, 对Archive集不断进行更新, 直至Archive集已满或MO-PSO算法达到最大迭代次数, 再从最终输出的Archive集中选出在MDR与FAR之间进行权衡的一组解\begin{document}$ ( {\hat \lambda , {{\hat K}_{\rm SD}}} ) $\end{document} , 从而避免仅凭经验选取参数使MDR/FAR较大的问题, 即
36 \begin{document}$ \begin{align} &( {\hat \lambda , {{\hat K}_{\rm SD}}} ) =\\ & \mathop {\arg \min }\limits_{( {\lambda , {K_{\rm SD}}} ) \in {R^n}} \{ {{\rm FAR}( {\lambda , {K_{\rm SD}}} ), } {\rm MDR}( {\lambda , {K_{\rm SD}}} ) \}. \end{align} $\end{document}
MO-PSO算法对SS-DEWMA图参数进行优化的具体步骤如下:
step 1:初始化粒子. 在参数\begin{document}$ ( {\lambda , {K_{\rm SD}}} ) $\end{document} 的取值范围内随机生成每个粒子的位置\begin{document}$ {\rm Pop}_0^z $\end{document} , 令速度\begin{document}$ v_0^z = 0 $\end{document} . 其中: 粒子数\begin{document}$ z = 1, 2, \ldots , N $\end{document} , \begin{document}$ N $\end{document} 为粒子总数; 下角标\begin{document}$ 0 $\end{document} 表示初值.令个体极值Pbest\begin{document}$ _0^z = {\rm Pop}_0^z $\end{document} .
step 2:计算FAR\begin{document}$ _p^z $\end{document} 和MDR\begin{document}$ _p^z $\end{document} . 将Pop\begin{document}$ _p^z $\end{document} 代入SS-DEWMA图求得\begin{document}$ \{ {l_k^z} \} $\end{document} 和UCL\begin{document}$ ^z $\end{document} , 并据式(32)和(33)计算FAR\begin{document}$ _p^z $\end{document} 和MDR\begin{document}$ _p^z $\end{document} . 下角标\begin{document}$ p $\end{document} 表示第\begin{document}$ p $\end{document} 次迭代.
step 3:生成/更新Archive集. 根据帕累托支配原则对每个粒子所对应的FAR\begin{document}$ _p^z $\end{document} 和MDR\begin{document}$ _p^z $\end{document} 进行评价, 若Archive集为空集, 则将非劣解拷贝到Archive集中; 若Archive集非空集, 只要该粒子优于或独立于原Archive集中的某个粒子, 则将该粒子插入Archive集.
step 4:判断Archive集. 如果Archive集中的粒子数超出了规定的大小, 即Archive集已满, 则对Archive集进行如下截断操作:
37 \begin{document}$ \begin{align} {\rm PN} = {\rm Int}\Big( {\dfrac{{| {{A_p}} | - {Z_{\rm ar}}}}{{| {{A_p}} |}} \times {\rm Grid}[q, 2] + 0.5} \Big). \end{align} $\end{document}
其中: \begin{document}$ | {{A_p}} | $\end{document} 代表Archive集包含的粒子数, \begin{document}$ {Z_{\rm ar}} $\end{document} 为Archive集大小, Grid\begin{document}$ [q, 2] $\end{document} 表示2维网格\begin{document}$ q $\end{document} 中包含的粒子数.
若Archive集未满, 则执行step 5;否则跳至step 7.
step 5:计算全局极值Gbest\begin{document}$ _i^z $\end{document} . 运用网格法计算Archive集的拥挤度, 在Archive集中选取拥挤度低的粒子\begin{document}$ Z $\end{document} 作为Gbest\begin{document}$ _i^z $\end{document} .
step 6:更新粒子位置和速度. 有
38 \begin{document}$ \begin{align} &{\rm Pop}_{p + 1}^z = {\rm Pop}_p^z + v_p^z, \end{align} $\end{document}
39 \begin{document}$ \begin{align} &v_{p + 1}^z = w \times v_p^z + {c_1} \times {r_1}( {{\rm Pbest}_p^z - {\rm Pop}_p^z} )+\\ &\quad \quad \quad \; {c_2} \times {r_2}( {{\rm Gbest}_p^z - {\rm Pop}_p^z} ). \end{align} $\end{document}
其中: \begin{document}$ w = \dfrac{2}{{\phi - 2 + \sqrt {{\phi ^2} - 4 \phi } }} $\end{document} 为惯性权重, \begin{document}$ \phi = {\phi _1} + {\phi _2} $\end{document} , \begin{document}$ \phi > 4 $\end{document} , 通常取\begin{document}$ \phi = 4.1 $\end{document} , \begin{document}$ {\phi _1} = {\phi _2} = 2.05 $\end{document} , \begin{document}$ {c_1} = w{\phi _1} $\end{document} , \begin{document}$ {c_2} = w{\phi _2} $\end{document} , \begin{document}$ {r_1} $\end{document} 和\begin{document}$ {r_2} $\end{document} 是0 \begin{document}$ \sim $\end{document} 1之间的随机数.
令\begin{document}$ p = p + 1 $\end{document} , 若\begin{document}$ i < M $\end{document} , 则跳至step 2, 其中\begin{document}$ M $\end{document} 为最大迭代次数; 否则跳至step 7.
step 7:输出Archive集, 并在Archive集中选出同时使MDR和FAR最小的一组参数作为最优参数\begin{document}$ ( {\hat \lambda , {{\hat K}_{\rm SD}}} ) $\end{document} .
上述算法中, 式(39)借鉴了文献[16 ]提出的带收缩因子的PSO, 以避免陷入局部最优问题. 另外, 如果Pop\begin{document}$ _{p{{ + }}1}^z $\end{document} 越过了边界范围, 即\begin{document}$ {\lambda ^z} $\end{document} 和/或\begin{document}$ K^z_{\rm SD} $\end{document} 超出了其取值范围, 则选取其相应的边界值, 再将\begin{document}$ v_{p{\rm{ + }}1}^z $\end{document} 符号取反, 重新计算式(38).
2.3
基于MO-SS-DEWMA图的在线故障检测
在利用CESCKF得到残差序列的基础上, 使用MO-SS-DEWMA图实现传感器故障在线检测的流程如图 2 所示.
2
基于MO-SS-DEWMA图的在线故障检测流程
仍使用图 1 仿真中的数据样本, 采用MO-SS-DEWMA图算法进行仿真, 得到的参数优化后的MDR和FAR三维图类似于图 1 , 但平滑参数和控制宽度的坐标取值范围变小, 分别为\begin{document}$ [ {0.25, 0.45} ] $\end{document} 和\begin{document}$ [ {3.5, 5} ] $\end{document} , 此时MDR和FAR在0附近, 小于1\begin{document}$ \% $\end{document} . 因此, 在此范围内取值可以较精确地检测出系统是否存在故障. 该算法并不仅限于用在输出传感器故障检测, 还适用于其他类型的故障检测. 本文仅对其在输出传感器故障检测方面的应用进行研究.
3
基于CESCKF的MO-SS-DEWMA图方法在传感器故障检测中的应用
本节以伺服电机驱动的连铸结晶器振动系统为对象进行故障检测仿真和实验研究. 在仿真中考虑了位移传感器发生偏移故障和方差故障的情况. 同时, 在连铸结晶器模拟振动台上, 对容易发生偏移故障的位移传感器进行实验验证.
3.1
系统模型及滤波估计
基于文献[17 ]给出的伺服电机驱动的连铸结晶器振动位移跟踪系统为对象进行仿真. 其数学模型如下:
40 \begin{document}$ \begin{align} &{{\dot x}_{{p}}} = 3\Big(\dfrac{{2π }}{{60{i_f}}}n \Big)\cos \theta, \\[7pt] &\dot \theta = \dfrac{{2π }}{{60{i_f}}}n, \\[7pt] &\dot n = \dfrac{1.5p\psi _f}{J}\dfrac{60}{2π }i_q- \dfrac{B}{J}n - \dfrac{60}{2π }\dfrac{T_L}{J}, \\[7pt] &{{\dot i}_{{q}}} = - \dfrac{{2π }}{{60}}pn{i_{{d}}} - \dfrac{{{R_{{s}}}}}{L}{i_{{q}}} - \dfrac{{p{\psi _{{f}}}}}{L}\dfrac{{2π }}{{60}}n + \dfrac{{{u_{{q}}}}}{L}, \\[7pt] &{{\dot i}_{{d}}} = - \dfrac{{{R_{{s}}}}}{L}{i_{{d}}} + \dfrac{{2π }}{{60}}pn{i_{{q}}} + \dfrac{{{u_{{d}}}}}{L}. \end{align} $\end{document}
其中: \begin{document}$ {x_p} $\end{document} 为结晶器的振动位移; \begin{document}$ \theta $\end{document} 为偏心轴转角; \begin{document}$ n $\end{document} 为伺服电机的转速; \begin{document}$ i $\end{document} 为减速器的减速比; 伺服电机为永磁同步电动机, 采用磁场定向矢量控制策略, \begin{document}$ {u_d} $\end{document} 、\begin{document}$ {u_q} $\end{document} 分别为定子电压\begin{document}$ d $\end{document} 、\begin{document}$ q $\end{document} 轴分量; \begin{document}$ {i_d} $\end{document} 、\begin{document}$ {i_q} $\end{document} 为定子电流\begin{document}$ d $\end{document} 、\begin{document}$ q $\end{document} 轴分量; \begin{document}$ L $\end{document} 为定子绕组等效电感; \begin{document}$ {R_s} $\end{document} 为定子等效电阻; \begin{document}$ p $\end{document} 为磁极对数; \begin{document}$ {\psi _{{f}}} $\end{document} 为转子永磁体磁链; \begin{document}$ J $\end{document} 为转子转动惯量; \begin{document}$ B $\end{document} 为粘性摩擦系数; \begin{document}$ {T_L} $\end{document} 为负载转矩. 在仿真中, 控制器设计采用文献[17 ]的控制方法, 系统参数如表 1 所示.
1
连铸结晶器振动系统参数
参数名称/单位
取值
参数名称/单位
取值
额定功率\begin{document}$ {P_{{N}}} $\end{document} /kW
20.4
额定转速\begin{document}$ {n_{{N}}} $\end{document} /(r /min)
1 500
额定电流\begin{document}$ {I_{{N}}} $\end{document} /A
45
定子电阻\begin{document}$ {R_{{s}}} $\end{document} /\begin{document}$ \Omega $\end{document}
0.14
定子等效电感\begin{document}$ L $\end{document} /mH
4.6
转子转动惯量\begin{document}$ J $\end{document} /(kg\begin{document}$ \cdot{\rm m}^2 $\end{document} )
0.054 7
转子永磁体磁链\begin{document}$ {\psi _{{f}}} $\end{document} /Wb
0.96
磁极对数\begin{document}$ p $\end{document}
3
粘性摩擦系数\begin{document}$ B $\end{document}
0.004
减速器减速比\begin{document}$ {{i_{{r}}}} $\end{document}
5.1
振动频率\begin{document}$ f $\end{document} /min
120
采样周期\begin{document}$ T $\end{document} /s
0.001
将\begin{document}$ {i_q} $\end{document} 作为控制输入, \begin{document}$ {x_p} $\end{document} 、\begin{document}$ \theta $\end{document} 和\begin{document}$ n $\end{document} 作为测量输出. 对式(40)的前3个状态方程, 利用欧拉离散化方法进行处理, 得到如下方程, 用于CESCKF状态估计:
41 \begin{document}$ \begin{align} \begin{cases} {x_p}_{, k} = {x_{p, }}_{k - 1} + T\Big( {3\Big(\dfrac{{2{{π }}}}{{60{i_{{r}}}}}{n_{k - 1}}\Big)\cos ( {{\theta _{k - 1}}} )}\Big), \\[9pt] {\theta _k} = {\theta _{k - 1}} + T\Big( {\dfrac{{2{{π }}}}{{60i}}{n_{k - 1}}} \Big), \\[9pt] {n_k} = {n_{k - 1}} + T\Big( {\dfrac{{1.5{p_n}{\psi _{{f}}}}}{J}\dfrac{{60}}{{2{{π }}}}( {{i_{{{q, }}}}_{k - 1}} )}-\\[9pt] \quad \quad\; \dfrac{B}{J}n_{k - 1} - \dfrac{60}{2π }\dfrac{T_L}{J} \Big), \\[7pt] {y_k} = [{{x_p}_{, k}}\; \; {{\theta _k}}\; \; {{n_k}}]^{\rm T}. \end{cases} \end{align} $\end{document}
其中: \begin{document}$ k $\end{document} 为采样步数, \begin{document}$ T $\end{document} 为采样周期, \begin{document}$ {y_k} $\end{document} 为传感器测量输出, 令状态\begin{document}$ {x_k} = {[ {{x_p}_{, k}\; \; {\theta _k}\; \; {n_k}} ]^{\rm{T}}} $\end{document} . 初始化参数: 采样周期\begin{document}$ T = 0.001 $\end{document} s, 过程白噪声\begin{document}$ {\upsilon _{k{{ - }}1}} $\end{document} 满足\begin{document}$ {\upsilon _k}\sim N( 0, {\varXi _{\upsilon , k - 1}} ) $\end{document} , \begin{document}$ {\varXi _{\upsilon , k - 1}} = {\rm diag}[ {0.01}\; \; {0.01}\; \; {0.01}] $\end{document} , 测量白噪声\begin{document}$ {\omega _k} $\end{document} 满足\begin{document}$ {\omega _k}\sim N( {0, {\varXi _{\omega , k}}} ) $\end{document} , \begin{document}$ {\varXi _{\omega , k}} = {\rm diag}[ {0.01}, \; {0.01}, \; {0.01}] $\end{document} . 状态初始值\begin{document}$ {x_0} = [ 0\; \; 0\; \; 0 ]^{\rm T} $\end{document} , 状态初始估计\begin{document}$ {\hat x_0} = [0\; \; 0\; \; 0]^{\rm T} $\end{document} , 状态协方差矩阵\begin{document}$ {P_0} = {\rm diag}[1\; \; 1\; \; 1] $\end{document} .
3.2
小波分解与多目标优化参数设置
采用小波分析方法对SS-DEWMA输出序列的多尺度分解过程为: 将SS-DEWMA输出序列\begin{document}$ \{ {{L_k}} \} $\end{document} 作为小波分析的原始信号, 通过一个低通滤波器卷积原始信号可以得到一个近似信号, 原始信号与近似信号之间的差值称为细节信号, 可以通过高通滤波器对原始信号进行卷积得到, 原始信号可以表示为最后一个近似信号与所有细节信号之和.
用Matlab实现过程为: 首先, 利用\begin{document}$ [C, L] = {\rm wavedec}( {\{ {{L_k}} \}, 3, '{\rm db}3'} ) $\end{document} 函数对\begin{document}$ \{ {{L_k}} \} $\end{document} 进行3尺度分解. 其中: \begin{document}$ C $\end{document} 由近似小波系数和各个尺度下的细节小波系数组成, \begin{document}$ L $\end{document} 为各小波系数的长度, db3表示对\begin{document}$ \{ {{L_k}} \} $\end{document} 分解所使用的小波. 其次, 分别利用\begin{document}$ A = {\rm appcoef}( {C, L, '{\rm db}3'} ) $\end{document} 和\begin{document}$ D = {\rm appcoef}( {C, L} ) $\end{document} 提取近似小波系数和细节小波系数. 再次, 使用\begin{document}$ {\rm thr} = {\rm thselect}( {D, '{\rm rigrsure}'} ) $\end{document} 函数利用Stein的无偏风险估计原理得到阈值. 最后, 使用wdencmp函数进行去噪和重构得到\begin{document}$ \{ {{l_k}} \} $\end{document} .
仿真中, MO-PSO算法的参数选取如表 2 所示.
2
MO-PSO参数设置
名称
取值
名称
取值
粒子总数\begin{document}$ N $\end{document}
100
迭代次数\begin{document}$ M $\end{document}
100
Archive集大小\begin{document}$ {Z_{\rm ar}} $\end{document} /A
100
惯性权重\begin{document}$ w $\end{document}
0.729 8
个体学习系数\begin{document}$ {c_1} $\end{document}
1.496 2
全局学习系数\begin{document}$ {c_2} $\end{document}
1.496 2
\begin{document}$ \lambda $\end{document} 的取值范围
\begin{document}$ ( {0, 1} ) $\end{document}
\begin{document}$ {{K_{{\rm{SD}}}}} $\end{document} 的取值范围
\begin{document}$ ( {0, 5} ) $\end{document}
3.3
传感器故障检测与对比仿真
MO-SS-DEWMA图监测数据时需要已知过程的均值和方差, 通常需要利用大量的数据样本得到均值和方差的近似值, 样本量越大, 越接近真实值; 同时, 数据监测需要满足及时性, 当监测的数据量较少时, 能够提高实际应用中的灵活性和多样性. 兼顾上述两方面要求, 本文选取1 000个残差向量样本数据作为MO-SS-DEWMA图要监控的数据量, 并假设故障发生在样本的[501, 1 000]采样处.
下面分别对位移传感器偏移故障、方差故障以及同时存在两种故障的情况进行仿真, 并与M-DEWMA图[12 ] 、多尺度优化指数加权移动平均(MS-OEWMA)图[5 ] 和SS-DEWMA图[13 ] 进行比较, 验证MO-SS-DEWMA图能否提高故障的检测精度. 仿真中, 4种检测方法的样本大小取值与故障参数均相同.
1) 偏移故障的情况.
由于连铸结晶器振动系统工作环境恶劣, 存在高温高粉尘等影响, 位移传感器容易发生偏移故障.
仿真中, 样本均方差\begin{document}$ \sigma = 0.04 $\end{document} , 首先, 假设偏移故障的偏移量为0.08 mm, 即\begin{document}$ a = 2, b = 1 $\end{document} , 仿真对比结果如表 3 所示.
3
传感器故障下的残差评价结果
残差评价方法
偏移量为2\begin{document}$ \sigma $\end{document} , 即0.08 mm
均方差为2\begin{document}$ \sigma $\end{document}
偏移量为2\begin{document}$ \sigma $\end{document} , 均方差为2\begin{document}$ \sigma $\end{document}
MDR\begin{document}$ /\% $\end{document}
FAR\begin{document}$ /\% $\end{document}
\begin{document}$ \lambda $\end{document}
\begin{document}$ {K_{\rm SD}} $\end{document}
MDR\begin{document}$ /\% $\end{document}
FAR\begin{document}$ /\% $\end{document}
\begin{document}$ \lambda $\end{document}
\begin{document}$ {K_{\rm SD}} $\end{document}
MDR \begin{document}$ /\% $\end{document}
FAR \begin{document}$ /\% $\end{document}
\begin{document}$ \lambda $\end{document}
\begin{document}$ {K_{\rm SD}} $\end{document}
M-DEWMA图
1.4
0
0.4
1.5
15
18.4
0.3
0.2
17
17.2
0.3
0.3
MS-OEWMA图
0
2.2
0.338 8
4.862 1
19.8
0
0.338 8
4.862 1
15.2
6.2
0.338 8
4.862 1
SS-DEWMA图
1.2
0
0.4
3
14
9.6
0.3
0.15
14.8
8.8
0.3
0.2
MO-SS-DEWMA图
0
0.2
0.324 7
3.470 3
8.4
5.2
0.321 4
0.141 5
7.8
6.2
0.303 0
0.215 8
由表 3 可见, 当故障的偏移量为0.08 mm时, 相比于现有方法, 本文所提出的MO-SS-DEWMA图方法都能准确地检测出故障信息, 其MDR为0, FAR值为0.2 \begin{document}$ \% $\end{document} . 若以(1-FAR-MDR)为检测精度, 则对于偏移量为2倍及以上样本均方差的偏移故障, MO-SS-DEWMA图方法的检测精度为99.8 \begin{document}$ \% $\end{document} , 说明本文所提出方法能够有效检测传感器偏移故障.
2) 方差故障的情况.
由于存在噪声干扰和电磁干扰等, 使得输出位移传感器测量值包含较大干扰信号, 从而导致测量不准确, 即发生位移传感器方差故障.
仿真中, 样本均方差仍为\begin{document}$ \sigma = 0.04 $\end{document} , 假设发生故障时的测量输出均方差为无故障时均方差的2倍(2\begin{document}$ \sigma $\end{document} ), 即\begin{document}$ a = 0, b = 2 $\end{document} , 其仿真对比结果如表 3 中间部分所示.
由表 3 可见, 相比于偏移故障, 对方差故障的准确检测相对较难, 漏报率都偏高. 相对而言, MO-SS-DEWMA图算法的检测精度为86.4 \begin{document}$ \% $\end{document} , 其他方法的漏报率均超过10 \begin{document}$ \% $\end{document} , 难以在实际中应用.
3) 偏移故障与方差故障同时出现的情况.
上述两组仿真是针对位移传感器偏移故障和方差故障单独进行的, 在本组仿真中, 考虑发生2\begin{document}$ \sigma $\end{document} 的偏移故障和均方差为2\begin{document}$ \sigma $\end{document} 的方差故障, 即\begin{document}$ a = 2, b = 2 $\end{document} , 其仿真对比结果如表 3 右侧部分所示.
由表 3 可见, 相对而言, MO-SS-DEWMA图算法的检测精度为86 \begin{document}$ \% $\end{document} , 其他方法的检测精度均低于80 \begin{document}$ \% $\end{document} .
由于在实际应用中无法预判发生偏移故障或方差故障, 采用同时存在偏移和方差故障的数据进行参数优化, 得到MO-SS-DEWMA图; 然后, 分别检测只有偏移故障、只有方差故障以及偏移故障和方差故障同时存在的情况. 仿真结果如表 4 所示.
4
MO-SS-DEWMA图对不同故障的评价结果
故障类型
参数值(0.303 0, 0.215 8)
MDR \begin{document}$ /\% $\end{document}
FAR \begin{document}$ /\% $\end{document}
只有偏移故障(\begin{document}$ a = 2, b = 1 $\end{document} )
0
1.4
只有方差故障(\begin{document}$ a = 0, b = 2 $\end{document} )
8.8
5.4
偏移故障和方差故障(\begin{document}$ a = 2, b = 2 $\end{document} )
7.8
6.2
从表 4 可以看出, 利用同时存在偏移和方差故障的数据得到的MO-SS-DEWMA图对于只有偏移故障的检测精度为98.6 \begin{document}$ \% $\end{document} , 对于只有方差故障的检测精度为85.8 \begin{document}$ \% $\end{document} , 与表 3 中的检测精度相近.
综上所述, 结晶器振动位移系统实际应用中, 可采用同时存在偏移和方差故障数据进行参数优化, 将得到的MO-SS-DEWMA图用来检测输出位移传感器是否存在故障.
3.4
传感器故障检测与对比实验
本节所采用的实验装置是基于西门子高性能运动控制器SIMOTION D425的连铸结晶器模拟振动台.在输出位移传感器数据[501, 1 000]中加入0.08 mm模拟传感器偏移故障. 通过采集各1 000组的系统输出和控制输入数据, 采用CESCKF获得残差数据, 之后, 使用与3.3节中的3)相同参数的MO-SS-DEWMA图和与上述相同的对比方法对残差数据进行检测, 结果如表 5 所示.
5
传感器偏移故障残差在模拟振动台上的评价结果
残差评价方法
MDR \begin{document}$ /\% $\end{document}
FAR \begin{document}$ /\% $\end{document}
\begin{document}$ \lambda $\end{document}
\begin{document}$ {K_{\rm SD}} $\end{document}
M-DEWMA图
0.4
2.2
0.3
0.3
MS-OEWMA图
0
1
0.338 8
4.862 1
SS-DEWMA图
0.6
0
0.3
0.2
MO-SS-DEWMA图
0.2
0.6
0.303 0
0.215 8
由表 5 可见, 当在模拟平台上对位移传感器偏移故障进行检测时, 相较于其他评价方法, 在优先考虑MDR的情况下, 采用同时存在偏移和方差故障数据得到的MO-SS-DEWMA图方法的检测精度最高, 该方法适用于对实际系统传感器的故障检测.