发明名称 一种卫星信号阻塞情况下的持续导航方法
摘要 一种卫星信号阻塞情况下的持续导航方法,所述卫星信号阻塞是指可见卫星数会不足4颗,实现步骤(1)建立导航解算状态方程,选择用户三轴方向的位置、速度和加速度作为状态量并建立用户运动模型,同时将用户运动模型和时钟模型组合在一起并离散化作为状态方程。本发明在少于4颗星的情况下分别针对高动态和低动态用户提出不同的持续导航方法解决方法。
申请公布号 CN103675880A 申请公布日期 2014.03.26
申请号 CN201310632105.4 申请日期 2013.11.29
申请人 航天恒星科技有限公司 发明人 刘欢;魏玉峤;褚鹏飞;侯春青
分类号 G01S19/49(2010.01)I 主分类号 G01S19/49(2010.01)I
代理机构 中国航天科技专利中心 11009 代理人 安丽
主权项 1.一种卫星信号阻塞情况下的持续导航方法,所述卫星信号阻塞是指可见卫星数会不足4颗,其特征在于实现步骤 如下: (1)建立导航解算状态方程,选择用户三轴方向的位置、速度和加速度作为状态量并建立用户运动模型,同时将用户运动模型和时钟模型组合在一起并离散化作为状态方程; X<sub>k</sub>=Φ<sub>k,k-1</sub>X<sub>k-1</sub>+W<sub>k-1</sub>                 (1) 式中 <img file="FDA0000426607200000011.GIF" wi="1528" he="182" /><img file="FDA0000426607200000012.GIF" wi="1412" he="114" /><img file="FDA0000426607200000013.GIF" wi="1450" he="172" />上式中各量的物理含义如下: <img file="FDA0000426607200000014.GIF" wi="939" he="101" />——系统的状态向量,包括载体在ECEF坐标系中三个方向上的位置、速度、加速度,以及系统的钟差、钟漂;x<sub>u</sub>——载体在ECEF坐标系中X轴方向的位置; v<sub>ux</sub>——载体在ECEF坐标系中X轴方向的速度; a<sub>ux</sub>——载体在ECEF坐标系中X轴方向的加速度; y<sub>u</sub>——载体在ECEF坐标系中Y轴方向的位置; v<sub>uy</sub>——载体在ECEF坐标系中Y轴方向的速度; a<sub>uy</sub>——载体在ECEF坐标系中Y轴方向的加速度; z<sub>u</sub>——载体在ECEF坐标系中Z轴方向的位置; v<sub>uz</sub>——载体在ECEF坐标系中Z轴方向的速度; a<sub>uz</sub>——载体在ECEF坐标系中Z轴方向的加速度; ct<sub>u</sub>——系统钟差; <img file="FDA0000426607200000021.GIF" wi="71" he="75" />——系统钟漂;Φ<sub>k,k-1</sub>——状态转移矩阵,表征系统状态由k-1时刻向k时刻状态转移的情况。状态转移矩阵Φ<sub>k,k-1</sub>由运动模型状态转移矩阵Φ<sub>u</sub>(k,k-1)和时钟模型转移矩阵Φ<sub>c</sub>(k,k-1)组成; Φ<sub>u</sub>(k,k-1)——运动模型状态转移矩阵,表征在状态转移过程中与系统运动模型相关的分量; Φ<sub>c</sub>(k,k-1)——时钟模型转移矩阵,表征在状态转移过程中与系统时钟模型相关的分量; W<sub>k-1</sub>——白噪声序列,表征在k-1时刻系统的过程噪声。白噪声序列W<sub>k-1</sub>由位置白噪声序列W<sub>u</sub>(k-1)和时钟白噪声序列W<sub>c</sub>(k-1)组成; W<sub>u</sub>(k-1)——位置白噪声序列; W<sub>c</sub>(k-1)——时钟白噪声序列; Q<sub>k-1</sub>——白噪声序列的均方差; (2)选择伪距和伪距率作为观测量,建立伪距方程和多普勒方程作为导航解算的量测方程; 伪距方程和多普勒方程分别表示如下: <img file="FDA0000426607200000022.GIF" wi="1714" he="141" /><img file="FDA0000426607200000023.GIF" wi="1820" he="126" />上式中各量的物理含义如下: ρ<sub>j</sub>——接收机测得的第j颗卫星的伪距观测量; D<sub>j</sub>——接收机测得的第j颗卫星的多普勒观测量; x<sub>j</sub>——第j颗卫星在ECEF坐标系中X轴方向的位置; y<sub>j</sub>——第j颗卫星在ECEF坐标系中Y轴方向的位置; z<sub>j</sub>——第j颗卫星在ECEF坐标系中Z轴方向的位置; v<sub>xj</sub>——第j颗卫星在ECEF坐标系中X轴方向的速度; v<sub>yj</sub>——第j颗卫星在ECEF坐标系中Y轴方向的速度; v<sub>zj</sub>——第j颗卫星在ECEF坐标系中Z轴方向的速度; <img file="FDA0000426607200000031.GIF" wi="812" he="121" />——第j颗卫星相对于接收机的向径;x<sub>u</sub>——接收机在ECEF坐标系中X轴方向的位置; y<sub>u</sub>——接收机在ECEF坐标系中Y轴方向的位置; z<sub>u</sub>——接收机在ECEF坐标系中Z轴方向的位置; v<sub>ux</sub>——接收机在ECEF坐标系中X轴方向的速度; v<sub>uy</sub>——接收机在ECEF坐标系中Y轴方向的速度; v<sub>uz</sub>——接收机在ECEF坐标系中Z轴方向的速度; ct<sub>u</sub>——接收机的钟差; <img file="FDA0000426607200000032.GIF" wi="71" he="75" />——接收机的钟漂;然后,根据公式(5)、(6)建立量测方程: Z<sub>k</sub>=H<sub>k</sub>X<sub>k</sub>+V<sub>k</sub>                                (7) 上式中各量的物理含义如下: X<sub>k</sub>——状态向量; Z<sub>k</sub>——观测向量; H<sub>k</sub>——观测矩阵; V<sub>k</sub>——观测噪声协方差矩阵; 设m为接收机观测到的卫星数目,则有: Z<sub>k</sub>=[ρ<sub>1</sub>,ρ<sub>2</sub>,…ρ<sub>m</sub>,D<sub>1</sub>,D<sub>2</sub>,…,D<sub>m</sub>]<sup>T</sup>                    (8) <img file="FDA0000426607200000041.GIF" wi="1443" he="157" /><img file="FDA0000426607200000042.GIF" wi="1462" he="339" /><img file="FDA0000426607200000043.GIF" wi="411" he="152" /><img file="FDA0000426607200000044.GIF" wi="1424" he="154" /><img file="FDA0000426607200000045.GIF" wi="408" he="140" /><img file="FDA0000426607200000046.GIF" wi="1524" he="351" /><img file="FDA0000426607200000047.GIF" wi="794" he="168" /><img file="FDA0000426607200000048.GIF" wi="804" he="168" /><img file="FDA0000426607200000049.GIF" wi="1486" he="178" /><img file="FDA00004266072000000410.GIF" wi="428" he="148" /><img file="FDA00004266072000000411.GIF" wi="435" he="154" /><img file="FDA00004266072000000412.GIF" wi="425" he="153" /><img file="FDA00004266072000000413.GIF" wi="1571" he="94" /><img file="FDA00004266072000000414.GIF" wi="1467" he="117" />上述各式中各量的物理含义如下: Z<sub>k</sub>——观测向量; H<sub>k</sub>——观测矩阵; H<sub>1</sub>——与伪距观测量有关的观测矩阵分量; H<sub>2</sub>——与多普勒观测量有关的观测矩阵分量; <img file="FDA0000426607200000051.GIF" wi="86" he="86" />——观测量ρ<sub>j</sub>对状态变量x<sub>u</sub>的偏导数,反映了系统测量值ρ<sub>j</sub>和系统状态x<sub>u</sub>之间的关系;<img file="FDA0000426607200000052.GIF" wi="86" he="89" />——观测量ρ<sub>j</sub>对状态变量y<sub>u</sub>的偏导数,反映了系统测量值ρ<sub>j</sub>和系统状态y<sub>u</sub>之间的关系;<img file="FDA0000426607200000053.GIF" wi="81" he="84" />——观测量ρ<sub>j</sub>对状态变量z<sub>u</sub>的偏导数,反映了系统测量值ρ<sub>j</sub>和系统状态z<sub>u</sub>之间的关系;<img file="FDA0000426607200000054.GIF" wi="86" he="79" />——观测量D<sub>j</sub>对状态变量x<sub>u</sub>的偏导数,反映了系统测量值D<sub>j</sub>和系统状态x<sub>u</sub>之间的关系;<img file="FDA0000426607200000055.GIF" wi="89" he="83" />——观测量D<sub>j</sub>对状态变量y<sub>u</sub>的偏导数,反映了系统测量值D<sub>j</sub>和系统状态y<sub>u</sub>之间的关系;<img file="FDA0000426607200000056.GIF" wi="83" he="92" />——观测量D<sub>j</sub>对状态变量z<sub>u</sub>的偏导数,反映了系统测量值D<sub>j</sub>和系统状态z<sub>u</sub>之间的关系;<img file="FDA0000426607200000057.GIF" wi="92" he="84" />——观测量D<sub>j</sub>对状态变量v<sub>ux</sub>的偏导数,反映了系统测量值D<sub>j</sub>和系统状态v<sub>ux</sub>之间的关系;<img file="FDA0000426607200000058.GIF" wi="94" he="88" />——观测量D<sub>j</sub>对状态变量v<sub>uy</sub>的偏导数,反映了系统测量值D<sub>j</sub>和系统状态v<sub>uy</sub>之间的关系;<img file="FDA0000426607200000059.GIF" wi="94" he="86" />——观测量D<sub>j</sub>对状态变量v<sub>uz</sub>的偏导数,反映了系统测量值D<sub>j</sub>和系统状态v<sub>uz</sub>之间的关系;ρ<sub>j</sub>——接收机测得的第j颗卫星的伪距观测量; D<sub>j</sub>——接收机测得的第j颗卫星的多普勒观测量; J<sub>j</sub>——卫星与接收机的相对速度与相对位置向量的点积; r<sub>j</sub>——第j颗卫星相对于接收机的向径; (3)实现滤波计算回路和增益计算回路 根据步骤(1)和步骤(2)中建立的状态方程和量测方程,将状态初始量输入状态方程实现状态一步预测,通过均方误差阵一步预测和滤波增益方程实现增益矩阵K的更新,增益矩阵K代入量测方程得到状态增量;状态的一步预测加上状态增量作为下一次状态一步预测的输入即实现了滤波计算回路;根据增益矩阵K计算新的均方误差矩阵P作为下一次均方误差一步预测的输入即实现了增益计算回路; (4)根据步骤(3)中滤波计算回路和增益计算回路中的状态方程选取运动模型,在低动态用户时,选取位置和速度模型作为运动模型,所述低动态用户是指汽车、轮船及地面手持设备;在高动态用户时,选取位置、速度和加速度模型,所述高动态用户是指导弹、飞机及卫星,所述加速度模型中的加速度信息使用惯性导航提供; (5)将上述步骤中建立的状态方程、量测方程、滤波计算回路执行流程、增益计算回路执行流程以及不同动态用户模型的选取,得到完整的扩展卡尔曼滤波方程,用于实现卫星阻塞情况下用户的持续导航算法;其中扩展卡尔曼滤波方程如下: <img file="FDA0000426607200000061.GIF" wi="994" he="109" /><img file="FDA0000426607200000062.GIF" wi="1117" he="90" /><img file="FDA0000426607200000063.GIF" wi="1079" he="105" />Z<sub>k</sub>=H<sub>k</sub>X<sub>k</sub>+V<sub>k</sub>+y<sub>k</sub>            (12) <img file="FDA0000426607200000064.GIF" wi="1107" he="111" /><img file="FDA0000426607200000071.GIF" wi="1277" he="92" />上式中各量的物理含义如下: <img file="FDA0000426607200000072.GIF" wi="106" he="92" />——k-1时刻系统状态变量的最优估计值(后验估计值);<img file="FDA0000426607200000073.GIF" wi="133" he="84" />——基于<img file="FDA0000426607200000074.GIF" wi="92" he="84" />推出的k时刻系统状态变量的先验估计值;Φ<sub>k,k-1</sub>——系统由k-1时刻到k时刻的状态转移矩阵; P<sub>k/k-1</sub>——k时刻系统状态变量先验估计误差的均方差矩阵; P<sub>k-1</sub>——k-1时刻系统状态变量后验估计误差的均方差矩阵; Q<sub>k-1</sub>——k-1时刻的过程噪声协方差矩阵; Γ<sub>k-1</sub>——k-1时刻的噪声输入矩阵; K<sub>k</sub>——k时刻的滤波增益矩阵; H<sub>k</sub>——系统在k时刻的观测矩阵; R<sub>k</sub>——k时刻的观测噪声协方差矩阵; Z<sub>k</sub>——k时刻的观测向量; X<sub>k</sub>——k时刻的状态向量; V<sub>k</sub>——k时刻的观测噪声协方差矩阵; <img file="FDA0000426607200000075.GIF" wi="76" he="84" />——k时刻系统状态变量的最优估计值(后验估计值);P<sub>k</sub>——k时刻系统状态变量后验估计误差的均方差矩阵; I——单位矩阵。 
地址 100086 北京市海淀区知春路82号院