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

专利名称：一种动态系统中自适应卡尔曼滤波方法
技术领域：
本发明主要涉及所有的卡尔曼(KALMAN)滤波的实际应用，尤其是对需要快速可靠地适应环境的动态系统的控制。
现有技术在解释此发明前，有必要先了解一下传统的卡尔曼(KALMAN)递推公式的原有技术及用于校正传感器系统PCF/FI90/00122(WO90/13794)及控制一个大的动态系统PCT/FI93/00192(WO 93/22625)的快速卡尔曼滤波方法(FKF)下面的马尔可夫(MARKOV)(有限记忆)过程是通过方程式(1)至(3)描述的。第一个方程是关于量测矢量Yt与处于时间点t1(t＝0，1，2...)上的状态矢量St的关系，这是一个线性化的量测(或观察)方程yt＝Htst+et(1)矩阵Ht是设计(雅可比)矩阵，它是从实际物理关系的偏导中产生出来的。第二个方程描述了整个系统时间进展，它作为一个线性化的系统(或状态)方程st＝Atst-1+Btut-1+at(2)矩阵At是状态转移(雅可比)矩阵，Bt是控制增益(雅可比)矩阵。方程式(2)描述了整个系统的当前状态St是如何从前一状态St-1，控制/外力ut-1及随机误差at推导出来的。当测量误差et和系统误差既不是自动的(即白噪声)也不是相互关联的而是由下列协方差矩阵给出的时Ret=Cov(et)=E(etet')]]>and (3)Rat=Cov(at)=E(atat')]]>这样著名的卡尔曼正向递推公式(4)至(6)给出我们当前状态St是最佳线性无偏估计(BLUE)t如下t＝Att-1+Btut-1+Kt{yt-Ht(Att-1+Btut-1)}(4)及其估计误差协方差矩阵如下Pt＝Cov(t)＝E{(t-st)(t-st)’}＝Pt-1-KtH’tPt-1(5)式中卡尔曼增益矩阵Kt由下列公式确定Kt=(AtPt-1At'+Rat)Ht'{Ht(AtPt-1At'+Rat)Ht'+Re}-1--(6)]]>此递推线性结果是局部最佳的，卡尔曼滤波器(KF)的稳定性要求其可观测性及可控制性条件必须得到满足(卡尔曼，1960)然而在方程(6)中经常要求对一个较大的矩阵求逆，矩阵行(及列)的数目n与量测矢量Yt中的元素一样多，这需要一个大的n值以满足可观测性及可控制性条件，这就是本发明及PCF/FI90/00122PCT/FI93/00192所要解决的问题。
下列状态方程修正式导出Att-1+Btut-1＝Ist+At(t-1-st-1)-at(7)并与量测方程(1)结合以得到所谓的增广模型ytAts^t-1+Btut-1=HtIst+etAt(s^t-1-st-1)-at--(8)]]>pi.e.zt＝Ztst+ ηt即状态参数可通过使用人们所熟知的如下回归分析法计算出。s^t=(Zt'Vt-1Zt)-1Zt'Vt-1zt----(9)]]>该结果在代数上等同于使用卡尔曼递推式，但在数值上不相等(参见Harvey，1981“Time Series Models”，Philip allanPublishers Ltd，Oxford，UK，(“时代”系列模型飞利蒲Allan出版社，牛津，)PP 101-119)。此时方程(9)中可逆变矩阵的维数等于状态矢量St中元素的数目(m)。Harvey的计算方法是快速卡尔曼滤波(FKF)的各种不同的变形的基础。
为了满足可观察性条件的需要，任何大的卡尔曼滤波器(KF)的初始化或临时排序可通过Lange(兰格)的高通滤波器而完成(Lange 1988)。它利用一个解析稀疏矩阵逆变公式来解决带有下列所谓标准的分块对角矩阵结构的回归模型。

这就是例如完整的风口测量相互比较实验的量测方程的矩阵表示，矢量b1，b2…bk为连续的位置座标(如气象气球的座标)，但它也可包括那些带有一个显著的时间或空间偏差的校正参数，矢量e为其他的在抽样时间内为常数的校正参数。
对于所有大的多重传感器系统，它们的设计矩阵Ht是稀疏的，这样人们可以一种或其他同类的矩阵分块方法来完成


式中c1表示t时刻校正参数bt，k表示在时间或空间中所有其他的状态参数A表示t时刻状态变化矩阵(分块对角矩阵)B表示t时刻状态独立作用Vt矩阵(分块对角矩阵)若矩阵分块是不明显的，人们可以通过使用一个特定的算法自动解决，这个算法把每一个稀疏的线性系统逆变为标准分块对角形式(参见用于分析及其它算法的矩阵块角形式和再安排，管理科学，18卷1号，1971年9月98-107页(Weil and kettler，1971“Rearranging Matrices to Block-angular Form forDecomposition(and other)Algorithms”，Management Science，Vo1.18，No.1，Semptember 1971，Pages 98-107.))然而，随机误差et的协方差矩阵可解会使原来简单的对角性变得不严格。
结果我们将面临如下这个难解回归分析问题一个用于空间量的增广模型，例如对于一组包含K个连续气球位置的气球跟踪实验数据

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

请注意后一个矩阵方程具有一个“嵌套”分块对角结构，有两种形式“校正”参数11Ct及Ct，这些参数的第一组Ct可以随时间的变化而变化，而这些参数的第二组Ct在长度为L的移动时间窗口范围内为常量(至少近似为常量)，后一种参数Ct使卡尔曼滤波过程能够自适应。在从(4)至(6)的传统卡尔曼递推式中求出后一种参数时产生导致一个可观测性问题，就计算原因而言，长度L必须是短的，但是对于PCF/FI90/00122的FKCF公式，采样尺度可能很大以至于根本不需要初始化(或排序)。
在解释PCT/FI93/00192方法之前，有必要先了解一些应用于实验性数据天气预测(NWP)系统中的卡尔曼滤波器理论的现有技术。因为以前他们也利用方程(1)量测方程yt＝Htst+et(线性化回归)式中状态矢量St描述的是t时刻的大气状态，现在St通常代表所有大气变化的格点值，例如不用压力水平上的位势高度(实际上它们与实际数值的偏差量可以通过某些方法进行估计)。
空气动力学是由一个人们熟知的偏微分方程(原始方程)决定的，通过利用例如NWP模型的切线逼近得出，下列方程(2)的表达式是用于计算在某一时间段上状态参数St随时间推移的变化(实际上，他们相对于参数空间轨迹的偏移量产生于非线性NWP模型)。
状态方程st＝Ast-1+But-1+at(离散的动态随机模型)四维的数据同化结果(

>)及NWP预测(

>)分别从卡尔曼滤波器系统得到s^t=st~+Kt(yt-Hts~t)]]>s~t=As^t-1+But-1--(12)]]>式中

>Qt＝Cov(at)＝Eata’t…(系统噪声)Rt＝Cov(et)＝Eete’t…(量测噪声)关键性的校正计算是用下列卡尔曼递推式来完成的Kt＝PtH’t(HtPtH’t+Rt)-1…(卡尔曼增益矩阵)Cov(t)＝Pt-KtHtPt…(计算精度)这里所需的用于卡尔曼增益矩阵计算的矩阵逆变极难用于实际NWP系统计算，因为数据固化系统必须在同一时刻处理几百万个数据元素，T.GaL-Chen博士在此问题上于1988年报告如下“希望大型并行超级计算机的开发(例如1000台台式计算机CRAYS协力工作)能够使计算更接近于最佳....)”见“重要回顾小组的报告--低对流层分布论文需要及技术”美国气象科学文摘71卷5号，1990年5月684页。
用PCT/FI93/00192中的方法从方程(8)推出增广模型算法ytAs^t-1+But-1=HtIst+etA(s^t-1st-1)-at]]>zt＝Ztst+ ηt即下列两组方程用于校正的目的s^t=(Zt'Vt-1Zt)-1Zt'Vt-1zt]]>...(最佳估计，由高斯-马尔可夫得出)={Ht'Rt-1Ht+Pt-1}-1(Ht'Rt-1yt+Pt-1s~t)--(13)]]>或

及Cov(t)＝E(t-st)(t-st)’＝(Z’tV-1Zt)-1(14)

其中

Pt=Cov(s~t)=ACov(s^t-1)A'+Qt--(15)]]>而代替Kt＝PtH’t(HtPtH’t+Rt)-1…(卡尔曼增益矩阵)PCT/FI93/00192中的FKF方法用Kt=Cov(s^t)Ht'Rt-1--(16)]]>
对于一个大的输入数据Yt矢量，增广模型算法优于传统卡尔曼递推式，因为当

是未知的时，卡尔曼增益矩阵Kt的计算要示对大的矩阵求逆，两种方法在代数上或统计学上是等价的，但在数值不等。
然而，增广的模型公式仍难于在数字上解决，首先，这是由于状态矢量St中包括用于对大气进行实际描述的矢量(＝m)格点数据，其次，对于一个实际的NWP系统，许多其它的状态参数必须被包含在状态矢量中。首先是有关观察系统的系统(校正)误差及所谓的小规模大气过程的物理参数表。
在PCT/FI93/00192中校正问题是通过使用去耦合状态的方法来解决的。通过进行下列分块而完成

and (17)At=At,1··At,kAt,cand,Bt=Bt,1··Bt,kBt,c]]>式中Ct表示t时刻的校正参数Ct表示格点k(k1，....k)处大气参数值A表示t时刻状态转移矩阵(子矩阵A1....Ak，Ac)B表示控制增益矩阵(子矩阵B1，...，Bk，Bc)结果 面临下列艰巨的回归分析问题

在任何时间点t上的递推快速卡尔曼滤波器(FKF)公式如下b^t,k={Xt,k'Vt,k-1Xt,k}-1Xt.k'Vt,k-1(yt,k-Gt,kc^t)fork=1,2,…k]]>c^t={Σk=0kGt,k'Rt,kGt,k}-1Σk=0kGt,k'Rt,kyt,k---(19)]]>式中设k＝1，2，....，k，Rt,k=Vt,k-1{I-Xt,k{Xt,k'Vt,k-1Xt,k}-1Xt,k'Vt,k-1}]]>Vt,k=Cov(et,k)Cov{Ak(s^t-1-st-1)-abt,k}]]>yt,k=yt,kAks^t-1+Bkut-1]]>Xt,k=[Xt,kI]]]>

且，即设k＝0Rt,0=Vt,0-1]]>Vt,0=Cov{Ac(s^t-1-st-1)-act}]]>yt，0＝Act-1+Bcut-1Gt，0＝I.从方程(20)得到下列数据同化精度Cov(s^t)=Cov(b^t,1,···,b^t,K,c^t)---(20)]]>

式中Ck={Xt,k'Vt,k-1Xt,k}-1]]>设k＝1，2，....，k，Dk={Xt,k'Vt,k-1Xt,k}-1Xt,k'Vt,k-1Gt,k]]>设k＝1，2，....，k，s={Σk=0kGt,k'Rt,kGt,k}-1]]>卡尔曼滤波器(KF)研究也曾被报道过，例如StephenE.Cohn和David F.Parrish(斯蒂芳E考恩及大卫·F帕力士)(1991)“二维卡尔曼滤波器预测误差协方差特点”美国气象科学每月天气回顾，119卷，1757至1785页。然而，对于四维(即空间及时间)理想的卡尔曼滤波器系统，在所有这些报告中都没有涉及到。这样就需要一个可靠的对状态参数误差协方差距阵进行估计和逆变的方法，如中部区域天气预报欧洲中心(ECMWF)的He；kklJarvinen博士所述“在气象学中，状态参数矢量St的维数(＝m)可能是100,000-10,000,000。这使得在实际中不可能准确处理，误差协方差矩阵，”见“作为一个派生问题的气候学数据同化”报告第43号(1995)气候系Helsink大学第10页。中部区域天气预报欧洲中心的Adrian Simmons博士确认说“在理论上卡尔曼滤波基本算法是很完善的，但是计算上的要求使得整个执行过程难以处理，见中部区域天气预报欧洲中心69号简报(1995年春)第12页。
从PCF/FI90/00122及PCT/FI93/00192中得知的快速卡尔曼滤波器(FKF)公式利用这样的假设方程(9)及(13)中的误差矩阵分别是分块对角的，参见FKF公式(19)式中这些对角块被表述如下Vt,k=Cov(et,k)Cov{Ak(s^t-1-st-1)-abt,k}]]>尤其是对于自适应卡尔曼滤波(及4维数据比较)情况，很显然连续状态参数矢量St-1，St-2，St-2的估计值是相互且自动相关的。
有必要将法斯特·卡尔曼滤波方法的原理应用于自适应卡尔曼滤波(AKF)，比起其它卡尔曼滤波方法，它具有同样或更快的运算速度、可靠性、准确度和低成本。在此将公开本发明的一个处理误差协方差的具体方法。
发明概要通过提供一种用于校正/调整实时或接近实时的动态系统的各种参数的自适应快速卡尔曼滤波法可以基本满足上述需要，测量误差和系统误差都是白化和部分正交化的。如详细说明中所描述这种FKF滤波运算与在可观察和可控制条件下的最佳卡尔曼滤波相接近，预测误差的方差和协方差提供了一种监视滤波稳定性的办法。
本发明最佳实施例我们重述一下线性量测(或观察)方程yt＝Htst+FtyCt+et(21)式中et代表白噪声，它既不与et-1，et-2.....相关，也不与St-1，St-2相关，且不与at，at-1，at-2相关，矩阵Ht是与前述从测量值Yt和状态参数St间物理关系的偏导中出的设计矩阵相同，见前面4页的矩阵分块(11)(以往对矩阵A和B的分块对角性假设不再有效)矩阵Fty描述了测量的系统误差与校正或“校正型”参数，状态矢量Ct之间的关系，该矢量lt在时间上是常数或变化极慢。矩阵Fty，Ft-1y，Ft-2y的列代表偏导数，类似于正弦波、矩形波的波形、“随机编码”函数等，以及根据已知物理关系、测量的系统误差的回归和自动回归(AR)确定的经验正交函数(EOF)。预测矢量Ct的元素将决定着红噪声系数的波幅，让我们看一下第5页上的一个类似的用于“自化”所观察到的“新的”测量序列的动态时间量增广模型。
同样我们重述一下线性系统(或状态)方程st＝(At+dAt)st-1+Btut-1+FtsCt+at(22)式中at表示白噪声，既与et，et-1，et-2不相关，也与t-1，t-2不相关，且与at-1，at-2不相关。矩阵At同样是状态转移矩阵，它来自于状态St与前一状态St-1之间的物理关系的偏导。矩阵Fts描述动态模型(如NWP)的系统误差与校正或“校正型”参数及矢量Ct的关系。其中G不随时间变化或随时间缓慢地变化。矩阵的列Fts，Ft-1s，Ft- 2s表示偏导数、如正弦波、矩形波这样的波型、随机编码函数等，及根据已知物理关系及模型的系统误差的回归和自动回归所确定的经验正交函数。估计矢量

的元素，将决定着红噪声系数的波幅。
矩阵dAt描述了动态(NWP)模型的系统状态转移误差是如何与当前(气候)条件成比例的。如果他们是未知的，但变化极慢，可以结合FKF法通过动态求均值来调整。下面就介绍这种调整方法。从系统方程(22)中得到结果并重述如下dAtst-1＝＝[s1I(mxm)，s2I(mxm)，…，smI(mxm)][da11，da21，…，dam1，da12，…，damm]’＝Mt-1rt(23)式中mt-1是一个由尺寸为m×m的m个对角线矩阵组成的矩阵，S1，S2....Sm是状态矢量St-1的m标量元素，rt是矩阵dAt的m×m个元素的列矢量。
请注意方程(23)颠倒了乘法顺序，这样可以象普通回归参数一样估计矩阵dAt的元素。
结果面临下面难解的回归分析问题一个动态时间量的增广模型(即用于对长度L动态抽样的偏差et和at的更新序列进行白化)

请注意上述矩阵方程具有一个嵌套分块对角结构，有三种不同形式的校正参数第一种Ct是嵌入每个时间段步骤t数据中。另两种形式由状态矢量Ct表示。第一组参数(即rt)用于校正状态转移矩阵的严重误差，第二组用于对测量误差和系统误差进行的白化和部分正交化(即用于对协方差矩阵分块对角化)最后两组参数在长的移动时间窗中或多或少保持为常量从而使卡尔曼滤波过程能够自适应。
同样请注意矩阵M不能象在方程(23)中那样取其最大尺寸(m×m2)。这是因为由于有太多的未知数，可观测性条件不能得到满足。因而矩阵M必须被有效地压缩以使其只表示与严重的转换误差相关的矩阵At中的那些元素。这种转化是通过利用所谓的最大相关法发现的。事实上偶然和缓慢位移的图形可以在状态参数矢量的空间中形成。这些一般是小范围的现象，而且它们不能被仅由模型方程导出的状态转移矩阵所充分描述，为了保证滤波稳定性，所有的矩阵dAt的估计元素必须在测量的动态平均中保持可观察，如通过方程(20)中监控它们的估计方差。
在时刻t，用于长度为L的时间窗口的法斯特卡尔曼滤波(FKF)公式如下s^t-l={Xt-l'Vt-l-1Xt-l}-1Xt-l'Vt-l-1(yt-l-Gt-lc^t)]]>设l＝0，1，2，....，L-1c^t={Σl=0LGt-l'Rt-lGt-l}-1Σl=0LGt-l'Rt-lyt-l--(25)]]>式中设l＝0，1，2，....，L-1，Rt-l=Vt-l-1{I-Xt-l{Xt-l'Vt-l-1Xt-l}-1Xt-l'Vt-l-1}]]>Vt-l=Cov(et-l)Cov{At-l(s^t-l-1-st-l-1)-at-l}]]>yt-l=yt-lAt-ls^t-l-1+Bt-lut-l-1]]>Xt-l=[Ht-lI]]]>

且，即设l＝L，Rt-L=Vt-L-1]]>Vt-L=Cov{Ac(C^t-1-Ct-1)-act]]>yt-L=AcC^t-1+Bcuct-1]]>Gt-L＝I.
为了最优化，有时有必要对滤波器的一些误差项进行具体化。如果这样做了恒方程(I)矩阵将从FKF公式中消失并必须被相应地代换。
此处的和在PCF/FI90/00122和PCT/FI93/00192中的FKF公式是基于假设误差协方差矩阵是分块对角的。因为运算上的限制禁止采用足够长的窗口，由于严重的可观察性和可控制性难题，要用传统的卡尔曼递推式(4)-(6)求得所有参数Ct，注定要失败。幸运的是，通过运用FKF公式，时间窗口可采用足够的长度以致于滤波的初始化或暂时排序变得完全多余了。
各种快速自适应卡尔曼滤波可以通过使用不同的佛罗本尼公式递归从大型线性回归方程(24)的标准方程系统中推导出ABCD-1=A-1+A-1BH-1CA-1-A-1BH-1-H-1CA-1H-1---(26)]]>式中H＝D-(A-1Bo)公式(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.一种通过自适应卡尔曼滤波用于大型传感系统的模型调整和参数校正的方法，这种大型传感系统的传感器输出单元根据外部事件输出信号，而要同步处理的传感器输出信号的数目超过50个，这种方法包括的步骤有a提供一种用于储存信息的数据的数据库单元，这些信息是—为一些所述传感器提供多个测试点传感器的输出信号值和为对应于所述测试点传感器所述外部事件提供多个数值，或从相邻的传感器提供所述输出信号值的同一时间序列用于对照；—伴随着用于所述模型和校正参数的数值及用于对应某一情形的外部事件的数值产生的所述传感器输出信号数值；—所述传感器的控制及对应于新的情况的外部事件的改变；b提供一种逻辑单元，以存取带有所述模型及校正参数的所述传感器信号输出值，所说逻辑单元有双路通道与所述数据库单元相连，按需要用Lange(兰格)高通滤波器计算未知模型的初始值及具有估计精度的校正参数；c将来自传感器的传感器输出信号提供给所述逻辑单元；d将所述控制和变化信息提供给数据库单元；e访问所述模型及状态转换矩阵的校正参数和元素的当前值，用FKF滤波器公式进行计算，该公式从前述逻辑单元Lange(兰格)的稀疏矩阵逆变公式得出，在所述逻辑单元中模型与校正参数的更新、外部事件的数值及其精度对应于所述的新情况；f通过监控所述的逻辑单元的精度估计并指出何时需要新的测试点、比较或系统重构，来控制所说的卡尔曼滤波的稳定性；g调整所述模型校正参数值使得可以实现持续地更新数据。
全文摘要
本发明是根据将Lange的快速卡尔曼滤波(FKF
文档编号G01D18/00GK1202240SQ96198347
公开日1998年12月16日 申请日期1996年11月15日 优先权日1995年11月15日
发明者安蒂·阿米·依玛里·兰格 申请人:安蒂·阿米·依玛里·兰格