CN102564486A - 一种传感器慢偏故障的校正方法 - Google Patents

Abstract

本发明提供一种传感器慢偏故障的校正方法，包含步骤a：对传感器的测量值建立状态方程和观测方程；步骤b：对所述状态方程和观测方程进行时间上的扩展；步骤c：利用扩展状态的卡尔曼滤波估计偏差，再将偏差值从测量值中减去。本发明能够实时估计出传感器的慢变偏差，从而能在传感器出现缓慢变化的偏差时，将该偏差估计出来，从传感器测量值中减去。

Description

一种传感器慢偏故障的校正方法

技术领域

本发明属于传感器故障诊断及校正技术领域，涉及一种可用于校正传感器慢偏故障的方法。

背景技术

传感器是一种检测装置，它能检测到被测量的信息并将其按照一定的规律变换成为相应的电信号或其他所需形式的信息输出，以满足信息的传输、处理、存储、显示、记录和控制等要求，是实现自动检测和自动控制的首要环节。随着信息处理技术以及微处理机和计算机技术的高速发展，传感器已广泛应用于航空航天、军事、医疗卫生、工农业生产以及环境保护监测等各个领域。尤其是近几十年来，随着材料科学领域的突破性进展，出现了越来越多新型的传感器，如光纤、超声波和微机械传感器等等。

作为信息采集系统的前端单元，传感器的准确性和稳定性是整个系统良好运行的先决条件。对于一些复杂的大型动态系统，如飞机，电站等，需要依赖准确可靠的传感器去获得测量信号来判断和监视系统状态，以进行有效控制。可以说，传感器的准确性和可靠性极大程度的决定了整个系统的性能、可靠性和安全性。

传感器的测量精度等信息通常由传感器生产厂家在产品出厂时给定或用户使用前通过其他方法标定。这两种标定精度通常是在一定条件下获得的，而实际应用中环境条件复杂多变，很难始终保持其与标定时的条件一致，当传感器的工作条件与它的标定条件不一致时，传感器精度往往不能真实反映传感器测量值的实际准确程度。很多传感器的工作环境非常恶劣，例如风力发电机状态监测系统中的传感器就工作在低温、大风或者高盐雾条件下，因此非常容易出现各种各样的故障。如果故障传感器的输出数据一直作为控制系统的输入时，可能会导致整个系统运行性能严重下降，甚至出现灾难性的事故。

为了保证测量的精度和可靠性，传统的方法主要是通过一个较高精度的传感器测量值作为基准，用来校正其余传感器测量值。这种方法的缺点在于，首先精度较高的传感器通常价格昂贵，使得成本大大增加；其次，即便在不需顾及成本的应用场合中，复杂的环境因素往往也会导致高精度的传感器精度下降，这样其基准作用也就无从谈起。另外一种实际中极为常用的方法是对系统中的传感器定期进行校准和维护，在定期校准带来大量的人力物力损耗的同时，人们还是无从得知两次校准之间传感器的输出是否准确，而校准周期的确定也是另外一个难题。因此怎样能够准确有效的在传感器出现测量偏差时，对这个偏差进行校正，具有十分重要的实用价值。

发明内容

本发明的目的在于克服现有技术的不足，提供一种传感器慢偏故障的校正方法，该方法能够实时估计出传感器的慢变偏差，从而能在传感器出现缓慢变化的偏差时，将该偏差估计出来，从传感器测量值中减去。这样就实现了测量偏差和测量真值的分离。

为达到上述目的，本发明采用如下技术方案：

一种传感器慢偏故障的校正方法，包含如下步骤：

步骤a：对传感器的测量值建立状态方程和观测方程；

步骤b：对所述状态方程和观测方程进行时间上的扩展；

步骤c：利用扩展状态的卡尔曼滤波估计偏差，再将偏差值从测量值中减去。

进一步地，所述步骤a可以按如下步骤实施：

步骤a1：将含有传感器测量偏差的线性动态系统的状态空间模型描述为：

x_t+1＝Fx_t+w_t    (1)

y_t＝Φx_t+b_t+v_t    (2)

b_t+1＝Ab_t+η_t    (3)

其中，下标t表示时间，x_t+1为N×1维的状态向量，F为状态转移矩阵，y_t为p×1维的观测向量，Φ为测量矩阵，b_t为M×1维的传感器测量偏差向量，A为b_t的状态转移矩阵，w_t为状态噪声，v_t为测量噪声，η_t为偏差状态噪声，且状态噪声w_t，测量噪声v_t及偏差状态噪声η_t均为相互独立的零均值的高斯白噪声，且满足：

E(w_t)＝0

$E\left(w_j{w}_t^T\right)=R{\mathrm{\δ}}_{\mathrm{jt}}---\left(4\right)$

E(v_t)＝0

$E\left(v_j{v}_t^T\right)=Q{\mathrm{\δ}}_{\mathrm{jt}}---\left(5\right)$

E(η_t)＝0

$E\left({\mathrm{\η}}_j{\mathrm{\η}}_t^T\right)=G{\mathrm{\δ}}_{\mathrm{jt}}---\left(6\right)$

式中R、Q和G分别是状态、测量和偏差状态噪声的方差；

步骤a2，将所述偏差向量扩展为状态向量，获得状态方程和观测方程：

所述状态方程为：

${\tilde{x}}_{t+1}=\tilde{F}{\tilde{x}}_t+{\tilde{w}}_t$

所述观测方程为：

$y_t=\widetilde{\mathrm{\Φ}}{\tilde{x}}_t+v_t$

其中， ${\tilde{x}}_{t+1}=\left[\begin{array}{c}{x}_{t+1}\\ b_{t+1}\end{array}\right],$ $\tilde{F}=\left[\begin{array}{cc}F& 0\\ 0& A\end{array}\right],$ ${\tilde{w}}_t=\left[\begin{array}{c}{w}_t\\ {\mathrm{\η}}_t\end{array}\right],$ $\widetilde{\mathrm{\Φ}}=\left[\begin{array}{cc}\mathrm{\Φ}& I\end{array}\right],$ I为单位矩阵。

进一步地，所述步骤b可以包含：

步骤b1，将状态变量进行时间上的扩展，即

$X_{t+1}=[x_{t+1}^1,x_t^1,...,x_{t-K+2}^1,x_{t+1}^2,x_t^2,...,x_{t-K+2}^2,...x_{t+1}^N,x_t^N,...,x_{t-K+2}^N]$

其中上角标为维数，

步骤b2，将传感器测量值进行时间上的扩展，即

$Y_{t+1}=[y_{t+1}^1,y_t^1,...,y_{t-K+2}^1,y_{t+1}^2,y_t^2,...,y_{t-K+2}^2,...y_{t+1}^p,y_t^p,...,y_{t-K+2}^p]$

其中上角标为维数，

步骤b3，使状态转移矩阵为：

${\stackrel{\‾}{F}}_{\mathrm{ppf}}=\left[\begin{array}{cccc}{F}_{\mathrm{ppf}}^{\left(1\right)}& 0& ...& 0\\ 0& F_{\mathrm{ppf}}^{\left(2\right)}& ...& 0\\ ...& ...& ...& ...\\ 0& ...& 0& F_{\mathrm{ppf}}^{\left(N\right)}\end{array}\right]$

$F_{\mathrm{ppf}}=\left[\begin{array}{cccc}h\left(0\right)& h\left(1\right)& ...& h(K-1)\\ 1& 0& ...& 0\\ ...& ...& ...& ...\\ 0& ...& 1& 0\end{array}\right]$

步骤b4，使测量矩阵为：

${\mathrm{\Φ}}_{\mathrm{ppf}}=\left[\begin{array}{ccccccc}\mathrm{\φ}\left(\mathrm{1,1}\right)& 0_{1\×(K-1)}& \mathrm{\φ}\left(\mathrm{1,2}\right)& 0_{1\×(K-1)}& ...& \mathrm{\φ}(1,N)& 0_{1\×(K-1)}\\ 0& \mathrm{\φ}\left(\mathrm{1,1}\right)& 0_{1\×(K-1)}& \mathrm{\φ}\left(\mathrm{1,2}\right)& ...& \mathrm{\φ}(1,N)& 0_{1\×(K-2)}\\ ...& ...& ...& ...& ...& ...& ...\\ 0& 0& ...& \mathrm{\φ}\left(\mathrm{1,1}\right)& 0& ...& \mathrm{\φ}(1,N)\\ ...& ...& ...& ...& ...& ...& ...\\ \mathrm{\φ}(p,1)& 0_{1\×(K-1)}& \mathrm{\φ}(p,2)& 0_{1\×(K-1)}& ...& \mathrm{\φ}(p,N)& 0_{1\×(K-1)}\\ ...& ...& ...& ...& ...& ...& ...\\ 0& 0& ...& \mathrm{\φ}(p,1)& 0& ...& \mathrm{\φ}(p,N)\end{array}\right]$

步骤b5，扩展后的状态方程和观测方程分别为：

$X_{t+1}={\stackrel{\‾}{F}}_{\mathrm{ppf}}{X}_t$

Y_t＝Φ_ppfX_t+V_t

其中，V_t为扩展后的测量噪声，且有V_t＝[v_t，…，v_t-K+1]^T。

进一步地，所述步骤c可以包含：

步骤c1，假设传感器的测量偏差与可测物理量c_t之间存在线性关系，即：

b_t≈α_tc_t

状态方程和观测方程分别为：

${\tilde{X}}_{t+1}={\tilde{F}}_c{\tilde{X}}_t$

$Y_t={\widetilde{\mathrm{\Φ}}}_c{\tilde{X}}_t+V_t$

其中

${\tilde{X}}_{t+1}={[x_{t+1},x_t,...,x_{t-K+2},{\mathrm{\α}}_{t+1},{\mathrm{\α}}_t,...,{\mathrm{\α}}_{t-K+2}]}^T$

Y_t＝[y_t，…，y_t-K+1]^T

${\tilde{F}}_c=\left[\begin{array}{cc}{F}_{\mathrm{ppf}}& 0\\ 0& F_{\mathrm{ppf}}\end{array}\right]$

${\widetilde{\mathrm{\Φ}}}_c=\left[\begin{array}{cccccccc}\mathrm{\Φ}& 0& ...& 0& c_t& 0& ...& 0\\ 0& \mathrm{\Φ}& ...& 0& 0& c_{t-1}& ...& 0\\ 0& ...& ...& ...& 0& ...& ...& ...\\ 0& ...& 0& \mathrm{\Φ}& 0& ...& 0& c_{t-K+1}\end{array}\right]$

c_t≠c_t-1≠…≠c_t-K+1.

步骤c2，假设t＝0时的状态

和方差P₀已知， ${\hat{\tilde{X}}}_{t+1}={[{\hat{x}}_{t+1},{\hat{x}}_t,...,{\hat{x}}_{t-K+2},{\widehat{\mathrm{\α}}}_{t+1},{\widehat{\mathrm{\α}}}_t,...,{\widehat{\mathrm{\α}}}_{t-K+2}]}^T,$ 对t＝1，2，…(t为自然数)，

对一阶信号，令

$h\left(k\right)=\frac{4K-6k-4}{K(K-1)}$

对二阶信号，令

$h\left(k\right)=\frac{{9K}^2+(-27-36k)K+{30k}^2+42k+18}{K^3-3K^2+2K}$

步骤c21，预测：

${\hat{\tilde{X}}}_{t+1}^{-}={\tilde{F}}_c{\hat{\tilde{X}}}_t$

$P_{t+1}^{-}={\tilde{F}}_c{P}_t{\tilde{F}}_c^T$

为t+1时刻扩展状态预测值，

为扩展的状态转移矩阵，

为t时刻的扩展状态估计值，

为t+1时刻状态估计方差预测值，

P_t为t时刻状态估计方差，

步骤c22，更新

$K_{t+1}=P_{t+1}^{-}{\tilde{H}}_c^T({\tilde{H}}_c{P}_{t+1}^{-}{\tilde{H}}_c^T+R_{t+1})$

K_t+1为t+1时刻新息加权值，

为扩展后的测量矩阵，R_t+1为测量方差阵，

${\hat{\tilde{X}}}_{t+1}={\hat{\tilde{X}}}_{t+1}^{-}+K_{t+1}(Y_{t+1}-{\tilde{H}}_c{\hat{\tilde{X}}}_{t+1}^{-})$

$P_{t+1}=(I-K_{t+1}{\tilde{H}}_c)P_{t+1}^{-}$

Y_t+1为扩展后的测量值，P_t+1为t+1时刻的状态估计方差阵；

步骤c23，计算出测量偏差

${\hat{b}}_{t+1}={\widehat{\mathrm{\α}}}_{t+1}{c}_{t+1}$

步骤c3，将偏差值从测量值中减去。

本发明的有益效果是：

通过对状态方程和测量方程进行扩展，并利用扩展状态的卡尔曼滤波估计偏差，能在进行系统状态估计的同时估计出测量偏差，进而可以从测量值中将其除去，因此估计出的系统状态与真实值非常接近，完全消除了传感器测量偏差对状态估计结果的影响。

附图说明

附图1所示为仿真中传感器的工作温度变化图。

附图2所示为系统的真实状态与传感器的测量值。

附图3所示为根据现有技术和根据本发明方法进行状态估计时，所得到的状态估计结果的区别。

附图4所示为根据现有技术和根据本发明方法估计测量偏差时，所得到的测量偏差的区别。

附图5所示为根据本发明方法进行测量偏差估计时，得到的测量偏差估计结果。

附图6所示为根据本发明方法进行测量偏差估计时，得到的测量偏差估计标准差。

具体实施方式

下面结合实例详细说明本发明的具体实施。

在一种实施方式，一个含有传感器测量偏差的线性动态系统的状态空间模型可以描述为：

x_t+1＝Fx_t+w_t    (1)

y_t＝Φx_t+b_t+v_t    (2)

b_t+1＝Ab_t+η_t    (3)

其中，下标t表示时间，x_t+1为N×1维的状态向量，F为状态转移矩阵，y_t为p×1维的观测向量，Φ为测量矩阵，b_t为M×1维的传感器测量偏差向量，A为b_t的状态转移矩阵，w_t为状态噪声，v_t为测量噪声，η_t为偏差状态噪声，且状态噪声w_t，测量噪声v_t及偏差状态噪声η_t均为相互独立的零均值的高斯白噪声，且满足：

E(w_t)＝0

$E\left(w_j{w}_t^T\right)=R{\mathrm{\δ}}_{\mathrm{jt}}---\left(4\right)$

E(v_t)＝0

$E\left(v_j{v}_t^T\right)=Q{\mathrm{\δ}}_{\mathrm{jt}}---\left(5\right)$

E(η_t)＝0

$E\left({\mathrm{\η}}_j{\mathrm{\η}}_t^T\right)=G{\mathrm{\δ}}_{\mathrm{jt}}---\left(6\right)$

式中R、Q和G分别是状态、测量和偏差状态噪声的方差。

将偏差向量扩展为状态向量时，式(1)-(3)可写为：

${\tilde{x}}_{t+1}=\tilde{F}{\tilde{x}}_t+{\tilde{w}}_t---\left(7\right)$

$y_t=\widetilde{\mathrm{\Φ}}{\tilde{x}}_t+v_t---\left(8\right)$

其中， ${\tilde{x}}_{t+1}=\left[\begin{array}{c}{x}_{t+1}\\ b_{t+1}\end{array}\right],$ $\tilde{F}=\left[\begin{array}{cc}F& 0\\ 0& A\end{array}\right],$ ${\tilde{w}}_t=\left[\begin{array}{c}{w}_t\\ {\mathrm{\η}}_t\end{array}\right],$ $\widetilde{\mathrm{\Φ}}=\left[\begin{array}{cc}\mathrm{\Φ}& I\end{array}\right],$ I为单位矩阵。

为了通过扩展状态的卡尔曼滤波估计出测量偏差，需要利用一种全新的建模方法对状态方程建模，即建立状态方程的多项式预测模型。

下面简述状态方程的一种多项式预测建模方法。

1、对状态变量扩展，

将状态变量进行时间上的扩展，即

$X_{t+1}=[x_{t+1}^1,x_t^1,...,x_{t-K+2}^1,x_{t+1}^2,x_t^2,...,x_{t-K+2}^2,...x_{t+1}^N,x_t^N,...,x_{t-K+2}^N]---\left(9\right)$

其中上角标为维数。

2、对传感器测量值扩展，

将传感器测量值进行时间上的扩展，即

$Y_{t+1}=[y_{t+1}^1,y_t^1,...,y_{t-K+2}^1,y_{t+1}^2,y_t^2,...,y_{t-K+2}^2,...y_{t+1}^p,y_t^p,...,y_{t-K+2}^p]---\left(10\right)$

其中上角标为维数。

3、使状态转移矩阵为，

${\stackrel{\‾}{F}}_{\mathrm{ppf}}=\left[\begin{array}{cccc}{F}_{\mathrm{ppf}}^{\left(1\right)}& 0& ...& 0\\ 0& F_{\mathrm{ppf}}^{\left(2\right)}& ...& 0\\ ...& ...& ...& ...\\ 0& ...& 0& F_{\mathrm{ppf}}^{\left(N\right)}\end{array}\right]---\left(11\right)$

$F_{\mathrm{ppf}}=\left[\begin{array}{cccc}h\left(0\right)& h\left(1\right)& ...& h(K-1)\\ 1& 0& ...& 0\\ ...& ...& ...& ...\\ 0& ...& 1& 0\end{array}\right]---\left(12\right)$

4、使测量矩阵为，

${\mathrm{\Φ}}_{\mathrm{ppf}}=\left[\begin{array}{ccccccc}\mathrm{\φ}\left(\mathrm{1,1}\right)& 0_{1\×(K-1)}& \mathrm{\φ}\left(\mathrm{1,2}\right)& 0_{1\×(K-1)}& ...& \mathrm{\φ}(1,N)& 0_{1\×(K-1)}\\ 0& \mathrm{\φ}\left(\mathrm{1,1}\right)& 0_{1\×(K-1)}& \mathrm{\φ}\left(\mathrm{1,2}\right)& ...& \mathrm{\φ}(1,N)& 0_{1\×(K-2)}\\ ...& ...& ...& ...& ...& ...& ...\\ 0& 0& ...& \mathrm{\φ}\left(\mathrm{1,1}\right)& 0& ...& \mathrm{\φ}(1,N)\\ ...& ...& ...& ...& ...& ...& ...\\ \mathrm{\φ}(p,1)& 0_{1\×(K-1)}& \mathrm{\φ}(p,2)& 0_{1\×(K-1)}& ...& \mathrm{\φ}(p,N)& 0_{1\×(K-1)}\\ ...& ...& ...& ...& ...& ...& ...\\ 0& 0& ...& \mathrm{\φ}(p,1)& 0& ...& \mathrm{\φ}(p,N)\end{array}\right]---\left(13\right)$

系统的状态方程和观测方程可写为：

$X_{t+1}={\stackrel{\‾}{F}}_{\mathrm{ppf}}{X}_t---\left(14\right)$

Y_t＝Φ_ppfX_t+V_t    (15)

其中，V_t为扩展后的测量噪声，且有V_t＝[v_t，…，v_t-K+1]^T。

在进行传感器测量偏差估计时，假设与传感器的测量偏差与一可测物理量c_t之间存在一近似线性关系，例如传感器的测量偏差与其工作温度存在这样的关系，即：

b_t≈α_tc_t    (16)

则式(7)和(8)可写为：

${\tilde{X}}_{t+1}={\tilde{F}}_c{\tilde{X}}_t---\left(17\right)$

$Y_t={\widetilde{\mathrm{\Φ}}}_c{\tilde{X}}_t+V_t---\left(18\right)$

其中

${\tilde{X}}_{t+1}={[x_{t+1},x_t,...,x_{t-K+2},{\mathrm{\α}}_{t+1},{\mathrm{\α}}_t,...,{\mathrm{\α}}_{t-K+2}]}^T---\left(10\right)$

Y_t＝[y_t，…，y_t-K+1]^T    (20)

${\tilde{F}}_c=\left[\begin{array}{cc}{F}_{\mathrm{ppf}}& 0\\ 0& F_{\mathrm{ppf}}\end{array}\right]---\left(21\right)$

${\widetilde{\mathrm{\Φ}}}_c=\left[\begin{array}{cccccccc}\mathrm{\Φ}& 0& ...& 0& c_t& 0& ...& 0\\ 0& \mathrm{\Φ}& ...& 0& 0& c_{t-1}& ...& 0\\ 0& ...& ...& ...& 0& ...& ...& ...\\ 0& ...& 0& \mathrm{\Φ}& 0& ...& 0& c_{t-K+1}\end{array}\right]---\left(22\right)$

当c_t≠c_t-1≠…≠c_t-K+1时，易知

和

满足控制系统的可观测性条件，即可以用扩展的卡尔曼滤波估计偏差。

综上，基于多项式预测模型的传感器偏差估计方法具体步骤如下：

假设t＝0时的状态

和方差P₀已知， ${\hat{\tilde{X}}}_{t+1}={[{\hat{x}}_{t+1},{\hat{x}}_t,...,{\hat{x}}_{t-K+2},{\widehat{\mathrm{\α}}}_{t+1},{\widehat{\mathrm{\α}}}_t,...,{\widehat{\mathrm{\α}}}_{t-K+2}]}^T,$ 对t＝1，2，…

通常对一阶信号，

$h\left(k\right)=\frac{4K-6k-4}{K(K-1)}---\left(23\right)$

对二阶信号，

$h\left(k\right)=\frac{{9K}^2+(-27-36k)K+{30k}^2+42k+18}{K^3-3K^2+2K}---\left(24\right)$

1、预测：

${\hat{\tilde{X}}}_{t+1}^{-}={\tilde{F}}_c{\hat{\tilde{X}}}_t---\left(25\right)$

$P_{t+1}^{-}={\tilde{F}}_c{P}_t{\tilde{F}}_c^T---\left(26\right)$

为t+1时刻扩展状态预测值，

为扩展的状态转移矩阵，

为t时刻的扩展状态估计值，

为t+1时刻状态估计方差预测值，

P_t为t时刻状态估计方差，

2、更新：

$K_{t+1}=P_{t+1}^{-}{\tilde{H}}_c^T({\tilde{H}}_c{P}_{t+1}^{-}{\tilde{H}}_c^T+R_{t+1})---\left(27\right)$

K_t+1为t+1时刻新息加权值，

为扩展后的测量矩阵，R_t+1为测量方差阵，

${\hat{\tilde{X}}}_{t+1}={\hat{\tilde{X}}}_{t+1}^{-}+K_{t+1}(Y_{t+1}-{\tilde{H}}_c{\hat{\tilde{X}}}_{t+1}^{-})---\left(28\right)$

$P_{t+1}=(I-K_{t+1}{\tilde{H}}_c)P_{t+1}^{-}---\left(29\right)$

Y_t+1为扩展后的测量值，P_t+1为t+1时刻的状态估计方差阵；

3、计算出测量偏差：

${\hat{b}}_{t+1}={\widehat{\mathrm{\α}}}_{t+1}{c}_{t+1}---\left(30\right)$

下面给出本发明的一个具体例子：

设一个传感器测量系统模型如式(31)(32)所描述：

x_t＝0.1t    (31)

y_t＝x_t+b_t+v_t    (32)

其中，x_t表示系统状态，y_t表示传感器测量输出，b_t表示传感器测量偏差，假设传感器测量偏差与工作温度的关系由式(33)表示，TE_t代表t时刻可测量的传感器工作温度。

b_t＝0.001t·TE_t    (33)

上述h(k)的计算公式是从一个复杂的方程推导出来的，该方程具有参数L和K，不同的L和K的取值，对应不同的表达形式，在公式23中，N＝1，L＝1；在公式24中，N＝1，L＝2。具体实施时，若要求出在不同L和K的取值下的h(k)的值，可参阅文献“P.Heinonen，Y.Neuvo.FIR-Median HybridFilters with Predictive FIR Substructures.IEEE Trans.on Acoustics，Speech andSignal Processing.1998，36(6)：892-899”

如果直接采用扩展状态的卡尔曼滤波算法，则测量偏差无法满足可观测性条件，因此无法获得准确的估计结果。而采用本发明提出的算法，通过控制传感器工作的温度并测量出其随时间变化的温度值，那么通过本发明的模型和方法就能够同时准确估计出系统状态和传感器的测量偏差，计算机仿真结果如图1-图6所示。

仿真中假设人为的使传感器的工作温度按照图1所示的规律变化，图2所示为系统的真实状态与传感器测量值，由于仿真中的测量矩阵为单位阵，传感器的测量值与真实状态在均值上是相等的，但从图2中可以明显看出，由于与温度有关的测量偏差的存在，传感器的测量输出随时间逐渐偏离真实值。

图3和图4所示为根据现有技术未校准测量偏差进行状态估计与根据本发明方法法进行状态估计的结果比较。从图3中可以看出，采用未经处理的测量数据直接进行状态估计，得到的估计结果由于受到测量偏差的影响，远远偏离了真实值；而采用本发明方法，如图4所示，由于可以同时估计出测量偏差，进而可以从测量值中将其除去，因此估计出的系统状态与真实值非常接近，完全消除了传感器测量偏差对状态估计结果的影响。

图5和图6所示为传感器测量偏差的估计结果。从图中可以看出，采用本发明算法，在进行系统状态估计的同时可以获得传感器测量偏差的准确估计，如果不采用本发明的多项式预测建模方法，在未知系统状态和测量偏差动态特性的情况下是无法同时估计出二者的。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员，在不脱离本发明构思的前提下，还可以做出若干改变、改进或润饰也应视为本发明的保护范围。

Claims (4)

1.一种传感器慢偏故障的校正方法，其特征在于，包含如下步骤：

步骤a：对传感器的测量值建立状态方程和观测方程；

步骤b：对所述状态方程和观测方程进行时间上的扩展；

步骤c：利用扩展状态的卡尔曼滤波估计偏差，再将偏差值从测量值中减去。

2.根据权利要求1所述的传感器慢偏故障的校正方法，其特征在于，所述步骤a包含：

步骤a1：将含有传感器测量偏差的线性动态系统的状态空间模型描述为：

x_t+1＝Fx_t+w_t    (1)

y_t＝Φx_t+b_t+v_t    (2)

b_t+1＝Ab_t+η_t    (3)

其中，下标t表示时间，x_t+1为N×1维的状态向量，F为状态转移矩阵，y_t为p×1维的观测向量，Φ为测量矩阵，b_t为M×1维的传感器测量偏差向量，A为b_t的状态转移矩阵，w_t为状态噪声，v_t为测量噪声，η_t为偏差状态噪声，且状态噪声w_t，测量噪声v_t及偏差状态噪声η_t均为相互独立的零均值的高斯白噪声，且满足：

E(w_t)＝0

$E\left(w_j{w}_t^T\right)=R{\mathrm{\δ}}_{\mathrm{jt}}---\left(4\right)$

E(v_t)＝0

$E\left(v_j{v}_t^T\right)=Q{\mathrm{\δ}}_{\mathrm{jt}}---\left(5\right)$

E(η_t)＝0

$E\left({\mathrm{\η}}_j{\mathrm{\η}}_t^T\right)=G{\mathrm{\δ}}_{\mathrm{jt}}---\left(6\right)$

式中R、Q和G分别是状态、测量和偏差状态噪声的方差；

步骤a2，将所述偏差向量扩展为状态向量，获得状态方程和观测方程：

所述状态方程为：

${\tilde{x}}_{t+1}=\tilde{F}{\tilde{x}}_t+{\tilde{w}}_t$

所述观测方程为：

$y_t=\widetilde{\mathrm{\Φ}}{\tilde{x}}_t+v_t$

其中， ${\tilde{x}}_{t+1}=\left[\begin{array}{c}{x}_{t+1}\\ b_{t+1}\end{array}\right],$ $\tilde{F}=\left[\begin{array}{cc}F& 0\\ 0& A\end{array}\right],$ ${\tilde{w}}_t=\left[\begin{array}{c}{w}_t\\ {\mathrm{\η}}_t\end{array}\right],$ $\widetilde{\mathrm{\Φ}}=\left[\begin{array}{cc}\mathrm{\Φ}& I\end{array}\right],$ I为单位矩阵。

3.根据权利要求1所述的传感器慢偏故障的校正方法，其特征在于，所述步骤b包含：

步骤b1，将状态变量进行时间上的扩展，即

$X_{t+1}=[x_{t+1}^1,x_t^1,...,x_{t-K+2}^1,x_{t+1}^2,x_t^2,...,x_{t-K+2}^2,...x_{t+1}^N,x_t^N,...,x_{t-K+2}^N]$

其中上角标为维数，

步骤b2，将传感器测量值进行时间上的扩展，即

$Y_{t+1}=[y_{t+1}^1,y_t^1,...,y_{t-K+2}^1,y_{t+1}^2,y_t^2,...,y_{t-K+2}^2,...y_{t+1}^p,y_t^p,...,y_{t-K+2}^p]$

其中上角标为维数，

步骤b3，使状态转移矩阵为：

${\stackrel{\‾}{F}}_{\mathrm{ppf}}=\left[\begin{array}{cccc}{F}_{\mathrm{ppf}}^{\left(1\right)}& 0& ...& 0\\ 0& F_{\mathrm{ppf}}^{\left(2\right)}& ...& 0\\ ...& ...& ...& ...\\ 0& ...& 0& F_{\mathrm{ppf}}^{\left(N\right)}\end{array}\right]$

$F_{\mathrm{ppf}}=\left[\begin{array}{cccc}h\left(0\right)& h\left(1\right)& ...& h(K-1)\\ 1& 0& ...& 0\\ ...& ...& ...& ...\\ 0& ...& 1& 0\end{array}\right]$

步骤b4，使测量矩阵为：

${\mathrm{\Φ}}_{\mathrm{ppf}}=\left[\begin{array}{ccccccc}\mathrm{\φ}\left(\mathrm{1,1}\right)& 0_{1\×(K-1)}& \mathrm{\φ}\left(\mathrm{1,2}\right)& 0_{1\×(K-1)}& ...& \mathrm{\φ}(1,N)& 0_{1\×(K-1)}\\ 0& \mathrm{\φ}\left(\mathrm{1,1}\right)& 0_{1\×(K-1)}& \mathrm{\φ}\left(\mathrm{1,2}\right)& ...& \mathrm{\φ}(1,N)& 0_{1\×(K-2)}\\ ...& ...& ...& ...& ...& ...& ...\\ 0& 0& ...& \mathrm{\φ}\left(\mathrm{1,1}\right)& 0& ...& \mathrm{\φ}(1,N)\\ ...& ...& ...& ...& ...& ...& ...\\ \mathrm{\φ}(p,1)& 0_{1\×(K-1)}& \mathrm{\φ}(p,2)& 0_{1\×(K-1)}& ...& \mathrm{\φ}(p,N)& 0_{1\×(K-1)}\\ ...& ...& ...& ...& ...& ...& ...\\ 0& 0& ...& \mathrm{\φ}(p,1)& 0& ...& \mathrm{\φ}(p,N)\end{array}\right]$

步骤b5，扩展后的状态方程和观测方程分别为：

$X_{t+1}={\stackrel{\‾}{F}}_{\mathrm{ppf}}{X}_t$

Y_t＝Φ_ppfX_t+V_t

其中，V_t为扩展后的测量噪声，且有V_t＝[v_t，…，v_t-K+1]^T。

4.根据权利要求1所述的传感器慢偏故障的校正方法，其特征在于，所述步骤c包含：

步骤c1，假设传感器的测量偏差与可测物理量c_t之间存在线性关系，即：

b_t≈α_tc_t

状态方程和观测方程分别为：

${\tilde{X}}_{t+1}={\tilde{F}}_c{\tilde{X}}_t$

$Y_t={\widetilde{\mathrm{\Φ}}}_c{\tilde{X}}_t+V_t$

其中

${\tilde{X}}_{t+1}={[x_{t+1},x_t,...,x_{t-K+2},{\mathrm{\α}}_{t+1},{\mathrm{\α}}_t,...,{\mathrm{\α}}_{t-K+2}]}^T$

Y_t＝[y_t，…，y_t-K+1]^T

${\tilde{F}}_c=\left[\begin{array}{cc}{F}_{\mathrm{ppf}}& 0\\ 0& F_{\mathrm{ppf}}\end{array}\right]$

${\widetilde{\mathrm{\Φ}}}_c=\left[\begin{array}{cccccccc}\mathrm{\Φ}& 0& ...& 0& c_t& 0& ...& 0\\ 0& \mathrm{\Φ}& ...& 0& 0& c_{t-1}& ...& 0\\ 0& ...& ...& ...& 0& ...& ...& ...\\ 0& ...& 0& \mathrm{\Φ}& 0& ...& 0& c_{t-K+1}\end{array}\right]$

c_t≠c_t-1≠…≠c_t-K+1.

步骤c2，假设t＝0时的状态

和方差P₀已知， ${\hat{\tilde{X}}}_{t+1}={[{\hat{x}}_{t+1},{\hat{x}}_t,...,{\hat{x}}_{t-K+2},{\widehat{\mathrm{\α}}}_{t+1},{\widehat{\mathrm{\α}}}_t,...,{\widehat{\mathrm{\α}}}_{t-K+2}]}^T,$ 对t＝1，2，…(t为自然数)，

对一阶信号，令

$h\left(k\right)=\frac{4K-6k-4}{K(K-1)}$

对二阶信号，令

$h\left(k\right)=\frac{{9K}^2+(-27-36k)K+{30k}^2+42k+18}{K^3-3K^2+2K}$

步骤c21，预测：

${\hat{\tilde{X}}}_{t+1}^{-}={\tilde{F}}_c{\hat{\tilde{X}}}_t$

$P_{t+1}^{-}={\tilde{F}}_c{P}_t{\tilde{F}}_c^T$

为t+1时刻扩展状态预测值，

为扩展的状态转移矩阵，

为t时刻的扩展状态估计值，

为t+1时刻状态估计方差预测值，

P_t为t时刻状态估计方差，

步骤c22，更新

$K_{t+1}=P_{t+1}^{-}{\tilde{H}}_c^T({\tilde{H}}_c{P}_{t+1}^{-}{\tilde{H}}_c^T+R_{t+1})$

K_t+1为t+1时刻新息加权值，

为扩展后的测量矩阵，R_t+1为测量方差阵，

${\hat{\tilde{X}}}_{t+1}={\hat{\tilde{X}}}_{t+1}^{-}+K_{t+1}(Y_{t+1}-{\tilde{H}}_c{\hat{\tilde{X}}}_{t+1}^{-})$

$P_{t+1}=(I-K_{t+1}{\tilde{H}}_c)P_{t+1}^{-}$

Y_t+1为扩展后的测量值，P_t+1为t+1时刻的状态估计方差阵；

步骤c23，计算出测量偏差

${\hat{b}}_{t+1}={\widehat{\mathrm{\α}}}_{t+1}{c}_{t+1}$

步骤c3，将偏差值从测量值中减去。

Images