CN106871905B - 一种非理想条件下高斯滤波替代框架组合导航方法 - Google Patents

一种非理想条件下高斯滤波替代框架组合导航方法 Download PDF

Info

Publication number
CN106871905B
CN106871905B CN201710120792.XA CN201710120792A CN106871905B CN 106871905 B CN106871905 B CN 106871905B CN 201710120792 A CN201710120792 A CN 201710120792A CN 106871905 B CN106871905 B CN 106871905B
Authority
CN
China
Prior art keywords
time
state
covariance
noise
value
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201710120792.XA
Other languages
English (en)
Other versions
CN106871905A (zh
Inventor
宋申民
谭立国
赵凯
张秀杰
吴骁航
司译文
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Harbin Institute of Technology
Original Assignee
Harbin Institute of Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Harbin Institute of Technology filed Critical Harbin Institute of Technology
Priority to CN201710120792.XA priority Critical patent/CN106871905B/zh
Publication of CN106871905A publication Critical patent/CN106871905A/zh
Application granted granted Critical
Publication of CN106871905B publication Critical patent/CN106871905B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C21/00—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00
    • G01C21/20—Instruments for performing navigational calculations
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/38—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system
    • G01S19/39—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system the satellite radio beacon positioning system transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/42—Determining position
    • G01S19/48—Determining position by combining or switching between position solutions derived from the satellite radio beacon positioning system and position solutions derived from a further system
    • G01S19/49—Determining position by combining or switching between position solutions derived from the satellite radio beacon positioning system and position solutions derived from a further system whereby the further system is an inertial position system, e.g. loosely-coupled

Abstract

本发明涉及一种非理想条件下高斯滤波替代框架组合导航方法，属于航天导航系统领域，为了解决现有的高斯滤波算法将噪声和状态分开处理，导致计算量大和估计精度不高的缺点，而提出一种非理想条件下高斯滤波替代框架组合导航方法，满足：k+1时刻的状态关于前k步量测值的一步后验概率密度函数是高斯的；k+1时刻的实际量测值关于前k步量测值的一步后验预测概率密度函数是高斯的；所述方法包括：建立状态在k+1时刻的预测均值和协方差表达式；建立状态在k+1时刻的后验均值和协方差表达式；基于CKF滤波算法对步骤一和步骤二得到的表达式进行求解，得到修正后的均值和方差；根据修正后的均值和方差修正航向，实现导航。本发明适用于航天飞行器导航系统。

Description

一种非理想条件下高斯滤波替代框架组合导航方法
技术领域
本发明涉及一种非理想条件下高斯滤波替代框架组合导航方法，属于航天导航系统领域。
背景技术
导航是引导飞行器、船舶或汽车等沿一定航线从一点运动到另一点的方法，其在军用领域和民用领域都有着广泛的应用。现代导航系统种类繁多，例如全球定位系统(GPS)、天文导航(CNS)、多普勒测速系统(Doppler)、惯性导航系统(INS)等[1]。由于GPS和INS两者之间具有较强的非相似性和互补性，将它们组合起来，便可以取长补短充分发挥各自的优势，同时克服GPS易受地形地物遮挡而导致定位中断和INS定位误差随时间而累计的缺陷。在GPS/INS组合导航系统中，由于系统本身元器件的不稳定性以及外部应用环境不确定因素的影响，导致系统噪声有时具有相关特性，比如当系统的输入源与传感器的输出且传感器测量值具有随机特性时，系统的过程噪声和量测噪声协方差将不为零。对于采用伪距的表达方式的紧耦合模型，并且由于实际应用中通信带宽的限制，使得组合导航系统可能同时具有非线性、时滞和噪声相关的特性，因此，设计更一般的算法是十分有必要。
下面简要介绍上述三类问题解决方案的发展现状并由此引出本文的算法。
针对常见的非线性系统，如制导或导航系统[2,3]，目标跟踪系统[1]等，在贝叶斯框架下依据概率密度非线性滤波算法可以分为两类。一类是高斯滤波算法[8]，如EKF[5]是以非线性函数的一阶泰勒展开形式来近似函数本身，对于强非线性函数近似能力较差，且存在计算Jacobi矩阵的问题，计算量大，不适合实时计算；UKF[6]算法是利用无迹变换来逼近状态后验概率密度函数，由于利用了样本点，因此不再需要计算Jacobi矩阵，但是在参数选择不当的时候容易造成状态误差协方差负定；Nrgaard M等提出的DDF算法[9]以斯特林多项式插值的方式去近似非线性函数，避免了陷入局部线性化的问题；CKF[7]算法是利用球径容积法则来实现对状态后验概率密度函数的近似，相比于UKF算法，减少了一个样本点但是数值稳定性却有了很大的提高。另一类是非高斯滤波算法，例如PF，但是计算量很大，不适合实时计算。
针对噪声相关问题，例如迎风航天器模型的过程噪声和量测噪声[10]，部分学者提出了解相关框架[11]，他们在状态方程中加入前一时刻的量测方程，将新构造出伪状态方程作为新系统的状态方程，通过选取合理的增益系数，达到新系统过程噪声和量测噪声解耦的目的，从而利用标准KF框架进行求解；另一部分学者采用两步预测代替一步预测的方式[12]，从而使得噪声之间相关性不再存在达到解决问题的目的；后来，经学者GuobinChang指出这两种方法是等价的[13]。文章[14]通过将噪声看做是状态增量联和状态一起进行估计，给出了一种新颖的解决噪声相关问题的方法，并从理论上证明了其与前两种方法的等价性，仿真说明了此方法对解决线性系统且噪声相关问题的有效性。
针对量测量不能即时获取以及过程噪声和量测噪声相关的情况，研究主要从两方面展开，一类是确定性时滞，一类是随机时滞。针对确定性时滞问题目前主要的解决方法有重计算法、Alexander法和状态增量法等[15]；对于随机时滞问题，一般以满足Bernoulli分布的随机序列来构造新的量测数据，Hermoso-Corazo A等[16]针对一步随机时滞和两步随机时滞状态估计问题，给出了其在EFK算法和UKF算法框架下的处理方法。Wang X等[17]利用对状态后验概率密度估计的两步预测思想，融合状态增量的方法，解决了带有随机时滞和相关噪声的状态估计问题。
虽然针对这些非理想情况已有丰硕的研究成果，但是对于多种非理想条件共存的系统研究成果甚少，或者给出的算法形式不够简洁或统一，因此需要一种新的方法，从设计算法的角度，在新的框架下给出此类问题的解决方案。
发明内容
本发明的目的是为了解决现有技术的替代框架下的高斯滤波算法将噪声和状态分开处理，导致估计精度不高，并且不具有一般性的缺点，而提出一种非理想条件下高斯滤波替代框架组合导航方法，所述方法基于非线性离散系统模型，并且满足如下条件：
k+1时刻的状态关于前k步量测值的一步后验概率密度函数是高斯的；
k+1时刻的实际量测值关于前k步量测值的一步后验预测概率密度函数是高斯的；
其特征在于，所述方法包括如下步骤：
步骤一：根据飞行器内的传感器测量值，获取飞行器在k时刻的状态值xk；
步骤二：建立状态在k+1时刻的预测均值
Figure BDA0001236972260000021
和协方差
Figure BDA0001236972260000022
表达式：
Figure BDA0001236972260000023
Figure BDA0001236972260000031
其中，
Figure BDA0001236972260000032
为k时刻对k+1时刻的预测估计值；
Figure BDA0001236972260000033
是k+1时刻的n维状态向量的协方差矩阵；Qk+1为k+1时刻的过程噪声ωk+1的协方差矩阵；Rk+1为k+1时刻的量测噪声νk+1的协方差矩阵；Sk+1为k+1时刻的过程噪声和量测噪声的协方差矩阵，且Sk+1≠0；
步骤三：建立状态在k+1时刻的后验均值
Figure BDA0001236972260000034
和协方差
Figure BDA0001236972260000035
表达式：
Figure BDA0001236972260000036
Figure BDA0001236972260000037
其中，
Figure BDA0001236972260000038
为增广后的状态向量的新息，
Figure BDA0001236972260000039
其中
Figure BDA00012369722600000310
为k+1时刻的新息nk+1与k+1时刻的预测均值
Figure BDA00012369722600000311
的互协方差；
步骤四：基于CKF滤波算法对步骤二和步骤三得到的表达式进行求解，得到修正后的均值和方差。
步骤五：根据所述修正后的均值和方差修正航向，实现导航。
本发明的有益效果为：本发明不同于以往噪声和状态分开处理的情况，本发明将过程噪声和量测噪声一起做为状态的一个扩张量进行滤波更新，将已知的噪声信息作为先验信息，在状态和量测的更新方程中用更新后的噪声的后验估计代替先验信息；通过仿真证明了其估计精度优于现有的滤波算法。同时，本发明提出的算法更具有一般性，选取合适的参数可以看出，一般高斯滤波算法是其在噪声不相关或不存在随机时滞情况的特例。
附图说明
图1为本发明的非理想条件下高斯滤波替代框架组合导航方法的流程图；
图2为EKF、CKF、UKF和CKF-CNRD的状态估计误差对比图；横坐标为时间，纵坐标为状态估计误差；
图3为EKF、CKF、UKF and CKF-CNRD的状态估计均方差对比图；横坐标为时间，纵坐标为状态估计均方差；
图4为CKF-CNRD的噪声估计均值效果图；横坐标为时间，纵坐标为噪声估计均值；
图5为CKF-CNRD噪声估计均方差效果图；横坐标为时间，纵坐标为噪声估计均方差；
图6为INS位置误差及KF算法、CKF-cnrd算法滤波误差对比图；
图7为INS速度误差及KF算法，CKF-cnrd算法滤波误差对比图；
图8为组合导航系统KF算法、CKF-cnrd算法位置误差对比图；
图9为组合导航系统KF算法、CKF-cnrd算法速度误差对比图。
具体实施方式
具体实施方式一：本实施方式的非理想条件下高斯滤波替代框架组合导航方法，是基于非线性离散系统模型的，非线性离散系统模型的表达式如下：
xk＝f(xk-1)+ωk-1 (1)
zk＝h(xk)+νk (2)
这里xk和zk分别为k时刻的n维状态和m维观测向量，f(·)和h(·)分别为非线性函数，过程噪声ωk-1和量测噪声νk分别为k-1时刻和k时刻的满足高斯分布的白噪声，其协方差阵满足
Figure BDA0001236972260000041
这里δkl为Kronecker delta函数；Sk≠0是过程噪声和量测噪声的协方差矩阵，并且这里假设初始状态x0是高斯的，其均值为
Figure BDA0001236972260000042
并且是独立于ωk和νk。
假设实际量测值zk在k＝1时的值总是可以获得的。然而，正如前面所讨论的，传感器的测量可能是随机延迟的。换句话说，实际测量输出yk可能是当前输出zk或者前面任意一个时刻的值，基于合理的假设，针对一步随机时滞的情况，我们以满足Bernoulli分布的随机序列来构造随机时滞量测数据，实际量测函数可以被表示为
Figure BDA0001236972260000043
这里，{γk；k＞1}表示一个独立的Bernoulli随机变量序列，取值为
p(γk＝1)＝E[γk]＝pk (5)
我们的目标是对系统(1)、(2)、(4)在公式(3)、(5)及前面讨论的基础上设计替代框架算法，也就在最小均方误差意义下，利用前k时刻已知的量，通过计算状态的前两阶矩，得到一步预测和滤波估计后的状态值，即
Figure BDA0001236972260000051
Figure BDA0001236972260000052
注2.1:上标″∧″和″～″分别表示估计值及真值与估计值的差，″k+1|k″及″k+1|k+1″分别表示由k时刻预测k+1时刻及由k+1时刻修正k+1时刻，文章后面的部分均采用此定义，不再说明。
注2.2:由于本文考虑一步随机时滞的情况，预测值和修正值均是基于实际量测值Yk的，这里Yk＝{y1,y2,…,yk}，而并非理论量测值zk，不再赘述。
基于文章[18]，我们给出如下两个假设条件，这也是为了实现本发明需要的必要条件：假设2.1.状态xk+1关于条件Yk的一步后验概率密度函数是高斯的，即
Figure BDA0001236972260000053
假设2.2.关于实际测量值yk+1在前k步测量值Yk已知的条件下的一步后验预测概率密度函数依旧是高斯的，即
Figure BDA0001236972260000054
实际量测值的一步预测均值和协方差分别为
上述两个假设条件可以总结为：k+1时刻的状态关于前k步量测值的一步后验概率密度函数是高斯的；k+1时刻的实际量测值关于前k步量测值的一步后验预测概率密度函数是高斯的。
基于上面的假设，再将带有一步随机测量时滞和相关噪声的非线性系统替代框架，即：
将方程(2)代入(4)得到
将方程(1)和(11)代入公式(6)和(10)中一步预测均值的定义可得
Figure BDA0001236972260000061
显然，这里存在和
Figure BDA0001236972260000063
即过程噪声和量测噪声相对于实际量测值的期望，由于Yk中含有噪声项υk且本文中过程噪声和量测噪声相关，即Yk与噪声ωk及υk相关，故其期望均不为0，因此在一步随机量测时滞存在的情况下，不能利用传统的计算方法来直接得到状态和量测一步预测的估计值。本文采用替代框架的处理思想，将ωk的均值0和协方差Qk以及υk的均值0和协方差Rk看做ωk和υk的一个先验信息，然后通过实际观测值yk进行更新，我们假设更新后的均值和方差分别为
Figure BDA0001236972260000064
及在滤波递推过程中，我们利用噪声的后验值而不是先验值进行计算。为了使得证明看起来更清晰，我们将证明过程中用到的前后估计的量都列写出来，以便形成一个完整的闭环。下面我们给出具体的计算步骤，分为两步来完成，即状态预测和状态修正。
本发明的方法包括如下步骤：
步骤一：根据飞行器内的传感器测量值，获取飞行器在k时刻的状态值xk；
步骤二(即状态预测)：建立状态在k+1时刻的预测均值
Figure BDA0001236972260000066
和协方差
Figure BDA0001236972260000067
表达式：
Figure BDA0001236972260000068
Figure BDA0001236972260000069
这里
Figure BDA00012369722600000610
Figure BDA00012369722600000611
其中，
Figure BDA00012369722600000612
为k时刻对k+1时刻的预测估计值；是k+1时刻的n维状态向量的协方差矩阵；Qk+1为k+1时刻的过程噪声ωk+1的协方差矩阵；Rk+1为k+1时刻的量测噪声νk+1的协方差矩阵；Sk+1为k+1时刻的过程噪声和量测噪声的协方差矩阵，且Sk+1≠0。
步骤三(即状态修正)：建立状态在k+1时刻的后验均值
Figure BDA0001236972260000071
和协方差
Figure BDA0001236972260000072
表达式：
Figure BDA0001236972260000073
Figure BDA0001236972260000074
其中，
Figure BDA0001236972260000075
为增广后的状态向量的新息，其中
Figure BDA0001236972260000077
为k+1时刻的新息nk+1与k+1时刻的预测均值
Figure BDA0001236972260000078
的互协方差。
步骤四：基于CKF滤波算法对步骤二和步骤三得到的表达式进行求解，得到修正后的均值和方差。
步骤五：根据所述修正后的均值和方差修正航向，实现导航。
具体实施方式二：本实施方式与具体实施方式一不同的是：
步骤二中建立状态在k+1时刻的预测均值
Figure BDA0001236972260000079
和协方差表达式的具体过程为：
令带有噪声的增广状态向量为
Figure BDA00012369722600000711
则k+1时刻状态变量x的一步预测均值和协方差阵分别为
Figure BDA00012369722600000712
由于ωk与xk中所含的噪声ωk-1不相关，故代入式(20)可得
Figure BDA00012369722600000715
这里
Figure BDA00012369722600000716
和均为已知信息，通过简单的数学推导,可以得到一步预测均值和方差
Figure BDA0001236972260000081
Figure BDA0001236972260000082
由于ωk+1和υk+1均与Yk不相关，显然有
p(ωk+1|Yk)＝p(ωk+1)＝N(ωk+1；0,Qk+1) (24)
p(υk+1|Yk)＝p(υk+1)＝N(υk+1；0,Rk+1) (25)
故由公式(11)、(24)、(25)，根据高斯分布的计算规则[14]，
Figure BDA0001236972260000083
在Yk已知的条件下，其与yk+1联和后验概率密度函数也是高斯的，即
Figure BDA0001236972260000084
由贝叶斯规则可知
Figure BDA0001236972260000085
重新整理公式(26)，可以得到
Figure BDA0001236972260000086
这里
Figure BDA0001236972260000087
进一步，p(lk+1|Yk+1)也是高斯的，可以依据公式(30)进行更新。
结合本文的假设及文献[18]，联和噪声的后验概率密度函数也是高斯的，满足：
Figure BDA0001236972260000091
以上即完成了对步骤二中公式(14)和(15)的证明。
其它步骤及参数与具体实施方式一相同。
具体实施方式三：本实施方式与具体实施方式一或二不同的是：
步骤三中，建立状态在k+1时刻的后验均值
Figure BDA0001236972260000092
和协方差
Figure BDA0001236972260000093
表达式的具体过程为：
考虑非线性系统(1)和(2)，并且假设k时刻的一步预测值已由步骤二给出，则k+1时刻的后验均值和协方差由下式给出：
Figure BDA0001236972260000095
这里
Figure BDA0001236972260000096
定义如下：
Figure BDA0001236972260000097
同步骤二中的讨论，可得
这里
Figure BDA0001236972260000099
和
Figure BDA00012369722600000910
分别由公式(43)和(44)给定。
关于公式(43)和(44)中的各个变量的定义，由前面的定义，实际观测的一步预测均值为
Figure BDA00012369722600000911
由新息的定义及实际观测值yk+1和其一步预测均值可得
Figure BDA00012369722600000912
状态和量测的一步预测协方差为
Figure BDA0001236972260000101
Figure BDA0001236972260000102
其中
Figure BDA0001236972260000104
Figure BDA0001236972260000107
Figure BDA0001236972260000108
其它步骤及参数与具体实施方式一或二相同。
具体实施方式四：本实施方式与具体实施方式一至三之一不同的是：
步骤四中，为了便于求解(例如使用matlab进行求解)，还需要一些步骤将步骤二和步骤三中的公式转化为可求解的形式，具体过程为：
①一步预测：
假设k+1时刻之前的所有信息已知，我们可以依据下列的数值格式进行迭代：
1)矩阵分解
Figure BDA0001236972260000111
2)容积点计算
Figure BDA0001236972260000112
这里ξi为容积点.
3)容积点传播
Xi,k+1|k＝f(Xi,k|k),i＝1,2,…,2n (47)
一步预测均值和协方差
Figure BDA0001236972260000113
Figure BDA0001236972260000114
②量测更新：
1)矩阵分解
Figure BDA0001236972260000116
2)容积点计算
Figure BDA0001236972260000117
Figure BDA0001236972260000118
这里n′＝2n+m，
Figure BDA0001236972260000119
和
Figure BDA00012369722600001110
分别表示状态，过程噪声和量测噪声分量。3)容积点传播
θi,k+1|k＝h(Xi,k+1|k),i＝1,2,…,2n (54)
Figure BDA00012369722600001111
可得
Figure BDA00012369722600001112
Figure BDA0001236972260000121
Figure BDA0001236972260000123
Figure BDA0001236972260000124
将式(56)-(60)代入(30)-(32),将所得公式代入修正表示式中，可以得到修正后的均值和方差。
具体实施方式五：本实施方式与具体实施方式一至四之一不同的是：
在步骤四执行完毕后，还包括数值仿真步骤，包括：
步骤A1：建立单变量非静态增长模型(UNGM)，来测试本发明的有效性，表达式为：
Figure BDA0001236972260000127
Figure BDA0001236972260000128
yk＝(1-γk)zk+γkzk-1
其中ωk和υk的协方差分别为Qk＝10和Rk＝3，并且互协方差为Sk＝1。
步骤A2：对于所述单变量非静态增长模型，执行100次独立的蒙特卡罗仿真以计算均方根误差RMSE；其中，k时刻的均方根误差RMSE为：
Figure BDA0001236972260000129
Figure BDA0001236972260000131
和
Figure BDA0001236972260000132
分别表示k时刻蒙特卡罗执行的真值和估计值。
步骤A3：根据步骤A2的计算结果得到统计结果。
图2和3表示EKF、CKF、UKF and CKF-CNRD之间的对比。一方面，因为EKF利用了一阶泰勒展开，其处理强非线性问题的能力很差，因此，我们可以看到EKF的估计值相当的差，相比来说，CKF和UKF的结果要更好一些。另一方面，由于测量值是随机延迟的，而EKF、CKF和UKF不具有处理延迟的能力，因此仿真中出现了许多尖点。从仿真图我们可以看到我们提出的CKF-CNRD算法完美解决了这两个问题。图4和5展示了ωk和υk的最小均方误差估计。我们可以看到噪声的均值不再为0，因为8％的过程噪声超出了范围[-2,2]并且8.5％的量测噪声超出了范围[-0.8,0.8]，这意味着我们做这样的一个噪声估计是很有必要的。从仿真中可以看出，噪声估计值的尖点出现的时刻总是相同的，这与推导过程中的表达式是相符的，因为噪声的估计值总是相差一个倍数。
其它步骤及参数与具体实施方式一至四之一相同。
本发明还可有其它多种实施例，在不背离本发明精神及其实质的情况下，本领域技术人员当可根据本发明作出各种相应的改变和变形，但这些相应的改变和变形都应属于本发明所附的权利要求的保护范围。
具体实施方式六：本实施方式与具体实施方式一至五之一不同的是：
在步骤四执行完毕后，还可以执行应用仿真步骤，具体如下：
利用伪距、伪距率的紧耦合表达方式对GPS/INS组合系统进行仿真，其中系统状态方程和量测方程如文献[19]所示，状态方程中的系数矩阵如文献[20]所示，系统方程中噪声矩阵惯导部分由陀螺仪马氏过程，陀螺仪白噪声和加速度计马氏过程组成，取值分别为(0.04/(57*3600))^2,(0.01/(57*3600))^2,(1e-3)^2，假设系统是一步随机时滞的，时滞概率取为0.3，利用Kalman滤波及本文提出的算法分别进行500s仿真，得到惯导系统及组合导航系统的位置误差和速度误差分别如图6至图9所示：
从仿真图6、7可以看出，单独的惯性导航系统存在一个误差越来越大的趋势，两种滤波算法均可减小误差且可得到较为精确的系统状态估计值。从图8、9可以看出，本文所提出的算法在KF算法性能较差的时候依然能够得到较为精确的目标位置和速度，反映出了本文提出算法在系统存在时滞情况的有效性。
综上，针对GPS/INS组合导航系统可能存在噪声相关和随机量测时滞的情况，本发明提出一种替代框架下的高斯滤波算法，不同于以往噪声和状态分开处理的情况，本发明将过程噪声和量测噪声一起作为状态的一个扩张量进行滤波更新，将已知的噪声信息作为先验信息，在状态和量测的更新方程中用更新后的噪声的后验估计代替先验信息，仿真证明了其估计精度优于现有的滤波算法。同时，本发明提出的方法更具有一般性，选取合适的参数可以看出，一般高斯滤波算法是其在噪声不相关或不存在随机时滞情况的特例。
<参考文献>
[1]黄晓瑞,崔平远,崔祜涛,GPS/INS组合导航系统自适应滤波算法与仿真研究，飞行力学19(2),2001,6,69-77.
[2]Seong Yun Cho,Byung Doo Kim,Adaptive IIR/FIR fusion filter and itsapplication to the INS/GPS integrated system,Automatica 44(8)(2008)2040–2047.
[3]A.Noureldin,A.El-Shafie,M.Bayoumi,GPS/INS integration utilizingdynamic neural networks for vehicular navigation,Inf.Fusion 12(1)(2011)48–57.
[4]C.G.Park,S.Y.Cho,Y.Jin,H.W.Park,G.I.Jee,J.G.Lee,Analysis ofmeasurement delay in INS/GPS integration with Kalman filtering,in:Proceedingsof the 3rd ASCC,Shanghai,China,Jul.4–7,2000.
[5]Yaakov B S,Li X,Thiagalingam K.Estimation with Applications toTracking and Navigation[M].New York:John Wiley and Sons,2001.
[6]Julier S J,Uhlmann J K.Unscented Filtering and NonlinearEstimation[J].Proceedings of the IEEE.2004,92(3):401-422.
[7]Arasaratnam I,Haykin S.Cubature Kalman Filters[J].IEEETransactions on Automatic Control.2009,54(6):1254-1269.
[8]Ito,K.,Xiong,K.:‘Gaussian filters for nonlinear filteringproblems’,IEEE Trans.Autom.Control,2000,45,(5),pp.910–927
[9]Nrgaard M,Poulsen N K,Ravn O.New Developments in State Estimationfor Nonlinear Systems[J].Automatica.2000,36(11):1627-1638.
[10]Simon,D.:‘Optimal state estimation:Kalman,H∞,and nonlinearapproaches’(Wiley-Interscience,2006)
[11]Feng,X.,Ge,Q.,Wen,C.:‘An optimal sequential filter for the linearsystem with correlated noises’.Control and Decision Conf.,CCDC’09,2009
[12]Wang,X.,Liang,Y.,Pan,Q.,&Yang,F.(2012).A Gaussian approximationrecursive filter for nonlinear systems with correlated noises.Automatica,48,2290–2297.
[13]Chang,G.:‘Comment on“a Gaussian Approximation Recursive Filterfor Nonlinear Systems with Correlated Noises”’,Automatica,2014,50,(2),pp.655–656
[14]Chang,G.:‘Alternative formulation of the Kalman filter forcorrelated process and observation noise’,IET Sci.Meas.Technol.,2014,8,(5),pp.310–318
[15]Gopalakrishnan A,Kaisare N S,Narasimhan S.Incorporating Delayedand Infrequent Measurements in Extended Kalman Filter Based Nonlinear StateEstimation[J].Journal of Process Control.2011,21(1):119-129.
[16]Hermoso-Carazo A,Linares-Perez J.Extended and Unscented FilteringAlgorithms Using One-Step Randomly Delayed Observations[J].AppliedMathematics and Computation.2007,190(2):1375-1393.
[17]Wang X,Liang Y,Pan Q,et al.Design and Implementation of GaussianFilter for Nonlinear System with Randomly Delayed Measurements and CorrelatedNoises[J].Applied Mathematics and Computation.2014,232:1011-1024.
[18]X.X.Wang,Y.Liang,Q.Pan,Zhao Chunhui,Gaussian filter for nonlinearsystems with one-step randomly delayed measure-ments,Automatica 49(4)(2013)976–986.
[19]Wang W,Liu Z,Xie R.Quadratic extended Kalman filter approach forGPS/INS integration[J].Aerospace science and technology,2006,10(8):709-713.
[20]袁信,俞济祥,陈哲.导航系统[M].北京航空工业出版社，1992。

Claims (4)

1.一种非理想条件下高斯滤波替代框架组合导航方法，所述方法基于非线性离散系统模型，并且满足如下条件：
k+1时刻的状态关于前k步量测值的一步后验概率密度函数是高斯的；
k+1时刻的实际量测值关于前k步量测值的一步后验预测概率密度函数是高斯的；
其特征在于，所述方法包括如下步骤：
步骤一：根据飞行器内的传感器测量值，获取飞行器在k时刻的状态值xk；
步骤二：根据所述状态值xk，建立状态在k+1时刻的预测均值
Figure FDA0002212528460000011
和协方差
Figure FDA0002212528460000012
表达式：
Figure FDA0002212528460000013
其中，
Figure FDA0002212528460000015
为k时刻对k+1时刻的预测估计值；
Figure FDA0002212528460000016
是k+1时刻的n维预测状态向量的协方差矩阵；Qk+1为k+1时刻的过程噪声ωk+1的协方差矩阵；Rk+1为k+1时刻的量测噪声νk+1的协方差矩阵；Sk+1为k+1时刻的过程噪声和量测噪声的协方差矩阵，且Sk+1≠0；
步骤三：建立状态在k+1时刻的后验均值和协方差
Figure FDA0002212528460000018
表达式：
Figure FDA0002212528460000019
其中，
Figure FDA00022125284600000111
为增广后的状态向量的新息，
Figure FDA00022125284600000112
其中为k+1时刻的新息nk+1与k+1时刻的预测均值
Figure FDA00022125284600000114
的互协方差；
步骤四：基于CKF滤波算法对步骤二和步骤三得到的表达式进行求解，得到修正后的均值和方差；
步骤五：根据所述修正后的均值和方差修正航向，实现导航。
2.根据权利要求1所述的方法，其特征在于，所述步骤二中建立状态在k+1时刻的预测均值
Figure FDA0002212528460000021
和协方差
Figure FDA0002212528460000022
表达式的具体过程为：
令带有噪声的增广状态向量为
则k+1时刻状态变量x的一步预测均值
Figure FDA0002212528460000024
和协方差矩阵
Figure FDA0002212528460000025
分别为
Figure FDA0002212528460000026
Figure FDA0002212528460000027
其中f(xk)为任意以xk为自变量的非线性函数；Yk为前k步量测值；
Figure FDA0002212528460000028
为由k时刻修正k时刻的过程噪声的估计值；
由于ωk与xk中所含的噪声ωk-1不相关，故ωkf(xk)T＝0，代入协方差矩阵
Figure FDA0002212528460000029
中，可得：
Figure FDA00022125284600000210
其中
Figure FDA00022125284600000211
为
Figure FDA00022125284600000212
的协方差矩阵，通过计算可得一步预测均值和协方差
Figure FDA00022125284600000213
Figure FDA00022125284600000214
由于ωk+1和υk+1均与Yk不相关，则有
p(ωk+1|Yk)＝p(ωk+1)＝N(ωk+1；0,Qk+1)
p(υk+1|Yk)＝p(υk+1)＝N(υk+1；0,Rk+1)
联合噪声
Figure FDA0002212528460000031
在Yk已知的条件下，其与yk+1联和后验概率密度函数也是高斯的，即
Figure FDA0002212528460000032
其中yk+1为k+1时刻的实际量测输出；
由贝叶斯规则可知
Figure FDA0002212528460000033
则有
Figure FDA0002212528460000034
其中，为k时刻对k+1时刻预测的实际量测输出估计值
Figure FDA0002212528460000036
的协方差矩阵，
Figure FDA0002212528460000037
为k时刻对k+1时刻的联合噪声的预测估计值，
Figure FDA0002212528460000038
为
Figure FDA0002212528460000039
与
Figure FDA00022125284600000310
的互协方差矩阵，
Figure FDA00022125284600000311
且p(lk+1|Yk+1)也是符合高斯分布的，则联合噪声的后验概率密度也是符合高斯分布的，满足：
Figure FDA00022125284600000312
其中，
Figure FDA00022125284600000313
为由k+1时刻修正k+1时刻联合噪声的估计值，
Figure FDA00022125284600000314
为的协方差矩阵；状态在k+1时刻的后验均值
Figure FDA00022125284600000316
和协方差
Figure FDA00022125284600000317
表达式为：
Figure FDA00022125284600000318
Figure FDA00022125284600000319
3.根据权利要求2所述的方法，其特征在于，所述步骤三中，建立状态在k+1时刻的后验均值
Figure FDA0002212528460000041
和协方差
Figure FDA0002212528460000042
表达式的具体过程为：
实际观测的一步预测均值为
Figure FDA0002212528460000043
由新息的定义及实际观测值yk+1和其一步预测均值可得：
γk+1为符合伯努利随机分布的序列的取值；pk+1＝p(γk+1＝1)；状态和量测的一步预测协方差为
Figure FDA0002212528460000045
Figure FDA0002212528460000046
4.根据权利要求1所述的方法，其特征在于，在所述步骤四执行完毕后，还包括数值仿真步骤，包括：
步骤A1：建立单变量非静态增长模型，表达式为：
Figure FDA0002212528460000048
yk＝(1-γk)zk+γkzk-1
其中ωk和υk的协方差分别为Qk＝10和Rk＝3，并且互协方差为Sk＝1；
步骤A2：对于所述单变量非静态增长模型，执行100次独立的蒙特卡罗仿真以计算均方根误差RMSE；其中，k时刻的均方根误差RMSE为：
Figure FDA0002212528460000049
和
Figure FDA00022125284600000411
分别表示k时刻蒙特卡罗执行的真值和估计值；
步骤A3：根据步骤A2的计算结果得到统计结果。
CN201710120792.XA 2017-03-02 2017-03-02 一种非理想条件下高斯滤波替代框架组合导航方法 Active CN106871905B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710120792.XA CN106871905B (zh) 2017-03-02 2017-03-02 一种非理想条件下高斯滤波替代框架组合导航方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710120792.XA CN106871905B (zh) 2017-03-02 2017-03-02 一种非理想条件下高斯滤波替代框架组合导航方法

Publications (2)

Publication Number Publication Date
CN106871905A CN106871905A (zh) 2017-06-20
CN106871905B true CN106871905B (zh) 2020-02-11

Family

ID=59169390

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710120792.XA Active CN106871905B (zh) 2017-03-02 2017-03-02 一种非理想条件下高斯滤波替代框架组合导航方法

Country Status (1)

Country Link
CN (1) CN106871905B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107765347B (zh) * 2017-06-29 2020-06-16 河海大学 一种高斯过程回归和粒子滤波的短期风速预测方法
CN109521454B (zh) * 2018-12-06 2020-06-09 中北大学 一种基于自学习容积卡尔曼滤波的gps/ins组合导航方法
CN111193528B (zh) * 2019-12-30 2021-09-17 哈尔滨工业大学 基于非理想条件下非线性网络系统的高斯滤波方法

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104121907A (zh) * 2014-07-30 2014-10-29 杭州电子科技大学 一种基于平方根容积卡尔曼滤波器的飞行器姿态估计方法
CN105737823A (zh) * 2016-02-01 2016-07-06 东南大学 基于五阶ckf的gps/sins/cns组合导航方法
CN105973238A (zh) * 2016-05-09 2016-09-28 郑州轻工业学院 一种基于范数约束容积卡尔曼滤波的飞行器姿态估计方法
CN106291645A (zh) * 2016-07-19 2017-01-04 东南大学 适于高维gnss/ins深耦合的容积卡尔曼滤波方法
CN106323291A (zh) * 2016-09-29 2017-01-11 安徽工程大学 一种基于平方根容积卡尔曼滤波的多机器人协同定位算法
CN106352876A (zh) * 2016-07-25 2017-01-25 北京航空航天大学 一种基于h∞和ckf混合滤波的机载分布式pos传递对准方法

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104121907A (zh) * 2014-07-30 2014-10-29 杭州电子科技大学 一种基于平方根容积卡尔曼滤波器的飞行器姿态估计方法
CN105737823A (zh) * 2016-02-01 2016-07-06 东南大学 基于五阶ckf的gps/sins/cns组合导航方法
CN105973238A (zh) * 2016-05-09 2016-09-28 郑州轻工业学院 一种基于范数约束容积卡尔曼滤波的飞行器姿态估计方法
CN106291645A (zh) * 2016-07-19 2017-01-04 东南大学 适于高维gnss/ins深耦合的容积卡尔曼滤波方法
CN106352876A (zh) * 2016-07-25 2017-01-25 北京航空航天大学 一种基于h∞和ckf混合滤波的机载分布式pos传递对准方法
CN106323291A (zh) * 2016-09-29 2017-01-11 安徽工程大学 一种基于平方根容积卡尔曼滤波的多机器人协同定位算法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
Alternative framework of the Gaussian filter for non-linear systems with synchronously correlated noises;Yu Han,etc;《IET Science, Measurement & Technology》;20160630;第10卷(第4期);第306-315页 *
具有随机时滞和异步相关噪声的非线性系统的高斯滤波器设计;于浛，等;《中国惯性技术学报》;20150430;第23卷(第2期);第238-247页 *
考虑随机量测时滞和同步相关噪声的改进高斯滤波算法;于浛，等;《控制理论与应用》;20160229;第33卷(第2期);第133-145页 *

Also Published As

Publication number Publication date
CN106871905A (zh) 2017-06-20

Similar Documents

Publication Publication Date Title
CN107255795B (zh) 基于ekf/efir混合滤波的室内移动机器人定位方法和装置
WO2018014602A1 (zh) 适于高维gnss/ins深耦合的容积卡尔曼滤波方法
CN108255791B (zh) 基于分布式传感器一致性的机动目标跟踪方法
CN104020480B (zh) 一种带自适应因子的交互式多模型ukf的卫星导航方法
CN106871905B (zh) 一种非理想条件下高斯滤波替代框架组合导航方法
CN104392136A (zh) 一种面向高动态非高斯模型鲁棒测量的高精度数据融合方法
CN111965618B (zh) 一种融合多普勒量测的转换量测跟踪方法及系统
Xu et al. A multi-sensor information fusion method based on factor graph for integrated navigation system
CN111193528B (zh) 基于非理想条件下非线性网络系统的高斯滤波方法
CN111964667A (zh) 一种基于粒子滤波算法的地磁_ins组合导航方法
CN113791240B (zh) 基于高阶滑膜跟踪微分器的加速度估计方法、系统、设备和介质
CN114139109A (zh) 一种目标跟踪方法、系统、设备、介质及数据处理终端
Shin et al. A new fusion formula and its application to continuous-time linear systems with multisensor environment
CN116734864B (zh) 一种常值观测偏差条件下航天器自主相对导航方法
KR101502721B1 (ko) 적응형 상호작용 다중모델 추정기를 이용한 정밀 위치정보 제공 방법 및 장치
Jwo et al. GPS navigation processing using the quaternion-based divided difference filter
CN111578931B (zh) 基于在线滚动时域估计的高动态飞行器自主姿态估计方法
CN111709438B (zh) 一种异构传感器信息融合方法
CN110989341B (zh) 一种约束辅助粒子滤波方法及目标跟踪方法
CN109115228B (zh) 一种基于加权最小二乘容积卡尔曼滤波的目标定位方法
CN114993341B (zh) 一种基于天基测量的运载火箭弹道估计方法及装置
CN115859626A (zh) 针对周期运动目标的自适应无迹卡尔曼滤波器设计方法
CN107966697B (zh) 一种基于渐进无迹卡尔曼的移动目标跟踪方法
CN110912535A (zh) 一种新型无先导卡尔曼滤波方法
Dubois et al. Performance evaluation of a moving horizon estimator for multi-rate sensor fusion with time-delayed measurements

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
CB03 Change of inventor or designer information
CB03 Change of inventor or designer information

Inventor after: Song Shenmin

Inventor after: Tan Liguo

Inventor after: Zhao Kai

Inventor after: Zhang Xiujie

Inventor after: Wu Xiaohang

Inventor after: Si Yiwen

Inventor before: Song Shenmin

Inventor before: Zhao Kai

Inventor before: Zhang Xiujie

Inventor before: Wu Xiaohang

Inventor before: Si Yiwen

GR01 Patent grant
GR01 Patent grant