CN103217172B - 一种卡尔曼滤波传感器信息融合的故障检测方法 - Google Patents

Abstract

本发明公开了一种卡尔曼滤波传感器信息融合的故障检测方法，包括以下几个步骤：步骤一：根据卡尔曼滤波理论，建立线性动态系统的状态方程和观测方程；步骤二：根据步骤一得出的观测方程，利用最小二乘方法获取状态估计和相应的均方误差阵、新息序列；步骤三：利用已知的新息序列，得到不同的渠道归一化的新息序列；并且组成m通道平行传感器的创新矩阵；步骤四：根据步骤三所得的创新矩阵，获取创新矩阵的谱范数和谱范数的均值；步骤五：对卡尔曼滤波传感器的信息融合进行故障检测；本发明采用数理统计和区间估计的方法，简化了复杂的计算，极大地提高了故障检测速度。

Description

一种卡尔曼滤波传感器信息融合的故障检测方法

技术领域

本发明涉及一种基于标准化创新矩阵的卡尔曼滤波传感器信息融合故障检测方法，属于卡尔曼滤波器信息融合的故障检测技术领域。

背景技术

智能化、高精度、高可靠性是未来飞行器对导航系统的要求。信息融合技术可满足未来飞行器对导航系统的这些要求。在导航和控制领域中，不可避免要对卡尔曼滤波器进行有效的测试。基于此目的已经产生了一些算法并且也产生了相应测试卡尔曼滤波测试的算法技术，这样不仅能够确保故障定位和检测，还可以进行估计修正。目前，有很多卡尔曼滤波器测试的算法，应用哲学算法可以对卡尔曼滤波器不同特征标志进行故障检测。尽管卡尔曼滤波检测算法种类繁多，但是迄今为止，多通道卡尔曼滤波器修正的监测和诊断的问题还没有得到解决。

在航空航天和海军导航系统中，单一的导航系统已经无法满足系统的要求，组合导航系统越来越受到高度的重视，GPS，DGPS，GLONASS和INS系统都通过卡尔曼滤波集成不同的组合。在多传感器融合的复杂系统中，即同时从多个来源的动态系统状态向量来获取信息，联邦或平行卡尔曼滤波器用来整合不同的导航系统够能到得到和满意的效果。Kalman滤波器信息融合。系统参数和状态多通道估计算法已成熟并用于估计一个动态系统的一个数学模型，以及多个测量通道的量测。在卡尔曼滤波器的这些模型中，有效数据的联合处理，可以改善状态向量的估计精度和数据处理的可靠性。这些算法对检测Kalman滤波器信息融合的应用而言，每个估计信道都需要它自己的“故障检测器”，这样就需要很大的计算量。所以，考虑到多信道估计程序需要大量的计算，采用先前的技术实现测试Kalman滤波器信息融合就不是一个简单的问题。因此，有必要开发用于测试的多通道估计程序的简单算法，以执行实时测试滤波器没有故障情况下参量变化的先验信息。

发明内容

本发明的目的是为了简化在GPS/INS/DVL组合导航系统中卡尔曼滤波器的传感器信息融合故障检测方法，提出一种卡尔曼滤波传感器信息融合的故障检测方法，本发明能够简单快速检测卡尔曼滤波器的传感器信息融合情况，在得到了标准化创新矩阵谱范数后，再对其进行数理统计和区间估计，从而对卡尔曼滤波器的传感器信息融合进行简单快速的故障检测。

一种卡尔曼滤波传感器信息融合的故障检测方法，包括以下几个步骤：

步骤一：根据卡尔曼滤波理论，建立线性动态系统的状态方程和观测方程；

步骤二：根据步骤一得出的观测方程，利用最小二乘方法获取状态估计和相应的均方误差阵、新息序列；

步骤三：利用已知的新息序列，得到不同的渠道归一化的新息序列；并且组成m通道平行传感器的创新矩阵；

步骤四：根据步骤三所得的创新矩阵，获取创新矩阵的谱范数和谱范数的均值；

步骤五：对卡尔曼滤波传感器的信息融合进行故障检测；

本发明的优点在于：

（1）本发明采用平行算法同时联立所有测量通道组合为矩阵，这样不需要对每个通道分别进行检测，降低了系统的复杂程度；

（2）本发明采用数理统计和区间估计的方法，简化了复杂的计算，极大地提高了故障检测速度。

附图说明

图1是GPS/INS/DVL组合导航系统的容错结构图；

图2是卡尔曼滤波器的传感器信息融合图；

图3是本发明的方法流程图。

具体实施方式

下面将结合附图和实施例对本发明作进一步的详细说明。

如图1所示，在GPS/INS/DVL组合导航系统中（即卫星/惯导/多普勒计程仪组合导航系统），认为惯导系统除了其固有的漂移之外不发生其它故障。当任意子系统发生故障时，对故障进行有效的检测和滤波处理。当故障不可消除时，将有故障传感器隔离，当故障消除后，自动进行系统恢复。如图2所示，在卡尔曼滤波器中有m个传感器通道，要同时对m个通道的传感器进行信息融合。本发明提出一种卡尔曼滤波传感器信息融合的故障检测方法。基于引进的标准规范化创新矩阵的光谱范数数学期望的统计。这种方法允许数学期望和实时的创新序列方差的同时测试，并且不需要先验信息在故障时刻变化值的统计特征。使用这种方法，可对信息融合的卡尔曼滤波器进行故障检测。

本发明是一种基于标准化创新矩阵的卡尔曼滤波传感器信息融合的故障检测方法，流程如图3所示，包括以下几个步骤：

步骤一：根据卡尔曼滤波理论，建立线性动态系统的状态方程和观测方程；

线性动态系统的状态方程为：

x(k+1)＝Φ(k+1,k)x(k)+G(k+1,k)ω(k)（1）

式中：x(k+1)是系统k+1时刻的n维状态向量，x(k)是系统n维状态向量，Φ(k+1,k)是系统n×n维转移矩阵；G(k+1,k)是n×n维系统噪声驱动阵，ω(k)是n维系统激励噪声序列。k表示卡尔曼滤波的更新时刻，状态x(k)是由m个测量通道所组成，所以x(k)的量测满足以下关系，第i个通道的观测方程：

z_i(k)＝H_i(k)x(k)+ν_i(k)（2）

式中：z_i(k)是第i个测量通道的n维测量向量；H_i(k)是系统第i个通道的n×n维的测量矩阵；ν_i(k)第i个通道的n维量测噪声序列。第i个通道的量测噪声ν(k)的均值是零方差为R_ii(k)，相关阵为各通道相互独立，ν_i(j)是第i个通道j时刻的量测噪声序列， $\mathrm{\δ}\left(\mathrm{kj}\right)=\left\{\begin{array}{cc}0& k\≠j\\ 1& k=j\end{array}\right..$

步骤二：根据步骤一得出的观测方程，利用最小二乘方法获取状态估计和相应的均方误差阵、新息序列；

系统的状态向量可以通过卡尔曼滤波传感器信息融合来估计，可以表示为：

$\hat{x}(k/k)=\hat{x}(k/k-1)+\underset{i=1}{\overset{m}{\mathrm{\Σ}}}P(k/k){H_i}^T\left(k\right)R_{\mathrm{ii}}^{-1}\left(k\right){\mathrm{\Δ}}_i\left(k\right)---\left(3\right)$

式中：是系统的状态估计，P(k/k)表示估计均方误差的相关矩阵，Δ_i(k)表示i通道的新息序列，是状态一步预测：

$\hat{x}(k/k-1)=\mathrm{\Φ}(k,k-1)\hat{x}(k-1)---\left(4\right)$

Δ_i(k)是i通道的新息序列：

${\mathrm{\Δ}}_i\left(k\right)=z_i\left(k\right)-H_i\left(k\right)\hat{x}(k/k-1)---\left(5\right)$

估计均方误差的相关矩阵：

$P^{-1}(k/k)=P^{-1}(k/k-1)+\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{H}_i^T\left(k\right)R_{11}^{-1}\left(k\right)H_i\left(k\right)---\left(6\right)$

一步预测均方误差的相关矩阵：

P(k/k-1)＝Φ(k,k-1)P(k-1/k-1)Φ^T(k,k-1)+G(k,k-1)Q(k-1)G^T(k,k-1)（7）

步骤三：利用已知的新息序列，得到不同的渠道归一化的新息序列；并且组成m通道平行传感器的创新矩阵；

最优估计是通过同时联立所有通道新息序列Δ_i(k)的权重总和定义的。符合方程式(3)-(7)的传感器信息融合滤波器为并行过滤器，正如(3)-(7)所示，估计状态向量平行算法要采用m个滤波器，每个滤波器的增益矩阵、滤波误差的相关矩阵等为已知估计。对于正常运作的卡尔曼滤波传感器信息融合,不同的通道归一化的新息序列为：

${\widetilde{\mathrm{\Δ}}}_i\left(k\right)={[H_i\left(k\right)P(k/k-1)H_i^T\left(k\right)+R_{\mathrm{ii}}\left(k\right)]}^{-1/2}{\mathrm{\Δ}}_i\left(k\right)---\left(8\right)$

其中，是第i个通道的归一化新息序列，并且它服从N(0,1)分布。故障能够导致测量信道特性、计算机失效、异常测量、对象噪声或测量统计特性的变化、卡尔曼滤波产生的真正过程和估计的轨迹的一个方差等突变，这样的一个突变可导致序列的先验特征的变化，所以有必要开发一个有效的方法对标准化新息序列的数学期望和方差同时联立检测。为此，建立两种假设：γ₀-Kalman滤波器工作在正常状态，γ₁-估计系统发生故障。

为此设定：A是一个n×m矩阵(n≥2,m≥2)且为m通道Kalman滤波器规范化创新矩阵，其中，各列为不同通道的规范化创新向量（不同的通道归一化创新的向量的维度是相等的），也就是说创新矩阵为：

$A\left(k\right)=[{\widetilde{\mathrm{\Δ}}}_1\left(k\right),{\widetilde{\mathrm{\Δ}}}_2\left(k\right),...{\widetilde{\mathrm{\Δ}}}_m\left(k\right)]---\left(9\right)$

其中，A(k)是k时刻的m通道平行传感器信息融合kalman滤波器的规范化创新矩阵，是各通道k时刻的归一化新息序列。对于假设检验γ₀和γ₁，矩阵A^T(k)A(k)服从白噪声分布，所以可以应用其特征值、最大值和最小值的特性。由于为随机矩阵特征值确定信心域(或间隔)有一定的困难，所以对于多维动态系统的故障检测问题的研究应用是非常复杂的。

步骤四：根据步骤三所得的创新矩阵，获取创新矩阵的谱范数和谱范数的均值；

在假设检验γ₀和γ₁研究中，要采用矩阵A(k)的谱范数，实数矩阵A(k)的谱范数||||₂的定义公式：

${\left|\right|A\left(k\right)\left|\right|}_2\≡\max \left\{{\left({\mathrm{\λ}}_i\right[A^T\left(k\right)A\left(k\right)\left]\right)}^{1/2}\right\}---\left(10\right)$

式中：λ_i[A^T(k)A(k)]是A^T(k)A(k)的特征值。

矩阵A^T(k)A(k)的特征值的平方根（即λ_i[A^T(k)A(k)]值）为矩阵A(k)的奇异值，所以，矩阵A(k)的范数等于它的最大奇异值，它是实数且为正定的。因为对于任意矩阵，确定奇异数和范数的计算相对于获得特征值简单，这说明卡尔曼滤波的归一化创新矩阵的谱范数的选择可以作为一个检测尺度。假设检验γ₀和γ₁，引进一个矩阵A(k)的谱范数数学期望的一维统计量。由于一个数学期望的估计需要使用其算术平均值，对于k时刻的数学期望值表达式可以写成（即A(k)的谱范数均值）：

$E\left\{{\left|\right|A\left(k\right)\left|\right|}_2\right\}\≈{\stackrel{\‾}{\left|\right|A\left(k\right)\left|\right|}}_2=\frac{1}{k}\underset{j=1}{\overset{k}{\mathrm{\Σ}}}{\left|\right|A\left(j\right)\left|\right|}_2---\left(11\right)$

这里借助Hansen的思想，他曾为一个由零均值和标准差为σ(a_ij∈N(0,1))的随机值所构成的随机矩阵A(k)∈R^n×m范数的数学期望提过许多的限定。假设r_i和a_j代表矩阵A的行和列，所以最大行-列范数为：

μ≡max[||r_i||₂,||a_j||₂]（12）

其中||r_i||和||a_j||是相应的向量范数，并且有下列不等式成立：

E{μ}≤E{||A||₂}≤[max(n,m)]^1/2E{μ}（13）

在实际的估算中公式（13）很难实现，很难估计E{μ}，所以把E{μ}用它的下界来替换：

$\mathrm{\σ}\sqrt{\max (n,m)}=\max \left[E\right\{{\left|\right|r_i\left|\right|}_2\},E\{{\left|\right|a_j\left|\right|}_2\left\}\right]\≤E\left\{\mathrm{\μ}\right\}---\left(14\right)$

所以公式（13）可以改为：

$\mathrm{\σ}\sqrt{\max (n,m)}\≤E\left\{{\left|\right|A\left|\right|}_2\right\}\≤f\left(\max (n,m)\right)\×\mathrm{\σ}\sqrt{\max (n,m)}---\left(15\right)$

式中f是定义的未知函数。并且有在n＝m→∞时函数f渐近为2且f的区间是（1,2），所以定义f的估计值为2，针对以上所述，E{||A(k)||₂}可以有以下定义：

$\mathrm{\σ}\sqrt{\max (n,m)}\≤E\left\{{\left|\right|A\left|\right|}_2\right\}\≤2\mathrm{\σ}\sqrt{\max (n,m)}---\left(16\right)$

表达式（16）表现了随机矩阵A标准差和它的二阶范数之间的关系。考虑到应用规范化创新矩阵来观测由分布的随机元素零均值和有限方差a_ij∈N(0,1)组成的传感器信息融合的卡尔曼滤波器故障，在本发明中不等式（16）明确表达了可以应用到解决诊断问题，可以通过卡尔曼滤波器的归一化新息矩阵的元素a_ij是否服从N(0,1)分布来规定，并要满足不等式（16）。不满足不等式（16）则说明元素a_ij的均值不在为零或者方差不再为单位方差或者{a_ij}不再服从白噪声。

所以根据是否满足不等式（16）来对该传感器信息融合卡尔曼滤波器进行检测，鉴于σ＝1和表达式（11）就可以得到谱范数均值符合以下关系形式：

$\sqrt{\max (n,m)}\≤{\stackrel{\‾}{\left|\right|A\left(k\right)\left|\right|}}_2\≤2\sqrt{\max (n,m)}---\left(17\right)$

步骤五：对卡尔曼滤波传感器的信息融合进行故障检测；

根据以上分析，测试传感器信息融合卡尔曼滤波器的故障检测问题，决策规则可假设为以下形式：

${\mathrm{\&gamma;}}_0:\sqrt{\max (n,m)}<{\stackrel{\&OverBar;}{\left|\right|A\left(k\right)\left|\right|}}_2<2\sqrt{\max (n,m)},$ ：无故障运行

${\mathrm{\&gamma;}}_1:{\stackrel{\&OverBar;}{\left|\right|A\left(k\right)\left|\right|}}_2\&le;\sqrt{\max (n,m)}$ 或 ${\stackrel{\&OverBar;}{\left|\right|A\left(k\right)\left|\right|}}_2\&GreaterEqual;2\sqrt{\max (n,m)},$ ：故障（18）

这种方法能够容易的确定m通道卡尔曼滤波器的归一化新息序列的谱范数的数学期望的界限，并且对所有通道的新息序列的数学期望和方差同时检查。

在GPS/INS/DVL组合导航系统中，采用上述方法后，就可简便的对系统进行故障诊断，不需要其它方法的辅助，若判断出GPS/INS子系统发生故障，则可判断GPS故障；若判断出INS/DVL子系统发生故障，则可判断DVL发生故障。

Claims (1)

1.一种卡尔曼滤波传感器信息融合的故障检测方法，包括以下几个步骤：

步骤一：建立线性动态系统的状态方程和观测方程；

线性动态系统的状态方程为：

x(k+1)＝Φ(k+1,k)x(k)+G(k+1,k)ω(k)(1)

式中：x(k+1)是系统k+1时刻的n维状态向量，x(k)是系统n维状态向量，Φ(k+1,k)是系统n×n维转移矩阵，G(k+1,k)是n×n维系统噪声驱动阵，ω(k)是n维系统激励噪声序列，k表示卡尔曼滤波的更新时刻，状态x(k)是由m个测量通道所组成，第i个通道的观测方程：

z_i(k)＝H_i(k)x(k)+ν_i(k)(2)

式中：z_i(k)是第i个测量通道的n维测量向量，H_i(k)是系统第i个通道的n×n维的测量矩阵，ν_i(k)第i个通道的n维量测噪声序列，第i个通道的量测噪声ν(k)的均值是零，方差为R_ii(k)，相关阵为E[ν_i(k)ν_i ^T(j)]＝R_ii(k)δ(kj)，各通道相互独立，ν_i(j)是第i个通道j时刻的量测噪声序列， $\mathrm{\&delta;}\left(kj\right)=\left\{\begin{array}{cc}0& k\&NotEqual;j\\ 1& k=j\end{array}\right.;$

步骤二：根据步骤一得出的观测方程，利用最小二乘方法获取状态估计和相应的均方误差阵、新息序列；

系统的状态向量通过卡尔曼滤波传感器信息融合来估计，表示为：

$\hat{x}(k/k)=\hat{x}(k/k-1)+\underset{i=1}{\overset{m}{\&Sigma;}}P(k/k)H_i^T\left(k\right)R_{ii}^{-1}\left(k\right){\mathrm{\&Delta;}}_i\left(k\right)---\left(3\right)$

式中：是系统的状态估计，P(k/k)表示估计均方误差的相关矩阵，Δ_i(k)表示i通道的新息序列，是状态一步预测：

$\hat{x}(k/k-1)=\mathrm{\&Phi;}(k,k-1)\hat{x}(k-1)---\left(4\right)$

${\mathrm{\&Delta;}}_i\left(k\right)=z_i\left(k\right)-H_i\left(k\right)\hat{x}(k/k-1)---\left(5\right)$

估计均方误差的相关矩阵：

$P^{-1}(k/k)=P^{-1}(k/k-1)+\underset{i=1}{\overset{m}{\&Sigma;}}{H}_i^T\left(k\right)R_{11}^{-1}\left(k\right)H_i\left(k\right)---\left(6\right)$

一步预测均方误差的相关矩阵：

P(k/k-1)＝Φ(k,k-1)P(k-1/k-1)Φ^T(k,k-1)+G(k,k-1)Q(k-1)G^T(k,k-1)(7)

步骤三：利用已知的新息序列，得到不同的渠道归一化的新息序列，并且组成m通道平行传感器的创新矩阵；

对于正常运作的卡尔曼滤波传感器信息融合，不同的通道归一化的新息序列为：

${\widetilde{\mathrm{\&Delta;}}}_i\left(k\right)={\&lsqb;H_i\left(k\right)P(k/k-1)H_i^T\left(k\right)+R_{ii}\left(k\right)\&rsqb;}^{-1/2}{\mathrm{\&Delta;}}_i\left(k\right)---\left(8\right)$

其中，是第i个通道的归一化新息序列，服从N(0,1)分布；

建立两种假设：γ₀-Kalman滤波器工作在正常状态，γ₁-估计系统发生故障；

设定A是一个n×m矩阵且为m通道Kalman滤波器规范化创新矩阵，n≥2,m≥2，其中，各列为不同通道的规范化创新向量，规范化创新矩阵为：

$A\left(k\right)=\&lsqb;{\widetilde{\mathrm{\&Delta;}}}_1\left(k\right),{\widetilde{\mathrm{\&Delta;}}}_2\left(k\right),...{\widetilde{\mathrm{\&Delta;}}}_m\left(k\right)\&rsqb;---\left(9\right)$

其中，A(k)是k时刻的m通道平行传感器信息融合kalman滤波器的规范化创新矩阵，是各通道k时刻的归一化新息序列；对于假设检验γ₀和γ₁，矩阵A^T(k)A(k)服从白噪声分布；

步骤四：根据步骤三所得的规范化创新矩阵，获取规范化创新矩阵的谱范数和谱范数的均值；

规范化创新矩阵A(k)的谱范数||||₂为：

||A(k)||₂≡max{(λ_i[A^T(k)A(k)])^1/2}(10)

式中：λ_i[A^T(k)A(k)]是A^T(k)A(k)的特征值；

A(k)的谱范数均值为：

$E\left\{\right|\left|A\left(k\right)\right|{|}_2\}\&ap;{\stackrel{\&OverBar;}{\left|\right|A\left(k\right)\left|\right|}}_2=\frac{1}{k}\underset{j=1}{\overset{k}{\&Sigma;}}|\left|A\left(j\right)\right|{|}_2---\left(11\right)$

假设r_i和a_j代表规范化创新矩阵A的行和列，最大行-列范数为：

μ≡max[||r_i||₂,||a_j||₂](12)

其中||r_i||和||a_j||是向量r_i和a_j的范数，并且有下列不等式成立：

E{μ}≤E{||A||₂}≤[max(n,m)]^1/2E{μ}(13)

把E{μ}用它的下界来替换：

$\mathrm{\&sigma;}\sqrt{max(n,m)}=max\&lsqb;E\left\{\right|\left|r_i\right|{|}_2\},E\{\left|\right|a_j\left|\right|2\}\&rsqb;\&le;E\{\mathrm{\&mu;}\}---\left(14\right)$

所以公式(13)改为：

$\mathrm{\&sigma;}\sqrt{max(n,m)}\&le;E\left\{\right|\left|A\right|{|}_2\}\&le;f\left(max\right(n,m\left)\right)\&times;\mathrm{\&sigma;}\sqrt{max(n,m)}---\left(15\right)$

式中f为未知函数，并且有在n＝m→∞时函数f渐近为2且f的区间是(1,2)，所以定义f的估计值为2，针对以上所述，E{||A(k)||₂}有以下定义：

$\mathrm{\&sigma;}\sqrt{max(n,m)}\&le;E\left\{\right|\left|A\right|{|}_2\}\&le;2\mathrm{\&sigma;}\sqrt{max(n,m)}---\left(16\right)$

表达式(16)表现了规范化创新矩阵A标准差和它的二阶范数之间的关系；

规范化创新矩阵A标准差和它的二阶范数之间的关系通过卡尔曼滤波器的归一化新息矩阵的元素a_ij是否服从N(0,1)分布来确定，并要满足不等式(16)；不满足不等式(16)则说明元素a_ij的均值不再为零或者方差不再为单位方差或者{a_ij}不再服从白噪声；

鉴于σ＝1和表达式(11)就得到谱范数均值符合以下关系形式：

$\sqrt{max(n,m)}\&le;{\stackrel{\&OverBar;}{\left|\right|A\left(k\right)\left|\right|}}_2\&le;2\sqrt{max(n,m)}---\left(17\right)$

步骤五：对卡尔曼滤波传感器的信息融合进行故障检测；

测试传感器信息融合卡尔曼滤波器的故障检测问题，决策方法为以下形式：

γ₀: $\sqrt{m\mathrm{\&alpha;}(n,m)}<{\stackrel{\&OverBar;}{\left|\right|A\left(k\right)\left|\right|}}_2<2\sqrt{max(n,m)},\&ForAll;k:$ 无故障运行

γ₁: ${\stackrel{\&OverBar;}{\left|\right|A\left(k\right)\left|\right|}}_2\&le;\sqrt{max(n,m)}$ 或 ${\stackrel{\&OverBar;}{\left|\right|A\left(k\right)\left|\right|}}_2\&GreaterEqual;2\sqrt{max(n,m)},\&Exists;k:$ 故障(18)。