发明名称 一种基于GPS/INS组合导航系统不同测量特性的自适应滤波方法
摘要 本发明提供一种基于GPS/INS组合导航系统不同测量特性的自适应滤波方法,该方法利用IMU测量得到的数据与INS的初始数据进行捷联惯导实时解算,建立GPS/INS速度位置组合滤波系统方程,并根据GPS/INS双系统测量互差分序列,统计计算GPS测量噪声协方差估计值,进行自适应的卡尔曼滤波解算。本发明实现了GPS测量噪声的实时跟踪,以及滤波增益矩阵的自适应调节,提高了组合导航系统定位精度。
申请公布号 CN102096086B 申请公布日期 2012.09.05
申请号 CN201010552746.5 申请日期 2010.11.22
申请人 北京航空航天大学 发明人 张海;常艳红;沈晓蓉;毛友泽;车欢;周启帆
分类号 G01S19/49(2010.01)I 主分类号 G01S19/49(2010.01)I
代理机构 北京永创新实专利事务所 11121 代理人 周长琪
主权项 1.一种基于GPS/INS组合导航系统不同测量特性的自适应滤波方法,包括:步骤一:利用IMU惯性测量单元测量得到沿载体轴相对于惯性空间的角速率和加速度分量信息,与INS惯性导航系统的初始位置、速度、姿态及初始姿态矩阵信息,进行捷联惯导实时解算,得到INS位置、速度、姿态的输出值;步骤二:根据INS的平台误差角方程、速度误差方程、位置误差方程以及陀螺、加速度计误差模型,并根据INS与GPS的位置差、速度差,分别建立GPS/INS组合导航系统卡尔曼滤波器的状态方程与观测方程;步骤三:通过连续系统离散化,建立离散型卡尔曼滤波器的递推方程;其特征在于:还包括下述步骤:步骤四:根据GPS数据采集频率,计算GPS/INS双系统测量互差分序列;(1)INS测量差分序列为:<maths num="0001"><![CDATA[<math><mrow><mi>&Delta;INS</mi><mrow><mo>(</mo><mi>k</mi><mo>)</mo></mrow><mo>=</mo><mfenced open='{' close=''><mtable><mtr><mtd><msub><mi>Z</mi><mi>INS</mi></msub><mrow><mo>(</mo><mi>k</mi><mo>)</mo></mrow><mo>-</mo><msub><mi>Z</mi><mi>INS</mi></msub><mrow><mo>(</mo><mi>k</mi><mo>-</mo><mn>1</mn><mo>)</mo></mrow></mtd><mtd><mi>k</mi><mo>&le;</mo><msub><mi>t</mi><mn>0</mn></msub></mtd></mtr><mtr><mtd><msub><mover><mi>Z</mi><mo>^</mo></mover><mi>INS</mi></msub><mrow><mo>(</mo><mi>k</mi><mo>)</mo></mrow><mo>-</mo><msub><mi>Z</mi><mrow><mi>GPS</mi><mo>/</mo><mi>INS</mi></mrow></msub><mrow><mo>(</mo><mi>k</mi><mo>-</mo><mn>1</mn><mo>)</mo></mrow></mtd><mtd><mi>k</mi><mo>></mo><msub><mi>t</mi><mn>0</mn></msub></mtd></mtr></mtable></mfenced></mrow></math>]]></maths>其中,k为组合导航系统滤波解算时刻,即GPS的数据采集时刻;t<sub>0</sub>为惯性可信时间阈值,根据惯导系统测量精度而确定;Z<sub>INS</sub>(k)为k时刻捷联惯导实时解算的位置、速度;Z<sub>GPS/INS</sub>(k-1)为k-1时刻的组合导航系统输出的位置、速度;<img file="FDA00001649846000012.GIF" wi="125" he="63" />通过如下方式得到:以Z<sub>GPS/INS</sub>(k-1)为k-1时刻的初值,采用k-1~k时间段内IMU测量数据,通过捷联惯导解算步骤,进行k-1~k时间段的惯性递推运算,得到k时刻的位置、速度,即为<img file="FDA00001649846000013.GIF" wi="152" he="65" />(2)GPS测量差分序列为:ΔGPS(k)=Z<sub>GPS</sub>(k)-Z<sub>GPS</sub>(k-1)    k=1,2,3,...其中,Z<sub>GPS</sub>(k)为k时刻GPS输出的位置、速度;(3)GPS/INS双系统测量互差分序列为C(k):C(k)=ΔINS(k)-ΔGPS(k)k=1,2,3,...;步骤五:设定时间阈值T,若组合导航系统滤波解算时刻k大于时间阈值T则直接进行步骤六;若组合导航系统滤波解算时刻k不大于时间阈值T,则对步骤四中得到的GPS/INS双系统测量互差分序列进行小样本统计的可信度判别:1)对GPS/INS双系统测量互差分序列,进行无重叠时间段的开窗统计,得到GPS/INS双系统测量的方差值序列D(k):<img file="FDA00001649846000021.GIF" wi="1064" he="250" />其中,k为组合导航系统滤波解算时刻,即GPS的数据采集时刻;N为无重叠统计窗的窗口宽度,C(i)为GPS/INS的双系统测量互差分序列在i时刻的值;2)对1)中得到的GPS/INS双系统测量的方差值序列D(k),进行小样本统计的可信度判别,该判别条件为:<maths num="0002"><![CDATA[<math><mrow><mi>max</mi><mrow><mo>(</mo><mo>|</mo><mi>D</mi><mrow><mo>(</mo><mi>j</mi><mo>)</mo></mrow><mo>-</mo><mfrac><mn>1</mn><mrow><mi>k</mi><mo>-</mo><mi>N</mi><mo>+</mo><mn>1</mn></mrow></mfrac><munderover><mi>&Sigma;</mi><mrow><mi>i</mi><mo>=</mo><mi>N</mi></mrow><mi>k</mi></munderover><mi>D</mi><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow><mo>|</mo><mo>)</mo></mrow><mo>&lt;</mo><mn>0.4</mn><mo>*</mo><mfrac><mn>1</mn><mrow><mi>k</mi><mo>-</mo><mi>N</mi><mo>+</mo><mn>1</mn></mrow></mfrac><munderover><mi>&Sigma;</mi><mrow><mi>i</mi><mo>=</mo><mi>N</mi></mrow><mi>k</mi></munderover><mi>D</mi><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mrow></math>]]></maths>k≥3N,j={k-2N,k-N,k}其中,max为求最大值的函数;若组合导航系统滤波解算时刻k小于或等于时间阈值T时已满足上述可信度判别条件,则执行步骤六,且进入步骤六时的GPS测量噪声协方差估计值为:<maths num="0003"><![CDATA[<math><mrow><msub><mover><mi>R</mi><mo>^</mo></mover><mi>k</mi></msub><mo>=</mo><mfrac><mn>1</mn><mn>2</mn></mfrac><mi>D</mi><mrow><mo>(</mo><mi>k</mi><mo>)</mo></mrow></mrow></math>]]></maths>k≤T若组合导航系统滤波解算时刻k等于时间阈值T时仍不满足所述的可信度判别条件,则执行步骤六,且进入步骤六时GPS测量噪声协方差估计值为:<maths num="0004"><![CDATA[<math><mrow><msub><mover><mi>R</mi><mo>^</mo></mover><mi>k</mi></msub><mo>=</mo><mfrac><mn>1</mn><mrow><mi>k</mi><mo>-</mo><mn>1</mn></mrow></mfrac><munderover><mi>&Sigma;</mi><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>k</mi></munderover><msup><mrow><mo>[</mo><mi>C</mi><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow><mo>-</mo><mfrac><mn>1</mn><mi>k</mi></mfrac><munderover><mi>&Sigma;</mi><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>k</mi></munderover><mi>C</mi><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow><mo>]</mo></mrow><mn>2</mn></msup><mo>/</mo><mn>2</mn></mrow></math>]]></maths>k=T;若组合导航系统滤波解算时刻k小于时间阈值T时不满足上述可信度判别条件,则等待下一时刻的统计结果;步骤六:根据GPS/INS双系统测量互差分序列,采用连续滑动窗口法计算GPS测量噪声协方差估计值,将其代入步骤三所述的卡尔曼滤波器的递推方程,进行自适应的卡尔曼滤波解算;①将步骤四中的GPS/INS的双系统测量互差分序列C作为统计样本,计算该统计样本在连续滑动窗内的均值和方差,计算公式依次为:<maths num="0005"><![CDATA[<math><mrow><mfenced open='{' close=''><mtable><mtr><mtd><mi>E</mi><mrow><mo>(</mo><mi>k</mi><mo>)</mo></mrow><mo>=</mo><mfrac><mn>1</mn><mi>M</mi></mfrac><munderover><mi>&Sigma;</mi><mrow><mi>i</mi><mo>=</mo><mi>k</mi><mo>-</mo><mi>M</mi><mo>+</mo><mn>1</mn></mrow><mi>k</mi></munderover><mi>C</mi><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow></mtd></mtr><mtr><mtd><msup><mi>&sigma;</mi><mn>2</mn></msup><mrow><mo>(</mo><mi>k</mi><mo>)</mo></mrow><mo>=</mo><mfrac><mn>1</mn><mrow><mi>M</mi><mo>-</mo><mn>1</mn></mrow></mfrac><munderover><mi>&Sigma;</mi><mrow><mi>i</mi><mo>=</mo><mi>k</mi><mo>-</mo><mi>M</mi><mo>+</mo><mn>1</mn></mrow><mi>k</mi></munderover><msup><mrow><mo>[</mo><mi>C</mi><mrow><mo>(</mo><mi>i</mi><mo>)</mo></mrow><mo>-</mo><mi>E</mi><mrow><mo>(</mo><mi>k</mi><mo>)</mo></mrow><mo>]</mo></mrow><mn>2</mn></msup></mtd></mtr></mtable></mfenced><mo>,</mo><mi>k</mi><mo>&GreaterEqual;</mo><mi>M</mi></mrow></math>]]></maths>其中,k为组合导航系统滤波解算时刻,即GPS的数据采集时刻;M为连续滑动窗的窗口宽度,C(i)为GPS/INS的双系统测量互差分序列在i时刻的值,E(k)为k时刻的均值,σ(k)为k时刻的标准差,σ<sup>2</sup>(k)为k时刻的方差;②针对连续滑动窗内的GPS/INS的双系统测量互差分序列,进行野值点判别,该判别条件为:|C(i)-E(k)|>α·σ(k)    k-M+1≤i≤k其中,α为可调剔除因子,若C(i)满足上述野值点判别条件,则将C(i)置为均值E(k),得到新的统计样本,再执行步骤①,重新计算该统计样本在连续滑动窗内的均值和方差,最后执行步骤③;否则直接执行步骤③;③基于统计样本在连续滑动窗内的均值和方差,计算GPS测量噪声协方差的初步估值<img file="FDA00001649846000032.GIF" wi="156" he="64" /><maths num="0006"><![CDATA[<math><mrow><msub><mover><mi>R</mi><mo>^</mo></mover><mi>GPS</mi></msub><mrow><mo>(</mo><mi>k</mi><mo>)</mo></mrow><mo>=</mo><mfrac><mn>1</mn><mn>2</mn></mfrac><msup><mi>&sigma;</mi><mn>2</mn></msup><mrow><mo>(</mo><mi>k</mi><mo>)</mo></mrow></mrow></math>]]></maths>④根据该GPS测量噪声协方差的初步估值<img file="FDA00001649846000034.GIF" wi="157" he="64" />计算GPS测量噪声协方差估计值,计算公式为:<maths num="0007"><![CDATA[<math><mrow><mfenced open='{' close=''><mtable><mtr><mtd><msub><mi>d</mi><mi>k</mi></msub><mo>=</mo><mrow><mo>(</mo><mn>1</mn><mo>-</mo><mi>b</mi><mo>)</mo></mrow><mo>/</mo><mrow><mo>(</mo><mn>1</mn><mo>-</mo><msup><mi>b</mi><mrow><mi>k</mi><mo>+</mo><mn>1</mn></mrow></msup><mo>)</mo></mrow></mtd></mtr><mtr><mtd><msub><mover><mi>R</mi><mo>^</mo></mover><mi>k</mi></msub><mo>=</mo><mrow><mo>(</mo><mn>1</mn><mo>-</mo><msub><mi>d</mi><mi>k</mi></msub><mo>)</mo></mrow><msub><mover><mi>R</mi><mo>^</mo></mover><mrow><mi>k</mi><mo>-</mo><mn>1</mn></mrow></msub><mo>+</mo><msub><mi>d</mi><mi>k</mi></msub><msub><mover><mi>R</mi><mo>^</mo></mover><mi>GPS</mi></msub><mrow><mo>(</mo><mi>k</mi><mo>)</mo></mrow></mtd></mtr></mtable></mfenced><mo>,</mo><mn>0</mn><mo>&lt;</mo><mi>b</mi><mo>&lt;</mo><mn>1</mn></mrow></math>]]></maths>其中,b为遗忘因子,0<b&lt;1;d<sub>k</sub>为k时刻的权系数,<img file="FDA00001649846000036.GIF" wi="42" he="57" />为k时刻的GPS测量噪声协方差估计值;<img file="FDA00001649846000037.GIF" wi="66" he="58" />为k-1时刻的GPS测量噪声协方差估计值;⑤将GPS测量噪声协方差估计值<img file="FDA00001649846000038.GIF" wi="42" he="57" />代入步骤三中卡尔曼滤波器的递推方程,实时调节滤波增益矩阵K<sub>k</sub>,进行自适应卡尔曼滤波解算:<maths num="0008"><![CDATA[<math><mfenced open='{' close=''><mtable><mtr><mtd><msub><mover><mi>X</mi><mo>^</mo></mover><mrow><mi>k</mi><mo>|</mo><mi>k</mi><mo>-</mo><mn>1</mn></mrow></msub><mo>=</mo><msub><mi>&Phi;</mi><mrow><mi>k</mi><mo>,</mo><mi>k</mi><mo>-</mo><mn>1</mn></mrow></msub><msub><mover><mi>X</mi><mo>^</mo></mover><mrow><mi>k</mi><mo>-</mo><mn>1</mn></mrow></msub></mtd></mtr><mtr><mtd><msub><mover><mi>X</mi><mo>^</mo></mover><mi>k</mi></msub><mo>=</mo><msub><mover><mi>X</mi><mo>^</mo></mover><mrow><mi>k</mi><mo>|</mo><mi>k</mi><mo>-</mo><mn>1</mn></mrow></msub><mo>+</mo><msub><mi>K</mi><mi>k</mi></msub><mrow><mo>(</mo><msub><mi>Z</mi><mi>k</mi></msub><mo>-</mo><msub><mi>H</mi><mi>k</mi></msub><msub><mover><mi>X</mi><mo>^</mo></mover><mrow><mi>k</mi><mo>|</mo><mi>k</mi><mo>-</mo><mn>1</mn></mrow></msub><mo>)</mo></mrow></mtd></mtr><mtr><mtd><msub><mi>K</mi><mi>k</mi></msub><mo>=</mo><msub><mi>P</mi><mrow><mi>k</mi><mo>|</mo><mi>k</mi><mo>-</mo><mn>1</mn></mrow></msub><msub><msup><mi>H</mi><mi>T</mi></msup><mi>k</mi></msub><msup><mrow><mo>(</mo><msub><mi>H</mi><mi>k</mi></msub><msub><mi>P</mi><mrow><mi>k</mi><mo>|</mo><mi>k</mi><mo>-</mo><mn>1</mn></mrow></msub><msub><msup><mi>H</mi><mi>T</mi></msup><mi>k</mi></msub><mo>+</mo><msub><mover><mi>R</mi><mo>^</mo></mover><mi>k</mi></msub><mo>)</mo></mrow><mrow><mo>-</mo><mn>1</mn></mrow></msup></mtd></mtr><mtr><mtd><msub><mi>P</mi><mrow><mi>k</mi><mo>|</mo><mi>k</mi><mo>-</mo><mn>1</mn></mrow></msub><mo>=</mo><msub><mi>&Phi;</mi><mrow><mi>k</mi><mo>,</mo><mi>k</mi><mo>-</mo><mn>1</mn></mrow></msub><msub><mi>P</mi><mrow><mi>k</mi><mo>-</mo><mn>1</mn></mrow></msub><msub><msup><mi>&Phi;</mi><mi>T</mi></msup><mrow><mi>k</mi><mo>,</mo><mi>k</mi><mo>-</mo><mn>1</mn></mrow></msub><mo>+</mo><msub><mi>&Gamma;</mi><mrow><mi>k</mi><mo>-</mo><mn>1</mn></mrow></msub><msub><mi>Q</mi><mrow><mi>k</mi><mo>-</mo><mn>1</mn></mrow></msub><msubsup><mi>&Gamma;</mi><mrow><mi>k</mi><mo>-</mo><mn>1</mn></mrow><mi>T</mi></msubsup></mtd></mtr><mtr><mtd><msub><mi>P</mi><mi>k</mi></msub><mo>=</mo><mrow><mo>(</mo><mi>I</mi><mo>-</mo><msub><mi>K</mi><mi>k</mi></msub><msub><mi>H</mi><mi>k</mi></msub><mo>)</mo></mrow><msub><mi>P</mi><mrow><mi>k</mi><mo>|</mo><mi>k</mi><mo>-</mo><mn>1</mn></mrow></msub><msup><mrow><mo>(</mo><mi>I</mi><mo>-</mo><msub><mi>K</mi><mi>k</mi></msub><msub><mi>H</mi><mi>k</mi></msub><mo>)</mo></mrow><mi>T</mi></msup><mo>+</mo><msub><mi>K</mi><mi>k</mi></msub><msub><mover><mi>R</mi><mo>^</mo></mover><mi>k</mi></msub><msub><msup><mi>K</mi><mi>T</mi></msup><mi>k</mi></msub></mtd></mtr></mtable></mfenced></math>]]></maths>其中,X<sub>k</sub>为k时刻的状态向量,即被估计向量,<img file="FDA00001649846000042.GIF" wi="76" he="58" />为状态X<sub>k-1</sub>的卡尔曼滤波估计值,<img file="FDA00001649846000043.GIF" wi="94" he="63" />为利用<img file="FDA00001649846000044.GIF" wi="75" he="57" />得到的对X<sub>k</sub>的一步预测,Φ<sub>k,k-1</sub>为k-1到k时刻的状态转移矩阵,Γ<sub>k-1</sub>为系统噪声矩阵,K<sub>k</sub>为k时刻的滤波增益矩阵,P<sub>k|k-1</sub>为一步预测均方误差矩阵,R<sub>k</sub>为观测噪声协方差矩阵,Q<sub>k-1</sub>为k-1时刻系统噪声协方差矩阵,P<sub>k</sub>为估计均方误差矩阵,H<sub>k</sub>为k时刻的量测矩阵,Z<sub>k</sub>为k时刻的观测向量,I为单位矩阵。
地址 100191 北京市海淀区学院路37号