CN1121093C - 一种动态系统中自适应卡尔曼滤波方法 - Google Patents

Abstract

本发明是根据将兰格的快速卡尔曼滤波(FKF^TM)原理，用于大型过程控制、预测或报警系统，在这种过程或系统中如果用其他计算方法的话不是太慢就是因为舍位误差而导致计算失败。本发明的方法使得在具有大的移动数据窗口的动态多参数系统的自适应卡尔曼滤波中采用快速卡尔曼滤波方法成为可能。

Description

一种动态系统中自适应卡尔曼滤波方法

技术领域

本发明主要涉及所有的卡尔曼(KALMAN)滤波的实际应用，尤其是对需要快速可靠地适应环境的动态系统的控制。

现有技术

在解释此发明前，有必要先了解一下传统的卡尔曼(KALMAN)递推公式的原有技术及用于校正传感器系统PCF/F190/00122(WO90/13794)及控制一个大的动态系统PCT/F193/00192(WO 93/22625)的快速卡尔曼滤波方法(FKF)

下面的马尔可夫(MARKOV)(有限记忆)过程是通过方程式(1)至(3)描述的。第一个方程是关于量测矢量Y_t与处于时间点t₁(t＝0，1，2...)上的状态矢量S_t的关系，这是一个线性化的量测(或观察)方程

y_t＝H_ts_t+e_t                                       (1)

矩阵H_t是设计(雅可比)矩阵，它是从实际物理关系的偏导中产生出来的。第二个方程描述了整个系统时间进展，它作为一个线性化的系统(或状态)方程

s_t＝A_ts_t+B_t+u_t-1+a_t                             (2)

矩阵A_t是状态转移(雅可比)矩阵，B_t是控制增益(雅可比)矩阵。方程式(2)描述了整个系统的当前状态S_t是如何从前一状态S_t-l，控制/强制u_t-l及随机误差a_t推导出来的。当测量误差e_t和系统误差既不是自动的(即白噪声)也不是相互关联的而是由下列协方差矩阵给出的时：

Q_t＝Cov(a_t)＝E(a_ta_t ^*)

              and                                (3)

R_t＝Cov(e_t)＝E(e_te_t ^*)

这样著名的卡尔曼正向递推公式(4)至(6)给出我们当前状态S_t是最佳线性无偏估计(BLUE) 如下： ${\hat{s}}_t=A_t{\hat{s}}_{t-1}+B_t{u}_{t-1}+K_t\{y_t-H_t(A_t{\hat{s}}_{t-1}+B_t{u}_{t-1})\}----\left(4\right)$ 及其估计误差协方差矩阵如下： ${\hat{P}}_t=\mathrm{Cov}\left({\hat{s}}_t\right)=E\{({\hat{s}}_t-s_t)({\hat{s}}_t-s_t)\text{'}\}=A_t{\hat{P}}_{t-1}{A}_t^{\text{'}}+Q_t-K_t{H}_T(A_t{\hat{P}}_{t-1}{A}_t^{\text{'}}+Q_t)---\left(5\right)$ 式中卡尔曼增益矩阵K_t由下列公式确定 $K_t=(A_t{\hat{P}}_{t-1}{A}_t^{\text{'}}+Q_t)H_t^{\text{'}}\{H_T(A_t{\hat{P}}_{t-1}{A}_t^{\text{'}}+Q_t)H_t^{\text{'}}+R_t{\}}^{-1}---\left(6\right)$

此递推线性结果是局部最佳的，卡尔曼滤波器(KF)的稳定性要求其可观测性及可控制性条件必须得到满足(卡尔曼，1960)然而在方程(6)中经常要求一个太大的矩阵逆变，矩阵行(及列)的数目n与测量矢量Y_t中的元素一样多，这需要一个大致n以满足可观测性及可控制性条件，这就是本发明及PCF/F190/C0122PCT/F193/00192所要解决的问题。

下列状态方程修正式导出 $A_t{\hat{s}}_{t-1}+B_t{u}_{t-1}={\mathrm{Is}}_t+A_t({\hat{s}}_{t-1}-s_{t-1})-a_t---\left(7\right)$ 并与测量方程(1)结合以得到所说的增量模型： $\left[\begin{array}{c}{y}_t\\ A_t{\hat{s}}_{t-1}+B_t{u}_{t-1}\end{array}\right]=\left[\begin{array}{c}{H}_t\\ I\end{array}\right]s_t\left[\begin{array}{c}{e}_t\\ A_t({\hat{s}}_{t-l}-s_{t-l})-a_t\end{array}\right]---\left(8\right)$ $i.e.z_t=z_t{s}_t+{\mathrm{\η}}_t$ 即：状态参数可通过使用人们所熟知的如下回归分析法计算出。 ${\hat{s}}_t={\left(Z_t^{\text{'}}{V}_t^{-1}{Z}_t\right)}^{-1}{Z}_t^{\text{'}}{V}_t^{-1}{z}_t---\left(9\right)$

该结果在代数上等同于使用卡尔曼递推式，但在数值上不相等(参见：Harvey，1981：“Time Series Models”，Philip allanPublishers Ltd，Oxford，UK，(“时代”系列模型飞利蒲Allan出版社，牛津，)PP 101-119)。此时方程(9)中可逆变矩阵的维数等于状态矢量St中元素的数目(m)。Harvey的计算方法是快速卡尔曼滤波(FKF)的各种不同的变形的基础。

为了满足可观察性条件的需要，任何大的卡尔曼滤波器(KF)的初始化或临时排序可通过Lange(兰格)的高通滤波器而完成(Lange 1988)。它利用一个解析稀疏矩阵逆变公式来解决带有下列所谓标准的分块对角矩阵结构的回归模型。 $\left[\begin{array}{c}{y}_1\\ y_2\\ .\\ .\\ y_K\end{array}\right]=\left[\begin{array}{cccccc}{X}_1& & & & & G_1\\ & X_2& & & & G_2\\ & & .& & & .\\ & & & .& & .\\ & & & & X_K& G_K\end{array}\right]\left[\begin{array}{c}{b}_1\\ .\\ .\\ b_K\\ c\end{array}\right]+\left[\begin{array}{c}{e}_1\\ e_2\\ .\\ .\\ e_K\end{array}\right]---\left(10\right)$

这就是例如完整的风口测量相互比较实验的量测方程的矩阵表示，矢量b₁，b₂...b_k为连续的位置座标(如气象气球的座标)，但它也可包括那些带有一个显著的时间或空间偏差的校正参数。矢量e为其他的在抽样时间内为常数的校正参数。

对于所有大的多重传感器系统，它们的设计矩阵H_t是稀疏的，这样人们可以一种或其他同类的矩阵分块方法来完成： $s_t=\left[\begin{array}{c}{b}_{t,1}\\ .\\ .\\ b_{t.K}\\ c_t\end{array}\right]y_t=\left[\begin{array}{c}{y}_{t,1}\\ y_{t,2}\\ .\\ .\\ y_{t,K}\end{array}\right]H_t=\left[\begin{array}{cccccc}{X}_{t,1}& & & & & G_{t,1}\\ & X_{t,2}& & & & G_{t,2}\\ & & .& & & .\\ & & & .& & .\\ & & & & X_{t,K}& G_{t,K}\end{array}\right]$ $A=\left[\begin{array}{ccccc}{A}_1& & & & \\ & .& & & \\ & & .& & \\ & & & A_K& \\ & & & & A_C\end{array}\right]$ and， $B=\left[\begin{array}{ccccc}{B}_1& & & & \\ & .& & & \\ & & .& & \\ & & & B_K& \\ & & & & B_C\end{array}\right]---\left(11\right)$ 式中c₁表示t时刻校正参数

b_t，k表示在时间或空间中所有其他的状态参数

A表示t时刻状态变化矩阵(分块对角矩阵)

B表示t时刻状态独立作用V_t矩阵(分块对角矩阵)

若矩阵分块是不明显的，人们可以通过使用一个特定的算法自动解决，这个算法把每一个稀疏的线性系统逆变为标准分块对角形式(参见：用于分析及其它算法的矩阵块角形式和再安排，管理科学，18卷1号，1971年9月98-107页(Weil and kettler，1971：“Rearranging Matrices to  Block-angular Form forDecomposition(and other)Algorithms”，Management Science，Vol.18，No.1，Semptember 1971，Pages 98-107.))然而，随机误差et的协方差矩阵可解会使原来简单的对角性变得不严格。

结果我们将面临如下这个难解回归分析问题：

一个用于空间量的增广模型，例如：对于一组包含K个连续气球位置的气球跟踪实验数据：

一个用于动态时间量的增广模型：(例如：用于“白化”一个所观测到的长度为L的动态采样的偏差et的“修正”序列)：

请注意后一个矩阵方程具有一个“嵌套”分块对角结构，有两种形式“校正”参数11C_t及C_t，这些参数的第一组C_t可以随时间的变化而变化，而这些参数的第二组Ct在长度为L的移动时间窗口范围内为常量(至少近似为常量)，后一种参数C_t使卡尔曼滤波过程能够自适应。在从(4)至(6)的传统卡尔曼递推式中求出后一种参数时产生导致一个可观测性问题，就计算原因而言，长度L必须是短的，但是对于PCF/FI90/00122的FKCF公式，采样尺度可能很大以至于根本不需要初始化(或排序)。

在解释PCT/FI193/00192方法之前，有必要先了解一些应用于实验性数据天气预测(NWP)系统中的卡尔曼滤波器理论的现有技术。因为以前他们也利用方程(1)：量测方程：y_t＝H_ts_t+e_t(线性化回归)

式中状态矢量S_t描述的是t时刻的大气状态，现在S_t通常代表所有大气变化的格点值，例如：不用压力水平上的位势高度(实际上它们与实际数值的偏差量可以通过某些方法进行估计)。

空气动力学是由一个人们熟知的偏微分方程(原始方程)决定的，通过利用例如NWP模型的切线逼近得出，下列方程(2)的表达式是用于计算在某一时间段上状态参数S_t随时间推移的变化(实际上，他们相对于参数空间轨迹的偏移量产生于非线性NWP模型)。状态方程：s_t＝As_t-1+Bu_t-1+a_t(离散的动态随机模型)

四维的数据同化结果 及NWP预测 分别从卡尔曼滤波器系统得到： ${\hat{s}}_t={\tilde{s}}_t+K_t(y_t-H_t{\tilde{s}}_t)$ ${\tilde{s}}_t=A{\hat{s}}_{t-1}+{\mathrm{Bu}}_{t-1}---\left(12\right)$ 式中 $P_t=\mathrm{Cov}\left({\tilde{s}}_t\right)=\mathrm{ACov}\left({\hat{s}}_{t-1}\right)A\text{'}+Q_t$ ...(预测精度)Q_t＝Cov(a_t)＝Ea_ta’_t             ...(系统噪声)R_t＝Cov(e_t)＝Ee_te’_t             ...(量测噪声)关键性的校正计算是用下列卡尔曼递推式来完成的K_t＝P_t H’_t(H_tP_tH’_t +R_t)^-1    ...(卡尔曼增益矩阵) $\mathrm{Cov}\left({\hat{s}}_t\right)=P_t-K_t{H}_t{P}_t$     ...(计算精度)

这里所需的用于卡尔曼增益矩阵计算的矩阵逆变极难用于实际NWP系统计算，因为数据固化系统必须在同一时刻处理几百万个数据元素，T.GaL-Chen博士在此问题上于1988年报告如下：“希望大型并行超级计算机的开发(例如：1000台台式计算机CRAYS协力工作)能够使计算更接近于最佳....)”见“重要回顾小组的报告一低对流层分布论文：需要及技术”美国气象科学文摘71卷5号，1990年5月684页。

用PCT/FI93/00192中的方法从方程(8)推出增广模型算法 $\left[\begin{array}{c}{y}_t\\ A_t{\hat{s}}_{t-1}+B{u}_{t-1}\end{array}\right]=\left[\begin{array}{c}{H}_t\\ I\end{array}\right]s_t+\left[\begin{array}{c}{e}_t\\ A_t({\hat{s}}_{t-l}-s_{t-l})-s_t\end{array}\right]$ $z_t=Z_t{s}_t+{\mathrm{\η}}_t$ 即下列两组方程用于校正的目的

               …(最佳估计，

                 由高斯一马尔可夫得出) ${\hat{s}}_t={\left(Z_t^{\text{'}}{V}_t^{-1}{Z}_t\right)}^{-1}{Z}_t^{\text{'}}{V}_t^{-1}{z}_t$ $={\{H_t^{\text{'}}{R}_t^{-1}{H}_t+P_t^{-1}\}}^{-1}(H_t^{\text{'}}{R}_t^{-1}{y}_t+P_t^{-1}{\tilde{s}}_t)---\left(13\right)$ 或： $={\tilde{s}}_t+K_t(y_t-H_t{\tilde{s}}_t)$ …(用另一种形式)及： $\mathrm{Cov}\left({\hat{s}}_t\right)=E({\hat{s}}_t-s_t)({\hat{s}}_t-s_t)\text{'}={\left(Z_t^{\text{'}}{V}^{-1}{Z}_T\right)}^{-1}---\left(14\right)$ $={\{H_t^{\text{'}}{R}_t^{-1}{H}_t+P_t^{-1}\}}^{-1}$ …(估计精度)其中 ${\tilde{s}}_t=A{\hat{s}}_{t-1}+{\mathrm{Bu}}_{t-1}$ …(NWP“预测”) $P_t=\mathrm{Cov}\left({\tilde{s}}_t\right)=\mathrm{ACov}\left({\hat{s}}_{t-1}\right)A\text{'}+Q_t---\left(15\right)$ 而代替K_t＝P_tH’_t(H_tP_tH’_t+R_t)^-1…(卡尔曼增益矩阵)PCT/F193/00192中的FKF方法用 $K_t=\mathrm{Cov}\left({\hat{s}}_t\right)H_t^{\text{'}}{R}_t^{-1}---\left(16\right)$ 对于一个大的输入数据Y_t矢量，增广模型算法优于传统卡尔曼递推式，因为当Cov 是未知的时，卡尔曼增益矩阵K_t的计算要示对大的矩阵求逆，两种方法在代数上或统计学上是等价的，但在数值不等。

然而，增广的模型公式仍难于在数字上解决，首先，这是由于状态矢量S_t中包括用于对大气进行实际描述的矢量(＝m)格点数据，其次，对于一个实际的NWP系统，许多其它的状态参数必须被包含在状态矢量中。首先是有关观察系统的系统(校正)误差及所谓的小规模大气过程的物理参数表。

在PCT/F193/00192中校正问题是通过使用去耦合状态的方法来解决的。通过进行下列分块而完成： $s_t=\left[\begin{array}{c}{b}_{t,1}\\ .\\ .\\ b_{t.K}\\ c_t\end{array}\right]y_t=\left[\begin{array}{c}{y}_{t,1}\\ y_{t,2}\\ .\\ .\\ y_{t,K}\end{array}\right]H_t=\left[\begin{array}{cccccc}{X}_{t,1}& & & & & G_{t,1}\\ & X_{t,2}& & & & G_{t,2}\\ & & .& & & .\\ & & & .& & .\\ & & & & X_{t,K}& G_{t,K}\end{array}\right]$

             and                                    (17) $A_t=\left[\begin{array}{c}{A}_{t,1}\\ .\\ .\\ A_{t,K}\\ A_{t,c}\end{array}\right]$ and， $B_t=\left[\begin{array}{c}{B}_{t,1}\\ .\\ .\\ B_{t,K}\\ B_{t,c}\end{array}\right]$ 式中：C_t表示t时刻的校正参数

  C_t表示格点k(k1，...k)处大气参数值

  A表示t时刻状态转移矩阵(子矩阵A₁....A_k，A_c)

  B表示控制增益矩阵(子矩阵B₁，...，B_k，B_c)结果  面临下列艰巨的回归分析问题在任何时间点t上的递推快速卡尔曼滤波器(FKF)公式如下： ${\hat{b}}_{t,k}={\left\{X_{t,K}^{\text{'}}{V}_{t,k}^{-1}{X}_{t,K}\right\}}^{-1}{X}_{t,K}^{\text{'}}{V}_{t,K}^{-1}(y_{t,k}-G_{t,k}{\hat{c}}_t)\mathrm{fork}=\mathrm{1,2},...,K$ ${\hat{c}}_t={\left\{\underset{k=0}{\overset{K}{\mathrm{\Σ}}}{G}_{t,k}^{\text{'}}{R}_{t,k}{G}_{t,k}\right\}}^{-1}\underset{k=0}{\overset{K}{\mathrm{\Σ}}}{G}_{t,k}^{\text{'}}{R}_{t,k}{y}_{t,k}---\left(19\right)$ 式中：设k＝1，2，....，k， $R_{t,k}=V_{t,k}^{-1}\{I-X_{t,k}{\left\{X_{t,k}^{\text{'}}{V}_{t,k}^{-1}{X}_{t,k}\right\}}^{-1}{X}_{t,k}^{\text{'}}{V}_{t,k}^{-1}\}$ $V_{t,k}=\left[\begin{array}{cc}\mathrm{Cov}\left(e_{t,k}\right)& \\ & \mathrm{Cov}\{A_k({\hat{s}}_{t-1}-s_{t-1})-a_{b_{t,k}}\}\end{array}\right]$ $y_{t,k}=\left[\begin{array}{c}{y}_{t,k}\\ A_k{\hat{s}}_{t-1}+B_k{u}_{t-1}\end{array}\right]$ $X_{t,k}=\left[\frac{X_{t,k}}{I}\right]$

且，即设k＝0 $R_{t,0}=V_{t,0}^{-1}$ $V_{t,0}=\mathrm{Cov}\{A_c({\hat{s}}_{t-1}-s_{t-1})-a_{c_t}\}$ $y_{t,0}=A_c{\hat{s}}_{t-1}+B_c{u}_{t-1}$ G_t，0＝I.从方程(20)得到下列数据同化精度 $\mathrm{Cov}\left({\hat{s}}_t\right)=\mathrm{Cov}({\hat{b}}_{t,1},...,{\hat{b}}_{t,K},{\hat{c}}_t)---\left(20\right)$ 式中 $C_k={\left\{X_{t,k}^{\text{'}}{V}_{t,k}^{-1}{X}_{t,k}\right\}}^{-1}$   设k＝1，2，....，k， $D_k={\left\{X_{t,k}^{\text{'}}{V}_{t,k}^{-1}{X}_{t,k}\right\}}^{-1}{X}_{t,k}^{\text{'}}{V}_{t,k}^{-1}{G}_{t,k}$   设k＝1，2，....，k， $s={\left\{\underset{k=0}{\overset{K}{\mathrm{\Σ}}}{G}_{t,k}^{\text{'}}{R}_{t,k}{G}_{t,k}\right\}}^{-1}$

卡尔曼滤波器(KF)研究也曾被报道过，例如：StephenE.Cohn和David F.Parrish(斯蒂芳E考恩及大卫·F帕力士)(1991)：“二维卡尔曼滤波器预测误差协方差特点”美国气象科学每月天气回顾，119卷，1757至1785页。然而，对于四维(即空间及时间)理想的卡尔曼滤波器系统，在所有这些报告中都没有涉及到。这样就需要一个可靠的对状态参数误差协方差距阵进行估计和逆变的方法，如：中部区域天气预报欧洲中心(ECMWF)的He；kk1Jarvinen博士所述：“在气象学中，状态参数矢量St的维数(＝m)可能是100,000-10,000,000。这使得在实际中不可能准确处理，误差协方差矩阵，”见“作为一个派生问题的气候学数据同化”报告第43号(1995)气候系Helsink大学第10页。中部区域天气预报欧洲中心的Adrian Simmons博士确认说：“在理论上卡尔曼滤波基本算法是很完善的，但是计算上的要求使得整个执行过程难以处理，见中部区域天气预报欧洲中心69号简报(1995年春)第12页。

从PCF/F190/00122及PCT/F193/00192中得知的快速卡尔曼滤波器(FKF)公式利用这样的假设：方程(9)及(13)中的误差矩阵分别是分块对角的，参见FKF公式(19)式中这些对角块被表述如下： $V_{t,k}=\left[\begin{array}{cc}\mathrm{Cov}\left(e_{t,k}\right)& \\ & \mathrm{Cov}\{A_k({\hat{s}}_{t-1}-s_{t-1})-a_{b_{t,k}}\}\end{array}\right]$ 尤其是对于自适应卡尔曼滤波(及4维数据比较)情况，很显然连续状态参数矢量S_t-1，S_t-2，S_t-2的估计值是相互且自动相关的。

有必要将快速卡尔曼滤波方法的原理应用于自适应卡尔曼滤波(AKF)，比起其它卡尔曼滤波方法，它具有同样或更快的运算速度、可靠性、准确度和低成本。在此将公开本发明的一个处理误差协方差的具体方法。

发明概要

通过提供一种用于校正/调整实时或接近实时的动态系统的各种参数的自适应快速卡尔曼滤波法可以基本满足上述需要，测量误差和系统误差都是白化和部分正交化的。如详细说明中所描述：这种FKF滤波运算与在可观察和可控制条件下的最佳卡尔曼滤波相接近，预测误差的方差和协方差提供了一种监视滤波稳定性的办法。

根据本发明第一方面的一种用于调整传感系统的模型和/或校正参数的方法，该传感系统配备有在现实中发生的外部事件的所述模型，其中：

·所述系统的传感器输出单元响应所述外部事件提供信号；

·N大于50，在此N是来自在任一情况下进行处理的所述传感器输出单元的数值的总数；以及

·所述方法使用自适应快速卡尔曼滤波(FKF)，其包括步骤a)到g)：

a)提供一种用于储存信息的数据的数据库单元，这些信息是：

-为一些所述传感器提供多个测试点传感器的输出信号值和为对应于所述测试点传感器输出信号值的所述外部事件提供多个数值，和/或从相邻的传感器提供所述输出信号值的同一时间序列用于对照；

-对应于在时间t的情况的所述传感器输出信号值，所述模型和/或校正参数的值，以及所述外部事件的值；以及

-对应于所述在时间t的情况的所述传感器的控制以及/或者外部事件的强制；

b)提供一种逻辑单元，以存取所述传感器信号输出值和所述模型和/或校正参数值，所说逻辑单元有双路通道与所述数据库单元相连，以及具有按需要用兰格高通滤波器计算未知模型和/或校正参数的初始值的能力；

c)将可用的来自所述传感器的所述传感器输出信号值提供给所述逻辑单元；

d)将所述传感器的控制和/或所述外部事件的所述强制方面的信息提供给所述数据库单元；

e)访问所述模型和/或校正参数以及状态转换矩阵的当前值，用快速卡尔曼滤波器公式计算所述的模型和/或校正参数的更新值： ${\hat{s}}_{t-1}={\left\{{X_{t-l}}^{\′}{X_{t-l}}^{-1}{X}_{t-l}\right\}}^{-1}{X_{t-l}}^{\′}{V_{t-l}}^{-1}(y_{t-l}-G_{t-l}{\hat{c}}_t)$ ${\hat{c}}_t={\left\{\underset{l=0}{\overset{L}{\mathrm{\Σ}}}{G}_{t-l}^{\′}{R}_{t-l}{G}_{t-l}\right\}}^{-1}\underset{l=0}{\overset{L}{\mathrm{\Σ}}}{G}_{t-l}^{\′}{R}_{t-l}{y}_{t-l}$ 在此，＝对于时间块t-l的模型和/或校正参数的滤波的值，

＝对所述时间t的情况共同的，那些模型和/或校正参数的滤波值， $y_{t-1}=\left[\begin{array}{c}{y}_{t-l}\\ A_{t-l}{\hat{s}}_{t-l-1}+B_{t-l}{u}_{t-l-1}\end{array}\right]$ y_t-l＝对于时间块t-1的所述传感器输出信号值的矢量， $A{\hat{s}}_{t-l-1}+{\mathrm{Bu}}_{t-l-1}$ ＝对于时间块t-l的所述预期的模型和/或校正参数的矢量，A＝状态转换矩阵，

＝对于前一时间块t-l-1的模型和/或校正参数的滤波值，B＝控制增益矩阵，u_t-l-1＝对于前一时间块t-l-1的所述传感器的控制或所述外部事件的强制， $X_{t-l}=\left[\frac{H_{t-l}}{I}\right]$ ＝对于状态参数s_t-l的增广的雅可比矩阵，H_t-l＝对于状态参数s_t-l的雅可比矩阵，I＝单位矩阵， $G_{t-l}=\left[\frac{F_{t-l}^y}{F_{t-l}^s{M}_{t-l-1}}\right]$ ＝对于状态参数c_t的增广的雅可比矩阵，M_t-l-1＝用于校正从时间块t-l-1到时间块t-l的状态转换误差的系数矩阵， $R_{t-l}={V_{t-l}}^{-1}\{I-X_{t-l}{\left\{{X_{t-l}}^{\′}{V_{t-l}}^{-1}{X}_{t-l}\right\}}^{-1}{X_{t-l}}^{\′}{V_{t-l}}^{-1}\},$ V_t-l＝在y_t-l中数据的所述增加的矢量的误差协方差，l＝将在时间t所述状态使用的数据的L个时间块上时间块索引，其中改进包括增广模型(24)的误差协方差矩阵的块对角线化：

它是对于在时间t的所述状态，在所述逻辑单元中，用快速卡尔曼滤波公式发展所述公因子F^y或F^s的任何一个或所述系数M所描述的；

f)通过监控所述的逻辑单元中模型和/或校正参数的所述更新值的精度估计，并指示需要传感器输出信号值、测试点数据、传感器比较或系统重构，控制所述自适应快速卡尔曼滤波的稳定性；

g)如果可以稳定地更新，调整所述模型和/或校正参数。

根据本发明的另一方面，在上述调整传感系统的模型和/或校正参数的方法中，在步骤e)中，该快速卡尔曼滤波器公式是通过应用佛罗本尼变换公式(26)： ${\left[\begin{array}{cc}A& B\\ C& D\end{array}\right]}^{-1}=\left[\begin{array}{cc}{A}^{-1}+A^{-1}{\mathrm{BH}}^{-1}{\mathrm{CA}}^{-1}& -A^{-1}{\mathrm{BH}}^{-1}\\ -H^{-1}{\mathrm{CA}}^{-1}& H^{-1}\end{array}\right]$ 其中H＝D-CA^-1B获得的，其中A、B、C、D表示当对时间t计算增广模型(8)： $\left[\begin{array}{c}{y}_t\\ A_t{\hat{s}}_{t-1}+B_t{u}_{t-1}\end{array}\right]=\left[\begin{array}{c}{H}_t\\ I\end{array}\right]s_t\left[\begin{array}{c}{e}_t\\ A_t({\hat{s}}_{t-l}-s_{t-l})-a_t\end{array}\right]$ 的参数s_t的估值s_t时，标准方程的不同的分块系数矩阵。

本发明最佳实施例

我们重述一下线性量测(或观察)方程 $y_t=H_t{s}_t+F_t^y{C}_t+e_t---\left(21\right)$

式中e_t代表白噪声，它既不与e_t-1，e_t-2.....相关，也不与S_t-1，S_t-2相关，且不与a_t，a_t-1，a_t-2相关，矩阵H_t是与前述从测量值Y_t和状态参数S_t间物理关系的偏导中出的设计矩阵相同，见前面4页的矩阵分块(11)(以往对矩阵A和B的分块对角性假设不再有效)矩阵F_t ^y描述了测量的系统误差与校正或“校正型”参数，状态矢量C_t之间的关系，该矢量1t在时间上是常数或变化极慢。矩阵F_t ^y，F_t-1 ^y，F_t-2 ^y的列代表偏导数，类似于正弦波、矩形波的波形、“随机编码”函数等，以及根据已知物理关系、测量的系统误差的回归和自动回归(AR)确定的经验正交函数(EOF)。预测矢量c_t的元素将决定着红噪声系数的波幅，让我们看一下第5页上的一个类似的用于“自化”所观察到的“新的”测量序列的动态时间量增广模型。同样我们重述一下线性系统(或状态)方程： $s_t=(A_t+d{A}_t)s_{t-1}+B_t{u}_{t-1}+F_t^s{c}_t+a_t---\left(22\right)$ 式中a_t表示白噪声，既与e_t，e_t-1，e_t-2不相关，也与 不相关，且与a_t-1，a_t-2不相关。矩阵A_t同样是状态转移矩阵，它来自于状态St与前一状态St-1之间的物理关系的偏导。矩阵F_t ^s描述动态模型(如NWP)的系统误差与校正或“校正型”参数及矢量C_t的关系。其中G不随时间变化或随时间缓慢地变化。矩阵的列F_t ^s，F_t-1 ^s，F_t-2 ^s表示偏导数、如正弦波、矩形波这样的波型、随机编码函数等，及根据已知物理关系及模型的系统误差的回归和自动回归所确定的经验正交函数。估计矢量

的元素，将决定着红噪声系数的波幅。

矩阵dA_t描述了动态(NWP)模型的系统状态转移误差是如何与当前(气候)条件成比例的。如果他们是未知的，但变化极慢，可以结合FKF法通过动态求均值来调整。下面就介绍这种调整方法。从系统方程(22)中得到结果并重述如下：dA_ts_t-1＝＝[s₁I_(mxm)，s₂I_(mxm)，…，s_mI_(mxm)][da₁₁，da₂₂，…，da_m1，da₁₂，…，da_mm]’＝M_t-1r_t                                                               (23)式中m_t-1是一个由尺寸为m×m的m个对角线矩阵组成的矩阵S₁，S₂....S_m是状态矢量S_t-1的m标量元素r_t是矩阵dAt的m×m个元素的列矢量

请注意方程(23)颠倒了乘法顺序，这样可以象普通回归参数一样估计矩阵dA_t的元素。

结果：面临下面难解的回归分析问题：

一个动态时间量的增广模型(即用于对长度L动态抽样的偏差e_t和a_t的更新序列进行白化)

请注意上述矩阵方程具有一个嵌套分块对角结构，有三种不同形式的校正参数。第一种形式C_t是嵌入每个时间段t数据中。另两种形式由状态矢量C_t表示。第一组参数用于对测量误差和系统误差白化和部分正交化(即有关误差协方差矩阵的分块对角化)第二组(即r_t)用于校正状态转移矩阵的严重误差，最后两组参数在长的移动时间窗中或多或少保持为常量从而使卡尔曼滤波过程能够自适应。

同样请注意矩阵M不能象在方程(23)中那样取其最大尺寸(m×m²)。这是因为由于有太多的未知数，可观测性条件不能得到满足。因而矩阵M必须被有效地压缩以使其只表示与严重的转换误差相关的矩阵A_t中的那些元素。这种转化是通过利用所谓的最大相关法发现的。事实上偶然和缓慢位移的图形可以在状态参数矢量的空间中形成。这些一般是小范围的现象，而且它们不能被仅由模型方程导出的状态转移矩阵所充分描述，为了保证滤波稳定性，所有的矩阵dA₁的估计元素必须在测量的动态平均中保持可观察，如通过方程(20)中监控它们的估计方差。

在时刻t，用于长度为L的时间窗口的快速卡尔曼滤波(FKF)公式如下： ${\hat{s}}_{t-1}={\left\{X_{t-l}^{\text{'}}{V_{t-l}^{-1}X}_{t-l}\right\}}^{-1}{X}_{t-l}^{\text{'}}{V}_{t-l}^{-1}(y_{t-l}-G_{t-l}{\hat{c}}_t)$ 设1＝0，1，2，....，L-1 ${\hat{c}}_t={\left\{\underset{l=0}{\overset{L}{\mathrm{\Σ}}}{G}_{t-l}^{\′}{R}_{t-l}{G}_{t-l}\right\}}^{-1}\underset{l=0}{\overset{L}{\mathrm{\Σ}}}{G}_{t-l}^{\′}{R}_{t-l}{y}_{t-l}---\left(25\right)$ 式中设1＝0，1，2，....，L-1， $R_{t-l}=V_{t-l}^{-1}\{I-X_{t-l}{\left\{X_{t-l}^{\text{'}}{V}_{t-l}^{-1}{X}_{t-l}\right\}}^{-1}{X}_{t-l}^{\text{'}}{V}_{t-l}^{-1}\}$ $V_{t-l}=\left[\begin{array}{cc}\mathrm{Cov}\left(e_{t-l}\right)& \\ & \mathrm{Cov}\{A_{t-l}({\hat{s}}_{t-l-1}-s_{t-l-1})-a_{t-l}\}\end{array}\right]$ $y_{t-l}=\left[\begin{array}{c}{y}_{t-l}\\ A_{t-l}{\hat{s}}_{t-l-1}+B_{t-l}{u}_{t-l-1}\end{array}\right]$ $X_{t-l}=\left[\frac{H_{t-l}}{I}\right]$ 且，即设1＝L， $R_{t-L}=V_{t-L}^{-1}$ $V_{t-L}=\mathrm{Cov}\{A_c({\hat{C}}_{t-1}-C_{t-1})-a_{c_t}\}$ $y_{t-L}=A_c{\hat{C}}_{t-1}+B_c{u}_{c_{t-1}}$ $G_{t-L}=I.$ 为了最优化，有时有必要对滤波器的一些误差项进行具体化。如果这样做了恒方程(I)矩阵将从FKF公式中消失并必须被相应地代换。

此处的和在PCF/F190/00122和PCT/F193/00192中的FKF公式是基于假设误差协方差矩阵是分块对角的。因为运算上的限制禁止采用足够长的窗口，由于严重的可观察性和可控制性难题，要用传统的卡尔曼递推式(4)-(6)求得所有参数C_t，注定要失败。幸运的是，通过运用FKF公式，时间窗口可采用足够的长度以致于滤波的初始化或暂时排序变得完全多余了。

各种快速自适应卡尔曼滤波可以通过使用不同的佛罗本尼公式递归从大型线性回归方程(24)的标准方程系统中推导出： ${\left[\begin{array}{cc}A& B\\ C& D\end{array}\right]}^{-1}=\left[\begin{array}{cc}{A}^{-1}+A^{-1}{\mathrm{BH}}^{-1}{\mathrm{CA}}^{-1}& -A^{-1}{\mathrm{BH}}^{-1}\\ -H^{-1}{\mathrm{CA}}^{-1}& H^{-1}\end{array}\right]---\left(26\right)$

式中H＝D-CA^-1B。公式(20)和(25)及其它任何从佛罗本尼公式(26)得出的广大F型的公式都是依据本发明的方法。

例如，对于逆变对称的带状对角矩阵有些有效的运算方法。许多天气预报的误差协方差矩阵是典型的带状对角式的。我们可以直接从方程系(8)中得出而无需将状态参数S合并成大型回归分析式(18)的观察块。它们的逆差协方差矩阵可被作为一个大块而逆变，佛罗本尼公式的一个递归应用，使得FKF公式与公式(25)极相似。

所有用于解决大型回归分析模型而逆变的矩阵通过运用部分分解运算方法，保持足够地小。本发明的最佳实施例如附图1所示并描述如下：

一个以笔记本电脑为基础的超级导航仪通过应用推广的快速卡尔曼滤波(FKF)方法完成卡尔曼滤波逻辑单元1的功能。所有的接收部分原理包括一个积分传感器、遥感、数据处理及传输系统3，比方说，一个国家气象/海洋机构和任意一个现成的全球定位系统接收器。运行在笔记本电脑上的数据数据库单元2中包括控制输入4中和不同子系统的执行方向的更新信息及辅助信息，如地图。基于所有这些输入，逻辑单元1提供了实时三维图形，通过用FKF递归方程系(24)展示正在发生的情况和通过用方程(15)预测不远的将来将发生的情况。当众所周知的最佳卡尔曼过滤稳定条件由观察系统(3)满足了，可靠的精度信息也就得到。这些误差的方差和协方差由方程(15)和(20)计算。集中数据处理系统3提供了每一时刻t状态转换矩阵A的预测值。

这些矩阵因些被局部调整轨迹以将观察到出现在空气/海洋环境的小规模变化考虑进去。(见例如，Cotton Thompson和Mielde，1994“工作站实时中规模预测”《美国气象学文摘》75卷，第3号1994，3第349-362页。

对于本领域的技术人员还可以对本发明作各种改动，但都没脱离本发明的精神实质。因而本发明的应用范围不限于所述的特定实施例，而是以权利要求书所规定的范围为准。

参考文献：

(1)卡尔曼RE(1960)：“一种线性滤波的新方法和预测问题”，基础工程学美国机械工程师学报，82：35-45。

(2)兰格A A(1982)：“甚低频奥米伽信号的多径传播”，IEEE PLANS’82-定位和导航会议录，1982年12月，302-309。

(3)兰格A A(1984)“测风设备的集成、校正和相互比较”，世界气象组织仪表和观测方法报告15。

(4)兰格A A(1988a)：“用于观测系统最优化校正的高通滤波器”，大系统的仿真和最优化，奥夏达斯 A J编著，牛津大学出版社/牛津大学出版部印刷所，1988，311-327。

(5)兰格A A(1988b)：“使用卫星辐射亮度测量来确定无线电探空仪偏差的方法”，世界气象组织仪表和观测方法报告33，201-206。

(6)兰格A A(1990)：“用于校正传感器系统的设备和方法”，专利合作条约国际专利公开，世界知识产权组织，国际局，WO90/13794，PCT/FI90/00122，1990年11月15日。

(7)兰格A A(1993)：“大动态系统中的快速卡尔曼滤波方法”，专利合作条约国际专利公开，世界知识产权组织，国际局，WO93/22625，PCT/FI93/00192，1993年11月11日。

(8)兰格A A(1993)：“基于表面的高空探测系统”，世界气象组织仪表和观测方法报告57，175-177。

Claims (2)

1、一种用于调整传感系统的模型和/或校正参数的方法，该传感系统配备有在现实中发生的外部事件的所述模型，其中：

·所述系统的传感器输出单元响应所述外部事件提供信号；

·N大于50，在此N是来自在任一情况下进行处理的所述传感器输出单元的数值的总数；以及

·所述方法使用自适应快速卡尔曼滤波(FKF)，其包括步骤a)到g)：

a)提供一种用于储存信息的数据的数据库单元，这些信息是：

-为一些所述传感器提供多个测试点传感器的输出信号值和为对应于所述测试点传感器输出信号值的所述外部事件提供多个数值，和/或从相邻的传感器提供所述输出信号值的同一时间序列用于对照；

-对应于在时间t的情况的所述传感器输出信号值，所述模型和/或校正参数的值，以及所述外部事件的值；以及

-对应于所述在时间t的情况的所述传感器的控制以及/或者外部事件的强制；

b)提供一种逻辑单元，以存取所述传感器信号输出值和所述模型和/或校正参数值，所说逻辑单元有双路通道与所述数据库单元相连，以及具有按需要用兰格高通滤波器计算未知模型和/或校正参数的初始值的能力；

c)将可用的来自所述传感器的所述传感器输出信号值提供给所述逻辑单元；

d)将所述传感器的控制和/或所述外部事件的所述强制方面的信息提供给所述数据库单元；

e)访问所述模型和/或校正参数以及状态转换矩阵的当前值，用快速卡尔曼滤波器公式计算所述的模型和/或校正参数的更新值： ${\hat{s}}_{t-1}={\left\{{X_{t-l}}^{\′}{X_{t-l}}^{-1}{X}_{t-l}\right\}}^{-1}{X_{t-l}}^{\′}{V_{t-l}}^{-1}(y_{t-l}-G_{t-l}{\hat{c}}_t)$ ${\hat{c}}_t={\left\{\underset{l=0}{\overset{L}{\mathrm{\Σ}}}{G}_{t-l}^{\′}{R}_{t-l}{G}_{t-l}\right\}}^{-1}\underset{l=0}{\overset{L}{\mathrm{\Σ}}}{G}_{t-l}^{\′}{R}_{t-l}{y}_{t-l}$ 式中设1＝0，1，2......，L-1，在此，＝对于时间块t-l的模型和/或校正参数的滤波的值，＝对所述时间t的情况共同的，那些模型和/或校正参数的滤波值， $y_{t-1}=\left[\begin{array}{c}{y}_{t-l}\\ A_{t-l}{\hat{s}}_{t-l-1}+B_{t-l}{u}_{t-l-1}\end{array}\right]$ y_t-l＝对于时间块t-1的所述传感器输出信号值的矢量， ${A\hat{s}}_{t-l-1}+{\mathrm{Bu}}_{t-l-1}$ ＝对于时间块t-l的所述预期的模型和/或校正参数的矢量，A＝状态转换矩阵，

＝对于前一时间块t-l-1的模型和/或校正参数的滤波值，B＝控制增益矩阵，u_t-l-1＝对于前一时间块t-l-1的所述传感器的控制或所述外部事件的强制， $X_{t-l}=\left[\frac{H_{t-l}}{I}\right]$ ＝对于状态参数s_t-l的增广的雅可比矩阵，H_t-l＝对于状态参数s_t-l的雅可比矩阵，I＝单位矩阵， $G_{t-l}=\left[\frac{F_{t-l}^y}{F_{t-l}^s{M}_{t-l-1}}\right]$ ＝对于状态参数c_t的增广的雅可比矩阵，M_t-l-1＝用于校正从时间块t-l-1到时间块t-l的状态转换误差的系数矩阵，R_t-l＝V_t-l ^-1{I-X_t-l{X_t-l′V_t-1 ^-1X_t-l}^-1X_t-l′V_t-l ^-1}，V_t-l＝在y_t-l中数据的所述增加的矢量的误差协方差，l＝将在时间t所述状态使用的数据的L个时间块上时间块索引，其中改进包括增广模型(24)的误差协方差矩阵的块对角线化：

它是对于在时间t的所述状态，在所述逻辑单元中，用快速卡尔曼滤波公式发展所述公因子F^y或F^s的任何一个或所述系数M所描述的；

f)通过监控所述的逻辑单元中模型和/或校正参数的所述更新值的精度估计，并指示需要传感器输出信号值、测试点数据、传感器比较或系统重构，控制所述自适应快速卡尔曼滤波的稳定性；

g)如果可以稳定地更新，调整所述模型和/或校正参数。

2.根据权利要求1所述的调整传感系统的模型和/或校正参数的方法，其特征在于在步骤e)中，该快速卡尔曼滤波器公式是通过应用佛罗本尼变换公式(26)： ${\left[\begin{array}{cc}A& B\\ C& D\end{array}\right]}^{-1}=\left[\begin{array}{cc}{A}^{-1}+A^{-1}{\mathrm{BH}}^{-1}{\mathrm{CA}}^{-1}& -A^{-1}{\mathrm{BH}}^{-1}\\ -H^{-1}{\mathrm{CA}}^{-1}& H^{-1}\end{array}\right]$ 其中H＝D-CA^-1B获得的，其中A、B、C、D表示当对时间t计算增广模型(8)： $\left[\begin{array}{c}{y}_t\\ A_t{\hat{s}}_{t-1}+B_t{u}_{t-1}\end{array}\right]=\left[\begin{array}{c}{H}_t\\ I\end{array}\right]s_t\left[\begin{array}{c}{e}_t\\ A_t({\hat{s}}_{t-l}-s_{t-l})-a_t\end{array}\right]$ 的参数s_t的估值 时，标准方程的不同的分块系数矩阵。

Images