一种适用于弹道导弹ins/cns/gnss组合导航系统的最优数据融合方法
阅读说明:本技术 一种适用于弹道导弹ins/cns/gnss组合导航系统的最优数据融合方法 (Optimal data fusion method suitable for ballistic missile INS/CNS/GNSS combined navigation system ) 是由 陈熙源 柳笛 刘晓 石春凤 马振 于 2019-09-24 设计创作,主要内容包括:本发明公开了一种适用于弹道导弹INS/CNS/GNSS组合导航系统的最优数据融合方法,该方法包括以下步骤:构造INS/CNS/GNSS组合导航系统模型;在广义高阶CKF的时间更新阶段和量测更新阶段分别引入自适应渐消因子和最大相关熵准则进行INS/CNS子系统和INS/GNSS子系统的局部状态估计;根据最小方差原理和容积准则融合INS/CNS子系统和INS/GNSS子系统的局部估计得到全局最优状态估计。本发明可以同时抑制过程建模误差和非高斯量测噪声对状态估计的影响,提高弹道导弹INS/CNS/GNSS组合导航的自适应性和鲁棒性,获得全局最优的状态估计值。(The invention discloses an optimal data fusion method suitable for a ballistic missile INS/CNS/GNSS combined navigation system, which comprises the following steps: constructing an INS/CNS/GNSS integrated navigation system model; respectively introducing a self-adaptive fading factor and a maximum correlation entropy criterion to carry out local state estimation of the INS/CNS subsystem and the INS/GNSS subsystem in a time updating stage and a measurement updating stage of the generalized high-order CKF; and according to the minimum variance principle and the volume criterion, local estimation of the INS/CNS subsystem and the INS/GNSS subsystem is fused to obtain global optimal state estimation. The method can simultaneously inhibit the influence of process modeling errors and non-Gaussian measurement noise on state estimation, improve the adaptivity and robustness of ballistic missile INS/CNS/GNSS combined navigation, and obtain the globally optimal state estimation value.)
技术领域
本发明涉及一种数据融合方法,具体涉及一种适用于弹道导弹INS/CNS/GNSS组合导航系统的最优数据融合方法。
背景技术
传统的基于EKF、UKF的联邦多传感器数据融合方法使用过程噪声协方差的上界而不是过程噪声协方差本身来消除局部状态估计之间的相关性,因此得到的全局状态估是次优的。此外,该类方法要求非线性系统模型必须达到足够的精度。然而,在弹道导弹飞行的过程中,一方面,其弹载的INS/CNS/GNSS组合导航系统很容易受到复杂环境的影响,进而导致INS/CNS/GNSS组合导航系统的量测信息受到非高斯噪声的影响,在这种情况下如果继续使用传统的基于高斯假设的非线性滤波方法将会使组合导航系统的导航精度严重下降;另一方面,由于真实的动态系统模型非常复杂,我们建立的系统过程模型只能是真实模型的理论近似,因而存在建模误差,这种建模误差也会导致组合导航系统的导航精度的下降。
所以,需要一个新的技术方案来解决上述问题。
发明内容
发明目的:为了克服现有技术中存在的组合导航系统的导航精度差的不足,提供一种适用于弹道导弹INS/CNS/GNSS组合导航系统的最优数据融合方法,其实现在过程模型含有不确定性以及量测噪声为非高斯噪声的情况下依然可以获取全局最优的状态估计。
技术方案:本发明提供一种适用于弹道导弹INS/CNS/GNSS组合导航系统的最优数据融合方法,包括以下步骤:
S1:构造INS/CNS/GNSS组合导航系统滤波模型;
S2:在广义高阶CKF的时间更新阶段和量测更新阶段分别引入自适应渐消因子和最大相关熵准则进行INS/CNS子系统和INS/GNSS子系统的局部状态估计;
S3:根据最小方差原理和容积准则融合INS/CNS子系统和INS/GNSS子系统的局部估计得到全局最优状态估计。
进一步的,所述步骤S1中构造INS/CNS/GNSS组合导航系统滤波模型包括以下步骤:
S1-1:设置INS/CNS/GNSS组合导航系统的状态向量为x=[φx,φy,φz,δvx,δvy,δvz,x,y,z,εx,εy,εz,Δx,Δy,Δz]T,其中φx,φy,φz表示发射点惯性系下姿态失准角,δvx,δvy,δvz表示发射点惯性系下的速度误差,x,y,z表示发射点惯性系下的位置误差,εx,εy,εz表示弹体坐标系下的陀螺常值漂移,Δx,Δy,Δz表示弹体坐标系下的加速度计常值偏置,T为转置符号;
S1-2:根据INS/CNS/GNSS组合导航系统的15维状态向量x建立系统的状态方程:
xk=f(xk-1)+vk-1
其中,f(·)为非线性系统函数,xk-1和xk分别表示k-1时刻和k时刻的状态向量,vk-1表示过程噪声,vk-1的协方差为
S1-3:分别建立INS/GNSS子系统和INS/CNS子系统的量测方程:
在INS/GNSS子系统中,将INS与GNSS输出的位置的差值和速度的差值作为量测信息,建立该子系统的量测方程:
z1,k=H1,kxk+ω1,k
其中,z1,k表示INS/GNSS子系统k时刻的量测向量,H1,k表示INS/GNSS子系统k时刻的测量矩阵,ω1,k表示INS/GNSS子系统k时刻的量测噪声,其方差是
在INS/CNS子系统中,将INS与CNS输出的姿态角的差值作为量测信息,建立该子系统的量测方程:
z2,k=H2,kxk+ω2,k
其中,z2,k表示INS/CNS子系统k时刻的量测向量,H2,k表示INS/CNS子系统k时刻的测量矩阵,ω2,k表示INS/CNS子系统k时刻的量测噪声,其方差为
进一步的,所述步骤S2具体包括以下步骤:
S2-1:因为INS/CNS子系统和INS/GNSS子系统采用同样的滤波过程进行局部状态估计,为了避免重复,在此仅具体说明INS/GNSS子系统的滤波过程。因此,将步骤S1-3中的z1,k、z2,k统一写成zk,H1,k和H2,k统一写成Hk,ω1,k和ω2,k统一写成ωk,R1,k和R2,k统一写成Rk,根据INS/GNSS子系统的局部状态估计要求设置核宽度γ的值,设置初始的状态向量、状态误差协方差、渐消因子分别为
τk(0)=1,S0|0=chol(P0|0),chol(·)表示乔列斯分解运算;S2-2:根据公式
计算出传播容积点其中i=1,...,2n2+1,n表示状态向量的维数,[·]i表示集合[·]的第i列;S2-3:预测k时刻状态向量
误差协方差阵和计算Sk|k-1=chol(Pk|k-1),其中S2-4:根据
和Sk|k-1产生新的容积点并预测k时刻的量测信息S2-5:根据INS/GNSS子系统的量测方程和步骤S2-3中的
和Pk|k-1构建如下方程:
其中,I表示单位矩阵,
Mp,k|k-1=chol(Pk|k-1)、Mr,k=chol(Rk)和在方程
的两边左乘得到:令则上述方程可重写为:Dk=BkXk+ek;S2-6:对量测噪声方差进行更新:
其中,diag(·)表示矩阵的对角化运算,m为量测信息的维度, di,k表示Dk的第i个元素,bi,k表示Bk的第i行元素;S2-7:计算状态信息和量测信息之间的互协方差矩阵Pxz,k|k-1,并令i=0;
S2-8:令
并设定卡方分布的值,其中θ表示卡方分布的自由度,α表示卡方分布的分位数,如果则计算 返回步骤S2-2继续执行下一滤波周期;S2-9:
如果
计算其中,i=(0,1,...),i=i+1,根据得到自适应渐消因子τk更新ak(i)的值,判断ak(i)是否满足若不满足,则返回步骤S2-9继续执行;若满足,则得到最优的自适应渐消因子τk的值,然后继续执行后续步骤;S2-10:使用自适应渐消因子τk更新预测的状态预测协方差:
计算卡尔曼增益:估计后验状态和后验协方差矩阵:返回步骤S2-2继续执行下一滤波周期。进一步的,所述步骤S3包括以下步骤:
S3-1:在分别对INS/GNSS子系统和INS/CNS子系统执行完滤波过程后便可获得相应的局部后验状态估计
为了便于区分,令INS/GNSS子系统获得局部后验状态估计为INS/CNS子系统获得局部后验状态估计为53-2:根据最小方差原理和容积准则得到全局最优状态估计:
其中令β=[β1,β2]T,P11为INS/GNSS子系统执行局部状态估计得到的Pk|k,P22为INS/CNS子系统执行局部状态估计得到的Pk|k,P12表示和估计误差的互协方差矩阵,P21表示和估计误差的互协方差矩阵。进一步的,P12和P21通过容积准则近似得到:
其中,
分别表示INS/GNSS子系统与INS/CNS子系统执行步骤S2-2得到的传播容积点和分别表示INS/GNSS子系统与INS/CNS子系统执行步骤S2-3得到的状态向量预测值,和分别表示INS/GNSS子系统与INS/CNS子系统执行步骤S2-4得到的新的容积点Xi,k|k-1,和分别表示INS/GNSS子系统与INS/CNS子系统执行步骤S2-4得到的量测信息预测值 和分别表示INS/GNSS子系统与INS/CNS子系统执行局部状态估计过程得到的滤波增益 I表示n维的单位矩阵。有益效果:本发明与现有技术相比,在广义高阶CKF的时间更新阶段和量测更新阶段分别引入自适应渐消因子和最大相关熵准则进行INS/GNSS子系统和INS/CNS子系统的局部状态估计,最后根据最小方差原理和容积准则融合INS/CNS子系统和INS/GNSS子系统的局部估计,使得弹道导弹INS/CNS/GNSS组合导航系统在系统模型存在过程建模误差和量测噪声为非高斯分布的情况下,该组合导航系统依然可以获得全局最优的状态估计,本发明可以同时抑制过程建模误差和非高斯量测噪声对状态估计的影响,提高弹道导弹INS/CNS/GNSS组合导航的自适应性和鲁棒性,从而保证了INS/CNS/GNSS组合导航系统的导航精度。
附图说明
图1为本发明方法的流程图。
具体实施方式
下面结合附图和具体实施例,进一步阐明本发明。
如图1所示,本发明提供一种适用于弹道导弹INS/CNS/GNSS组合导航系统的最优数据融合方法,包括以下步骤:
S1:构造INS/CNS/GNSS组合导航系统滤波模型:
S1-1:设置INS/CNS/GNSS组合导航系统的状态向量为:
x=[φx,φy,φz,δvx,δvy,δvz,x,y,z,εx,εy,εz,Δx,Δy,Δz]T (1)
其中φx,φy,φz表示发射点惯性系下姿态失准角,δvx,δvy,δvz表示发射点惯性系下的速度误差,x,y,z表示发射点惯性系下的位置误差,εx,εy,εz表示弹体坐标系下的陀螺常值漂移,Δx,Δy,Δz表示弹体坐标系下的加速度计常值偏置,T为转置符号。
S1-2:根据INS/CNS/GNSS组合导航系统的15维状态向量x建立系统的状态方程:
xk=f(xk-1)+vk-1 (2)
其中,f(·)为非线性系统函数,xk-1和xk分别表示k-1时刻和k时刻的状态向量,vk-1表示过程噪声,vk-1的协方差为
S1-3:分别建立INS/GNSS子系统和INS/CNS子系统的量测方程:
在INS/GNSS子系统中,将INS与GNSS输出的位置和速度分别做差,并将其作为量测信息,建立该子系统的量测方程:
z1,k=H1,kxk+ω1,k (3)
其中,z1,k表示INS/GNSS子系统k时刻的量测向量,H1,k表示INS/GNSS子系统k时刻的测量矩阵,ω1,k表示INS/GNSS子系统k时刻的量测噪声,其方差是
在INS/CNS子系统中,该子系统的量测方程将INS与CNS输出的姿态角的差值作为量测信息,建立该子系统的量测方程:
z2,k=H2,kxk+ω2,k (4)
其中,z2,k表示INS/CNS子系统k时刻的量测向量,H2,k表示INS/CNS子系统k时刻的测量矩阵,ω2,k表示INS/CNS子系统k时刻的量测噪声,其方差为
S2:在广义高阶CKF的时间更新阶段和量测更新阶段分别引入自适应渐消因子和最大相关熵准则进行INS/CNS子系统和INS/GNSS子系统的局部状态估计,其具体过程如下:
S2-1:因为INS/GNSS子系统和INS/CNS子系统分别通过局部滤波器1和局部滤波器2采用同样的滤波过程进行局部状态估计,为了避免重复,在此仅具体说明INS/GNSS子系统的滤波过程:,因此,本实施例将步骤S1中的z1,k、z2,k统一写成zk,H1,k和H2,k统一写成Hk,ω1,k和ω2,k统一写成ωk,R1,k和R2,k统一写成Rk。根据INS/GNSS子系统的局部状态估计要求设置核宽度γ的值,设置初始的状态向量、状态误差协方差、渐消因子分别为
τk(0)=1,S0|0=chol(P0|0),chol(·)表示乔列斯分解运算。S2-2:根据公式(5)计算出传播容积点
其中i=1,...,2n2+1,n表示状态向量的维数,[·]i表示集合[·]的第i列。
S2-3:预测k时刻状态向量和误差协方差阵:
其中
S2-4:根据
阳Sk|k-1产生新的容积点:
其中,Sk|k-1=chol(Pk|k-1)。
S2-5:预测k时刻的量测信息:
S2-6:根据INS/GNSS子系统的量测方程和和Pk|k-1构建如下方程:
其中,I表示单位矩阵,Mp,k|k-1=chol(Pk|k-1)、Mr,k=chol(Rk)和
在方程
的两边左乘得到:
令
则公式(11)可重写为:Dk=BkXk+ek (12)
S2-7:对量测噪声方差进行更新:
其中,diag(·)表示矩阵的对角化运算,m为量测信息的维度,
di,k表示Dk的第i个元素,bi,k表示Bk的第i行元素。S2-8:计算状态信息和量测信息之间的互协方差矩阵Pxz,k|k-1:
令i=0,
并设定卡方分布的值,其中θ表示卡方分布的自由度,α表示卡方分布的分位数;如果
则计算:
Sk|k=chol(Pk|k) (19)
继续执行k=k+1,并返回步骤S2-2中的公式(5)所在位置开始执行下一滤波周期;
如果
计算:
根据得到自适应渐消因子τk更新ak(i)的值,判断ak(i)是否满足
若不满足,i=i+1,则返回公式(20)更新渐消因子τk的值;若满足,则得到最优的自适应渐消因子τk的值,然后继续执行下列步骤。S2-9:使用得到的自适应渐消因子τk更新预测的状态预测协方差:
计算卡尔曼增益:
S2-10:估计后验状态和后验协方差矩阵:
Sk|k=chol(Pk|k) (25)
执行k=k+1,并返回步骤S2-2中的公式(5)所在位置开始执行下一滤波周期。
S3:根据最小方差原理和容积准则融合INS/GNSS子系统和INS/CNS子系统的局部估计得到全局最优状态估计:
S3-1:在分别对INS/GNSS子系统和INS/CNS子系统执行完以上滤波过程后便可获得相应的局部后验状态估计
为了便于区分,令INS/GNSS子系统获得局部后验状态估计为INS/CNS子系统获得局部后验状态估计为S3-2:根据最小方差原理和容积准则得到全局最优状态估计:
其中令β=[β1,β2]T,
P11为INS/GNSS子系统执行局部状态估计得到的Pk|k,P22为INS/CNS子系统执行局部状态估计得到的Pk|k,P12表示和估计误差的互协方差矩阵,P21表示和估计误差的互协方差矩阵,P12和P21通过容积准则近似得到:
其中,分别表示INS/GNSS子系统与INS/CNS子系统执行步骤S2的过程中得到的传播容积点
和分别表示INS/GNSS子系统与INS/CNS子系统执行步骤S2的过程中得到的状态向量预测值,和分别表示INS/GNSS子系统与INS/CNS子系统执行步骤S2的过程中得到的新的容积点Xi,k|k-1,和分别表示INS/GNSS子系统与INS/CNS子系统执行步骤S2的过程中得到的量测信息预测值 和分别表示INS/GNSS子系统与INS/CNS子系统执行步骤S2的过程中得到的滤波增益 I表示n维的单位矩阵。