技术领域本发明属于组合导航的技术领域，具体涉及一种基于强跟踪的相关熵扩展卡尔曼滤波的组合导航方法，可用于提高飞行器、导弹、舰船、地面车辆的SINS/GNSS组合导航系统的精度。
背景技术：
高精度的组合导航一直以来是导航领域关注的关键问题。为了获得更高的导航精度，学者们和工业界人士致力于寻找具有鲁棒性和高精度的滤波方法，用于对载体的位置、速度、姿态等状态进行估计。捷联惯性导航系统(Serialinertialnavigationsystem，SINS)主要基于加速度二次积分的推算过程，自主性好、导航信息量大，但长时间工作导航精度变差。全球导航卫星系统(GlobalNavigationSatelliteSystem，GNSS)是以导航卫星为基础的无线电导航系统，具有较高的定位精度和测速精度，但是机载GNSS接收机接收信号时受载体飞行姿态的影响，地面GNSS接收机受所在地域的影响。这两种导航系统单独使用时，均存在不同的限制，难以达到理想的导航精度。而SINS/GNSS组合导航是组合上述两种导航设备的特点，利用导航滤波算法和计算机技术对多种导航信息进行综合处理，显著地提高了导航系统性能，特别是导航精度。针对SINS/GNSS组合导航系统而言，其组合滤波的非线性状态方程包括捷联惯导力学编排方程和姿态误差方程，故采用的导航滤波方法宜选取为非线性滤波方法。经典的扩展卡尔曼滤波方法通过时间更新和量测更新得到对状态进行估计，但当模型不精确或者噪声统计信息是非高斯时，滤波效果会很差，甚至会导致滤波过程发散。基于最大相关熵准则，有学者提出了相应的相关熵卡尔曼滤波算法，用于解决滤波过程中噪声的非高斯统计特性问题，取得了不错的效果。但是其未考虑模型的不确定性，针对这种模型的不确定性问题，基于强跟踪理论将渐消因子引入到相关熵卡尔曼滤波算法中，提出了一种基于强跟踪的相关熵扩展卡尔曼滤波的SINS/GNSS组合导航方法，能有效地解决SINS/GNSS组合导航系统的非线性问题，以及噪声统计特性的非高斯问题，提高了SINS/GNSS组合导航的精度。
技术实现要素：
针对SINS/GNSS组合导航系统模型的不确定性和噪声的非高斯特性，本发明提供了一种基于强跟踪的相关熵扩展卡尔曼滤波的组合导航方法，基于最大相关熵准则建立了相关熵卡尔曼滤波方法，将渐消因子引入到相关熵扩展卡尔曼滤波方法，提高滤波精度。该方法利用强跟踪理论的特点，最大限度的消除模型的不确定性对SINS/GNSS组合导航的精度的影响，提高导航滤波精度，增强导航滤波稳定性。为了达到上述目的，本发明提出了一种基于强跟踪的相关熵扩展卡尔曼滤波的组合导航方法，包括以下五个步骤：步骤一：建立SINS/GNSS组合导航系统的动力学方程在东北天地理坐标系下，SINS/GNSS组合导航系统的动力学方程为：φ·E=(ωiesinL+vERN+htanL)φN-(ωiecosL+vERN+h)φU-δvERM+h+vN(RM+h)2δh+ϵbE+wgE---(1)]]>φ·N=-(ωiesinL+vERN+htanL)φE-vNRM+hφU+δvERN+h-vE(RN+h)2δh-ωieδLsin+ϵbN+wgN---(2)]]>φ·U=(ωiecosL+vERN+h)φE+vNRM+hφN+tanLRN+hδvE-vEtanL(RM+h)2δh+(ωiecosL+vEsec2LRN+h)δL+ϵbU+wgU---(3)]]>v·E=(2ωiesinL+vERN+htanL)vN-(2ωiecosL+vERN+h)vU+fE+fNφU-fUφN---(4)]]>v·N=-(2ωiesinL+vERN+htanL)vE-vNRM+hvU+fEφU+fN-fUφE---(4)]]>v·U=(2ωiecosL+vERN+htanL)vE+vN2RM+h-g-fEφN+fNφE+fU---(5)]]>L·=vNRM+h---(6)]]>λ·=vE(RN+h)cosL---(7)]]>h·=vU---(8)]]>其中，φE,φN,φU分别是东、北、天的高度，vE,vN,vU分别是东、北、天的速度，L,λ,h分别为地理纬度、经度和高度；fE,fN,fU分别为东北天三个轴向的加速度计输出的比力测量值，δvE,δvN,δL,δh分别是惯导系统的速度、位置输出减去GNSS对应输出的近似值；RM为RM＝Re(1-2e+3esin2L)，RN为RN＝Re(1-esin2L)，其中Re＝6378137.0m为地球长半径，e＝1/298.257是地球扁率，ωie＝7.292×10-5rad/s是地球自转角速度，g＝9.78m/s2是重力加速度；地球陀螺东北天常值漂移εbE,εbN,εbU分别为：ϵ·bE=0,ϵ·bN=0,ϵ·bU=0---(9)]]>其中，wgE,wgN,wgU分别是陀螺东、北、天的噪声；状态向量为x＝[φE,φN,φU,vE,vN,vU,L,λ,h,εbE,εbN,εbU]T，则其动力学方程为x·=f(x)+w---(10)]]>式中，w为系统的噪声；步骤二：建立SINS/GNSS组合导航系统的量测方程选取GNSS输出的载体速度信息vEG,vNG,vUG和位置信息LG,λG,hG作为量测值组成量测方程：y＝[vEG,vNG,vUG,LG,λG,hG]T(11)则相应的量测方程为y＝Hx+v(12)其中，量测方程系数矩阵H为H＝[I6×6,O6×6]，I为单位矩阵，O为零矩阵，v为量测噪声。步骤三：离散化SINS/GNSS组合导航的动力学方程和量测方程将式(10)的动力学方程和式(12)的量测方程进行离散化，得到非线性的状态方程和线性的量测方程：xk＝f(xk-1)+wk-1(13)yk＝Hxk+vk(14)其中，xk为状态向量，f(·)为非线性离散函数，wk为系统噪声向量；yk为量测向量；vk为量测噪声向量，k是指第k步，代表tk时刻；系统噪声向量wk和量测噪声向量vk满足：Cov[wk,vj]=E[wkvjT]=0---(15)]]>步骤四：对非线性离散后的动力学方程进行线性化首先，将式(13)中的非线性离散函数f(xk-1)围绕估计值展开成泰勒级数，并略去二阶以上的项，得线性化之后的状态方程：xk＝Φk/k-1xk-1+uk-1+wk-1(16)其中，Φk/k-1=∂f(xk-1)∂xk-1|xk-1=x^k-1---(17)]]>uk-1=f(x^k-1)-∂f(xk-1)∂xk-1|xk-1=x^k-1·x^k-1---(18)]]>其中，表示状态x的估计值；步骤五：对线性化的动力学方程采用基于强跟踪的相关熵扩展卡尔曼滤波算法，输出SINS/GNSS组合导航的信息1.对由式(16)和式(14)组成的非线性离散系统的状态和误差方差阵Pk进行初始化和时间更新；对状态和误差方差阵分别赋予初值x^0=E[x0]---(19)]]>P0=E[(x0-x^0)(x0-x^0)T]---(20)]]>其中，k取值为k＝1,2,3,…,N，N由滤波时间和采样周期决定；2.第k-1步tk-1时刻的状态估计值为和误差方差阵为Pk-1，对第k步tk时刻的状态向量xk和误差方差阵Pk进行时间更新分别得到状态的一步预测x^k-=f(x^k-1)---(21)]]>和相应的一步预测的误差方差阵：Pk/k-1=λkΦk/k-1Pk-1Φk/k-1T+Qk-1---(22)]]>其中，“－”表示先验估计，Qk-1表示系统的噪声方差阵，λk为渐消因子，其计算过程如下：λk=λ,λ>11,λ≤1---(23)]]>其中，λ=tr(Nk)tr(Mk)---(24)]]>Nk＝Sk(Vk-HQk-1HT)-Rk(25)Mk=SkHΦk/k-1Pk-1Φk/k-1THT---(26)]]>其中，tr(·)为求矩阵迹的算子，Rk表示量测噪声的方差阵；Vk为实际输出残差序列的协方差阵，由下式估算：Vk=γ1γ1T,k=0ρVk-1+γk-1γk-1T1+ρ,k≥1---(27)]]>其中，为输出残差序列，ρ为遗忘因子，且满足0<ρ≤1；4.计算相关信息熵系数Sk：Sk=Gσ(||yk-Hx^k-||Rk-1)Gσ(||x^k--Φk/k-1x^k-1||Pk/k-1-1)---(28)]]>其中，Gσ(·)表示带宽为σ的高斯核函数；5.计算滤波的增益矩阵Kk：Kk=(Pk/k-1-1-SkHTRk-1H)-1SkHTRk-1---(29)]]>6.计算状态估计x^k=x^k-+Kk(yk-Hx^k-)---(30)]]>7.计算状态估计的误差方差阵Pk：Pk=(I-KkH)Pk/k-1(I-KkH)T+KkRkKkT---(31)]]>通过以上7步循环可得到载体三维的位置(L,λ,h)和速度(vE,vN,vU)的估计值。本发明的有益效果是：本发明将强跟踪理论中的渐消因子引入到相关熵扩展卡尔曼滤波算法中，有效地解决了非线性SINS/GNSS组合导航系统的模型不确定性和噪声统计特性的非高斯问题，从而提高了导航精度，增强了导航过程的稳定性。附图说明图1为本发明的流程图。图2为本发明中的基于强跟踪的相关熵扩展卡尔曼滤波算法的原理图。具体实施方式下面结合附图和实施例对本发明作详细说明。一种基于强跟踪的相关熵扩展卡尔曼滤波的组合导航方法，其流程图如图1所示，包括以下五个步骤:步骤一、建立SINS/GNSS组合导航系统的动力学方程在东北天地理坐标系下，SINS/GNSS组合导航系统的动力学方程为φ·E=(ωiesinL+vERN+htanL)φN-(ωiecosL+vERN+h)φU-δvERM+h+vN(RM+h)2δh+ϵbE+wgE---(1)]]>φ·N=-(ωiesinL+vERN+htanL)φE-vNRM+hφU+δvERN+h-vE(RN+h)2δh-ωieδLsin+ϵbN+wgN---(2)]]>φ·U=(ωiecosL+vERN+h)φE+vNRM+hφN+tanLRN+hδvE-vEtanL(RM+h)2δh+(ωiecosL+vEsec2LRN+h)δL+ϵbU+wgU---(3)]]>v·E=(2ωiesinL+vERN+htanL)vN-(2ωiecosL+vERN+h)vU+fE+fNφU-fUφN---(4)]]>v·N=-(2ωiesinL+vERN+htanL)vE-vNRM+hvU+fEφU+fN-fUφE---(4)]]>v·U=(2ωiecosL+vERN+htanL)vE+vN2RM+h-g-fEφN+fNφE+fU---(6)]]>L·=vNRM+h---(7)]]>λ·=vE(RN+h)cosL---(8)]]>h·=vU---(9)]]>式中，φE,φN,φU分别是东、北、天的高度，vE,vN,vU分别是东、北、天的速度，L,λ,h分别为地理纬度、经度和高度。fE,fN,fU分别为三个轴向的加速度计输出的比力测量值，δvE,δvN,δL,δh分别是惯导系统的速度、位置输出减去GNSS对应输出的近似值。RM为RM＝Re(1-2e+3esin2L)，RN为RN＝Re(1-esin2L)，Re＝6378137.0m为地球长半径，e＝1/298.257是地球扁率，ωie＝7.292×10-5rad/s是地球自转角速度，g＝9.78m/s2是重力加速度。地球陀螺常值漂移εbE,εbN,εbU定义为ϵ·bE=0,ϵ·bN=0,ϵ·bU=0---(10)]]>式中，wgE,wgN,wgU分别是陀螺东、北、天的噪声。取状态向量为x＝[φE,φN,φU,vE,vN,vU,L,λ,h,εbE,εbN,εbU]T，则其动力学方程为x·=f(x)+w---(11)]]>式中，w为系统噪声。步骤二、建立SINS/GNSS组合导航系统的量测方程选取GNSS输出的载体速度信息vEG,vNG,vUG和位置信息LG,λG,hG作为量测值组成量测方程：y＝[vEG,vNG,vUG,LG,λG,hG]T(12)则相应的量测方程为y＝Hx+v(13)式中，量测方程系数矩阵H为H＝[I6×6,O6×6]，I为单位矩阵,O为零矩阵，v为量测噪声。步骤三、离散化SINS/GNSS组合导航的动力学方程和量测方程离散化SINS/GNSS组合导航的动力学方程和量测方程，建立SINS/GNSS组合导航的非线性离散系统。将上述式(10)的动力学方程和式(12)的量测方程进行离散化，得到非线性的状态方程和线性的量测方程：xk＝f(xk-1)+wk-1(14)yk＝Hxk+vk(15)式中，xk为状态向量，f(·)为非线性离散函数，wk为系统噪声向量；yk为量测向量，vk为量测噪声向量，k是指第k步，代表tk时刻。系统噪声向量wk和量测噪声向量vk满足Cov[wk,vj]=E[wkvjT]=0---(16)]]>步骤四、对非线性离散后的动力学方程进行线性化首先，将式(13)中的非线性离散函数f(xk-1)围绕估计值展开成泰勒级数，并略去二阶以上项，得线性化之后的状态方程：xk＝Φk/k-1xk-1+uk-1+wk-1(17)式中，Φk/k-1=∂f(xk-1)∂xk-1|xk-1=x^k-1---(18)]]>uk-1=f(x^k-1)-∂f(xk-1)∂xk-1|xk-1=x^k-1·x^k-1---(19)]]>其中，表示状态x的估计值。利用泰勒级数对非线性离散函数进行泰勒级数展开时，只取前两项。以上四个步骤为前期准备过程，下面一个步骤是基于强跟踪的相关熵扩展卡尔曼滤波方法的滤波过程。步骤五、对线性化的动力学方程采用基于强跟踪的相关熵扩展卡尔曼滤波算法，输出SINS/GNSS组合导航的信息本发明的基于强跟踪的相关熵扩展卡尔曼滤波算法的原理图参见图2。1.对由式(16)和式(14)组成的非线性离散系统的状态和误差方差阵Pk进行初始化和时间更新：首先，对状态和误差方差阵分别赋予初值：x^0=E[x0]---(20)]]>P0=E[(x0-x^0)(x0-x^0)T]---(21)]]>其中，k取值为k＝1,2,3,…,N，N由滤波时间和采样周期决定。比如当滤波时间为100秒，采样周期为1赫兹时，N＝100/1＝100。2.设第k-1步tk-1时刻的状态估计值为和误差方差阵为Pk-1，在此基础之上对第k步即tk时刻的状态向量xk和误差方差阵Pk进行时间更新分别得到状态的一步预测x^k-=f(x^k-1)---(22)]]>和相应的一步预测的误差方差阵Pk/k-1=λkΦk/k-1Pk-1Φk/k-1T+Qk-1---(23)]]>式中，上标“－”表示先验估计，Qk-1表示系统的噪声方差阵；λk为渐消因子，其计算过程如下：λk=λ,λ>11,λ≤1---(24)]]>式中，λ=tr(Nk)tr(Mk)---(25)]]>Nk＝Sk(Vk-HQk-1HT)-Rk(26)Mk=SkHΦk/k-1Pk-1Φk/k-1THT---(27)]]>其中，tr(·)为求矩阵迹的算子，Rk表示量测噪声的方差阵。Vk为实际输出残差序列的协方差阵，一般由下式估算：Vk=γ1γ1T,k=0ρVk-1+γk-1γk-1T1+ρ,k≥1---(28)]]>其中，为输出残差序列，ρ为遗忘因子，满足0<ρ≤1，一般来说ρ取为0.95。4.计算相关信息熵系数Sk：Sk=Gσ(||yk-Hx^k-||Rk-1)Gσ(||x^k--Φk/k-1x^k-1||Pk/k-1-1)---(29)]]>其中，范数表示表示带宽为σ的高斯核函数，即Gσ(||yk-Hx^k-||Rk-)=exp{||yk-Hx^k-||Rk-122σ2