联邦滤波器在船舶导航中的应用.doc

联邦滤波器在船舶导航中的应用1 准备知识So federal filter with its design flexibility and small computation, fault tolerance is good wait for a characteristic to achieve high precision navigation provides an effective way. 联邦滤波器由多个子滤波器和一个主滤波器组成，是一种具有两级数据处理的分散化滤波方法。子滤波器的工作状态与其所受干扰情况密切相关。当干扰大到足以影响滤波器正常工作时，一般就认为子滤波器发生了故障。为了描述子滤波器的工作状态，定义如下新息序列：Ui,k=Zi,k-Hi,kXi,k/(k-1) 当子滤波器工作存在干扰时，观测方程为：Zi,k=Hi,kXi,k+Vi,k 新息序列也随之变为：ui,k=Hi,kXi,k-Xi,kk-1 +Vi,k 2 工作原理2.1 集中式卡尔曼滤波概述卡尔曼滤波是解决状态最优估计问题的一种常用方法。设随机动态系统的数学模型和有关随机向量的统计性质如下： 式中。是系统状态向量；是系统观测向量；是系统噪声向量；是观测噪声向量；是系统状态转移矩阵；是系统噪声矩阵；是系统观测矩阵。假定系统噪声和观测噪声是不相关的零均值高斯白噪声，初始状态向量是m维高斯随机向量，且： 随机系统的状态估计问题，就是根据选定的估计准则和获得的量测信息，对系统状态进行估计。其中状态估计确定了被估计的随机状态的转移过程，估计准则确定了状态估计最优性的含义，卡尔曼滤波的估计准则是：估计量是的无偏估计，即；同时也是的最小方差估计，即：。根据这两个准则推导出系统的完整的滤波算法。下面给出最终的结果：滤波方程： 预测方程： 增益矩阵： 预测误差协方差阵： 波误差协方差阵： 初值为 ,。卡尔曼滤波算法是一组递推算法。计算最优滤波值，只需要即时的观测值，无需存储以前的观测数据。利用上一时刻的最优估计值，由预测方程求得下一时刻的一步预测值，；当前测量值和之差，即新息，与相乘对进行修正，得出k+1时刻的最优估计值。如此反复递推运算，但这个计算顺序不是唯一的。在该计算中，我们把预测下一时刻的和称为滤波时间的更新，把得到和的过程称为滤波的测量更新。2.2 联邦滤波基本理论卡尔曼滤波及联邦滤波一个重要的应用是高精度船舶导航系统。现在可供运载体使用的导航装备很多，选择余地很大。由于采用各种原理实现导航设备的增加使量测信息增多，这对提高船舶导航系统精度提出了很高的要求2.2.1 联邦滤波器结构联邦滤波器是一种信息融合技术，其结构可通过图1来说明，它由若干个子滤波器和一个主滤波器组成，是一种具有两级数据处理的特殊的分散化滤波方法，其特殊性在于联邦滤波器采用的是信息分配原理。图中公共参考系统的输出一方面直接给主滤波器，另一方面给各子滤波器作为公共状态变量。各子滤波器的局部估计值（公共状态）及其协方差阵送入主滤波器和主滤波器的估计值一起进行融合以得到全局最优估计和相应的协方差阵。参考系统子系统1子系统2子系统N局部滤波器1局部滤波器2局部滤波器N主滤波器时间更新最优组合图1 联邦滤波器的一般结2.2.2 联邦滤波器算法设、表示联邦滤波器的最优估计值和方差，Xi、Pi表示第i个子滤波器的估计值和方差（i=1，2.n），、表示主滤波器的估计值和方差。全局状态方程为： 其中：是系统的状态转移矩阵，是系统的噪声矩阵。设是零均值的白噪声序列，其协方差阵为。第i个子系统的量测方程为： 其中：Zi(k)是第i个量测传感器的量测值，Hi（k）是第i个量测传感器的量测阵，Vi（k）是传感器的观测噪声阵，是独立于W(k)的零均值白噪声序列，其协方差阵为R（k）。信息分配过程：信息分配就是在各子滤波器和主滤波器之间分配系统的信息。系统的过程信息和按如下的信息分配原则在各子滤波器和主滤波器之间进行分配：， 其中是信息分配因子，并满足信息分配原理： 在系统噪声分配过程中，采用了方程上界技术，得到的子滤波器的滤波结果是次优的。物理意义上是由于过程噪声方差被上界取代，则每个滤波器更少的依靠时间更新后的状态值，更多的依靠测量更新。信息的时间更新：时间更新过程在各子滤波器和主滤波器之间独立进行，各子滤波器和主滤波器的滤波算法为： （13）信息的测量，由于主滤波器没有测量值，所以主滤波器没有测量更新。测量更新只在各个局部子滤波器中进行，测量更新通过下式起作用： 信息融合：联邦滤波器的核心算法是将各个局部滤波器的局部估计信息按下式进行融合，以得到全局的最优估计。 通过以上的信息分配、添加（时间更新和测量更新）和信息融合过程，在局部滤波器中由于方差上界技术引起的信息丢失，在融合过程中这种非最优性被重新合成，得到全局最优解。3 船舶航行姿态系统中联邦滤波器设计3.1 船舶航行系统中联邦滤波器结构在多传感器组合导航系统中，联邦滤波器能够利用信息分配原理实现多传感器信息的最优估计，同时使整个系统具有一定的容错能力。考虑到各个子系统的导航精度都较高，在船舶航行系统中，各子系统容易出现故障，为确保组合系统具有高容错性能，本文采用的是如图3-1所示的联邦滤波器无复位结构。子滤波器选择GPS（全球定位系统）、SINS（惯性导航系统）、TAN（地形辅助性导航系统）和北斗双星子系统。GPS子系统量测信息与SINS系统组成子滤波器1；北斗双星子系统量测信息与SINS 系统组成子滤波器2 ，TAN子系统量测信息与SINS 系统组成子滤波器3 ；主滤波器的输出是所有子系统与SINS 系统的量测所确定的全局最优估计。GPSBDTAN子滤波器1子滤波器2子滤波器3主滤波器时间更新信息融合日图2 联邦滤波器结构图中,分别为三个子滤波器的状态变量和状态协方差阵；为主滤波器状态变量和主滤波器状态协方差阵。系统选取的全局状态变量为 式中：为纬度弧长；为经度弧长；为海流东向速度分量；为海流北向速度分量；为船舶相对水的速度；为船舶航向。动态方程为：其中：为独立的高斯白噪声；海流用一阶马尔可夫过程来表示；为相关时间常数。全局状态方程状态转移矩阵 图3 联邦滤波器流程图设动态系统为: 观测方程为：3.2无重置估计仿真结果： 图4 无重置仿真状态X1仿真曲线 图5 无重置仿真状态X2仿真曲线图6 无重置仿真状态X3仿真曲线有重置估计仿真上图为400次滤波的状态误差。蓝色为滤波前X曲线，红色为滤波后X曲线，绿色为主滤波器估计值曲线。无重置结构的滤波器受噪声的影响大，系统容错性能差。3.3有重置估计仿真仿真的信息分配系数为 图7 有重置仿真状态X1仿真曲线 图8 有重置仿真状态X1仿真曲线图9 有重置仿真状态X3仿真曲线上图为400次滤波的状态误差。蓝色为滤波前X曲线，红色为滤波后X曲线，绿色为主滤波器估计值曲线。 很明显看出有重置的联邦滤波器能够更好地滤除噪声干扰，子滤波器的得到的估计值相似性非常高，可以更好地降低某一个子滤波器因为故障或者抗干扰能力差对系统造成的影响，极大提高了系统的容错性，极大地提高了系统的稳定性。附件无重置滤波仿真程序：clearN=400;w(1)=0;w=randn(1,N)x1(1)=0;x2(1)=0;x3(1)=0;a=1;for k=2:N;x1(k)=x1(k-1)+w(k-1);x2(k)=0.4*x2(k-1)+w(k-1);x3(k)=0.9*x3(k-1)+w(k-1);end V=randn(1,N); q1=std(V); Rvv=q1.2; q3=std(w); Rww=q3.2; a=1;b=0.4;c=0.9;d=0.2; Y1=d*x1+V; Y2=d*x2+V;Y3=d*x3+V;a1(1)=0; a2(1)=0;a3(1)=0;s1(1)= 0; s2(1)= 0; s3(1)= 0; e(1)=1f(1)=1 for t=2:N; p1(t)=a.2*a1(t-1)+Rww; b1(t)=d*p1(t)/(d.2*p1(t)+Rvv); s1(t)=a*s1(t-1)+b1 (t)*(Y1(t)-a*d*s1(t-1); a1(t)=p1(t)-d*b1(t)*p1(t); p2(t)=b.2*a2(t-1)+Rww; b2(t)=d*p1(t)/(d.2*p2(t)+Rvv); s2(t)=b*s2(t-1)+b2(t)*(Y2(t)-b*d*s2(t-1); a2(t)=p2(t)-d*b2(t)*p1(t); p3(t)=c.2*a3(t-1)+Rww; b3(t)=d*p3(t)/(d.2*p1(t)+Rvv); s3(t)=c*s3(t-1)+b3(t)*(Y3(t)-c*d*s3(t-1); a3(t)=p3(t)-d*b3(t)*p1(t); e(t)=1/(1/a1(t)+1/a2(t)+1/a3(t);f(t)=e(t)*(s1(t)/a1(t)+s2(t)/a2(t)+s3(t)/a3(t)end t=1:Nfigure(1);plot(t,s1,r,t,x1,b,t,f,g),xlabel(时间),ylabel(状态x1)figure(2);plot(t,s2,r,t,x2,b,t,f,g),xlabel(时间),ylabel(状态x2)figure(3);plot(t,s3,r,t,x3,b,t,f,g),xlabel(时间),ylabel(状态x3)无重置滤波仿真程序：clearN=400;w(1)=0;w=randn(1,N)x1(1)=0;x2(1)=0;x3(1)=0;a=1;for k=2:N;x1(k)=x1(k-1)+w(k-1);x2(k)=0.4*x2(k-1)+w(k-1);x3(k)=0.9*x3(k-1)+w(k-1);end V=randn(1,N); q1=std(V); Rvv=q1.2; q3=std(w); Rww=q3.2; a=1;b=0.4;c=0.9;d=0.2; Y1=d*x1+V; Y2=d*x2+V;Y3=d*x3+V;a1(1)=0; a2(1)=0;a3(1)=0;s1(1)=0; s2(1)=0; s3(1)=0; e(1)=1f(1)=1 for t=2:N; p1(t)=a.2*3*e(t-1)+3*Rww;a1(t)=1/(1/p1(t)+ d.2/Rvv);s1(t)=a1(t)*(a*f(t-1)/p1(t)+d*Y1(t)/Rvv); p2(t)=b.2*3*e(t-1)+3*Rww;a2(t)=1/(1/p2(t)+ d.2/Rvv);s2(t)=a2(t)*(b*f(t-1)/p2(t)+d*Y2(t)/Rvv); p3(t)=c.2*3*e(t-1)+3*Rww;a3(t)=1/(1/p3(t)+ d.2/Rvv);s3(t)=a3(t)*(c*f(t-1)/p3(t)+d*Y3(t)/Rvv); e(t)=1/(1/a1(t)+1/a2(t)+