基于矩阵分解的卡尔曼滤波技术分析及应用

1、-word资料-基于矩阵分解的卡尔曼滤波技术分析及应用【摘要】本文简要介绍了卡尔曼滤波研究的发展历程,重点对卡尔曼滤波及其在改善数值稳定性，提高计算效率等数值方面的研究与发展进行了综述，对Q-R分解，U-D分解，奇异值分解(SVD)等在卡尔曼滤波的应用进行了介绍。最后给出了一种基于Q-R矩阵分解的自适应滤波方法，仿真验证了其有效性。引言I960年美籍科学家卡尔曼(R.E.Kalman)在系统状态空间模型的基础上提出了著名的线性卡尔曼滤波器，它在线性的前提假设下是一个线性无偏、最小方差估计器，从而可以为线性滤波问题提供精确解析解。自该技术被提出以来，它已成为控制、信号处理与通信等领域最基本最重要

2、的计算方法和工具之一,并已成功地应用到航空、航天、电力系统及社会经济等不同领域。随着微型计算机的普及应用,对卡尔曼滤波的数值稳定性、计算效率、实用性和有效性的要求越来越高.为此,人们在如何改善卡尔曼滤波的计算复杂性和数值稳定性方面作了大量的探索工作,各种基于平方根滤波与平滑,U-D分解滤波与平滑,奇异值分解滤波与平滑,状态与偏差分离滤波以及并行与分散滤波等方法得到不断发展.本文给出了矩阵分解的一些基础知识，并着重从卡尔曼滤波数值计算方法入手,对现有的常规卡尔曼滤波、基于矩阵的因式分解滤波的数值计算方法进行了较系统的介绍和分析,并在第四章给出了一种基于Q-R矩阵分解的自适应滤波算法。常规卡尔曼滤

3、波协方差卡尔曼滤波考虑如下线性离散系统x=Ax+rw(2.1.1)k+1kkkkz=Cx+vkkkk(2.1.2)式中xeRn是状态向量，zeRm是量测向量，weRp是系统噪声向量，veRm是量kkkk测噪声向量假设系统噪声和量测噪声是互不相关的零均值高斯白噪声,方差阵分别为Qk，kR，则协方差卡尔曼滤波方程为：kx=Ax(2.1.3)kik-ik-ik-ip=apAt+rqrtkIk-1k-1k-1k-1k-1k-1k-1(2.1.4)ZYZYZYx=x+Kz一Cxkk|k-1kkkk|k-1(2.1.5)p=I-KCpkkkk|k-1(2.1.6)K=PCtCPCt+R-1k|k|k-1k

4、kk|k-1kk理论分析和实际应用均证明上述滤波公式是数值不稳定的,其原因是由于计算机有限字长的限制,计算中舍入误差和截断误差的累积、传递会使协方差阵P失去对称正定性,因此,Joseph提k出一种所谓“稳定化”卡尔曼滤波，其目的是减小滤波算法对计算舍入误差的灵敏性，保证P的对k称正定性,以提高滤波的数值稳定性,防止发散.其滤波阵公式，只是将(2.1.6)式改写为如下形式即可:P=I-KCPI-KCt+KRKtkkkk|k-1kkkkk但该算法由于所需计算量和存储量较大,而且并不一定很奏效,因而应用并不广泛.信息滤波为了解决在某些没有有关初始状态信息和先验知识可供采用情况下的滤波,Fraser提

5、出了信息滤波，即用协方差阵P的逆P-1来代替P的递推计算，这种算法对测量更新比kkk较有效,但时间更新所需计算量较大.推广卡尔曼滤波器推广卡尔曼滤波(EKF)是一种应用最广泛的非线性系统滤波方法。EKF与线性卡尔曼滤波公式完全类似,只是上述滤波公式中Ak，匚，和Ck要在由非线性函数的偏导计算得到，不能象线性滤波那样可事先离线计算增益和协方差阵，但EKF与常规卡尔曼滤波一样,数值稳定性差,初值不易确定.为了改善上述常规滤波算法的数值稳定性,并提高计算效率，自七十年代以来，人们提出了平方根滤波、UD分解滤波、奇异值分解滤波等一系列数值鲁棒的滤波算法.基于矩阵因式分解的滤波方法3.1预备知识定理31

6、1设A是实正定对称矩阵，则存在唯一正线下三角矩阵S,使得A=SSt(3.1.1)定理3.1.2Householder变换设ueCn，且uHu=1，则|H(u)=E(u,u;2)=E一2uuhn(3.1.2)称为初等酉阵，或Householder变换。定理313Cholesky分解设AeC阿是正定Hermite矩阵，用L表示单位下三角n矩阵，D是对角矩阵，则有A=(LD1/2)(LD1/2)H(3.1.3)定理314QR分解设AeR，则a可唯地分解为nA=QR(3.1.4)定理315奇异值分解设AeCmxn，cQ,是A的r个奇异值，J则存在m阶r12r酉矩阵和n阶酉矩阵V，使得DOA=U()VO

7、O(3.1.5)其中D=diag(8,8,)，且15l=c(i=1,2,r)。12rii3.2平方根协方差滤波(SRCF)首先提出平方根滤波思想的是Potter,他把戶上按Cholesky方法分解为下三角阵sk，即令P=SST，在滤波递推计算中用Sk的传递计算代替P的计算，由公式（3.1.1）可知，从而kkkkk保证了P的对称正定性.Potter的算法经美国阿波罗登月舱的实际应用，证明是很成功的.随k后,Potter的算法被推广来解决存在着系统噪声和量测量为向量的情形。Schmidt给出了向量量测既可以同时处理，也可以序列处理的一种处理过程噪声的方法.为了提高平方根滤波的计算效率,Carlso

8、n注意到传递阵通常是块上三角阵的特点，给出了一种量测更新和时间更新均为上三角阵形式的快速平方根滤波，减少了计算量.上述平方根滤波均把时间更新和量测更新按常规分成两个分离的过程，其算法的关键是通过利用正交变换获得上三角阵的平方根矩阵.为了减小计算量，人们对如何构造正交变换的问题给予了很大的注意，常用的正交变换方法是Householder变换，即公式（3.1.2）、修正的Gram-Sehmidt正交化法及Givens变换等。1975年,MorfKailath在总结上述平方根滤波基础上，把时间更新和量测更新两个过程结合起来，给出了一种量测和时间更新的联合更新方程，从而仅需一个正交变换，即完成滤波计算

9、，且无需计算滤波增益阵.3.2平方根信息滤波与平方根协方差滤波相对应，信息滤波的平方根滤波方法也得到人们的极大重视和研究。DyerMeReynolds基于Householder变换利用动态规划理论研究出一种平方根信息滤波（SRIF）,与SRCF类似,SRIF把信息矩阵P-1定义为平方根阵形式，即定义P-i=S-1S（-i）t由S-1kkkkk的递推计算来代替P-1的计算.Bierman利用“数据方程”法给出一种结构较简单的SRIF,并给k出有色噪声情况的滤波公式,该算法需要计算状态转移矩阵的逆,即要求状态转移矩阵是非奇异的，针对这一问题，给出一种对状态转移阵奇异仍适用的SRIF,Bierma在

10、此基础上，把SRIF应用于具有时间延迟系统的滤波，并把SRIF推广到大规模互联系统的情形，大大减少了计算量和存储量正如SRCF那样,将量测更新和时间更新结合起来,可以容易的求得联合SRIF更新方程,Paige-Saunders基于把卡尔曼滤波转换为最小二乘估计的思想，提出一种联合量测更新和时间更新的SRIF方案。由于SRIF在某些情况下，如对于多量测量系统,比常规卡尔曼滤波有更高的计算效率、更好的数值稳定性和精度，因此，在轨道确定、飞行状态估计和多传感器跟踪与辨识等方面得到应用。U-D分解滤波上述SRCF和SRIF，般来讲由于存在矩阵的求逆运算和平方根计算，所需计算量较常规卡尔曼滤波要大，因而

11、限制了在工程中的应用。Bierman在研究和应用SRIF及Carlson序列滤波的基础上，于19751977年间，提出了一套计算效率高、数值稳定的称之为“U-D分解”滤波的算法。该算法把协方差阵分解为单位上三角阵U和对角阵D,即有Pk=UkDkUT,UD1/2kkkk相当于协方差平方根阵S,即公式（3.1.3）。U-D分解滤波既具有平方根滤波的优点，即始终能保证协方差阵的正定性，同时避免了Carlson等平方根滤波算法中平方根的计算，因而具有与常规卡尔曼滤波相当的计算量，是上述滤波算法中效率最高的一种算法，并且在实际应用问题中，结合实际问题的特点,U-D分解算法计算效率还更高，因而，近年来在轨

12、道确定、目标跟踪和飞行状态估计及神经网络学习算法等方面得到广泛应用和发展,但该算法由于量测更新采用序列处理,对于量测量较多的系统,计算效率受到一定影响,对有色噪声的处理不如SRIF,SRCF方便。3.4基于奇异值分解(SVD)的滤波方法奇异值分解由于具有很强的数值鲁棒性和可靠性,广泛应用于最小二乘问题、病态方程组求解及广义逆计算等场合,并在控制、通讯与信号处理等领域越来越受到人们的极大重视在滤波问题中,也已得到应用。Oshman基于协方差阵的频谱分解,以SVD为计算工具,提出称之为v-Lambda滤波的方法,他把协方差阵分解成P=VVT形式，其中V是矩阵P的特征向量矩阵,人为对角元是P阵奇异值

13、的对角矩阵。首先给出量测更新为信息滤波模式,时间更新为协方差滤波模式的vLambda滤波方法。随后,Oshman又给出量测更新和时间更新方程均为信息滤波模式的滤波公式，该算法不需计算滤波增益，由于利用SVD，使得状态估计算法鲁棒性较之平方根滤波、U-D分解滤波更好，但该算法由于进行一步滤波迭代计算，需一次正交变换,两次奇异值分解，所以，其缺点是计算量较大，但其优异的数值鲁棒性以及随着奇异值分解并行处理的实现而随之带来计算时间的减少，使得此算法将成为种极富吸引力的滤波方法.另外利用协方差阵的对称正定性,给出一种类似U-D分解形式、计算量较小的基于SVD的滤波算法,本文又给出种基于SVD的推广卡尔

14、曼滤波算法，并应用于飞行状态估计问题，随后又提出种基于SVD的递推最小二乘辨识新方法，与递推最小二乘、基于U-D分解的递推最小二乘法相比,不仅收敛速度快、数值稳定性和辨识精度高，而且能得到系统参数的无偏估计。一种基于Q-R矩阵分解的自适应滤波算法4.1机动目标加速度模型及自适应卡尔曼算法令状态向量X为：X二x1,x2,x3T,x1,x1,x2,x3分别表示目标的位置、速度和加速度。假设目标当前加速度服从非零均值矿的一阶马尔可夫过程，即：x3=u+a(t)a(t)=-aa(t)+w(t)(4.1.2)目标的时域状态方程为：X二AX+Ba+Lw(t)4.1.3)式中：0一0_o1、B=0、L=0a

15、a101A=0000、w(t)为白噪声序列根据上述时域方程,经过离散化处理后,得到系统的离散状态方程和观测方程为:X(k+1)=(k+1,k)X(k)+B(k)a+W(k)4.1.4)Y(k)=HX(k)+V(k)4.1.5)式中：1T(k+1,k)=01(1+aT+eaT)a2(1eaT)a00eaTB(k)=aT21eaT+)2a1eaTTa1eaT(k+1)TW(k)=J1+a(k+1)T5)+ea(k+1)Ts)/a21一ea(k+1)t5)/aw(5)d5kTea(k+1)T5)H=100W(k)、V(k)均为白噪声序列。基于状态方程(4.1.4)和观测方程(4.1.5)，依据标准的

16、卡尔曼滤波方程，得到当前机动加速度统计模型的自适应卡尔曼滤波算法。1.时间更新X=O去+B(k)a(4.1.6)(k/k1)(k,k1)(k1/k1)(4.1.7)b=x3(k/k-1)82=(4-兀)(B-X3(k/k)2/K(4.1.8)bmaxp(k/k-1)二(k,k-1)P(k-1/k-1)t(k,k-1)+Q(4.1.9)wK(k)二P(k/k-1)Ht(HP(k/k-1)Ht+R)-1(4.1.10)其中，amaX表示最大可能的加速度。max2.测量更新X=X+K(k)Y(k)-hX(k/k-1)(4.1.11)(k/k)(k/k-1)P(k/k)=P(k/k1)K(k)HP(k

17、/k1)=IK(k)HP(k/k1)(4.1.12)4.2引入Q-R矩阵分解的自适应平方根滤波算法在机动“当前”统计自适应滤波算法中,在计算协方差阵时,存在矩阵相减的运算,由于计算或其他参数不匹配的影响,有可能导致协方差阵出现不对称或负定的情况。协方差阵是一个对称的非负定的矩阵,利用矩阵分解的技术，在基于平方根矩阵分解的基础上，引入Q-R矩阵分解,构造出协方差平方根自适应滤波算法。令P(k/k)二S(k/k)ST(k/k)P(k/k-1)二S(k/k-1)ST(k/k-1)Q=UUtw其中S(k/k)、S(k/kD、u为下三角矩阵。又P(k/k-1)=O(k,k-1)P(k-1/k-1)t(k

18、,k-1)+Q=w(k,k-1)S(k-1/k-1)St(k-1/k-1)t(k,k-1)+UUt=(k,k-1)S(k-1/k-1)U(k,k-1)S(k-1/k-1)Ut=S(k-1/k-1)ST(k-1/k-1)对S(k-1/k-1)进行Q-R分解得到：(k/k-1)=S(k/k-1)Q(k/k-1)(4.2.1)因此有P(k/k-1)=(k-1/k-1)t(k-1/k-1)=S(k/k-l)Q(k/k-1)Qt(k/k-1)St(k/k-1)=S(k/k-1)ST(k/k-1)经过推导，得到下述的基于Q-R分解的自适应平方根卡尔曼滤波算法。1.时间更新(k/k-1)=S(k/k-1)U

19、(4.2.2)进行Q-R分解得到(k/k-1)=S(k/k-1)Q(k/k-1)(4.2.3)F=St(k/k1)QHtka=FFt+R-1kkk1+-aRr=k-k1+aRk(4.2.4)(4.2.5)(4.2.6)K(k)=aS(k/k-1)Fkk(4.2.7)文(k/k-1)=(k/k-1)J?(k-1/k-1)+Ba(4.2.8)a=X3(k/k-1)(4.2.9)52=(4兀)(a-X3(k/k)2/兀amax(4.2.10)P(k/k-1)=S(k/k-1)ST(k/k-1)(4.2.11)2.测量更新S(k/k)=S(k/k-1)I-arFFT(4.2.12)kkkkX?=X?+

20、K(k)Y(k)-HX?(k/k-1)(4.2.13)(k/k)(k/k-1)P(k/k)=S(k/k)ST(k/k)(4.2.14)4.3仿真与结论为了验证算法的有效性，进行了一系列的仿真，并在一维的情况下，研究了基于Q-R矩阵分解的当前统计自适应滤波算法对常加速度目标运动的跟踪特性。仿真中，假设观测噪声方差(距离观测误差)与目标距离的平方成正比，即观测噪声为:V(k)=(0-x(k)+c)w(k)式中B为相对误差系数,c为固定观测误差,w(k)为均值为零,方差为1的标准正态分布函数因此,观测噪声方差为:R(k)=(0 x(k)+c)2E(w(k)w(k)仿真中所选参数为:c=100m,B=0.01,目标加速度a=10m/s2