CN115112126B - 一种gnss/ins组合导航系统保护级反演方法 - Google Patents

一种gnss/ins组合导航系统保护级反演方法 Download PDF

Info

Publication number
CN115112126B
CN115112126B CN202211048910.8A CN202211048910A CN115112126B CN 115112126 B CN115112126 B CN 115112126B CN 202211048910 A CN202211048910 A CN 202211048910A CN 115112126 B CN115112126 B CN 115112126B
Authority
CN
China
Prior art keywords
fault
calculating
gnss
ins
protection level
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
CN202211048910.8A
Other languages
English (en)
Other versions
CN115112126A (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.)
Jiaoxin Beidou Technology Co ltd
Jiaoxin Beidou Beijing Information Technology Co ltd
Original Assignee
Jiaoxin Beidou Technology Co ltd
Jiaoxin Beidou Beijing Information Technology Co ltd
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 Jiaoxin Beidou Technology Co ltd, Jiaoxin Beidou Beijing Information Technology Co ltd filed Critical Jiaoxin Beidou Technology Co ltd
Priority to CN202211048910.8A priority Critical patent/CN115112126B/zh
Publication of CN115112126A publication Critical patent/CN115112126A/zh
Application granted granted Critical
Publication of CN115112126B publication Critical patent/CN115112126B/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/10—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration
    • G01C21/12—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration executed aboard the object being navigated; Dead reckoning
    • G01C21/16—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration executed aboard the object being navigated; Dead reckoning by integrating acceleration or speed, i.e. inertial navigation
    • G01C21/165—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration executed aboard the object being navigated; Dead reckoning by integrating acceleration or speed, i.e. inertial navigation combined with non-inertial navigation instruments
    • 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
    • G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C25/00—Manufacturing, calibrating, cleaning, or repairing instruments or devices referred to in the other groups of this subclass
    • 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
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Remote Sensing (AREA)
  • Radar, Positioning & Navigation (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Data Mining & Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Automation & Control Theory (AREA)
  • General Engineering & Computer Science (AREA)
  • Software Systems (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Computing Systems (AREA)
  • Manufacturing & Machinery (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Evolutionary Biology (AREA)
  • Operations Research (AREA)
  • Probability & Statistics with Applications (AREA)
  • Navigation (AREA)
  • Position Fixing By Use Of Radio Waves (AREA)

Abstract

本发明公开了一种GNSS/INS组合导航系统保护级反演方法，属于导航技术领域，包括以下步骤：1、参数初始化；2、分别针对GNSS和INS进行故障检测和排除；3、利用故障检测和排除后的GNSS和INS测量信息进行组合导航解算；4、针对四种故障模式ḠĪ、ḠI、GSĪ和GSI进行完好性风险需求分配；5、计算ḠĪ故障模式保护级；6、计算ḠI故障模式保护级；7、计算GSĪ故障模式保护级；8、计算GSI故障模式保护级；9、计算最终保护级。本发明综合考虑了组合导航系统中的GNSS和INS故障风险，在进行保护级反演时能够兼顾GNSS和INS故障风险，有助于实现更为严谨和可靠的完好性监测。

Description

一种GNSS/INS组合导航系统保护级反演方法
技术领域
本发明涉及GNSS/INS组合导航技术领域，具体涉及一种GNSS/INS组合导航保护级反演方法。
背景技术
随着社会的发展，尤其是科学技术的进步，大大促进了社会生产力的发展，使得人们的生活得到了极大地改善与提高，特别是以导航系统为代表的信息化的广泛应用，为社会的进步注入了强劲的动力。
完好性是描述导航系统性能的重要指标之一，是对导航系统所提供信息正确性的测量，也包括系统在无法用于导航时向用户发出告警的能力。完好性风险为导航系统未检测到故障但定位误差超出指定告警门限的概率。保护级为指定完好性风险需求下用户位置误差的安全边界，因此保护级反演是完好性监测中不可或缺的环节。
现有完好性监测技术主要针对全球导航卫星系统(Global NavigationSatellite System,GNSS)设计，如专利申请号为CN202111461447.5公开了一种 APNT服务的定位和完好性监测方法及系统，该方法包括：确定目标场景下的定位精度需求；当所述定位精度需求为高精度定位时，采用组合定位算法确定航空器的位置，并采用多解分离方式对组合定位进行完好性监测；当所述定位精度需求为低精度定位时，判断航空器是否为高空用户；若否，则采用基于LDACS 的高空用户与低空用户的空对空定位算法确定航空器的位置，并采用基于最小二乘残差法对空对空定位进行完好性监测。
专利申请号为201811372441.9的公开了一种GNSS增强系统的完好性监测系统，其特征在于，包括：第一级监测模块和第二级监测模块，所述第一级监测模块与所述第二级监测模块连接；所述第一级监测模块包括播发前完好性监测模块，所述第二级监测模块包括播发后完好性监测模块。本发明还提供一种 GNSS增强系统的完好性监测方法，其特征在于，包括：T0时刻对即将播发的改正数进行播发前完好性监测的步骤；对播发前完好性监测结果进行电文编码播发的步骤；在T0+TP时刻对已播发的改正参数进行播发后完好性监测的步骤；对播发后完好性监测结果进行电文编码并播发的步骤。
由于GNSS单点定位和组合导航在原理上的差异，此类完好性监测中的保护级反演方法无法直接应用于GNSS/INS组合导航系统，进而无法保证组合导航定位结果的可靠性。而现有针对GNSS/INS组合导航的完好性监测技术往往以惯性导航系统(InertialNavigation System,INS)辅助GNSS(如专利申请号为 CN201810088308.4公开了一种组合导航慢变斜坡故障完好性监测方法。该方法针对传统组合导航系统中噪声级慢变斜坡故障，通过卫星导航与惯导组合导航扩展卡尔曼滤波的过程量测矩阵以及系统状态协方差矩阵，构建递推式指数型完好性检测统计量，满足导航作业中虚警率和漏警率需求的前提下，设置递推时间间隔并进行检测统计量与检测阈值比较，最终完成完好性监测并及时向用户发出告警信息)，在保护级反演过程中未考虑INS故障的可能性。但在实际应用场景中，低成本INS故障所引发的完好性风险不容忽视，此时若只考虑 GNSS故障，将难以保障完好性监测中保护级反演环节的严谨性。
因此，提供一种GNSS/INS组合导航系统保护级反演方法，是一个值得研究的问题。
发明内容
为了解决现有技术在保护级反演过程中难以顾及INS故障风险的问题，本发明提供了一种GNSS/INS组合导航系统保护级反演方法；该方法在进行保护级反演过程中能够兼顾GNSS和INS故障风险，从而实现更为严谨和可靠的完好性监测。
根据本申请的一方面，提供了一种组合导航系统保护级反演方法，该组合导航系统包括全球导航卫星系统和惯性导航系统，所述方法包括以下步骤：
步骤1、参数初始化，包括全球导航卫星系统卫星先验故障概率、惯性导航系统先验故障概率、完好性风险需求值；
步骤2、分别针对全球导航卫星系统和惯性导航系统进行故障检测和排除；
步骤3、利用故障检测和排除后的全球导航卫星系统和惯性导航系统测量信息进行组合导航解算，解算方法为扩展卡尔曼滤波算法；
步骤4、针对第一故障模式、第二故障模式、第三故障模式、第四故障模式进行完好性风险需求分配，其中分配给第一故障模式的完好性风险需求为第一完好性风险需求值，分配给第二故障模式的完好性风险需求为第二完好性风险需求值，分配给第三故障模式的完好性风险需求为第三完好性风险需求值，分配给第四故障模式的完好性风险需求为第四完好性风险需求值，第一故障模式为全球导航卫星系统无故障且惯性导航系统无故障；第二故障模式为全球导航卫星系统无故障且惯性导航系统故障；第三故障模式为全球导航卫星系统单星故障且惯性导航系统无故障；第四故障模式为全球导航卫星系统单星故障且惯性导航系统故障；
步骤5、计算第一故障模式的保护级：
步骤6、计算第二故障模式的保护级：
步骤7、计算第三故障模式的保护级：
步骤8、计算第四故障模式的保护级：
步骤9、计算最终保护级。
进一步地，第一完好性风险需求值、第二完好性风险需求值、第三完好性风险需求值、第四完好性风险需求值之和小于或等于完好性风险需求值。
进一步地，步骤5包括：
步骤5-1、根据第一故障模式的先验发生概率计算第一故障模式的双侧分位数，其中，第一故障模式的先验发生概率根据全球导航卫星系统卫星先验故障概率、惯性导航系统先验故障概率、当前时刻参与滤波解算的可见卫星数计算；
步骤5-2、根据第一故障模式的双侧分位数计算第一故障模式的保护级。
进一步地，步骤6包括：
步骤6-1、根据滤波器的增益矩阵、滤波器的观测矩阵计算残差传播矩阵；
步骤6-2、根据第二故障模式的先验发生概率计算第二故障模式的双侧分位数，其中，第二故障模式的先验发生概率根据全球导航卫星系统卫星先验故障概率、惯性导航系统先验故障概率、当前时刻参与滤波解算的可见卫星数计算；
步骤6-3、根据残差传播矩阵、滤波器的信息矢量、第二故障模式的双侧分位数计算第二故障模式的保护级。
进一步地，步骤7包括：
步骤7-1、根据残差传播矩阵、可见星观测噪声标准差、滤波器的增益矩阵、滤波器的观测矩阵计算每颗卫星对应的斜率；
步骤7-2、计算最小可检测偏差，包括：
根据第三完好性风险需求值、第三故障模式的先验发生概率、非中心化卡方分布概率累积分布函数、全球导航卫星系统故障检测时的检验统计量计算最小可检测偏差，其中，第三故障模式的先验发生概率根据全球导航卫星系统卫星先验故障概率、惯性导航系统先验故障概率、当前时刻参与滤波解算的可见卫星数计算；
步骤7-3、根据第三故障模式的先验发生概率计算第三故障模式的双侧分位数；
步骤7-4、根据卫星对应的斜率、最小可检测偏差、第三故障模式的双侧分位数计算第三故障模式的保护级。
进一步地，步骤8包括：
步骤8-1、根据第四故障模式的先验发生概率计算第四故障模式的双侧分位数，第四故障模式的先验发生概率根据全球导航卫星系统卫星先验故障概率、惯性导航系统先验故障概率、当前时刻参与滤波解算的可见卫星数计算；
步骤8-2、计算第四故障模式的保护级，包括：根据卫星对应的斜率、最小可检测偏差、第四故障模式的双侧分位数计算第四故障模式的保护级。
进一步地，步骤9包括：将第一故障模式的保护级、第二故障模式的保护级、第三故障模式的保护级、第四故障模式的保护级中绝对值的最大值作为最终保护级。
根据本申请的另一方面，提供了一种GNSS/INS组合导航系统保护级反演方法，包括以下步骤：
步骤1、参数初始化，包括GNSS卫星先验故障概率Psat、INS先验故障概率Pins、完好性风险需求值Ireq以及连续性风险需求值Creq；
步骤2、分别针对GNSS和INS进行故障检测和排除；
步骤3、利用故障检测和排除后的GNSS和INS测量信息进行组合导航解算，解算方法为扩展卡尔曼滤波算法；
步骤4、针对四种故障模式
Figure GDA0003896769470000051
和GsI进行完好性风险需求分配，其中分配给GNSS无故障且INS无故障
Figure GDA0003896769470000061
的完好性风险需求为
Figure GDA0003896769470000062
分配给 GNSS无故障且INS故障模式
Figure GDA0003896769470000063
的完好性风险需求为
Figure GDA0003896769470000064
分配给GNSS单星故障且INS无故障
Figure GDA0003896769470000065
的完好性风险需求为
Figure GDA0003896769470000066
分配给GNSS单星故障且 INS故障GsI的完好性风险需求为Ir(GsI)；
步骤5、计算
Figure GDA0003896769470000067
故障模式保护级：
步骤6、计算
Figure GDA0003896769470000068
故障模式保护级：
步骤7、计算
Figure GDA0003896769470000069
故障模式保护级：
步骤8、计算GsI故障模式保护级：
步骤9、计算最终保护级。
进一步地，步骤4中的四种故障模式下的完好性风险需求应满足：
Figure GDA00038967694700000610
进一步地，步骤5包括步骤5-1、计算
Figure GDA00038967694700000611
的双侧分位数
Figure GDA00038967694700000612
其表达式为：
Figure GDA00038967694700000613
其中：
Figure GDA00038967694700000614
表达
Figure GDA00038967694700000615
的先验发生概率，n为当前时刻参与滤波解算的可见卫星数；Q-1(·)表示标准正态分布概率累积分布函数的逆函数；
步骤5-2、计算
Figure GDA00038967694700000616
故障模式的保护级，其表达式为：
Figure GDA00038967694700000617
其中σq为q方向滤波估计误差标准差，其值为滤波估计协方差矩阵Pk第q 个对角线元素的平方根。
进一步地，步骤6包括步骤6-1、计算残差传播矩阵Bk，其表达式为：
Figure GDA00038967694700000618
其中，Ik表示单位矩阵，Kk表示滤波器的增益矩阵，Hk表示滤波器的观测矩阵；
步骤6-2、计算
Figure GDA0003896769470000071
的双侧分位数
Figure GDA0003896769470000072
其表达式为：
Figure GDA0003896769470000073
其中，
Figure GDA0003896769470000074
表达
Figure GDA0003896769470000075
的先验发生概率；
步骤6-3、计算计算
Figure GDA0003896769470000076
故障模式的保护级，其表达式为：
Figure GDA0003896769470000077
其中：γk为滤波器的新息矢量；(·)q表示对应矢量的第q个分量；
Figure GDA0003896769470000078
表示矩阵(Bk+Kk)Rk(Bk+Kk)T第q行对角线元素的平方根，Rk测量噪声协方差矩阵。
进一步地，步骤7包括步骤7-1、计算每颗卫星对应的斜率Slopei,q，其表达式为：
Figure GDA0003896769470000079
其中：tq,i表示矩阵Tk中第q行第i列的元素；si,i为矩阵S中第i行第i列的元素；σi为第i颗可见星观测噪声标准差；矩阵Tk和S的表达式分别为：
Tk＝Bk+Kk
Figure GDA00038967694700000710
步骤7-2、计算最小可检测偏差λa
最小可检测偏差λa通过求解以下非线性方程获得：
Figure GDA00038967694700000711
其中：Td表示进行GNSS故障检测时的检验统计量；
Figure GDA00038967694700000712
表示表示自由度为n-4，非中心化参数为λa的非中心化卡方分布概率累积分布函数在Td处的函数值；
Figure GDA00038967694700000713
表示
Figure GDA00038967694700000714
的先验发生概率；
Td的计算公式为：
Figure GDA0003896769470000081
Figure GDA0003896769470000082
的计算公式为：
Figure GDA0003896769470000083
步骤7-3、计算
Figure GDA0003896769470000084
的双侧分位数
Figure GDA0003896769470000085
其表达式为：
Figure GDA0003896769470000086
步骤7-4、计算
Figure GDA0003896769470000087
故障模式保护级，其表达式为：
Figure GDA0003896769470000088
进一步地，步骤8包括步骤8-1、计算GsI的双侧分位数
Figure GDA0003896769470000089
其表达式为：
Figure GDA00038967694700000810
其中：P(GSI)＝nPinsPsat(1-Psat)n-1表达GsI的先验发生概率；
步骤8-2、计算GsI故障模式的保护级，其表达式为：
Figure GDA00038967694700000811
进一步地，步骤9计算最终保护级，其表达式为：
Figure GDA00038967694700000812
垂向保护级VPL为：
VPL＝PLu,final
水平保护级HPL为：
Figure GDA00038967694700000813
其中，e,n,u分别表示东向，北向和垂向在卡尔曼滤波器状态向量中的索引。
根据本申请的又一方面，提供了一种GNSS/INS组合导航系统保护级反演方法，包括以下步骤：
步骤1、参数初始化，包括GNSS卫星先验故障概率Psat、INS先验故障概率Pins、完好性风险需求值Ireq以及连续性风险需求值Creq；
步骤2、分别针对GNSS和INS进行故障检测和排除；
步骤3、利用故障检测和排除后的GNSS和INS测量信息进行组合导航解算，解算方法为扩展卡尔曼滤波算法；
步骤4、针对四种故障模式
Figure GDA0003896769470000091
和GsI进行完好性风险需求分配，其中分配给GNSS无故障且INS无故障
Figure GDA0003896769470000092
的完好性风险需求为
Figure GDA0003896769470000093
分配给 GNSS无故障且INS故障模式
Figure GDA0003896769470000094
的完好性风险需求为
Figure GDA0003896769470000095
分配给GNSS单星故障且INS无故障
Figure GDA0003896769470000096
的完好性风险需求为
Figure GDA0003896769470000097
分配给GNSS单星故障且 INS故障GsI的完好性风险需求为Ir(GsI)
四种故障模式下的完好性风险需求应满足：
Figure GDA0003896769470000098
步骤5、计算
Figure GDA0003896769470000099
故障模式保护级：
步骤5-1、计算
Figure GDA00038967694700000910
的双侧分位数
Figure GDA00038967694700000911
其表达式为：
Figure GDA00038967694700000912
其中：
Figure GDA00038967694700000913
表达
Figure GDA00038967694700000914
的先验发生概率，n为当前时刻参与滤波解算的可见卫星数；Q-1(·)表示标准正态分布概率累积分布函数的逆函数。
步骤5-2、计算
Figure GDA00038967694700000915
故障模式的保护级，其表达式为：
Figure GDA00038967694700000916
其中σq为q方向滤波估计误差标准差，其值为滤波估计协方差矩阵Pk第q 个对角线元素的平方根。
步骤6、计算
Figure GDA00038967694700000917
故障模式保护级：
步骤6-1、计算残差传播矩阵Bk，其表达式为：
Figure GDA0003896769470000101
其中，Ik表示单位矩阵，Kk表示滤波器的增益矩阵，Hk表示滤波器的观测矩阵。
步骤6-2、计算
Figure GDA0003896769470000102
的双侧分位数
Figure GDA0003896769470000103
其表达式为：
Figure GDA0003896769470000104
其中，
Figure GDA0003896769470000105
表达
Figure GDA0003896769470000106
的先验发生概率。
步骤6-3、计算计算
Figure GDA0003896769470000107
故障模式的保护级，其表达式为：
Figure GDA0003896769470000108
其中：γk为滤波器的新息矢量；(·)q表示对应矢量的第q个分量；σBKq表示矩阵(Bk+Kk)Rk(Bk+Kk)T第q行对角线元素的平方根，Rk表示测量噪声协方差矩阵。
步骤7、计算
Figure GDA00038967694700001010
故障模式保护级：
步骤7-1、计算每颗卫星对应的斜率Slopei,q，其表达式为：
Figure GDA00038967694700001011
其中：tq,i表示矩阵Tk中第q行第i列的元素；si,i为矩阵S中第i行第i列的元素；σi为第i颗可见星观测噪声标准差。矩阵Tk和S的表达式分别为：
Tk＝Bk+Kk
Figure GDA00038967694700001012
步骤7-2、计算最小可检测偏差λa
最小可检测偏差λa通过求解以下非线性方程获得：
Figure GDA00038967694700001013
其中：Td表示进行GNSS故障检测时的检验统计量；
Figure GDA0003896769470000111
表示表示自由度为n-4，非中心化参数为λa的非中心化卡方分布概率累积分布函数在Td处的函数值；
Figure GDA0003896769470000112
表示
Figure GDA0003896769470000113
的先验发生概率。
Td的计算公式为：
Figure GDA0003896769470000114
Figure GDA0003896769470000115
的计算公式为：
Figure GDA0003896769470000116
步骤7-3、计算
Figure GDA0003896769470000117
的双侧分位数
Figure GDA0003896769470000118
其表达式为：
Figure GDA0003896769470000119
步骤7-4、计算
Figure GDA00038967694700001110
故障模式保护级，其表达式为：
Figure GDA00038967694700001111
步骤8、计算GsI故障模式保护级：
步骤8-1、计算GsI的双侧分位数
Figure GDA00038967694700001112
其表达式为：
Figure GDA00038967694700001113
其中：P(GSI)＝nPinsPsat(1-Psat)n-1表达GsI的先验发生概率。
步骤8-2、计算GsI故障模式的保护级，其表达式为：
Figure GDA00038967694700001114
步骤9、计算最终保护级，其表达式为：
Figure GDA00038967694700001115
垂向保护级VPL为：
VPL＝PLu,final
水平保护级HPL为：
Figure GDA0003896769470000121
其中，e,n,u分别表示东向，北向和垂向在卡尔曼滤波器状态向量中的索引。
积极有益效果：本发明将INS故障纳入完好性监测中，基于GNSS无故障且INS无故障、GNSS无故障且INS故障、GNSS单星故障且INS无故障和GNSS 单星故障且INS故障四种故障模式在扩展卡尔曼滤波器前提下的故障传播特性，分别设计了相应的保护级反演方案，并以四种故障模式中各保护级绝对值的最大值作为GNSS/INS组合导航系统的最终保护级，从而实现了针对GNSS/INS 组合导航系统并兼顾GNSS卫星和INS故障的保护级反演方法，有助于实现更为严谨和可靠的完好性监测。
附图说明
图1为本发明优选实施例中GNSS/INS组合导航系统保护级反演方法流程图；
图2为本发明优选实施例中INS和GNSS无故障条件下组合导航定位误差和
Figure GDA0003896769470000122
故障模式保护级曲线；
图3为本发明优选实施例中INS故障GNSS无故障条件下组合导航定位误差和
Figure GDA0003896769470000123
故障模式保护级曲线；
图4为本发明优选实施例中GNSS单星故障INS无故障条件下组合导航定位误差和
Figure GDA0003896769470000124
故障模式保护级曲线；
图5显为本发明优选实施例中GNSS单星故障且INS故障条件下组合导航定位误差和GsI故障模式保护级曲线；
图6为本发明优选实施例中组合导航定位误差和最终保护级曲线。
具体实施方式
下面通过附图和实施例对本发明进一步详细说明，通过这些说明，本发明的特点和优点将变得更为清楚明确。然而，本发明并不受限于以下所公开的示范性实施例；可以通过不同形式来对其加以实现。说明书的实质仅仅是帮助相关领域技术人员综合理解本发明的具体细节。
实施例1
本实施例针对使用GNSS/INS组合导航系统作为导航定位手段的无人机进行保护级反演。无人机的轨迹以及导航信息通过仿真生成。仿真过程中，无人机的飞行时间为900s，GNSS星座为GPS，INS的输出频率为40Hz，GPS的输出频率为1Hz。仿真过程中的传感器误差参数如表1所示：
表1传感器误差参数
Figure GDA0003896769470000131
如图1所示，本发明的一种GNSS/INS组合导航自主完好性监测方法。包括以下步骤：
步骤1、参数初始化：
本实施例中各参数分别初始化为：
GNSS卫星先验故障概率Psat＝1×10-5；
INS先验故障概率Pins＝1×10-3；
完好性风险需求值Ireq＝1×10-7
连续性风险需求值Creq＝4×10-6；
步骤2、故障检测和排除；
本实施例中针对GNSS卫星使用最小二乘残差接收机自主完好性监测算法进行故障检测与排除，针对INS使用基于滤波估计误差的故障检测方法进行故障检测。
步骤3、组合导航解算
利用故障检测和排除后的GNSS和INS测量信息进行组合导航解算，解算方法为扩展卡尔曼滤波算法；
步骤4、完好性风险需求分配
针对四种故障模式
Figure GDA0003896769470000141
和GsI进行完好性风险需求分配，其中分配给GNSS无故障且INS无故障
Figure GDA0003896769470000142
的完好性风险需求为
Figure GDA0003896769470000143
分配给GNSS无故障且INS故障模式
Figure GDA0003896769470000144
的完好性风险需求为
Figure GDA0003896769470000145
分配给GNSS单星故障且 INS无故障
Figure GDA0003896769470000146
的完好性风险需求为
Figure GDA0003896769470000147
分配给GNSS单星故障且INS故障GsI的完好性风险需求为Ir(GsI)
四种故障模式下的完好性风险需求应满足：
Figure GDA0003896769470000148
本实施例中，
Figure GDA0003896769470000149
步骤5、计算
Figure GDA00038967694700001410
故障模式保护级：
步骤5-1、计算
Figure GDA00038967694700001411
的双侧分位数
Figure GDA00038967694700001412
其表达式为：
Figure GDA00038967694700001413
其中：
Figure GDA00038967694700001414
表达
Figure GDA00038967694700001415
的先验发生概率，n为当前时刻参与滤波解算的可见卫星数；Q-1(·)表示标准正态分布概率累积分布函数的逆函数。
步骤5-2、计算
Figure GDA00038967694700001416
故障模式的保护级，其表达式为：
Figure GDA00038967694700001417
其中σq为q方向滤波估计误差标准差，其值为滤波估计协方差矩阵Pk第q 个对角线元素的平方根。
图2给出了INS和GNSS无故障条件下组合导航定位误差和
Figure GDA0003896769470000151
故障模式保护级曲线，可见本方法提供的保护级可以有效包络INS和GNSS无故障条件下水平和垂向的定位误差。
步骤6、计算
Figure GDA0003896769470000152
故障模式保护级：
步骤6-1、计算残差传播矩阵Bk，其表达式为：
Figure GDA0003896769470000153
其中，Ik表示单位矩阵，Kk表示滤波器的增益矩阵，Hk表示滤波器的观测矩阵。
步骤6-2、计算
Figure GDA0003896769470000154
的双侧分位数
Figure GDA0003896769470000155
其表达式为：
Figure GDA0003896769470000156
其中，
Figure GDA0003896769470000157
表达
Figure GDA0003896769470000158
的先验发生概率。
步骤6-3、计算计算
Figure GDA0003896769470000159
故障模式的保护级，其表达式为：
Figure GDA00038967694700001510
其中：γk为滤波器的新息矢量；(·)q表示对应矢量的第q个分量；σBKq表示矩阵(Bk+Kk)Rk(Bk+Kk)T第q行对角线元素的平方根，Rk测量噪声协方差矩阵。
图3给出了INS故障GNSS无故障条件下组合导航定位误差和
Figure GDA00038967694700001511
故障模式保护级曲线，其中INS故障条件为：在200-250s历元为垂向加速度计加入 20×10-4g的常值故障偏差。可见本方法提供的保护级可以有效包络INS故障 GNSS无故障条件下水平和垂向的定位误差。
步骤7、计算
Figure GDA00038967694700001512
故障模式保护级：
步骤7-1、计算每颗卫星对应的斜率Slopei,q，其表达式为：
Figure GDA0003896769470000161
其中：tq,i表示矩阵Tk中第q行第i列的元素；si,i为矩阵S中第i行第i列的元素；σi为第i颗可见星观测噪声标准差。矩阵Tk和S的表达式分别为：
Tk＝Bk+Kk
Figure GDA0003896769470000162
步骤7-2、计算最小可检测偏差λa
最小可检测偏差λa通过求解以下非线性方程获得：
Figure GDA0003896769470000163
其中：Td表示进行GNSS故障检测时的检验统计量；
Figure GDA0003896769470000164
表示表示自由度为n-4，非中心化参数为λa的非中心化卡方分布概率累积分布函数在Td处的函数值；
Figure GDA0003896769470000165
表示
Figure GDA0003896769470000166
的先验发生概率。
Td的计算公式为：
Figure GDA0003896769470000167
Figure GDA0003896769470000168
的计算公式为：
Figure GDA0003896769470000169
步骤7-3、计算
Figure GDA00038967694700001610
的双侧分位数
Figure GDA00038967694700001611
其表达式为：
Figure GDA00038967694700001612
步骤7-4、计算
Figure GDA00038967694700001613
故障模式保护级，其表达式为：
Figure GDA00038967694700001614
图4给出了GNSS单星故障INS无故障条件下组合导航定位误差和
Figure GDA00038967694700001615
故障模式保护级曲线，其中GNSS单星故障条件为：在200s-250s历元为G19分别注入10σ0的伪距故障偏差。可见本方法提供的保护级可以有效包络GNSS单星故障且INS无故障条件下水平和垂向的定位误差。
步骤8、计算GsI故障模式保护级：
步骤8-1、计算GsI的双侧分位数
Figure GDA0003896769470000171
其表达式为：
Figure GDA0003896769470000172
其中：P(GSI)＝nPinsPsat(1-Psat)n-1表达GsI的先验发生概率。
步骤8-2、计算GsI故障模式的保护级，其表达式为：
Figure GDA0003896769470000173
图5给出了GNSS单星故障且INS故障条件下组合导航定位误差和GsI故障模式保护级曲线，其中GNSS单星故障条件为：在200s-250s历元为G19分别注入20σ0的伪距故障偏差；INS故障条件为：在200-250s历元为垂向加速度计加入20×10-4g的常值故障偏差。可见本方法提供的保护级可以有效包络GNSS单星故障且INS故障条件下水平和垂向的定位误差。
步骤9、计算最终保护级，其表达式为：
Figure GDA0003896769470000174
垂向保护级VPL为：
VPL＝PLu,final
水平保护级HPL为：
Figure GDA0003896769470000175
其中，e,n,u分别表示东向，北向和垂向在卡尔曼滤波器状态向量中的索引。
图5给出了组合导航定位误差和最终保护级曲线，其中GNSS单星故障条件为：在500s-550s历元为G19号卫星的伪距注入10σ0故障偏差；INS故障条件为：在200-250s历元为垂向加速度计加入20×10-4g的常值故障偏差。可见本方法提供的保护级可以有效包络GNSS/INS组合导航系统水平和垂向的定位误差。
综合以上优选实施例中的GNSS/INS组合导航系统保护级反演过程，本发明所提供的一种GNSS/INS组合导航系统保护级反演方法，能够综合考虑GNSS 和INS的故障风险，针对每种故障模式反演出能够有效包络定位误差的保护级，并形成组合导航系统的最终保护级，从而实现更为严谨和可靠的完好性监测。
以上所述的具体描述，对发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (2)

1.一种GNSS/INS组合导航系统保护级反演方法，其特征在于，包括以下步骤：
步骤1、参数初始化，包括GNSS卫星先验故障概率Psat、INS先验故障概率Pins、完好性风险需求值Ireq以及连续性风险需求值Creq；
步骤2、分别针对GNSS和INS进行故障检测和排除；
步骤3、利用故障检测和排除后的GNSS和INS测量信息进行组合导航解算，解算方法为扩展卡尔曼滤波算法；
步骤4、针对四种故障模式
Figure FDA0003896769460000011
和GsI进行完好性风险需求分配，其中分配给GNSS无故障且INS无故障
Figure FDA0003896769460000012
的完好性风险需求为
Figure FDA0003896769460000013
分配给GNSS无故障且INS故障模式
Figure FDA0003896769460000014
的完好性风险需求为
Figure FDA0003896769460000015
分配给GNSS单星故障且INS无故障
Figure FDA0003896769460000016
的完好性风险需求为
Figure FDA0003896769460000017
分配给GNSS单星故障且INS故障GsI的完好性风险需求为Ir(GsI)；
步骤5、计算
Figure FDA0003896769460000018
故障模式保护级：
步骤6、计算
Figure FDA0003896769460000019
故障模式保护级：
步骤7、计算
Figure FDA00038967694600000110
故障模式保护级：
步骤8、计算GsI故障模式保护级：
步骤9、计算最终保护级；
所述的步骤5包括步骤5-1、计算
Figure FDA00038967694600000111
的双侧分位数
Figure FDA00038967694600000112
其表达式为：
Figure FDA00038967694600000113
其中：
Figure FDA00038967694600000114
表达
Figure FDA00038967694600000115
的先验发生概率，n为当前时刻参与滤波解算的可见卫星数；Q-1(·)表示标准正态分布概率累积分布函数的逆函数；
步骤5-2、计算
Figure FDA00038967694600000116
故障模式的保护级，其表达式为：
Figure FDA0003896769460000021
其中σq为q方向滤波估计误差标准差，其值为滤波估计协方差矩阵Pk第q个对角线元素的平方根；
所述的步骤6包括步骤6-1、计算残差传播矩阵Bk，其表达式为：
Figure FDA0003896769460000022
其中，Ik表示单位矩阵，Kk表示滤波器的增益矩阵，Hk表示滤波器的观测矩阵；
步骤6-2、计算
Figure FDA0003896769460000023
的双侧分位数
Figure FDA0003896769460000024
其表达式为：
Figure FDA0003896769460000025
其中，
Figure FDA0003896769460000026
表达
Figure DEST_PATH_FDA0003896769460000023
的先验发生概率；
步骤6-3、计算
Figure FDA0003896769460000028
故障模式的保护级，其表达式为：
Figure FDA0003896769460000029
其中：γk为滤波器的新息矢量；(·)q表示对应矢量的第q个分量；σBKq表示矩阵(Bk+Kk)Rk(Bk+Kk)T第q行对角线元素的平方根，Rk为测量噪声协方差矩阵；
所述的步骤7包括步骤7-1、计算每颗卫星对应的斜率Slopei,q，其表达式为：
Figure FDA00038967694600000210
其中：tq,i表示矩阵Tk中第q行第i列的元素；si,i为矩阵S中第i行第i列的元素；σi为第i颗可见星观测噪声标准差；矩阵Tk和S的表达式分别为：
Tk＝Bk+Kk
Figure FDA00038967694600000211
步骤7-2、计算最小可检测偏差λa
最小可检测偏差λa通过求解以下非线性方程获得：
Figure FDA0003896769460000031
其中：Td表示进行GNSS故障检测时的检验统计量；
Figure FDA0003896769460000032
表示自由度为n-4，非中心化参数为λa的非中心化卡方分布概率累积分布函数在Td处的函数值；
Figure FDA0003896769460000033
表示
Figure FDA0003896769460000034
的先验发生概率；
Td的计算公式为：
Figure FDA0003896769460000035
Figure FDA0003896769460000036
的计算公式为：
Figure FDA0003896769460000037
步骤7-3、计算
Figure FDA0003896769460000038
的双侧分位数
Figure DEST_PATH_IMAGE001
，其表达式为：
Figure FDA00038967694600000310
步骤7-4、计算
Figure FDA00038967694600000311
故障模式保护级，其表达式为：
Figure FDA00038967694600000312
所述的步骤8包括步骤8-1、计算GsI的双侧分位数
Figure FDA00038967694600000313
其表达式为：
Figure FDA00038967694600000314
其中：P(GSI)＝nPinsPsat(1-Psat)n-1表达GsI的先验发生概率；
步骤8-2、计算GsI故障模式的保护级，其表达式为：
Figure FDA00038967694600000315
所述的步骤9计算最终保护级，其表达式为：
Figure FDA00038967694600000316
垂向保护级VPL为：
VPL＝PLu,final
水平保护级HPL为：
Figure FDA0003896769460000041
其中，e,n,u分别表示东向，北向和垂向在卡尔曼滤波器状态向量中的索引。
2.根据权利要求1 所述的一种GNSS/INS组合导航系统保护级反演方法，其特征在于，所述的步骤4中的四种故障模式下的完好性风险需求应满足：
Figure FDA0003896769460000042
CN202211048910.8A 2022-08-30 2022-08-30 一种gnss/ins组合导航系统保护级反演方法 Active CN115112126B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211048910.8A CN115112126B (zh) 2022-08-30 2022-08-30 一种gnss/ins组合导航系统保护级反演方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211048910.8A CN115112126B (zh) 2022-08-30 2022-08-30 一种gnss/ins组合导航系统保护级反演方法

Publications (2)

Publication Number Publication Date
CN115112126A CN115112126A (zh) 2022-09-27
CN115112126B true CN115112126B (zh) 2022-11-18

Family

ID=83336238

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211048910.8A Active CN115112126B (zh) 2022-08-30 2022-08-30 一种gnss/ins组合导航系统保护级反演方法

Country Status (1)

Country Link
CN (1) CN115112126B (zh)

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109900300A (zh) * 2019-03-27 2019-06-18 北京航空航天大学 一种用于无人机的组合导航完好性监测系统
CN110133689A (zh) * 2019-05-24 2019-08-16 中国科学院国家授时中心 自适应用户自主完好性监测方法
CN114235007A (zh) * 2021-12-02 2022-03-25 北京航空航天大学 一种apnt服务的定位和完好性监测方法及系统
CN114545454A (zh) * 2022-02-15 2022-05-27 南京航空航天大学 一种面向自动驾驶的融合导航系统完好性监测方法
CN114721017A (zh) * 2022-03-04 2022-07-08 北京理工大学 一种gnss/ins组合导航自主完好性监测方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
FR2866171B1 (fr) * 2004-02-06 2006-06-30 Thales Sa Procede automatique de transmission des alertes de surveillance bord vers le sol

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109900300A (zh) * 2019-03-27 2019-06-18 北京航空航天大学 一种用于无人机的组合导航完好性监测系统
CN110133689A (zh) * 2019-05-24 2019-08-16 中国科学院国家授时中心 自适应用户自主完好性监测方法
CN114235007A (zh) * 2021-12-02 2022-03-25 北京航空航天大学 一种apnt服务的定位和完好性监测方法及系统
CN114545454A (zh) * 2022-02-15 2022-05-27 南京航空航天大学 一种面向自动驾驶的融合导航系统完好性监测方法
CN114721017A (zh) * 2022-03-04 2022-07-08 北京理工大学 一种gnss/ins组合导航自主完好性监测方法

Also Published As

Publication number Publication date
CN115112126A (zh) 2022-09-27

Similar Documents

Publication Publication Date Title
US10018729B2 (en) Selected aspects of advanced receiver autonomous integrity monitoring application to kalman filter based navigation filter
US8902105B2 (en) Method and apparatus for determining an integrity indicating parameter indicating the integrity of positioning information determined in a gobal positioning system
US9146322B2 (en) Hybrid system and device for calculating a position and for monitoring its integrity
CN110007317B (zh) 一种选星优化的高级接收机自主完好性监测方法
CN109100748B (zh) 一种基于低轨星座的导航完好性监测系统及方法
KR101206364B1 (ko) 다중 기준국 환경에서 이상위성의 판단방법 및 이를 이용한 판단장치
EP2081043A2 (en) Navigation system with apparatus for detecting accuracy failures
CN110879407B (zh) 一种基于完好性风险模型的卫星导航观测量新息检测方法
WO2021202004A2 (en) System and method for reconverging gnss position estimates
US20090182495A1 (en) Navigation system with apparatus for detecting accuracy failures
Dovis et al. Recent advancement on the use of global navigation satellite system-based positioning for intelligent transport systems [guest editorial]
US6298316B1 (en) Failure detection system
Han et al. GNSS/IMU tightly coupled scheme with weighting and FDE for rail applications
CN108507590B (zh) 定速评估方法及系统、车载终端
CN112198533B (zh) 一种多假设下的地基增强系统完好性评估系统及方法
CN115235463B (zh) 一种gnss/ins组合导航系统完好性风险需求分配方法
CN112179347B (zh) 一种基于光谱红移误差观测方程的组合导航方法
CN115112126B (zh) 一种gnss/ins组合导航系统保护级反演方法
EP3693761A1 (en) Alternate uncertainty limits in the presence of a detected satellite fault
KR101040053B1 (ko) 의사거리 측정잡음 감시기반 위성전파항법시스템 위성시계 고장검출 및 고장위성 식별 방법
KR101040054B1 (ko) 수신기 시계오차 감시기반 위성전파항법시스템 위성시계 고장검출 및 고장위성 식별 방법
CN115291253A (zh) 一种基于残差检测的车辆定位完好性监测方法及系统
KR101290085B1 (ko) 다중 기준국 환경에서 대류층 지연 변칙현상 모니터링 방법 및 이를 이용한 시스템
CN116859417B (zh) 用于北斗ppp-rtk/mems的完好性监测方法
Lee et al. Optimal continuity allocation for a tightly-coupled KF-based GNSS/IMU navigation system with redundant IMUs

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
GR01 Patent grant
GR01 Patent grant