容积卡尔曼滤波

2024-07-06

容积卡尔曼滤波（精选8篇）

容积卡尔曼滤波第1篇

目前, 对卫星轨道进行实时预报的方法有利用初轨直接积分、扩展卡尔曼滤波 (extended Kalman filter, EKF) 方法[2]和无迹卡尔曼滤波 (unscented Kalman filter, UKF) 方法。利用初轨和动力学模型直接积分的预报精度受初始轨道的精度影响非常大。在实时轨道预报中, 初轨精度通常不高, 所以很少直接积分进行预报。EKF是一种较常用的非线性系统滤波方法[3], 通过截取系统函数的一阶乃至二阶项, 将非线性模型线性化, 然后利用卡尔曼滤波进行计算[4]。但是, 对于含有多个摄动力的卫星运动模型来说, EKF的应用存在两个较大的问题。一是卫星运动模型非线性化程度高, 对系统进行线性化时舍去的高阶项不能完全忽略, 线性化后模型精度会明显降低;二是由于动力学参数较多, 造成状态向量维数较大, EKF计算过程中必然出现的Jocobin矩阵很复杂, 对其进行处理会造成计算量的增加, 大大降低计算速度。

容积卡尔曼滤波器 (cubature Kalman filter, CKF) 是Arasaratnam等人于2009年提出的一种新的滤波方法, 其核心是利用带权重的采样点近似概率分布, 通过球面积分方法求随机变量与函数乘积的积分, 进而估计变量的均值和方差。这样做既不需要对函数进行线性化, 也不需要计算Jacobin矩阵。所以, 理论上来说, CKF具备更优的非线性逼近性能、数值精度以及滤波稳定性, 且实现简单[5]。

1 卫星轨道短周期预报的一般过程

轨道预报通常是利用观测到的卫星轨道数据, 首先通过参数估计进行初始定轨, 然后利用初轨和力学模型进行计算得到预报的轨道值。参数估计是其中的关键环节, 估计的精度直接关系到最终结果的误差大小。参数估计有两大类方法, 一类是以最小二乘法为代表的批处理, 另一类是以卡尔曼滤波为代表的序贯处理[6]。实时轨道预报中批处理方法因计算量较大而很少采用, 序贯处理则采用递推的方式, 依次处理数据, 效率更高。

序贯处理的一般过程如图1所示。首先, 由已知的观测数据[Z1, Z2, …, Zk]T进行滤波, 估计动力学参数Lk和卫星状态Xk, 同时预报下一时刻的卫星状态Xk+1和观测值Zk+1, 然后读取下一时刻的观测数据Zk+1, 通过估计观测值Zk+1和实际观测值Zk+1的偏差计算增益K和协方差阵P, 然后更新预测值Xk+1, 以此作为这一时刻的卫星状态重复前述步骤, 形成一个递推的轨道预报过程。

卫星动力学模型是轨道预报的重要组成部分, 摄动力和参数的估计精度对最终预报精度有较大影响[7]。在此选取其中最重要的几个摄动力, 地球非球型引力摄动J2项、日月引力、大气阻力和太阳光压辐射[8], 具体动力学模型参见文献[9]、[10]。

2 基于CKF的卫星轨道短周期预报

2.1 CKF原理简介

应用Bayes滤波来解决实际问题时, 经常会遇到这样的积分难题。

式 (1) 中, f (x) 是任意函数, Y是随机变量, Rn为积分域。通常来讲, 这个积分是不能够得到准确结果的。CKF算法的重要突破就在于, 利用带权重的点w (x) 模拟随机变量的分布, 然后基于求积原理的用球面积分方法进行积分, 从而得到近似的概率分布[11]。具体来说, 就是对于积分

式 (2) 中, w (x) 是权重函数。选择一组具有权重值的采样点集 (ξi, ωi) 来近似积分。

对于如下的离散时间非线性系统

式 (4) 中, xk为系统状态向量, zk观测向量, wk和vk为相互独立的白噪声, 其统计特性满足

利用CKF算法进行状态预测可以分为选取采样点、时间更新、量测更新三个步骤[9]。

2.2 基于CKF的卫星轨道短周期预报

以广播星历计算的卫星状态作为观测量, 故, 将其离散化以后就可以利用CKF方法进行预报, 具体来说分为以下步骤

1) 起步将摄动参数L设为0, 则k时刻的状态向量Xk=[ZkT, 0]T, 通过滤波更新状态向量得到估计值Xk^。

2) CKF根据cubature变换得到2n个具有相同权值的cubature点来近似积分:

式 (7) 中, m表示容积点总数, m=2n, n为系统的状态维数。使用[1]表示对单位向量e的元素进行全排列和改变元素符号产生的点集, 称为完整全对称点集。[1]j表示点集中[1]j的第j个点[11]。

3) 预测卫星状态量^xk+1|k和误差方差阵Pk+1|k。已知k时刻的状态滤波值^xk|k, 考虑卫星动力学模型, 利用8阶龙格-库塔方法积分得k到k+1时刻的状态转移矩阵。

将误差方差阵Pk|k进行乔列斯基分解

由得到的容积点ξi和对应点权值Wi按照以下公式计算k+1时刻的预报状态量和误差方差阵Pk+1|k

整个过程中, 权值的确定和根据动力学模型进行积分是关键。

4) 预报观测量。仿真中采用广播星历拟合出的卫星坐标作为观测量, 由预报状态量^xk+1|k按照式 (12) 和式 (13) 计算经观测方程传播后的容积点Zi*, k+1|k

再由式 (15) 计算观测值估计^zk+1|k

5) 计算增益。根据式 (16) 和式 (17) 分别计算轨道的观测量方差阵、观测量与状态量的协方差阵, 然后求得增益Wk。

6) 更新预测值^xk+1|k+1和预测方差阵Pk+1|k+1

7) 重复以上过程, 可以依次进行预报。利用拉格朗日内插公式可以计算得到非整步长点的卫星状态。

3 仿真结果分析

以GPS星座PRN02号卫星为研究对象, 从2012年11月18日0时整开始的卫星广播星历文件brdc3230.12 n中读取轨道信息, 计算得到24 h内的卫星轨道状态作为观测值, 其误差为2 m左右。分别采用EKF和CKF算法进行轨道预报, EKF算法实现参见文献[3]。预报时间为24 h, 步长15min。同时, IGS发布的精密轨道精度达到2.5cm[10], 可作为参考计算两种预报方法的误差, 精密轨道信息从2012年11月18日的卫星轨道的数据从igr17150.sp3文件中读取。

仿真计算的结果如图2和表1所示。从图2中可以看到, CKF预报结果的误差明显小于EKF。而且, 预报时间越长, 差距越大, 说明CKF比EKF更加稳定。

从表1的数据中可以看到, 两种算法在径向和切向上的误差都远大于法向误差, 三个方向上EKF都较CKF误差更大。

在计算效率方面, 用上述两种方法进行仿真, 在同一台电脑上分别运行20次, EKF和CKF所用平均时间分别为196.5 s和145.8 s, 表明CKF在计算速度上也优于于EKF。

4 结论

CKF有效地避免了EKF中截去高阶项所带来的误差, 同时, 计算中不需要计算复杂的Jacobin矩阵, 节省了计算时间, 也提高了滤波稳定性。仿真中, 由于动力学模型的不完整和积分方法的选择等问题, 两种方法的预报误差都比较大, 但是排除其他因素后, CKF在精度、速度、稳定性三个方面的性能都好于EKF。所以, 基于CKF算法的卫星轨道实时预报方法能够有效改善轨道预报的精度和速度。

参考文献

[1] 周善石, 胡小工, 吴斌.区域监测网精密定轨与轨道预报精度分析.中国科学:物理学力学天文学, 2010; (6) :800—808

[2] 郑作亚, 陈永奇, 卢秀山.UKF算法及其在GPS卫星轨道短期预报中的应用.武汉大学学报 (信息科学版) , 2008; (3) :249 —252

[3] Matzuka B, Aoi M, Attarian A, et al.Nonlinear filtering methodologies for parameter estimation.Department of Mathematics, North Carolina State University.2012:1

[4] 姜伟南.基于UKF的低轨卫星实时定轨方法研究.长沙:国防科技大学, 2007

[5] 穆静, 蔡远利.迭代容积卡尔曼滤波算法及其应用.系统工程与电子技术, 2011;33 (7) :1454—1457, 1509

[6] 李敏.多模GNSS融合精密定轨理论及其应用研究.武汉:武汉大学, 2011:63

[7] Xu Guo Chang.Orbits.Berlin, Heidelberg:Springer-Verlag Berlin Heidelberg, 2008:39

[8] 李智, 徐冬梅, 董绪荣.GPS卫星受力分析与轨道确定.指挥技术学院学报, 1999;10 (2) :95—100

[9] 房媛媛.Cubature-Kalman滤波及性能分析.声学技术, 2011;30 (6) :89—90

[10] 刘伟平, 郝金明, 钱龙, 等.GPS轨道实时预报对比研究.大地测量与地球动力学, 2012; (1) :94—96

基于卡尔曼滤波的土壤水分同化试验第2篇

基于集合卡尔曼滤波的土壤水分同化试验

集合卡尔曼滤波是由大气数据同化发展的新的顺序同化算法,它利用蒙特卡罗方法计算背景场的误差协方差矩阵,克服了卡尔曼滤波需要线性化的模型算子和观测算子的.难点.我们发展了一个基于集合卡尔曼滤波和简单生物圈模型(SiB2,Simple Biosphere Model)的单点陆面数据同化方案.利用7月6日至8月9日青藏高原GAME-Tibet实验区MS3608站点的观测数据进行了同化试验.结果表明,利用集合卡尔曼滤波的数据同化方法可以明显地提高表层、根区、深层土壤水分的估算精度.

作者：黄春林李新 HUANG Chun-lin LI Xin 作者单位：中国科学院,寒区旱区环境与工程研究所,甘肃,兰州,730000 刊名：高原气象 ISTIC PKU英文刊名：PLATEAU METEOROLOGY 年，卷(期)：2006 25(4) 分类号：P33 关键词：陆面数据同化系统集合卡尔曼滤波简单生物圈模型土壤水分

单雷达航迹滤波与卡尔曼滤波算法第3篇

单雷达数据处理是空管自动化系统中的重要功能, 它接收原始的雷达数据进行点迹跟踪、航迹处理, 然后把得到的航迹信息进行分发, 供显示终端进行单路雷达航迹显示, 以及多雷达数据处理模块进行航迹融合。正确实现目标/航迹的关联、精确地对航迹进行滤波估计是提高融合效果的关键。本文介绍了单雷达航迹状态滤波原理, 并分析了目前广泛用于航迹滤波的卡尔曼滤波算法。

1航迹状态滤波

在实际的单雷达数据跟踪过程中, 数据相关与状态估计几乎是同时进行的。准确的预测航迹状态是数据相关的前提, 正确的航迹滤波是实现平滑、稳定航迹的关键, 它是单雷达数据处理中的难点。

航迹预测是在本次航迹滤波数值的基础上按照目标运动模型来估计目标未来的状态。滤波 (后验估计值) 用来估计目标当前的运动参数, 如位置、速率和加速度等。它是由本次观测值和预测估值 (先验估计值) 合并处理, 形成的新的航迹状态参数。对那些没有接收到观测的预测航迹来说, 上一次的预测估计就成了滤波估计。要预测航迹未来状态, 就要定义数学模型来描述与预测有关的物理现象, 把下一时刻的目标状态表示为当前时刻目标状态的函数, 当前的目标状态包括目标位置、目标速度、加速度, 其中加速度可以看作是具有随机特性的扰动输入, 目标运动方程如下:

$x_{k + 1}^{-} = A x_{k} + G a_{k}$

式中:xk是第k次雷达扫描时目标的滤波状态矢量;x-k+1表示根据第k次滤波状态矢量计算出的第k+1次的状态预测;A是转移矩阵;G是噪声增益矩阵;ak是目标的加速度。

一般, 常用的航迹滤波和预测算法有卡尔曼滤波和常系数滤波算法。卡尔曼滤波是根据最小均方误差准则建立起来的估计方法, 它基于包含嗓声的观测值对未知系统状态或参数进行估计, 观测噪声和扰动输入都看作符合高斯过程的白噪声。常用的常系数滤波算法有α-β滤波算法, 它是简化的卡尔曼滤波算法, 其特点为滤波系数是常数而不是动态计算出来的, 因此算法简单运行速度快。卡尔曼滤波有高的准确性, 而常系数滤波算法则有计算上的优势。

这两种算法都可以采用递归调用。过去接收到的数据包含在当前的估计计算中, 因此所有的历史数据都被使用。当有新的观测到来时, 产生“新量”即观测和预测的差, 这种“新量”经过最佳加权去更新目标状态估计, 然后利用动态系统模型和观测模型产生一个最佳的预测, 即在下一次观测到来前作出预测。图1给出了滤波和预测的过程。

图中的增益是用来提供滤波的状态估计和获取下一次扫描时间的预测估计。这些增益可能是固定的, 如常系数滤波算法, 也可能是动态计算出来的值, 如卡尔曼滤波算法。

2卡尔曼滤波算法

卡尔曼滤波算法是根据最小均方误差准则建立起来的估计方法。它用线性递推的方法, 对多个测量数

据和多个信号参数进行处理, 给出无偏差的最小均方误差估计。

卡尔曼滤波器用于估计离散时间过程的状态变量x∈Rn。这个离散时间过程由以下离散随机差分方程描述:

$x_{k} = A x_{k - 1} + B u_{k - 1} + w_{k - 1} (1)$

n×n阶增益矩阵A将上一时刻k-1的状态线性映射到当前时刻k的状态。实际中A可能随时间变化, 但在这儿假设为常数。n×l阶矩阵B代表可选的控制输入u∈Rl的增益。

定义观测变量z∈Rn, 得到量测方程:

$z_{k} = Η x_{k} + v_{k} (2)$

式中:n×m阶矩阵H表示状态变量xk对测量变量zk的增益。实际中H可能随时间变化, 但在这儿假设为常数。

随机信号wk和vk分别表示过程噪声和观测噪声。假设它们为相互独立, 正态分布的白色噪声, 均值为0, 方差分别为Q和R。

定义先验估计误差的协方差P-k和后验估计误差的协方差Pk:

$\begin{array}{l} Ρ_{k}^{-} = E ⌊ (x_{k} - \hat{x}_{k}^{-}) (x_{k} - \hat{x}_{k}^{-})^{Τ} ⌋ (3) \\ Ρ_{k} = E ⌊ (x_{k} - \hat{x}_{k}) (x_{k} - \hat{x}_{k})^{Τ} ⌋ (4) \end{array}$

卡尔曼滤波器中, 先验估计 $\hat{x}_{k}^{-}$ 和加权的测量变量zk及其预测 $Η \hat{x}_{k}^{-}$ 之差的线性组合构成了后验状态估计 $\hat{x}$ k:

$\hat{x}_{k} = \hat{x}_{k}^{-} + Κ (z_{k} - Η \hat{x}_{k}^{-}) (5)$

式 (5) 中测量变量及其预测之差 $(z_{k} - Η \hat{x}_{k}^{-})$ 称为测量过程的新量或残差, 它的大小反映预测值的准确度, 它的方差反映了这种估计的偏离程度。残差为0表明二者完全吻合。n×m阶矩阵K叫做残余的卡尔曼增益或混合因数, 作用是使后验估计误差协方差Pk最小, 即使得后验估计值 (滤波值) $\hat{x}$ k最接近真实值xk。

卡尔曼增益K的计算公式为:

$\begin{array}{l} Κ_{k} = Ρ_{k}^{-} Η^{Τ} (Η Ρ_{k}^{-} Η^{Τ} + R)^{- 1} = \frac{Ρ_{k}^{-} Η^{Τ}}{Η Ρ_{k}^{-} Η^{Τ} + R} \\ (6) \end{array}$

由式 (6) 可以看出, 观测噪声协方差R越小, 新量的增益K越大;先验估计误差协方差P-k越小, 新量的增益K越小。这可以这样理解:随着测量噪声协方差R趋于零, 测量变量zk的权重越来越大, 而zk的预测 $Η \hat{x}_{k}^{-}$ 的权重越来越小;随着先验估计误差协方差P-k趋于零, 测量变量zk的权重越来越小, 而zk的预测 $Η \hat{x}_{k}^{-}$ 的权重越来越大。

卡尔曼滤波器用反馈控制的方法估计过程状态:滤波器估计过程某一时刻的状态, 然后以含噪声的测量变量的方式获得反馈。因此卡尔曼滤波器可分为两个部分:时间更新方程和测量更新方程。其计算过程的工作原理如图2所示。

时间更新方程负责及时向前推算当前状态变量和误差协方差估计的值, 以便为下一个时间状态构造先验估计。测量更新方程负责反馈, 也就是说, 它将先验估计和新的测量变量结合以构造改进的后验估计。

3结束语

卡尔曼滤波算法既适用于暂态过程也适用于稳态过程, 从跟踪精度上来说它是最好的;从计算时间来说, 它的计算时间虽然较其它滤波方法都长, 但也算是在接受范围之内的一种比较简便的方法。因此在目前空管自动化系统的雷达处理模块中广泛应用, 起到关键的作用。

摘要：单雷达数据处理在民航空管自动化系统中起到基础性的作用, 其航迹质量直接关系着自动化系统的准确性和可靠性。航迹状态滤波是单雷达数据处理中的难点, 是实现平滑、稳定航迹的关键。本文以航迹状态滤波算法为切入点, 介绍了广泛应用于自动化系统中的卡尔曼滤波算法。

关键词：单雷达,航迹状态滤波,卡尔曼滤波算法

参考文献

[1]WELCH G, BISHOP G.An introduction to Kalman filtering[EB/OL].UNC-Chapel Hill, 2006-07-24.http://www.cs.unc.edu/～welch/kalman/.

[2]何友, 王国宏, 等.多传感器信息融合及应用[M].北京:电子工业出版社, 2007:19-25.

容积卡尔曼滤波第4篇

卡尔曼滤波是20世纪60年代发展起来的一种现代滤波方法, 它的一个重要作用在于系统的状态估计。当噪声是正态分布时, 这种滤波给出了状态的最小方差估计, 当不是正态分布情况时, 这种滤波给出了状态的线性最小方差估计。

1 滤波算法的发展

约200年前, 高斯 (Guass) 提出了参数估计最小二乘法, 它利用一系列观测数据Z (0) , Z (1) , ...Z (k) , 使它与估计值之间的误差平方和最小, 最小二乘法使用简单, 在一些简单的估计问题中, 仍得到广泛应用[2]。1940年左右, 维纳和柯尔莫格洛夫首次用统计方法研究随机系统, 提出了基于最小均方误差准则下的最优线性滤波一维纳滤波。维纳滤波理论奠定了用统计方法研究随机控制问题的基础。1960年左右, 卡尔曼 (Kalman) 和布西 (Bucy) 在维纳滤波的基础上提出了最优线性递推滤波——卡尔曼滤波。卡尔曼滤波实际上是一种数据处理的递推算法, 它克服了维纳滤波需要整段数据的缺点, 根据前一时刻的估计和现在的数据可以由递推方程得到新时刻的估计, 使存储量和计算量大大减少。它在线性系统中的应用已经很成熟, 对于非线性系统可采用近似方法求解。由于卡尔曼滤波的上述优点, 它在航天技术、导航、通信和经济管理等方面得到了广泛应用[3]。

2 扩展卡尔曼滤波器

2.1 两相静止坐标系下的感应电机的数学状态空间模型

在两相静止α-β坐标系下, 感应电机的4阶数学模型为:

感应电机的转速机械方程为:

φ (X) 和 (4) 式中的系数定义为式中, K为扭转弹性力矩系数;Tγ为转子时间常数;np为极对数, (2) 中Lm为各相间的互感;Ω为转速;is为定子电流, is=[isα, isβ]T;ψγ为转子磁链, ψγ=[ψαγ, ψβγ]T;σ为漏感系数;Ls为定子自感, (4) 中Lm为定子自感;J为转动惯量;Lγ为转子自感;TL为负载力矩;Rs为定子电阻;Rγ为转子电阻。

2.2 扩展卡尔曼滤波器设计模型描述

模型描述为在两相静止α-β坐标系下, 为了利用卡尔曼滤波算法对电机转子转速进行预测, 需加入转子转速ωγ作为新增的状态变量。在采样周期很短的条件下, 可认为ωγ的变化为零, 即ωγ·=0。这样两相静止坐标系下的状态空间模型扩为五阶, 选择定子电流iαs, iβs, 转子磁链ψαγ, ψβγ和转子转速ωγ为状态变量, 感应电机状态空间模型表示[4]:

其中:

式中, Rs为定子电阻;Lm为互感;Ls为定子绕组自感;Lγ为转子绕组自感;ωγ为电机角速度;σ为漏磁系数,

可见α-β两相静止坐标系下, 新的感应电机状态空间模型为一个输入输出由定子电压、电流组成, 系统矩阵由定子和转子参数构成的五阶非线性状态方程。

3 仿真试验

本文用RMxprt软件通过输入相关的电机参数来建立电机结构模型[5], 但是RMxprt软件只能建立单相和三相感应电机的模型, 无法建立多相感应电机的模型。因此, 只能先利用RMxprt建立三相感应电机的模型, 然后再通过改变电机绕组的分布方式, 把三相60°相带感应电机改为六相30°相带的感应电机。图1给出了建立的六相电机的模型图。该电机定子绕组的材料为铜, 转子的导条材料为硅青铜, 定子和转子的冲片的材料为D23, 电机定子的槽数是36, 电机的每极每相槽数为3。图2为在RMxprt软件建立电机模型的基础上, 利用Maxwell软件建立的六相感应电机网格剖分图。

用Maxwell软件进行网格剖分可以大大降低分析人员的工作量, 但是为了进一步提高软件的分析精度, 还需要根据不同的电机部件适当地增加网格剖分密度, 电机模型的网格剖分的各部件的节点数如表1所示。由于电机的气隙、定子齿轭和转子齿轭是电机磁路的组成部分, 因此适当增加节点数目可以得到较高的计算精度。通过求解得到的六相感应电机的计算结果如表2所示。

感应电机A相电流的波形如图3所示, 可见将时间的步长定为0.000 2 s, 计算0.4 s, 仿真得到的电机转矩已经稳定。此时可以得到用扩展卡尔曼滤波器观测到的反电动势波形如图4所示, 由图可以看出电机的反电动势的幅值略低于300 V, 即低于电机额定电压311 V (幅值) 。图5为用扩展卡尔曼滤波器法观测到的六相电机输出转矩特性曲线图, 从图中可以看出在电机进入稳态后, 电机的转矩大小稳定在35 N·m左右。

图6为采用卡尔曼滤波器法观测到的六相电机输出转矩特性曲线图。扩展卡尔曼滤波器对电机低速时的观测效果更加令人满意。

4 结语

本文将卡尔曼滤波用于估计感应电机非线性系统的状态, 设计了扩展卡尔曼滤波器。扩展卡尔曼滤波利用泰勒展开式截断的方式将非线性系统线性化, 线性化以后在形式上同卡尔曼滤波无太多差别。通过对建立的六相感应电机的模型进行电机电磁场的有限元分析, 证明了采用扩展卡尔曼滤波器法能拥有更好的动态跟踪性能。一般来说, 扩展卡尔曼滤波不是最优的, 实际上可以把它作为一种限制复杂性的滤波器。它被限定成具有与线性滤波器类似的结构形式。由于使用了线性逼近, 这种滤波器有可能发散, 在使用过程中尤其要注意这点。

摘要：提出了一种扩展卡尔曼滤波器的设计, 将卡尔曼滤波算法推广至非线性感应电机系统的参数辨识。扩展卡尔曼滤波器采用工程实际中普遍采用的泰勒展开式截断的方法, 对非线性方程进行线性化处理。在仿真中将这种算法应用于六相感应电机模型中, 仿真结果证明扩展卡尔曼滤波器获得了很好的动态跟踪性能以及低速状态下良好的观测效果。

关键词：感应电机,卡尔曼滤波器,扩展卡尔曼滤波器

参考文献

[1]Kailath T.An Innovations Approach to Least-Squares Estimation:Part I.Linear Filtering in Additive White Noise[J].IEEE Trans.Autom.Contr, 1968 (13) :646～655

[2]Kailath T.The Innovations Approach to Detection and Estimaion Theory[J].Proc.IEEE, 1970 (58) :680～695

[3]张贤达.现代信号处理[M].清华大学出版社, 2002

[4]Bose, B.K Power Electronics and AC Drives[M].Prentice-Hall, 1986

基于卡尔曼滤波的需求预第5篇

随着市场竞争的激烈化,各大企业越来越注重组织货源的有效性。为降低库存成本,合理组织货源,订单预测逐渐被提到日程上来,但由于市场的复杂性,多变性,以及各种干扰因素的影响,订单所表现的并非市场的真正需求,其只能是社会需求的一个反映。对此,本文在已知订单的基础上,尝试运用卡尔曼滤波原理,滤掉市场上干扰因素对订单的影响,对社会需求做出最佳估计。期望获得比较贴近市场真实需求的预测结果。为下阶段的货源组织及销售计划制定提供更为合理的依据。

1 卡尔曼滤波需求预测模型

卡尔曼滤波源于Rudolf Emil Kalman在1960年发表的论文《A New Approach to Linear Filtering and Prediction Problems》(线性滤波与预测问题的新方法),是采用由状态方程和观测方程组成的线性随机系统的状态空间模型来描述滤波器,并利用状态方程的递推性,按线性无偏最小均方误差估计准则,采用一套递推算法对滤波器的状态变量作最佳估计,从而求得滤掉噪声的有用信号的最佳估计[1],卡尔曼滤波已广泛应用于各个领域,例如交通方面就有部分学者利用该原理进行短时交通流预测,并且取得了较好的效果。

根据卡尔曼滤波原理我们知道,离散卡尔曼滤波有五个核心方程[2,3,4],对应离散卡尔曼滤波的五个核心方程,即可以得到基于卡尔曼滤波的需求预测方程组:

参数定义:X(k|k-1):k时刻销量的第一估计值;X(k-1|k-1):k-1时刻销量的最优估计值;U(k-1):过程误差(在本文中通过分析上个时刻的估计误差得到);P(k|k-1):k时刻销量第一估计值误差的协方差;P(k-1|k-1):k-1时刻销量最优估计值误差的协方差;A、B:过程参数;Kg(k):卡尔曼增益;Q:过程噪音协方差;R:测量噪音协方差;Q、R会随着时刻的变化而变化,但在实际中一般令其为常数。H:测量参数。A是通过对历史数据的分析,得出的经验数,为方便计算过程参数B和测量参数H均为单位1。

2 模型实现

利用Matlab工具对预测程序进行编写、运行现将算法实现过程简述如下:

1)读取外部数据;

2)对历史数据进行分析,确定经验数A的大小;

3)引入公式(1),根据k-1时刻的销量最优估计值X(k-1|k-1)确定k时刻的销量第一估计值X(k|k-1);

4)引入公式(2),根据前次得到的滤波误差方差P(k-1|k-1),计算预测k时刻的误差方差P(k|k-1);

5)引入公式(3),计算k时刻卡尔曼增益Kg(k);

6)引入公式(4),根据k时刻的订单数据对k时刻的销量最优估计值做出预测;

7)引入公式(5),更新k时刻的滤波误差方差P(k|k);

8)通过计算k-1时刻订单与已统计出的社会真实需求之间的误差,确定k-1时刻的过程误差U(k-1);

9)当k等于1时,X(1[1]),P(1[1]),U(1)没有先验数据,但由于卡尔曼滤波器在递推过程中不断用新信息对状态估计进行修正,所以卡尔曼滤波器是渐近稳定的,也就是说当滤波时间充分长时,状态初值X(1[1])的值对未来估计的影响将衰减至近于零,初始协方差P(1[1])对未来滤波最优估计值误差协方差的影响也将衰减至近于零.因此,滤波的初始条件用近似的算法确定。在本文中令X(1[1])为历史同期的数据,P(1[1])为单位1,U(1)为0;

10)至此一次滤波循环结束,当前时刻由k变为k+1,返回第二步进入下一轮滤波循环;

3 预测实例与结果分析

本文利用云南某烟草公司2006年某品牌卷烟的销售数据和2007年该品牌卷烟的订单数据对模型进行实例验证,期望滤掉2007年订单数据中的干扰因素,以便获得2007年该品牌卷烟的社会真实需求,该实例以半月为周期,共计24个时间点。

由图1可得当Q=0,R=0时,即过程噪音协方差和测量噪音协方差均为0时,此时不存在干扰因素,故预测值即是订单值,两者的曲线除去初始点外全部相同。但是实际中Q,R不可能为0,干扰因素必然存在,因此在图1中订单数据(预测值)与社会真实需求之间有较大的差距。当令Q=1,R=1时,通过图2并对比图1可以发现,利用卡尔曼滤波模型对订单数据进行除噪后得出的预测数据开始偏离订单数据而向社会真实需求的方向靠拢,虽然仍有一些差距,但是较之图1已经有了部分改善。下面就从定量的角度进一步说明该模型的可行性。

4 误差分析

为了能评价和比较仿真实验结果,本文使用如下几个误差指标对模型结果进行误差分析:

上述指标中marerr为检测值分数的预测误差,maxrer为惩罚大的预测误差,因而提供更好的模型评价参数,EC表示预测值与实测值的拟合度,一般拟合度越接近1表示预测结果与实际拟合的越好。利用模型得出的结果及上述四个误差指标,可以得到误差分析表,如下表示。

在进行预测时因为初始数据具有很大的随意性,因此在误差分析时舍去k=1时对应的数据。仅仅对剩余的23个数据进行误差分析。

方案一为不对干扰因素进行过滤,直接以订单作为社会需求时,该需求与社会真实需求之间的各误差指标。

方案二则利用了上述卡尔曼滤波模型对订单数据进行处理,得出了基于卡尔曼滤波的最优社会真实需求,继而得出了预测数据与社会真实需求之间的各误差指标。

通过定量对比我们可以发现,利用卡尔曼滤波原理对订单数据处理后,除平均相对误差外其余各误差指标均优于未处理前的订单数据。

5 总结及展望

以往人们大多对订单进行预测,将订单作为了社会的真实需求而忽视了订单与社会真实需求的不同。本文首次尝试在订单数据的基础上以卡尔曼滤波模型为工具对社会真实需求进行了预测,并取得了一些成效,获得了较之订单更贴近社会真实需求的数据。

虽然本文在某些方面取得一些成果,但是仍有诸多不足需要进一步分析研究:

1)在本文中为方便计算过程参数B和测量参数H均取单位1,这与实际不相符,如何对历史数据进行分析,确定较合适的过程参数B和测量参数H是下一步要研究的工作;

2)在本文中过程噪音协方差Q和测量噪音协方差R均取某一常数,这与实际也不相符,其应随着时刻的不同而变化,如何建立过程噪音协方差Q和测量噪音协方差R与时刻k之间的函数关系也需要进一步的研究。

参考文献

[1]杨兆升,朱中.基于卡尔曼滤波理论的交通流量实时预测模型.中国公路学报,1999;13(3):63—67

[2]王明俊,刘开生,张起森,等.高速公路路段费用的卡尔曼滤波预测.系统工程,2001;19(4):76—79

[3]Okutaki I.Dynamic prediction of traffic volume through kalman filte-ring theory.Transp Res,1984;18(2):1—11

基于卡尔曼滤波的背景更新算法第6篇

1 基于卡尔曼滤波的背景更新算法

1.1 背景更新算法

通过背景差分法进行运动目标检测,虽然可以提取出完整的运动目标区域,但是随着时间的推移,光照和阴影会产生变化,导致最初确定的背景图像与真实背景存在一些误差,最终检测结果当中会出现很多扰动目标,严重影响了目标检测的效果。为了减小因为光照和阴影变化产生的影响,常用的办法是不断地进行背景更新。

1.2 卡尔曼滤波

卡尔曼滤波是Rudof Emil Kalman提出的一种高效率的递归滤波器,源自于贝叶斯估计理论,是一种建立在线性最小方差估计的基础上的一种算法。以估计为基础,将状态方程的概念引入,把信号过程视为白噪声作用下的一个线性系统的输出,建立系统的状态模型,通过状态方程表达系统的输入输出关系,最后通过系统状态转移方程,对输入输出状态进行时间和观测更新,从而估计出目标的状态,卡尔曼滤波算法不需要所有的观测数据,只根据当前帧的最优估计值和这一帧的实际观测值来估计下一帧背景的实际值,由此利用卡尔曼滤波算法可以对变化的背景进行实时估计。

1.3基于卡尔曼滤波的背景更新算法实现

通过卡尔曼滤波的时域低通滤波对视频检测中的背景图像进行更新可以被看做是一个多帧降噪的过程,其原理如图1所示。

当视频图像序列通过这个时域低通滤波器时,由于图像的快速变化过程中,图像序列的缓变部分被分离出来。其时域递归低通滤波器公式如下:

其中,B(i,j,k) 是当前检测的第K帧背景图像,B(i,j,k+1) 是下一帧也就是K+1帧的背景图像,I(i,j,k) 表示在K帧时,坐标(i,j)图像序列的特征值,M(i,j,k) 是当前运动目标图像被二值化以后的结果,背景因子是β ,运动因子是α是,Th是检测阈值,avg1是差值图像的平均值,avg2是当前第k帧图像的平均值。M(i,j,k) 为:如果|I(i,j,k )-B(i,j,k)| >Th ,则M(i,j,k)=1;否则M(i,j,k)=0。δ1、δ2是相应图像的标准方差。

根据卡尔曼滤波理论的原理,这个公式其实是关于一种背景的递归预测和更新,首先通过预测得出一个值,然后通过这个值和相关项来计算出新背景的更新值,但这种背景预测有一个大前提,运动目标的变化要远大于背景的变化,,α和β是一组常量,决定着背景提取是否自适应,一般取值需要足够小,取0到1之间,这样可以有效的吧运动目标从背景图像序列当中分离出来。在实际的应用当中,时间因素还需要被考虑进去,α和β的取值在实际中不需要通过公式来进行严格计算,完全可以根据经验来选取一些合适的值。该文中取α=0.048,β=0.13。而且在背景预测模型当中,初始帧的背景信息一般都是未知的,只能通过图像传感器来获取。所以当k=0时,取I(i,j,k)=B(i,j,k) ,即M(i,j,k)=0。M(i,j,k) 是已提取的经过二值化的运动目标轮廓,因为M(i,j,k)=0 ,这表示在运动目标检测初始化的时刻即k=0时,图像序列当中没有出现运动的目标,在这一时刻运动目标轮廓的值为0。

2实验结果与分析

我们在实验室拍摄一组运动目标图像序列,基于背景差分法通过matlab仿真对其进行运动目标检测实验,在实验过程中加入光照等条件,不断地进行背景的更新,最后提取出运动目标轮廓。

实验结果如图2所示,第一张图为实时更新的背景图像,第二张图是第168帧的目标检测图像,第三张图是对提取出的运动目标进行二值化和滤波处理后的检测结果,由于对已提取出的前景目标用外接矩形框采取了标记,所以其目标区域背景变化较小, 对非标记的目标区域进行实时的更新,可以有效地减小误差值,提高运动目标检测的精确度。

3 总结

基于卡尔曼粒子滤波的目标跟踪算法第7篇

根据跟踪目标的不同,目标跟踪算法通常可以分为3类[2],即基于偏微分方程的目标跟踪算法、基于Mean Shift的目标跟踪算法和基于滤波理论的目标跟踪方法。文中主要对基于滤波理论的目标跟踪方法进行讨论,这种算法实质是将序列图像的目标跟踪问题建模为目标状态的概率密度函数在时间序列上的传播,通常包括先验密度、观测密度以及后验密度3大要素。

1 基于滤波理论的3种算法描述

1.1 Kalman滤波算法描述

Kalman滤波器[3]是由卡尔曼(R.E.Kalman)于1960年提出的一种递推估计器。Kalman滤波是一种在时域内采用递归滤波对系统状态进行最小均方误差估计的方法,根据系统以前的状态对下一个状态作最优估计,预测时具有无偏、稳定和最优的特点[4]。

Kalman滤波使用要求比较严格,必须满足3个条件[5]:动态系统是线性的;观测是状态的线性组合;噪声是附加的,而且服从正态分布。当这些要求有几个不满足时,Kalman滤波尽管仍然可以使用,但却不能显示其优越性了。

1.2 粒子滤波算法描述

粒子滤波[6]的思想基于蒙特卡洛方法(Monte Carlo Methods),它是利用粒子集来表示概率,也就是通过随机抽取的加权粒子来代替状态的后验概率分布,它可以用于任何形式的状态空间模型上。当随机抽取的粒子数量很大时,结果也就无限接近于状态变量实际的后验密度。

在粒子滤波中,针对许多实际情况,无法直接从后验密度中采样粒子。因此,通常采用与后验密度相近的某一提议密度π(xk|xk-1,y1∶k),从该提议密度中采样新的粒子。在每一离散时刻,新的粒子从提议密度中进行采样,并且计算权重以补偿提议密度和后验密度的差异。将权重归一化之后,后验密度由一组粒子进行近似。

1.3 Kalman粒子滤波算法

由于标准粒子滤波算法选择先验概率密度作为重要密度函数,这种选取方法只有在测量精度要求不高的场合才能取得较好的结果。重要性权重的方差随着时间而随机递增,使得粒子的权重主要集中在少数粒子上,甚至在经过多步递归后,可能只有一个粒子有非零权值,而其他粒子的值很小,几乎可以忽略不计,从而使得大量的计算工作都被浪费在用来更新那些对p(xk|z1∶k)的估计不起任何作用的粒子上,结果粒子集无法表达真实的后验概率分布,这就是粒子滤波算法的退化问题[7]。粒子退化是粒子滤波算法中不可避免的现象,解决方法是选取好的重要性密度函数和粒子重采样。

基于以上原因,文中提出基于卡尔曼粒子滤波器的算法,卡尔曼滤波器选取高斯分布作为提议分布,采用最优提议分布线性化的方法可以得到较优的提议分布。对于每一个粒子,一个单独的卡尔曼滤波器用来产生和传播高斯提议分布。用卡尔曼算法可以得到一个更好的重要密度函数,能够使先验分布朝着高似然度区域移动,在一定程度上避免了粒子退化问题,提高了粒子滤波算法的估计精度。

2 卡尔曼粒子滤波在跟踪中的应用

基于卡尔曼粒子滤波器的跟踪算法流程图如图1所示。

(1)初始化。前文中已经提到,粒子滤波器的基本思想是蒙特罗卡模拟,其中状态的后验密度由一组有权重的粒子 ${(\tilde{x}_{k - 1}^{(i)}, W_{k - 1}^{(i)}), i = 1, \dots, Ν}$ 来近似。从先验密度p(x $_{0}^{(i)}$ )中采样,得到初始(令k=1)带权重的粒子 ${(\tilde{x}_{0}^{(i)}, Ρ_{0}^{(i)}), 1 / Ν, i = 1, \dots, Ν}$ ,其中k表示第k帧, $\tilde{x}_{k - 1}^{(i)}$ 为对系统状态x $_{k - 1}^{i}$ 的估计,W $_{k - 1}^{(i)}$ 为相应的权重,p $_{0}^{(i)}$ 是第i个粒子的协方差矩阵。

(2)采样阶段。对于i=1,…,N有

1)卡尔曼预测。

$\overset{⌢}{x}_{k / k - 1}^{(i)} = A \tilde{x}_{k - 1}^{(i)} + (Ι_{n} - A) \bar{x}$ (1)

$\overset{⌢}{Ρ}_{k}^{(i)} = A Ρ_{k - 1}^{(i)} A^{Τ} + Q$ (2)

其中,A是n×n维状态转移矩阵,由 $\tilde{x}_{k + 1 / k} = A \tilde{x}_{k / k}$ 得到;P为协方差矩阵;Q为控制点向量矩阵。

2)基于预测状态x $_{k / k - 1}^{(i)}$ 进行观测得到观测新息v $_{k}^{(i)}$

3)应用新息矩阵H计算卡尔曼更新

S $_{k}^{(i)}$ =HP $_{k}^{(i)}$ HT+σ2Im (3)

$Κ_{k}^{(i)} = \overset{⌢}{Ρ}_{k}^{(i)} Η^{Τ} (S_{k}^{i})^{- 1}$ (4)

$\overset{⌢}{x}_{k | k}^{(i)} = \overset{⌢}{x}_{k | k - 1}^{(i)} + Κ_{k}^{(i)} v_{k}^{(i)}$ (5)

$Ρ_{k}^{(i)} = \overset{⌢}{Ρ}_{k}^{(i)} - Κ_{k}^{(i)} S_{k}^{(i)} (Κ_{k}^{(i)})^{Τ}$ (6)

其中,S为函数空间,H为新息矩阵,

$Η = [\begin{matrix} B (S)^{Τ} \\ B (S)^{Τ} \end{matrix}], B$

为B样条基函数,σ2为方差,K为卡尔曼增益。

4)从下面的提议密度中采样得到粒子x $_{k}^{(i)}$

$π (x_{k}^{(i)} | \tilde{x}_{k - 1}^{(i)}, y_{1 ∶ k}) = Ν (\overset{⌢}{x}_{k | k}^{(i)}, Ρ_{k}^{(i)})$ (7)

其中,N为提议密度矩阵。

(3)计算新的权重。

(4)输出阶段。

(5)重采样阶段。

(6)k=k+1,若继续跟踪,转到步骤(2);否则,算法结束。

在卡尔曼粒子滤波器中,在采样阶段每个独立的粒子按照高斯过程进化,这个过程同时考虑了状态转移密度和观测密度。高斯提议密度表征了随机粒子的一阶和二阶统计特性,从该提议密度中进行采样得到新粒子与父粒子具有相似的分布,因此继承了父粒子的协方差。每个粒子的状态和协方差可以看作后验密度中的一个高斯分量密度,并按照卡尔曼滤波器进化。

在粒子滤波器中,随着粒子数目的增加,算法的性能逐渐改善,但计算量也变得越来越大。因此,在实际的目标跟踪中,需要对算法的性能和计算效率进行分析,找到临界点。下面就通过临界粒子数目对3种不同的跟踪器进行性能和效率比较。

当粒子数为500时,卡尔曼滤波和粒子滤波都取得了较好的结果。对于这两种跟踪器,逐渐减少其粒子数,直到得到临界粒子数目。采用一个经典的非线性例子[8]进行仿真,比较结果如表1和图2所示。

状态方程

xk=0.5xk-1+25xk-1/(1+x $_{k - 1}^{2}$ )+8 cos 1.2k-1+rk (8)

观测方程

pk=x $_{k - 1}^{2}$ /20+qk (9)

其中,rk和qk分别代表系统噪声和观测噪声,且都是均值为零的高斯白噪声。

由表1可知,粒子滤波和卡尔曼粒子滤波两种算法在各自的临界数目下跟踪误差大致相同,而卡尔曼粒子滤波临界粒子数远少于其他两种滤波,同时每帧跟踪时间也大幅减少。

由图2(a)可知,此时的卡尔曼滤波算法已经不能得到准确结果,而卡尔曼粒子滤波的结果明显更接近于真实值。图2(b)中卡尔曼粒子滤波的估计值误差远低于其他两种跟踪算法,该算法的执行效率和跟踪误差都优于卡尔曼滤波器和粒子滤波器。图3显示卡尔曼粒子滤波器的估计值都在真实值的95%置信区间内。综上可知,基于卡尔曼粒子滤波器的算法能够满足本系统的实时性和鲁棒性的要求。然而,卡尔曼粒子滤波器在硬件系统中实现起来比其他两种算法更复杂。

3 实验结果分析

实验素材选取于一段事先拍摄好的视频。由图4的跟踪结果可以看出,检测运动目标的标识符为白色“十”字符,而该 “十”字符中心点基本能够与目标的中心点重合,达到了预期的目标。由此可知,本文采用的基于卡尔曼粒子滤波的跟踪算法能较好地实现运动目标的跟踪,有效达到预期目的,有良好的实际应用价值和推广价值。

4 结束语

主要研究了基于滤波理论的目标跟踪方法,首先简单介绍了目前常用的基本算法,并对基于滤波理论的3种算法进行了比较分析,最后在Matlab平台上进行了仿真实验,实验结果表明Kalman 粒子滤波算法能够很好地对目标进行跟踪,并在实时性、鲁棒性等方面较为突出,比其他算法更有效。

参考文献

[1]刘卫光,李广鑫.一种通用的视频目标跟踪系统设计[J].计算机技术与发展,2009,19(10):110-113.

[2]NIXON M S,HARRIS C J.Statistical gait description viatemporal moments[C].Proceedings of Proc 4th IEEE South-west Symposium on Image Analysis and Interpretation,2000:291-295.

[3]王利循.无迹Kalman滤波在IMU和GPS组合导航系统中的应用研究[D].南昌:南昌大学,2010.

[4]程娟.复杂背景下运动目标识别算法研究[D].武汉:武汉理工大学,2008.

[5]ZHOU Feng,MENG Xiuyun.Adaptive kalman filter of transferalignment with un-modeled wing flexure of aircraft[J].Jour-nal of Beijing Institute of Technology,2008,17(4):434-438.

[6]朱娟.蒙特卡洛滤波算法在目标跟踪中的应用[D].长春:中国科学院研究生院,2010.

[7]LI P H,ZHANG T W,ARTHUR E C P.Visual contour track-ing based on particle filters[J].Image and Vision Compu-ting,2003(21):111-123.

容积卡尔曼滤波第8篇

1 无源定位原理

通常卫星导航系统需要至少观测4颗卫星,列写4个独立导航方程,才能完成三维位置和用户钟差解算。但北斗双星系统仅有两颗工作星,无法得到4个独立方程。因此,要想使双星系统从有源方式转换为无源方式,根本途径是增加导航源[2],方法之一是增加高精度时钟作为时间源。其基本思想是用户终端配备高精度铷原子钟。初始时刻进行用户钟与中心站钟同步,并设用户钟稳定度足够高。在载体工作期间其频率漂移足够小,把钟差引起的测量误差降低到最小。由测得的卫星伪距和数字高程信息,即可计算用户的空间位置。

在机载环境下,为了保证北斗接收机始终能够连续定位,不出现长时间数据丢失,在进行无源定位时,导航接收机在接收到信号后,利用伪距通过最小二乘算法进行定位。根据接收卫星的接收情况,有两种解算方法:一种是收到3颗卫星时的定位情况;另一种是收到两颗卫星时的定位情况。定位解算方程可以分别描述下:

(1)双星无源定位观测方程组

${\begin{cases} ρ^{i} = \sqrt{(x - x_{i})^{2} + (y - y_{i})^{2} + (z - z_{i})^{2}} + C (Δ t + Δ t_{t}^{i} + Δ t_{s}^{i}), i = 1, 2 \\ \frac{x^{2} + y^{2}}{(R_{e} + h)^{2}} + \frac{z^{2}}{(R_{p} + h)^{2}} = 1 \end{cases} (1)$

其中,Re为地球赤道半径;Rp为地球极半径;根据WGS-84模型,Re=6 378 137 m,Rp=Re(1-f);f为地球的椭圆度,f=1 298.257;h为用户的海拔高度。

伪距测量值ρi在铷钟经初始时间校正量Δt+Δt $_{s}^{i}$ 和大气层时延校正量Δt $_{t}^{i}$ 修正后得到用户至卫星的“真实”距离,利用最小二乘法可以解算用户位置,但无法计算钟差,定位精度会有一定损失。

(2)三星无源定位观测方程组

${\begin{cases} ρ^{i} = \sqrt{(x - x_{i})^{2} + (y - y_{i})^{2} + (z - z_{i})^{2}} + C (Δ t + Δ t_{t}^{i} + Δ t_{s}^{i}), i = 1, 2, 3 \\ \frac{x^{2} + y^{2}}{(R_{e} + h)^{2}} + \frac{z^{2}}{(R_{p} + h)^{2}} = 1 \end{cases} (2)$

利用最小二乘法解算3星无源定位观测方程组,可以得到用户位置和钟差,消除了铷原子频标频率漂移造成的伪距测量误差,提高了定位精度。实际应用中,东星或西星工作异常的概率很小,备份星在大多数时间内都可以作为独立的资源,因此,3星无源定位成为本系统的主要工作模式。

2 试飞数据分析

2.1 机载数据处理

根据试飞实际情况,对载机不作较大改动的情况下,在飞机前设备舱加装北斗接收机及天线,天线位置位于前设备舱盖中心线偏左。载机上同时装有GPS接收机和飞行参数记录仪。对数据进行分析,首先利用时间相关算法确定机载飞行参数记录仪中飞行参数和北斗无源系统记录数据的同步位置,其次判断北斗无源定位系统的工作状态。最后,根据机载记录设备输出的GPS高度和北斗无源设备记录的伪距值进行组合导航解算,输出计算结果。

2.2 卡尔曼滤波方程建立与解算步骤

为消除动态北斗定位数据中的随机误差,需要将真实的状态从各种随机干扰中的估计出来。在动态数据处理中,采用卡尔曼滤波方法来提取所需要的数据,卡尔曼滤波方程建立过程[3,4,5]。

建立离散型卡尔曼滤波方程如下

${\begin{cases} 状态方程 X_{k} = Φ_{k, k - 1} X_{k - 1} + Γ_{k, k - 1} W_{k - 1} \\ 观测方程 Y_{k} = Η X_{k} + V_{k} \end{cases} (3)$

其中,Φk,k-1,Γk-1矩阵分别是系统的状态转移矩阵,系统噪声矩阵。Wk-1,Vk分别是方差阵为Qk-1,Rk的白噪声系列。

对式(3)进行卡尔曼滤波,得

状态预测 $\hat{X}_{k, k - 1} = Φ_{k, k - 1} \hat{X}_{k - 1}$

状态滤波 $\hat{X}_{k} = \hat{X}_{k, k - 1} + Κ_{k} (Ζ_{k} - Η \hat{X}_{k, k - 1})$

滤波增益Kk=Pk,k-1HT(HPk,k-1HT+Rk)-1

预测协方差Pk,k-1=Φk,k-1Pk-1ΦTk,k-1+Γk,k-1Qk-1ΓTk,k-1

滤波协方差Pk=(I-KkH)Pk,k-1

其中,Φ,Γ,H,根据实际的情况确定为。其中,τ为数据采样时间间隔,GPS为每秒10次,即0.1 s;BD为每秒8次,即0.125 s。

航迹的初始条件设定:飞机初始位置为北纬40°,东经99°。伪距测量的误差涉及到卫星时钟误差、卫星星历误差、电离层延迟误差,对流层延迟误差、多经效应误差和接收机误差等。上述误差经过修改及补偿后可大大降低伪距测量误差,假设伪距测量噪声是标准差为100 m的白噪声。

3 北斗定位精度分析

以差分GPS输出数据为实际测量基准数据,设第i个测量点,差分GPS输出的基准数据与北斗输出的定位数据经纬度误差值为(Δλi,Δφi),单位(′),Ti为系统的采样间隔,径向长度误差为 $Δ L_{i} = 1 852 \sqrt{Δ φ_{i}^{2} + (Δ λ_{i} \cos φ_{i})^{2}}$ ,径向长度误差率为 $Δ R_{i} = \frac{Δ L_{i}}{Τ_{i}}$ 。北斗无源定位方式下,选取考核点±5分钟导航时间内,有效的测量点数量n累计≥340个,其径向误差率标准差计算公式为

$\sqrt{\frac{\sum_{i = 1}^{n} Δ R_{i}^{2}}{n}}$

大误差资料剔除原则:按照数据统计理论,误差值超过标准差值的3倍以上,则认为是大误差数据,在分析误差产生原因之后予以剔除。

以下试验结果是以铷钟型无源接收机接收数据为例。

3.1 飞行轨迹

图1和图2给出了飞机在地球地理坐标下的实际飞行的二维和三维飞行轨迹,飞机在东经101°～101.2°,北纬40.72°～40.83°;东经99.48°～99.52°,北纬40.79°～40.82°等区域处于大坡度转弯状态,北斗接收机天线受到遮挡,在对应时间段内有丢失卫星信号现象发生。根据本文第一部分无源定位原理可知,只要能够接收到两颗北斗卫星信号,北斗接收机就会有定位数据输出,定位解算根据式(1)进行,但定位精度会有一定损失。飞机在大部分的飞行区域内,北斗接收机可以接收3颗定位卫星信号,定位解算根据式(2)进行。在整个飞行过程中,由于3颗北斗卫星位置在赤道上空,北斗接收机接收信号会受到飞机飞行姿态的影响,不同飞行阶段飞机的姿态不同。下面针对不同的飞行阶段,进行分段处理。

3.2 试飞数据的分段处理

分析表明,在平飞阶段,数据较好,在拐弯(由东向转向北向)时,北斗接收机天线受到遮挡(接收不到东星信号),接收卫星信号不好,没有处理。

飞行资料分段分析结果如表1所示。

总共采样点数26 263,其中有效点数24 496,数据有效点数为93.27%,无源定位精度在3σ条件下为92.842 1 m。

(1)起飞爬升阶段数据分析以第一数据段为例,如图3所示。

(2)平飞阶段数据分析以第二数据段为例,如图4所示。

对飞行轨迹分析可以得出如下结论:

(1)飞行过程中,北斗和GPS系统的定位曲线重合较好,两者位置相差的平均距离为92.842 1 m,具有较好的位置精度,说明设计实现的北斗接收机和滤波算法,在飞机巡航状态下可以满足空中一般用户需要。

(2)在北斗接收机丢失信号的时间段内,位置数据开始有发散的趋势,但飞机的大坡度转弯时间不长,只要接收到两颗卫星信号的情况下,北斗接收机仍然有数据输出,没有影响到飞行轨迹曲线的连续性。

(3)在基于北斗无源定位系统的卡尔曼滤波算法中,能够根据北斗三星的有效位置信息进行估计和补偿,提高北斗导航系统的精度,具有较好的工程应用价值。

(4)如要进一步提高北斗三星无源系统定位精度,则需要进行两方面的深入研究。需要改进完善定位解算方法,进行自适应算法的研究。需要对气压式高度表辅助定位时,引起系统误差及其修正算法进一步研究。

4 结束语