CN115408652A - A method for determining the protection level of EKF integrated navigation based on truncation error estimation - Google Patents

A method for determining the protection level of EKF integrated navigation based on truncation error estimation Download PDF

Info

Publication number
CN115408652A
CN115408652A CN202211341296.4A CN202211341296A CN115408652A CN 115408652 A CN115408652 A CN 115408652A CN 202211341296 A CN202211341296 A CN 202211341296A CN 115408652 A CN115408652 A CN 115408652A
Authority
CN
China
Prior art keywords
error
truncation
protection level
truncation error
pseudo
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.)
Granted
Application number
CN202211341296.4A
Other languages
Chinese (zh)
Other versions
CN115408652B (en
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.)
Beihang University
Original Assignee
Beihang University
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 Beihang University filed Critical Beihang University
Priority to CN202211341296.4A priority Critical patent/CN115408652B/en
Publication of CN115408652A publication Critical patent/CN115408652A/en
Application granted granted Critical
Publication of CN115408652B publication Critical patent/CN115408652B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • 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/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • 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/183—Compensation of inertial measurements, e.g. for temperature effects
    • 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/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/03—Cooperating elements; Interaction or communication between different cooperating elements or between cooperating elements and receivers
    • G01S19/07—Cooperating elements; Interaction or communication between different cooperating elements or between cooperating elements and receivers providing data for correcting measured positioning data, e.g. DGPS [differential GPS] or ionosphere corrections
    • 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/393—Trajectory determination or predictive tracking, e.g. Kalman filtering
    • 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/396—Determining accuracy or reliability of position or pseudorange measurements
    • 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/45—Determining position by combining measurements of signals from the satellite radio beacon positioning system with a supplementary measurement
    • G01S19/47—Determining position by combining measurements of signals from the satellite radio beacon positioning system with a supplementary measurement the supplementary measurement being an inertial measurement, e.g. tightly coupled inertial
    • 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

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Mathematical Optimization (AREA)
  • Computational Mathematics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Data Mining & Analysis (AREA)
  • Mathematical Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • General Engineering & Computer Science (AREA)
  • Software Systems (AREA)
  • Databases & Information Systems (AREA)
  • Algebra (AREA)
  • Operations Research (AREA)
  • Computing Systems (AREA)
  • Automation & Control Theory (AREA)
  • Position Fixing By Use Of Radio Waves (AREA)

Abstract

The invention relates to an EKF (extended Kalman filter) combined navigation protection level determination method based on truncation error estimation; establishing a GNSS/INS tight combination navigation system state space model considering truncation errors, and estimating truncation error values and covariance thereof; modifying the model with the estimated truncation error and its covariance; calculating the statistical characteristic of the corrected positioning error by considering the residual amount of the truncation error to obtain a protection level expression comprising the truncation error; acquiring actual pseudo-range, pseudo-range rate observed quantity samples and corresponding residual truncation error sample data; respectively enveloping pseudo range, pseudo range rate observation error and residual truncation error sample data based on an extreme value theory to obtain corresponding amplification factors; and substituting the amplification factor into the protection level expression to obtain the protection level of the GNSS/INS tight combination navigation system. The invention combines more accurate protection level and more accurate error tail probability distribution, and ensures the reliability of the protection level.

Description

一种基于截断误差估计的EKF组合导航保护级确定方法A method for determining the protection level of EKF integrated navigation based on truncation error estimation

技术领域technical field

本发明属于组合导航技术领域，具体涉及一种基于截断误差估计的EKF组合导航保护级确定方法。The invention belongs to the technical field of integrated navigation, and in particular relates to an EKF integrated navigation protection level determination method based on truncation error estimation.

背景技术Background technique

单导航系统已无法满足无人机在复杂环境中精准、可信、可靠飞行的需求，GNSS/INS组合导航具有GNSS全天候长期稳定导航和INS短时精度高的优点，可作为复杂环境中无人机飞行的导航手段。为保证飞行安全，GNSS/INS组合导航必须满足给定的完好性需求，可实时计算指定风险概率下定位误差的置信上限，即保护级，当保护级超过告警限时，说明此时的导航系统不可用，需要及时向用户报警。A single navigation system can no longer meet the precise, credible, and reliable flight requirements of UAVs in complex environments. GNSS/INS integrated navigation has the advantages of GNSS all-weather long-term stable navigation and INS short-term high accuracy, and can be used as an unmanned aerial vehicle in complex environments. means of navigation for aircraft. In order to ensure flight safety, GNSS/INS integrated navigation must meet the given integrity requirements, and can calculate the upper confidence limit of positioning error under the specified risk probability in real time, that is, the protection level. When the protection level exceeds the warning limit, it means that the navigation system at this time cannot It is necessary to notify the user in time.

采用EKF（ExtendedKalmanFilter，扩展卡尔曼滤波器）解算非线性问题时，需要利用泰勒级数对非线性系统进行线性化处理，但线性化需要舍弃泰勒展开式中高阶分量。GNSS/INS组合导航系统完好性监测的核心是对各因素导致的定位误差的统计特性进行估计，从而计算保护级。在将EKF方法运用于GNSS/INS组合导航时，人为产生了舍去泰勒展开式中高阶项的计算误差，称为截断误差。When using EKF (Extended Kalman Filter, Extended Kalman Filter) to solve nonlinear problems, it is necessary to use Taylor series to linearize the nonlinear system, but the linearization needs to discard the high-order components in the Taylor expansion. The core of GNSS/INS integrated navigation system integrity monitoring is to estimate the statistical characteristics of positioning errors caused by various factors, so as to calculate the protection level. When the EKF method is applied to GNSS/INS integrated navigation, the calculation error of discarding the high-order items in the Taylor expansion is artificially generated, which is called the truncation error.

在系统非线性较强的时候，截断误差使计算出的定位误差协方差不准确，导致保护级计算不够准确，无法满足完好性需求。因此需在滤波过程中把截断误差考虑进来。When the nonlinearity of the system is strong, the truncation error makes the calculated positioning error covariance inaccurate, which leads to the inaccurate calculation of the protection level and cannot meet the integrity requirements. Therefore, the truncation error needs to be taken into account in the filtering process.

另外，目前组合导航保护级的相关研究中，均假设误差服从零均值高斯分布，而实际中如伪距、伪距率误差是非高斯非零均值分布的，导致计算的保护级不能包络实际误差尾部，可能会造成潜在的完好性风险。In addition, in the current research on the protection level of integrated navigation, it is assumed that the error obeys the zero-mean Gaussian distribution, but in practice, such as pseudo-range and pseudo-range rate errors are non-Gaussian and non-zero-mean distribution, resulting in the calculated protection level cannot envelop the actual error Tail, may pose a potential integrity risk.

发明内容Contents of the invention

鉴于上述的分析，本发明旨在公开了一种基于截断误差估计的EKF组合导航保护级确定方法，用于解决EKF线性化过程中产生的截断误差使保护级计算不够准确的问题。In view of the above analysis, the present invention aims to disclose a method for determining the protection level of EKF integrated navigation based on truncation error estimation, which is used to solve the problem that the truncation error generated in the EKF linearization process makes the protection level calculation not accurate enough.

本发明公开了一种基于截断误差估计的EKF组合导航保护级计算方法，包括以下步骤：The invention discloses an EKF integrated navigation protection level calculation method based on truncation error estimation, comprising the following steps:

建立考虑截断误差的GNSS/INS紧组合导航系统状态空间模型，估计截断误差值和其协方差；Establish the state space model of GNSS/INS compact integrated navigation system considering the truncation error, and estimate the truncation error value and its covariance;

用估计的截断误差和其协方差修正所述模型；并考虑截断误差残余量，计算修正后的定位误差统计特性，得到包括截断误差的保护级表达式；Using the estimated truncation error and its covariance to modify the model; and considering the residual amount of truncation error, calculating the statistical characteristics of the corrected positioning error to obtain a protection level expression including the truncation error;

获取实际的伪距、伪距率观测量样本以及对应的残余截断误差样本数据；Obtain the actual pseudorange, pseudorange rate observation samples and corresponding residual truncation error sample data;

基于极值理论分别包络伪距、伪距率观测误差和残余截断误差样本数据，得到对应的放大因子；Based on the extreme value theory, the sample data of pseudorange, pseudorange rate observation error and residual truncation error are respectively enveloped to obtain the corresponding amplification factor;

将放大因子代入所述保护级表达式得出GNSS/INS紧组合导航系统的保护级。The protection level of the GNSS/INS compact integrated navigation system is obtained by substituting the amplification factor into the protection level expression.

进一步地，所述考虑截断误差的GNSS/INS紧组合导航系统状态空间模型为：Further, the state space model of the GNSS/INS compact integrated navigation system considering the truncation error is:

Figure 849845DEST_PATH_IMAGE001
；
Figure 849845DEST_PATH_IMAGE001
;

其中，k为时刻，

Figure 100698DEST_PATH_IMAGE002
是n维的状态向量；
Figure 66249DEST_PATH_IMAGE003
是m维的量测向量；
Figure 402552DEST_PATH_IMAGE004
为n阶的状态一步转移矩阵，
Figure 721538DEST_PATH_IMAGE005
为m×n阶的量测矩阵；
Figure 967843DEST_PATH_IMAGE006
为截断误差向量；
Figure 714082DEST_PATH_IMAGE007
为n维的系统噪声向量，
Figure 272102DEST_PATH_IMAGE008
为m维的量测噪声向量。Among them, k is the moment,
Figure 100698DEST_PATH_IMAGE002
is an n -dimensional state vector;
Figure 66249DEST_PATH_IMAGE003
is the m -dimensional measurement vector;
Figure 402552DEST_PATH_IMAGE004
is the n -order state one-step transition matrix,
Figure 721538DEST_PATH_IMAGE005
is a measurement matrix of order m × n ;
Figure 967843DEST_PATH_IMAGE006
is the truncation error vector;
Figure 714082DEST_PATH_IMAGE007
is the n -dimensional system noise vector,
Figure 272102DEST_PATH_IMAGE008
is the m -dimensional measurement noise vector.

进一步地，GNSS/INS紧组合导航系统状态空间模型的状态向量为17维的误差状态向量：Furthermore, the state vector of the GNSS/INS compact integrated navigation system state space model is a 17-dimensional error state vector:

Figure 394779DEST_PATH_IMAGE009
Figure 394779DEST_PATH_IMAGE009

其中，

Figure 980744DEST_PATH_IMAGE010
是三维速度误差；
Figure 163463DEST_PATH_IMAGE011
是三维姿态角误差；
Figure 474359DEST_PATH_IMAGE012
是三维位置误差；
Figure 10514DEST_PATH_IMAGE013
是三维加速度计误差；
Figure 90465DEST_PATH_IMAGE014
是三维陀螺仪误差；
Figure 444086DEST_PATH_IMAGE015
为一维接收机始终偏差；
Figure 366911DEST_PATH_IMAGE016
为一维接收机时钟漂移。in,
Figure 980744DEST_PATH_IMAGE010
is the three-dimensional velocity error;
Figure 163463DEST_PATH_IMAGE011
is the three-dimensional attitude angle error;
Figure 474359DEST_PATH_IMAGE012
is the three-dimensional position error;
Figure 10514DEST_PATH_IMAGE013
is the three-dimensional accelerometer error;
Figure 90465DEST_PATH_IMAGE014
is the three-dimensional gyroscope error;
Figure 444086DEST_PATH_IMAGE015
is the constant bias of the one-dimensional receiver;
Figure 366911DEST_PATH_IMAGE016
is the one-dimensional receiver clock drift.

进一步地，GNSS/INS紧组合导航系统状态空间模型的量测量为GNSS观测伪距

Figure 831391DEST_PATH_IMAGE017
、伪距率
Figure 500270DEST_PATH_IMAGE018
与INS预测的伪距
Figure 165737DEST_PATH_IMAGE019
、伪距率
Figure 185646DEST_PATH_IMAGE020
之差；Furthermore, the quantity measurement of the state space model of the GNSS/INS compact integrated navigation system is the GNSS observation pseudorange
Figure 831391DEST_PATH_IMAGE017
, pseudorange rate
Figure 500270DEST_PATH_IMAGE018
Pseudoranges from INS predictions
Figure 165737DEST_PATH_IMAGE019
, pseudorange rate
Figure 185646DEST_PATH_IMAGE020
Difference;

量测噪声向量

Figure 188237DEST_PATH_IMAGE021
包括GNSS观测伪距和伪距率噪声向量。measurement noise vector
Figure 188237DEST_PATH_IMAGE021
Includes GNSS observed pseudoranges and pseudorange rate noise vectors.

进一步地，采用NPF方法，估计截断误差值和其协方差，具体包括：Further, the NPF method is used to estimate the truncation error value and its covariance, including:

1）获取当前时刻观测量与上一时刻观测量的估计的残差；1) Obtain the estimated residual of the observed quantity at the current moment and the observed quantity at the previous moment;

残差

Figure 193098DEST_PATH_IMAGE022
；residual
Figure 193098DEST_PATH_IMAGE022
;

其中，

Figure 154101DEST_PATH_IMAGE023
为k时刻第i颗星的伪距观测量，
Figure 864568DEST_PATH_IMAGE024
为k时刻第i颗星的伪距率观测量；
Figure 77374DEST_PATH_IMAGE025
和
Figure 455266DEST_PATH_IMAGE026
为k-1时刻第i颗星的观测量的估计值，等于函数
Figure 587170DEST_PATH_IMAGE027
在k-1时刻估计的定位解处泰勒展开后省略高阶项的值；i=1,…,p，p为可见星数目；in,
Figure 154101DEST_PATH_IMAGE023
is the pseudo-range observation of the i -th star at time k ,
Figure 864568DEST_PATH_IMAGE024
is the pseudorange rate observation of the i -th star at time k ;
Figure 77374DEST_PATH_IMAGE025
and
Figure 455266DEST_PATH_IMAGE026
is the estimated value of the observed quantity of the i -th star at time k -1, which is equal to the function
Figure 587170DEST_PATH_IMAGE027
The value of the high-order item is omitted after Taylor expansion at the estimated positioning solution at time k -1; i= 1,..., p , p is the number of visible stars;

2）利用残差估计截断误差的估计量

Figure 971884DEST_PATH_IMAGE028
和其协方差
Figure 785119DEST_PATH_IMAGE029
：2) Use the residual to estimate the estimator of the truncation error
Figure 971884DEST_PATH_IMAGE028
and its covariance
Figure 785119DEST_PATH_IMAGE029
:

估计量

Figure 814255DEST_PATH_IMAGE030
；Estimator
Figure 814255DEST_PATH_IMAGE030
;

协方差

Figure 992427DEST_PATH_IMAGE031
；Covariance
Figure 992427DEST_PATH_IMAGE031
;

其中，

Figure 739803DEST_PATH_IMAGE032
为半正定的权矩阵，
Figure 356729DEST_PATH_IMAGE033
为量测噪声矩阵。in,
Figure 739803DEST_PATH_IMAGE032
is a positive semi-definite weight matrix,
Figure 356729DEST_PATH_IMAGE033
is the measurement noise matrix.

进一步地，用截断误差的估计值和其协方差修正后的GNSS/INS紧组合导航系统状态空间模型为：Furthermore, the state space model of the GNSS/INS compact integrated navigation system corrected by the estimated value of the truncation error and its covariance is:

Figure 132049DEST_PATH_IMAGE034
；
Figure 132049DEST_PATH_IMAGE034
;

式中，

Figure 809018DEST_PATH_IMAGE035
为截断误差经NPF方法估计后截断误差残余量，
Figure 778111DEST_PATH_IMAGE036
。In the formula,
Figure 809018DEST_PATH_IMAGE035
is the truncation error residual after the truncation error is estimated by the NPF method,
Figure 778111DEST_PATH_IMAGE036
.

进一步地，包括截断误差的保护级表达式的确定过程包括：Further, the process of determining the protection level expression including the truncation error includes:

1）根据卡尔曼滤波的误差传播，确定出定位误差；1) According to the error propagation of Kalman filter, the positioning error is determined;

2）进行完好性的假设分配，得到H0假设下，GNSS/INS紧组合导航系统垂直定位误差表达式；2) Carry out the hypothesis assignment of integrity, and obtain the vertical positioning error expression of the GNSS/INS compact integrated navigation system under the H0 assumption;

3）根据垂直定位误差表达式，确定出垂直定位误差分布特性；3) According to the vertical positioning error expression, the distribution characteristics of the vertical positioning error are determined;

4）根据误差分布特性，得到H0假设下包括截断误差的保护级表达式。4) According to the error distribution characteristics, the protection level expression including the truncation error under the H0 assumption is obtained.

进一步地，H0假设下，所述垂直定位误差表达式：Further, under the H0 assumption, the vertical positioning error expression:

Figure 870832DEST_PATH_IMAGE037
；
Figure 870832DEST_PATH_IMAGE037
;

其中，

Figure 77823DEST_PATH_IMAGE038
；
Figure 722431DEST_PATH_IMAGE039
为卡尔曼增益；v为误差状态向量X中与高度信息所对应的行号；in,
Figure 77823DEST_PATH_IMAGE038
;
Figure 722431DEST_PATH_IMAGE039
is the Kalman gain; v is the line number corresponding to the height information in the error state vector X ;

Figure 303454DEST_PATH_IMAGE040
为
Figure 262182DEST_PATH_IMAGE041
矩阵中与高度信息对应的行中的元素，
Figure 385996DEST_PATH_IMAGE042
为
Figure 404768DEST_PATH_IMAGE043
矩阵中与高度信息对应的行中的元素，
Figure 489398DEST_PATH_IMAGE044
为
Figure 782976DEST_PATH_IMAGE045
矩阵中与高度信息对应的行中的元素；
Figure 303454DEST_PATH_IMAGE040
for
Figure 262182DEST_PATH_IMAGE041
the element in the row corresponding to the height information in the matrix,
Figure 385996DEST_PATH_IMAGE042
for
Figure 404768DEST_PATH_IMAGE043
the element in the row corresponding to the height information in the matrix,
Figure 489398DEST_PATH_IMAGE044
for
Figure 782976DEST_PATH_IMAGE045
elements in the row corresponding to the height information in the matrix;

垂直定位误差

Figure 698980DEST_PATH_IMAGE046
服从高斯分布特性；vertical positioning error
Figure 698980DEST_PATH_IMAGE046
Obey the characteristics of Gaussian distribution;

Figure 577068DEST_PATH_IMAGE047
；
Figure 577068DEST_PATH_IMAGE047
;

H0假设下，包括截断误差的保护级表达式：Under the H0 assumption, the protection level expression including the truncation error is:

Figure 8050DEST_PATH_IMAGE048
；
Figure 8050DEST_PATH_IMAGE048
;

其中

Figure 308581DEST_PATH_IMAGE049
为无故障漏检系数。in
Figure 308581DEST_PATH_IMAGE049
is the error-free missed detection coefficient.

进一步地，基于极值理论包络伪距样本数据得到对应的放大因子的方法，包括：Further, the method of obtaining the corresponding amplification factor based on the extreme value theory envelope pseudorange sample data includes:

1）将获取的实际伪距样本取绝对值后从小到大排列，并计算出样本的标准差；1) Arrange the obtained actual pseudorange samples in ascending order after taking the absolute value, and calculate the standard deviation of the samples;

2）用平均超出量函数图法选取出一个伪距样本作为阈值；2) Select a pseudo-range sample as the threshold using the average excess amount function graph method;

3）用Bootstrap法将伪距样本有放回地重采样B次，得到B组样本；3) Use the Bootstrap method to resample the pseudorange samples B times with replacement to obtain group B samples;

3）根据极值理论，对每组样本进行参数估计，确定出参数置信限值；3) According to the extreme value theory, estimate the parameters of each group of samples, and determine the parameter confidence limits;

4）将参数置信限值带入样本分布，根据组合导航伪距误差的概率分布的完好性确定出伪距误差放大因子。4) The parameter confidence limit is brought into the sample distribution, and the pseudo-range error amplification factor is determined according to the integrity of the probability distribution of the integrated navigation pseudo-range error.

进一步地，将放大因子代入所述保护级表达式得出GNSS/INS紧组合导航系统的保护级

Figure 16774DEST_PATH_IMAGE050
；其中，Further, the protection level of the GNSS/INS tight integrated navigation system is obtained by substituting the amplification factor into the protection level expression
Figure 16774DEST_PATH_IMAGE050
;in,

Figure 174086DEST_PATH_IMAGE051
；
Figure 174086DEST_PATH_IMAGE051
;

其中，

Figure 295626DEST_PATH_IMAGE052
为第i颗可见星的伪距误差放大因子，
Figure 321219DEST_PATH_IMAGE053
为第i颗可见星的伪距率误差放大因子，
Figure 211815DEST_PATH_IMAGE054
为第i颗可见星的伪距残余截断误差放大因子，
Figure 805607DEST_PATH_IMAGE055
为第i颗可见星的伪距率残余截断误差放大因子。in,
Figure 295626DEST_PATH_IMAGE052
is the pseudorange error amplification factor of the i -th visible star,
Figure 321219DEST_PATH_IMAGE053
is the pseudorange rate error amplification factor of the i -th visible star,
Figure 211815DEST_PATH_IMAGE054
is the pseudorange residual truncation error amplification factor of the i -th visible star,
Figure 805607DEST_PATH_IMAGE055
is the pseudorange rate residual truncation error amplification factor of the i -th visible star.

本发明可实现以下有益效果之一：The present invention can realize one of the following beneficial effects:

本发明的基于截断误差估计的EKF组合导航保护级确定方法，针对EKF线性化过程中产生的截断误差使保护级计算不够准确的问题，通过估计出的截断误差的大小和其协方差，减小了截断误差的不确定性，修正定位误差及其分布，从而得到更准确的保护级，使其满足给定的完好性风险；通过采用极值理论包络实际采样的伪距、伪距率等误差、截断残余误差分布，定位出更准确的误差尾部概率分布。结合更准确的保护级和更准确的误差尾部概率分布，保证了保护级的可靠性。The EKF integrated navigation protection level determination method based on truncation error estimation of the present invention aims at the problem that the truncation error generated in the EKF linearization process makes the calculation of the protection level inaccurate, through the estimated truncation error size and its covariance, reduce The uncertainty of the truncation error is corrected, and the positioning error and its distribution are corrected to obtain a more accurate protection level to meet the given integrity risk; by using the extreme value theory to envelop the actual sampling pseudo-range, pseudo-range rate, etc. Error, truncate the residual error distribution, and locate a more accurate error tail probability distribution. Combined with a more accurate protection level and a more accurate error tail probability distribution, the reliability of the protection level is guaranteed.

附图说明Description of drawings

附图仅用于示出具体实施例的目的，而并不认为是对本发明的限制，在整个附图中，相同的参考符号表示相同的部件；The accompanying drawings are only for the purpose of illustrating specific embodiments, and are not considered to limit the present invention. Throughout the accompanying drawings, the same reference symbols represent the same components;

图1为本发明实施例中的基于截断误差估计的EKF组合导航保护级确定方法流程图。Fig. 1 is a flowchart of a method for determining protection levels of EKF integrated navigation based on truncation error estimation in an embodiment of the present invention.

具体实施方式Detailed ways

下面结合附图来具体描述本发明的优选实施例，其中，附图构成本申请一部分，并与本发明的实施例一起用于阐释本发明的原理。Preferred embodiments of the present invention will be specifically described below in conjunction with the accompanying drawings, wherein the accompanying drawings constitute a part of the application and are used together with the embodiments of the present invention to explain the principles of the present invention.

本发明的一个实施例公开一种基于截断误差估计的EKF组合导航保护级确定方法，如图1所示，包括以下步骤：One embodiment of the present invention discloses a method for determining the protection level of EKF integrated navigation based on truncation error estimation, as shown in Figure 1, comprising the following steps:

步骤S1、建立考虑截断误差的GNSS/INS紧组合导航系统状态空间模型，估计截断误差值和其协方差；Step S1, establishing a GNSS/INS compact integrated navigation system state space model considering the truncation error, and estimating the truncation error value and its covariance;

步骤S2、用估计的截断误差和其协方差修正所述模型；并考虑截断误差残余量，计算修正后的定位误差统计特性，得到包括截断误差的保护级表达式；Step S2, using the estimated truncation error and its covariance to modify the model; and considering the truncation error residual, calculating the statistical characteristics of the corrected positioning error, and obtaining a protection level expression including the truncation error;

步骤S3、获取实际的伪距、伪距率观测量样本以及对应的残余截断误差样本数据；Step S3, obtaining actual pseudoranges, pseudorange rate observation samples and corresponding residual truncation error sample data;

所述残余截断误差样本数据包括伪距残余截断误差样本数据和伪距率残余截断误差样本数据；The residual truncation error sample data includes pseudorange residual truncation error sample data and pseudorange rate residual truncation error sample data;

步骤S4、基于极值理论分别包络伪距、伪距率观测误差和残余截断误差数据，得到对应的放大因子；Step S4, respectively enveloping the pseudorange, pseudorange rate observation error and residual truncation error data based on the extremum theory to obtain the corresponding amplification factor;

步骤S5、将放大因子代入所述保护级表达式得出GNSS/INS紧组合导航系统的保护级。Step S5, substituting the amplification factor into the protection level expression to obtain the protection level of the GNSS/INS compact integrated navigation system.

在步骤S1中，建立的考虑截断误差的GNSS/INS紧组合导航系统状态空间模型为：In step S1, the established state-space model of the GNSS/INS compact integrated navigation system considering the truncation error is:

Figure 820968DEST_PATH_IMAGE056
；
Figure 820968DEST_PATH_IMAGE056
;

其中，k为时刻，

Figure 728881DEST_PATH_IMAGE057
是n维的状态向量；
Figure 536300DEST_PATH_IMAGE003
是m维的量测向量；
Figure 238677DEST_PATH_IMAGE058
和
Figure 492066DEST_PATH_IMAGE059
是已知的系统结构参数，分别称为n阶的状态一步转移矩阵、m×n阶的量测矩阵；
Figure 734828DEST_PATH_IMAGE060
为截断误差向量；
Figure 600016DEST_PATH_IMAGE061
是n维的系统噪声向量，
Figure 410977DEST_PATH_IMAGE062
是m维的量测噪声向量，两者都是零均值的高斯白噪声向量序列（服从正态分布），且它们之间互不相关。Among them, k is the moment,
Figure 728881DEST_PATH_IMAGE057
is an n -dimensional state vector;
Figure 536300DEST_PATH_IMAGE003
is the m -dimensional measurement vector;
Figure 238677DEST_PATH_IMAGE058
and
Figure 492066DEST_PATH_IMAGE059
are the known system structure parameters, which are respectively called the n -order state one-step transition matrix and the m × n -order measurement matrix;
Figure 734828DEST_PATH_IMAGE060
is the truncation error vector;
Figure 600016DEST_PATH_IMAGE061
is the n -dimensional system noise vector,
Figure 410977DEST_PATH_IMAGE062
is the m -dimensional measurement noise vector, both of which are zero-mean Gaussian white noise vector sequences (subject to normal distribution), and they are not correlated with each other.

即：which is:

Figure 525564DEST_PATH_IMAGE063
。
Figure 525564DEST_PATH_IMAGE063
.

Figure 509700DEST_PATH_IMAGE064
为系统噪声矩阵；
Figure 416345DEST_PATH_IMAGE065
为量测噪声矩阵。
Figure 509700DEST_PATH_IMAGE064
is the system noise matrix;
Figure 416345DEST_PATH_IMAGE065
is the measurement noise matrix.

优选的，在本实施例中，GNSS/INS紧组合导航系统状态空间模型的状态向量为17维的误差状态向量：Preferably, in this embodiment, the state vector of the GNSS/INS compact integrated navigation system state space model is a 17-dimensional error state vector:

Figure 257262DEST_PATH_IMAGE066
Figure 257262DEST_PATH_IMAGE066

其中，in,

Figure 62407DEST_PATH_IMAGE067
是在ENU导航坐标系下的三维速度误差（东、北、高）；
Figure 62407DEST_PATH_IMAGE067
is the three-dimensional velocity error (east, north, height) in the ENU navigation coordinate system;

Figure 787918DEST_PATH_IMAGE068
是三维姿态角误差（俯仰角、横滚角、航向角）；
Figure 787918DEST_PATH_IMAGE068
is the three-dimensional attitude angle error (pitch angle, roll angle, heading angle);

Figure 362119DEST_PATH_IMAGE069
是在大地坐标系下的三维位置误差（纬度、经度、高度）；
Figure 362119DEST_PATH_IMAGE069
is the three-dimensional position error (latitude, longitude, height) in the geodetic coordinate system;

Figure 373937DEST_PATH_IMAGE070
是三维加速度计误差（含零偏）；
Figure 373937DEST_PATH_IMAGE070
is the three-dimensional accelerometer error (including zero bias);

Figure 89214DEST_PATH_IMAGE071
是三维陀螺仪误差（含零偏）；
Figure 89214DEST_PATH_IMAGE071
is the three-dimensional gyroscope error (including zero bias);

Figure 680733DEST_PATH_IMAGE072
为接收机始终偏差；
Figure 680733DEST_PATH_IMAGE072
is the receiver always bias;

Figure 171757DEST_PATH_IMAGE073
为接收机时钟漂移。
Figure 171757DEST_PATH_IMAGE073
is the receiver clock drift.

优选的，在本实施例中，GNSS/INS紧组合导航系统状态空间模型的量测量为GNSS观测伪距

Figure 557739DEST_PATH_IMAGE074
、伪距率
Figure 9580DEST_PATH_IMAGE075
与INS预测的伪距
Figure 201527DEST_PATH_IMAGE076
、伪距率
Figure 484741DEST_PATH_IMAGE077
之差；Preferably, in this embodiment, the quantity measurement of the state space model of the GNSS/INS compact integrated navigation system is the GNSS observation pseudorange
Figure 557739DEST_PATH_IMAGE074
, pseudorange rate
Figure 9580DEST_PATH_IMAGE075
Pseudoranges from INS predictions
Figure 201527DEST_PATH_IMAGE076
, pseudorange rate
Figure 484741DEST_PATH_IMAGE077
Difference;

量测噪声向量

Figure 228575DEST_PATH_IMAGE021
包括GNSS观测伪距和伪距率噪声向量。measurement noise vector
Figure 228575DEST_PATH_IMAGE021
Includes GNSS observed pseudoranges and pseudorange rate noise vectors.

在步骤S1中，采样NPF（非线性预测滤波）方法，估计截断误差的大小和其协方差。In step S1, the NPF (non-linear predictive filtering) method is sampled to estimate the magnitude of the truncation error and its covariance.

具体包括：Specifically include:

1）获取当前时刻观测量与上一时刻观测量的估计的残差；1) Obtain the estimated residual of the observed quantity at the current moment and the observed quantity at the previous moment;

残差

Figure 292346DEST_PATH_IMAGE078
；residual
Figure 292346DEST_PATH_IMAGE078
;

其中，

Figure 960087DEST_PATH_IMAGE023
为k时刻第i颗星的伪距观测量，
Figure 301070DEST_PATH_IMAGE024
为k时刻第i颗星的伪距率观测量；
Figure 91171DEST_PATH_IMAGE025
和
Figure 314342DEST_PATH_IMAGE026
为k-1时刻第i颗星的观测量的估计值，等于函数
Figure 739770DEST_PATH_IMAGE079
在k-1时刻估计的定位解处泰勒展开后省略高阶项的值；i=1,…,p，p为可见星数目。in,
Figure 960087DEST_PATH_IMAGE023
is the pseudo-range observation of the i -th star at time k ,
Figure 301070DEST_PATH_IMAGE024
is the pseudorange rate observation of the i -th star at time k ;
Figure 91171DEST_PATH_IMAGE025
and
Figure 314342DEST_PATH_IMAGE026
is the estimated value of the observed quantity of the i -th star at time k -1, which is equal to the function
Figure 739770DEST_PATH_IMAGE079
The value of the high-order term is omitted after Taylor expansion at the estimated positioning solution at time k -1; i= 1,…, p , p is the number of visible stars.

2）利用残差估计截断误差的估计量

Figure 997576DEST_PATH_IMAGE080
和其协方差
Figure 958578DEST_PATH_IMAGE081
：2) Use the residual to estimate the estimator of the truncation error
Figure 997576DEST_PATH_IMAGE080
and its covariance
Figure 958578DEST_PATH_IMAGE081
:

估计量

Figure 934625DEST_PATH_IMAGE082
；Estimator
Figure 934625DEST_PATH_IMAGE082
;

协方差

Figure 881852DEST_PATH_IMAGE083
；Covariance
Figure 881852DEST_PATH_IMAGE083
;

其中，

Figure 56482DEST_PATH_IMAGE084
为半正定的权矩阵上标“-1”代表矩阵的逆，上标“T”代表矩阵的转置。in,
Figure 56482DEST_PATH_IMAGE084
The superscript "-1" of the semi-positive definite weight matrix represents the inverse of the matrix, and the superscript " T " represents the transposition of the matrix.

在步骤S2中，用截断误差的估计值和其协方差修正原模型的过程包括：In step S2, the process of modifying the original model with the estimated value of the truncation error and its covariance includes:

1）将当前时刻的残差估计截断误差的估计量

Figure 391648DEST_PATH_IMAGE085
值代入，修正EKF的一步预测方程；1) Estimate the truncation error of the residual error at the current moment
Figure 391648DEST_PATH_IMAGE085
Substituting the value to modify the one-step prediction equation of EKF;

Figure 510783DEST_PATH_IMAGE086
；
Figure 510783DEST_PATH_IMAGE086
;

2）得到修正后的一步预测的协方差阵；2) Obtain the covariance matrix of the revised one-step forecast;

Figure 651914DEST_PATH_IMAGE087
；
Figure 651914DEST_PATH_IMAGE087
;

3）修正模型后截断误差还会留有一定的残余，考虑这些残余截断误差

Figure 884312DEST_PATH_IMAGE088
，假设
Figure 62484DEST_PATH_IMAGE089
服从零均值高斯分布
Figure 13122DEST_PATH_IMAGE090
，
Figure 426786DEST_PATH_IMAGE091
；重写状态方程得到用截断误差的估计值和其协方差修正后的GNSS/INS紧组合导航系统状态空间模型为：3) After correcting the model, the truncation error will still leave a certain amount of residue. Consider these residual truncation errors
Figure 884312DEST_PATH_IMAGE088
, assuming
Figure 62484DEST_PATH_IMAGE089
Follow a Gaussian distribution with zero mean
Figure 13122DEST_PATH_IMAGE090
,
Figure 426786DEST_PATH_IMAGE091
; Rewrite the state equation to obtain the state-space model of the GNSS/INS compact integrated navigation system corrected by the estimated value of the truncation error and its covariance:

Figure 202106DEST_PATH_IMAGE092
。
Figure 202106DEST_PATH_IMAGE092
.

在步骤S2中，考虑截断误差残余量，对定位误差的统计特性进行修正后，确定出包括截断误差的保护级表达式的过程包括：In step S2, the process of determining the protection level expression including the truncation error includes:

1）根据卡尔曼滤波的误差传播，确定出定位误差；1) According to the error propagation of Kalman filter, the positioning error is determined;

其中，卡尔曼滤波的误差传播过程包括：Among them, the error propagation process of Kalman filtering includes:

依据Kalman滤波公式，有状态后验估计值：According to the Kalman filter formula, there is a state posterior estimate:

Figure 613496DEST_PATH_IMAGE093
；
Figure 613496DEST_PATH_IMAGE093
;

将量测方程代入上式：Substitute the measurement equation into the above formula:

Figure 848168DEST_PATH_IMAGE094
；
Figure 848168DEST_PATH_IMAGE094
;

令

Figure 940889DEST_PATH_IMAGE095
，两边同时减去
Figure 882300DEST_PATH_IMAGE096
，并代入预测方程
Figure 792488DEST_PATH_IMAGE097
进行推导：make
Figure 940889DEST_PATH_IMAGE095
, subtracting both sides
Figure 882300DEST_PATH_IMAGE096
, and substitute into the prediction equation
Figure 792488DEST_PATH_IMAGE097
Do the derivation:

Figure 452139DEST_PATH_IMAGE098
；
Figure 452139DEST_PATH_IMAGE098
;

得到k时刻状态误差

Figure 332239DEST_PATH_IMAGE099
，其中定位误差为状态误差中的位置分量，下面都将定位误差简单写成
Figure 456053DEST_PATH_IMAGE099
。Get the state error at time k
Figure 332239DEST_PATH_IMAGE099
, where the positioning error is the position component in the state error, and the positioning error is simply written as
Figure 456053DEST_PATH_IMAGE099
.

由定位误差表达式

Figure 474825DEST_PATH_IMAGE100
可知，k时刻定位误差表示为k-1时刻状态误差、k-1时刻的过程噪声、k时刻的截断误差估计量、k时刻的截断误差残余、k时刻的GNSS量测噪声五部分的递归函数。By the positioning error expression
Figure 474825DEST_PATH_IMAGE100
It can be seen that the positioning error at time k is expressed as a recursive function of five parts: state error at time k -1, process noise at time k-1, truncation error estimate at time k , truncation error residual at time k , and GNSS measurement noise at time k .

2）进行完好性的假设分配，得到H0假设下，GNSS/INS紧组合导航系统垂直定位误差表达式；2) Carry out the hypothesis assignment of integrity, and obtain the vertical positioning error expression of the GNSS/INS compact integrated navigation system under the H0 assumption;

具体的，在两类假设：Specifically, there are two types of assumptions:

H0假设：无卫星故障发生；H0 assumption: no satellite failure occurs;

H1假设：卫星发生故障的情况。Hypothesis H1: The situation where the satellite fails.

完好性在两类假设间分配：Integrity is distributed between two types of assumptions:

H0假设下系统完好性由保护级完好性保证。H1假设下通过完好性故障监测算法检测卫星故障，实现故障监测完好性。Under the H0 assumption, the system integrity is guaranteed by the protection level integrity. Under the assumption of H1, satellite faults are detected through the integrity fault monitoring algorithm, and the integrity of fault monitoring is realized.

其中，保护级是定位误差的置信上限，

Figure 293876DEST_PATH_IMAGE101
是无故障漏检概率。Among them, the protection level is the upper confidence limit of the positioning error,
Figure 293876DEST_PATH_IMAGE101
is the probability of failure-free missed detection.

具体的，在H0假设下，所述垂直定位误差表达式可表示为：Specifically, under the assumption of H0, the expression of the vertical positioning error can be expressed as:

Figure 56296DEST_PATH_IMAGE102
；
Figure 56296DEST_PATH_IMAGE102
;

其中，

Figure 769037DEST_PATH_IMAGE103
，
Figure 647125DEST_PATH_IMAGE104
为卡尔曼增益；
Figure 281369DEST_PATH_IMAGE105
为
Figure 644217DEST_PATH_IMAGE106
矩阵中与大地坐标系下高度信息对应的行中的元素，
Figure 86831DEST_PATH_IMAGE107
为
Figure 447405DEST_PATH_IMAGE108
与大地坐标系下高度信息对应的行中的元素，
Figure 365682DEST_PATH_IMAGE044
为
Figure 391276DEST_PATH_IMAGE045
与大地坐标系下高度信息对应的行中的元素；in,
Figure 769037DEST_PATH_IMAGE103
,
Figure 647125DEST_PATH_IMAGE104
Gain for Kalman;
Figure 281369DEST_PATH_IMAGE105
for
Figure 644217DEST_PATH_IMAGE106
The elements in the row corresponding to the height information in the geodetic coordinate system in the matrix,
Figure 86831DEST_PATH_IMAGE107
for
Figure 447405DEST_PATH_IMAGE108
The elements in the row corresponding to the height information in the geodetic coordinate system,
Figure 365682DEST_PATH_IMAGE044
for
Figure 391276DEST_PATH_IMAGE045
The elements in the row corresponding to the height information in the geodetic coordinate system;

由于，在本实施例中，建立的17维的误差状态向量中高度误差为第9维；因此，

Figure 281872DEST_PATH_IMAGE109
为
Figure 875664DEST_PATH_IMAGE110
矩阵中第9行元素；
Figure 891025DEST_PATH_IMAGE111
为
Figure 533359DEST_PATH_IMAGE112
中第9行元素；
Figure 340778DEST_PATH_IMAGE044
为
Figure 454272DEST_PATH_IMAGE045
中第9行元素。Because, in the present embodiment, the height error is the 9th dimension in the 17-dimensional error state vector established; therefore,
Figure 281872DEST_PATH_IMAGE109
for
Figure 875664DEST_PATH_IMAGE110
The 9th row element in the matrix;
Figure 891025DEST_PATH_IMAGE111
for
Figure 533359DEST_PATH_IMAGE112
element in line 9;
Figure 340778DEST_PATH_IMAGE044
for
Figure 454272DEST_PATH_IMAGE045
element in line 9.

3）根据垂直定位误差表达式，确定出垂直定位误差分布特性；3) According to the vertical positioning error expression, the distribution characteristics of the vertical positioning error are determined;

Figure 81562DEST_PATH_IMAGE113
、
Figure 527587DEST_PATH_IMAGE114
、
Figure 330458DEST_PATH_IMAGE115
、
Figure 203736DEST_PATH_IMAGE116
、
Figure 52743DEST_PATH_IMAGE117
中各分量相互独立，垂直定位误差
Figure 489410DEST_PATH_IMAGE118
服从高斯分布特性：
Figure 81562DEST_PATH_IMAGE113
,
Figure 527587DEST_PATH_IMAGE114
,
Figure 330458DEST_PATH_IMAGE115
,
Figure 203736DEST_PATH_IMAGE116
,
Figure 52743DEST_PATH_IMAGE117
Each component in is independent of each other, and the vertical positioning error
Figure 489410DEST_PATH_IMAGE118
Obey the Gaussian distribution characteristics:

Figure 209104DEST_PATH_IMAGE119
；
Figure 209104DEST_PATH_IMAGE119
;

4）根据误差分布特性，得到H0假设下包括截断误差的保护级表达式：4) According to the error distribution characteristics, the protection level expression including the truncation error under the H0 assumption is obtained:

Figure 315601DEST_PATH_IMAGE120
；
Figure 315601DEST_PATH_IMAGE120
;

其中，

Figure 527270DEST_PATH_IMAGE121
为无故障漏检系数。
Figure 315098DEST_PATH_IMAGE122
，函数定义为：in,
Figure 527270DEST_PATH_IMAGE121
is the error-free missed detection coefficient.
Figure 315098DEST_PATH_IMAGE122
, the function is defined as:

Figure 951615DEST_PATH_IMAGE123
；
Figure 951615DEST_PATH_IMAGE123
;

该保护级考虑了非线性系统线性化过程产生的截断误差对定位误差的影响，从而使保护级更加严谨、准确。This protection level considers the impact of the truncation error generated by the linearization process of the nonlinear system on the positioning error, so that the protection level is more rigorous and accurate.

在步骤S3中，伪距、伪距率观测量样本和残余截断误差样本数据从实际接收机的观测数据中获取。In step S3, the pseudorange, pseudorange rate observation sample data and residual truncation error sample data are obtained from the actual receiver observation data.

为保证样本数据的独立性，数据采集的时间跨度尽量大、接收机位置分散，例如在北京、上海、拉萨3地对同一星座进行1个月的观测，按30s的时间间隔对全部的伪距、伪距率误差样本进行采样，收集每一颗卫星的伪距、伪距率误差样本各约40000个。In order to ensure the independence of the sample data, the time span of data collection should be as large as possible, and the locations of the receivers should be dispersed. For example, the observation of the same constellation in Beijing, Shanghai, and Lhasa for one month should be carried out at a time interval of 30s for all pseudoranges. , Pseudo-range rate error samples are sampled, and about 40,000 pseudo-range and pseudo-range rate error samples are collected for each satellite.

其中，残余截断误差样本根据实际定位结果计算获取。将每一时刻GNSS/INS组合导航的组合定位解和惯导定位解代入EKF观测方程，用原始观测量减去泰勒展开保留的一阶项，得到该时刻的残余截断误差。收集至少40000个样本。Wherein, the residual truncation error samples are calculated and obtained according to the actual positioning results. Substitute the combined positioning solution of GNSS/INS integrated navigation and inertial navigation positioning solution into the EKF observation equation at each moment, and subtract the first-order item retained by Taylor expansion from the original observations to obtain the residual truncation error at that moment. Collect at least 40000 samples.

在步骤S4中，采用相同的过程，基于极值理论分别包络伪距误差、伪距率误差和残余截断误差，得到伪距误差放大因子、伪距率误差放大因子和残余截断误差放大因子。In step S4, the same process is adopted to respectively envelope the pseudorange error, pseudorange rate error and residual truncation error based on the extremum theory to obtain the pseudorange error amplification factor, pseudorange rate error amplification factor and residual truncation error amplification factor.

其中，残余截断误差放大因子包括伪距残余截断误差放大因子和伪距率残余截断误差放大因子。Wherein, the residual truncation error amplification factor includes a pseudorange residual truncation error amplification factor and a pseudorange rate residual truncation error amplification factor.

具体的，基于极值理论包络伪距误差得到伪距误差放大因子的过程，包括：Specifically, the process of obtaining the pseudorange error amplification factor based on the envelope pseudorange error of the extreme value theory includes:

1）将获取的实际伪距样本取绝对值后从小到大排列，并计算出样本的标准差

Figure 166696DEST_PATH_IMAGE124
；1) Arrange the obtained actual pseudorange samples from small to large after taking the absolute value, and calculate the standard deviation of the samples
Figure 166696DEST_PATH_IMAGE124
;

将获取的实际伪距样本取绝对值后从小到大排列后，构成样本序列

Figure 147553DEST_PATH_IMAGE125
；N为获取的样本总数。并计算出样本的标准差
Figure 270229DEST_PATH_IMAGE124
。Take the absolute value of the obtained actual pseudorange samples and arrange them from small to large to form a sample sequence
Figure 147553DEST_PATH_IMAGE125
; N is the total number of samples obtained. and calculate the standard deviation of the sample
Figure 270229DEST_PATH_IMAGE124
.

2）用平均超出量函数图法选取出一个伪距样本作为阈值；2) Select a pseudo-range sample as the threshold using the average excess amount function graph method;

具体的，平均超出量函数为：

Figure 964516DEST_PATH_IMAGE126
，即对样本中大于阈值T的超出量求均值，得到阈值取T时的点
Figure 288181DEST_PATH_IMAGE127
，不断地从小到大以等间距小步长取T值，并在二维坐标上画出这些点连接的曲线，阈值可选取为曲线中接近线性且斜率为正的部分中的某一个T值。Specifically, the average excess function is:
Figure 964516DEST_PATH_IMAGE126
, that is, the average value of the excess amount greater than the threshold T in the sample is obtained, and the point when the threshold is T is obtained
Figure 288181DEST_PATH_IMAGE127
, continuously take T values from small to large with equal intervals and small steps, and draw the curve connecting these points on the two-dimensional coordinates. The threshold can be selected as a certain T value in the part of the curve that is close to linear and has a positive slope .

3）用Bootstrap法将伪距样本有放回地重采样B次，得到B组样本

Figure 599077DEST_PATH_IMAGE128
；
Figure 728707DEST_PATH_IMAGE129
，
Figure 464450DEST_PATH_IMAGE130
；3) Use the Bootstrap method to resample the pseudorange samples B times with replacement, and obtain the B group samples
Figure 599077DEST_PATH_IMAGE128
;
Figure 728707DEST_PATH_IMAGE129
,
Figure 464450DEST_PATH_IMAGE130
;

3）根据极值理论，对每组样本进行参数估计，确定出参数置信限值；3) According to the extreme value theory, estimate the parameters of each group of samples, and determine the parameter confidence limits;

超出量为每组样本中大于阈值T的样本与T的差，即

Figure 818071DEST_PATH_IMAGE131
，
Figure 819525DEST_PATH_IMAGE132
，
Figure 690529DEST_PATH_IMAGE133
，
Figure 93829DEST_PATH_IMAGE134
为每组样本中大于阈值T的样本的个数，
Figure 883930DEST_PATH_IMAGE135
的分布
Figure 795517DEST_PATH_IMAGE136
可以用样本的分布
Figure 532529DEST_PATH_IMAGE137
表示：The excess amount is the difference between the samples greater than the threshold T in each group of samples and T , that is
Figure 818071DEST_PATH_IMAGE131
,
Figure 819525DEST_PATH_IMAGE132
,
Figure 690529DEST_PATH_IMAGE133
,
Figure 93829DEST_PATH_IMAGE134
is the number of samples greater than the threshold T in each group of samples,
Figure 883930DEST_PATH_IMAGE135
Distribution
Figure 795517DEST_PATH_IMAGE136
The distribution of samples that can be used
Figure 532529DEST_PATH_IMAGE137
express:

Figure 852652DEST_PATH_IMAGE138
；
Figure 852652DEST_PATH_IMAGE138
;

因此可得样本分布

Figure 751337DEST_PATH_IMAGE139
。Therefore, the sample distribution
Figure 751337DEST_PATH_IMAGE139
.

令每组样本中大于阈值样本的概率为

Figure 399488DEST_PATH_IMAGE140
，
Figure 736928DEST_PATH_IMAGE141
。Let the probability of samples greater than the threshold in each group of samples be
Figure 399488DEST_PATH_IMAGE140
,
Figure 736928DEST_PATH_IMAGE141
.

极值理论指出，

Figure 849240DEST_PATH_IMAGE142
的极大值分布服从Gumbel分布，则超出量
Figure 371358DEST_PATH_IMAGE143
服从I型广义Pareto分布：Extreme value theory states that,
Figure 849240DEST_PATH_IMAGE142
The distribution of the maximum value of obeys the Gumbel distribution, then the excess
Figure 371358DEST_PATH_IMAGE143
Obey the type I generalized Pareto distribution:

Figure 569121DEST_PATH_IMAGE144
；
Figure 569121DEST_PATH_IMAGE144
;

为了得到样本分布，则需要估计两个参数

Figure 444673DEST_PATH_IMAGE145
和
Figure 349175DEST_PATH_IMAGE146
。In order to obtain the sample distribution, two parameters need to be estimated
Figure 444673DEST_PATH_IMAGE145
and
Figure 349175DEST_PATH_IMAGE146
.

每个重采样样本中样本超出阈值的概率

Figure 855243DEST_PATH_IMAGE147
，同理可得总体样本中样本超出阈值的概率
Figure 602619DEST_PATH_IMAGE148
。记自助统计量
Figure 642381DEST_PATH_IMAGE149
。The probability of a sample exceeding the threshold in each resampled sample
Figure 855243DEST_PATH_IMAGE147
, in the same way, the probability of the sample exceeding the threshold in the overall sample can be obtained
Figure 602619DEST_PATH_IMAGE148
. Record self-service statistics
Figure 642381DEST_PATH_IMAGE149
.

Figure 729286DEST_PATH_IMAGE146
可进行极大似然估计得到。
Figure 468572DEST_PATH_IMAGE150
的概率密度函数为：
Figure 729286DEST_PATH_IMAGE146
It can be estimated by maximum likelihood.
Figure 468572DEST_PATH_IMAGE150
The probability density function of is:

Figure 640927DEST_PATH_IMAGE151
；
Figure 640927DEST_PATH_IMAGE151
;

Figure 733648DEST_PATH_IMAGE152
；
Figure 733648DEST_PATH_IMAGE152
;

I型广义Pareto分布的对数似然函数为：The log-likelihood function of the type I generalized Pareto distribution is:

Figure 737376DEST_PATH_IMAGE153
；
Figure 737376DEST_PATH_IMAGE153
;

令

Figure 585246DEST_PATH_IMAGE154
，解得：
Figure 431849DEST_PATH_IMAGE155
；make
Figure 585246DEST_PATH_IMAGE154
,Solutions have to:
Figure 431849DEST_PATH_IMAGE155
;

即第b组重采样样本的

Figure 390577DEST_PATH_IMAGE156
最大似然估计值为
Figure 514391DEST_PATH_IMAGE157
。同理可得总体样本的参数估计
Figure 939687DEST_PATH_IMAGE158
。记自助统计量
Figure 86635DEST_PATH_IMAGE159
。That is, the b group of resampled samples
Figure 390577DEST_PATH_IMAGE156
The maximum likelihood estimate is
Figure 514391DEST_PATH_IMAGE157
. In the same way, the parameter estimation of the population sample can be obtained
Figure 939687DEST_PATH_IMAGE158
. Record self-service statistics
Figure 86635DEST_PATH_IMAGE159
.

这样可得B组重采样样本估计的参数值

Figure 911371DEST_PATH_IMAGE160
和
Figure 250211DEST_PATH_IMAGE161
。In this way, the parameter values estimated by group B resampling samples can be obtained
Figure 911371DEST_PATH_IMAGE160
and
Figure 250211DEST_PATH_IMAGE161
.

用Bootstrap法，依据规定的置信度

Figure 439884DEST_PATH_IMAGE162
，确定参数置信限值
Figure 870865DEST_PATH_IMAGE163
和
Figure 171397DEST_PATH_IMAGE164
。Using the Bootstrap method, according to the specified confidence
Figure 439884DEST_PATH_IMAGE162
, to determine the parameter confidence limits
Figure 870865DEST_PATH_IMAGE163
and
Figure 171397DEST_PATH_IMAGE164
.

4）将参数置信限值带入样本分布，根据组合导航伪距误差的概率分布的完好性确定出伪距误差放大因子。4) The parameter confidence limit is brought into the sample distribution, and the pseudo-range error amplification factor is determined according to the integrity of the probability distribution of the integrated navigation pseudo-range error.

将

Figure 879590DEST_PATH_IMAGE165
和
Figure 302481DEST_PATH_IMAGE166
值代入样本分布：Will
Figure 879590DEST_PATH_IMAGE165
and
Figure 302481DEST_PATH_IMAGE166
Substitute the values into the sample distribution:

Figure 158441DEST_PATH_IMAGE167
；
Figure 158441DEST_PATH_IMAGE167
;

得到GNSS/INS组合导航伪距误差的概率分布：Get the probability distribution of GNSS/INS integrated navigation pseudorange error:

Figure 449614DEST_PATH_IMAGE168
；
Figure 449614DEST_PATH_IMAGE168
;

并令

Figure 340210DEST_PATH_IMAGE169
，
Figure 934002DEST_PATH_IMAGE170
为指定的完好性风险。可解得
Figure 949363DEST_PATH_IMAGE171
的分位数
Figure 591697DEST_PATH_IMAGE172
：and order
Figure 340210DEST_PATH_IMAGE169
,
Figure 934002DEST_PATH_IMAGE170
is the specified integrity risk. Solvable
Figure 949363DEST_PATH_IMAGE171
Quantile of
Figure 591697DEST_PATH_IMAGE172
:

Figure 399116DEST_PATH_IMAGE173
；
Figure 399116DEST_PATH_IMAGE173
;

确定放大因子

Figure 524329DEST_PATH_IMAGE174
：Determine the magnification factor
Figure 524329DEST_PATH_IMAGE174
:

Figure 354882DEST_PATH_IMAGE175
。
Figure 354882DEST_PATH_IMAGE175
.

对伪距误差、伪距率误差和残余截断误差分别进行误差包络，得到伪距误差放大因子

Figure 863223DEST_PATH_IMAGE176
，伪距率放大因子
Figure 462832DEST_PATH_IMAGE177
，伪距残余截断误差放大因子
Figure 8214DEST_PATH_IMAGE178
，伪距率残余截断误差放大因子
Figure 388380DEST_PATH_IMAGE179
。The pseudorange error, pseudorange rate error and residual truncation error are respectively enveloped to obtain the pseudorange error amplification factor
Figure 863223DEST_PATH_IMAGE176
, Pseudorange rate magnification factor
Figure 462832DEST_PATH_IMAGE177
, pseudorange residual truncation error amplification factor
Figure 8214DEST_PATH_IMAGE178
, pseudorange rate residual truncation error amplification factor
Figure 388380DEST_PATH_IMAGE179
.

在步骤S5中，计算保护级中，In step S5, in calculating the protection level,

由Kalman滤波量测噪声，量测噪声向量

Figure 638095DEST_PATH_IMAGE180
中每一分量的方差为k时刻的量测噪声矩阵
Figure 279161DEST_PATH_IMAGE181
对角线值，其中，第i颗可见星的伪距量测噪声方差等于k时刻的量测噪声矩阵
Figure 323341DEST_PATH_IMAGE182
第2i-1行第2i-1列值，即
Figure 925223DEST_PATH_IMAGE183
，第i颗可见星的伪距率量测噪声方差等于k时刻的量测噪声矩阵第2i行第2i列值，即
Figure 650734DEST_PATH_IMAGE184
。第i颗可见星伪距的残余截断误差方差等于k时刻的截断误差矩阵
Figure 224935DEST_PATH_IMAGE185
第2i-1行第2i-1列值，即
Figure 502332DEST_PATH_IMAGE186
；第i颗可见星伪距率的残余截断误差方差等于k时刻的截断误差矩阵
Figure 217609DEST_PATH_IMAGE185
第2i行第2i列值，即
Figure 543549DEST_PATH_IMAGE187
。The measurement noise is filtered by Kalman, the measurement noise vector
Figure 638095DEST_PATH_IMAGE180
The variance of each component in is the measurement noise matrix at time k
Figure 279161DEST_PATH_IMAGE181
Diagonal values, where the pseudorange measurement noise variance of the i -th visible star is equal to the measurement noise matrix at time k
Figure 323341DEST_PATH_IMAGE182
Row 2 i -1, column 2 i -1 value, ie
Figure 925223DEST_PATH_IMAGE183
, the measurement noise variance of the pseudorange rate of the i-th visible star is equal to the value of the measurement noise matrix at the second i row and the second i column of the measurement noise matrix at time k , that is
Figure 650734DEST_PATH_IMAGE184
. The residual truncation error variance of the i -th visible star pseudorange is equal to the truncation error matrix at time k
Figure 224935DEST_PATH_IMAGE185
Row 2 i -1, column 2 i -1 value, ie
Figure 502332DEST_PATH_IMAGE186
; The residual truncation error variance of the i -th visible star pseudorange rate is equal to the truncation error matrix at time k
Figure 217609DEST_PATH_IMAGE185
The value of row 2i , column 2i , i.e.
Figure 543549DEST_PATH_IMAGE187
.

根据伪距误差放大因子，伪距率放大因子和残余截断误差放大因子；放大伪距误差方差、伪距率误差方差和残余截断误差方差后，伪距量测噪声、伪距率量测噪声和残余截断误差的分布为：According to the pseudorange error amplification factor, pseudorange rate amplification factor and residual truncation error amplification factor; after amplifying the pseudorange error variance, pseudorange rate error variance and residual truncation error variance, the pseudorange measurement noise, pseudorange rate measurement noise and The distribution of the residual truncation error is:

Figure 34573DEST_PATH_IMAGE188
；
Figure 34573DEST_PATH_IMAGE188
;

Figure 358238DEST_PATH_IMAGE189
；
Figure 358238DEST_PATH_IMAGE189
;

Figure 872396DEST_PATH_IMAGE190
；
Figure 872396DEST_PATH_IMAGE190
;

Figure 64343DEST_PATH_IMAGE191
；
Figure 64343DEST_PATH_IMAGE191
;

Figure 534507DEST_PATH_IMAGE192
为第i颗可见星的伪距误差放大因子，
Figure 91390DEST_PATH_IMAGE053
为第i颗可见星的伪距率误差放大因子，
Figure 155161DEST_PATH_IMAGE054
为第i颗可见星的伪距残余截断误差放大因子，
Figure 88482DEST_PATH_IMAGE055
为第i颗可见星的伪距率残余截断误差放大因子。
Figure 534507DEST_PATH_IMAGE192
is the pseudorange error amplification factor of the i -th visible star,
Figure 91390DEST_PATH_IMAGE053
is the pseudorange rate error amplification factor of the i -th visible star,
Figure 155161DEST_PATH_IMAGE054
is the pseudorange residual truncation error amplification factor of the i -th visible star,
Figure 88482DEST_PATH_IMAGE055
is the pseudorange rate residual truncation error amplification factor of the i -th visible star.

将各方差代入计算，得最终保护级：Substituting the variances into the calculation, the final protection level is obtained:

Figure 163886DEST_PATH_IMAGE193
；
Figure 163886DEST_PATH_IMAGE193
;

其中，in,

Figure 953987DEST_PATH_IMAGE051
。
Figure 953987DEST_PATH_IMAGE051
.

在本实施例中，建立的17维的误差状态向量中高度误差为第9维；对应的最终保护级

Figure 442737DEST_PATH_IMAGE194
；即v=9。In this embodiment, the altitude error in the established 17-dimensional error state vector is the ninth dimension; the corresponding final protection level
Figure 442737DEST_PATH_IMAGE194
; that is, v= 9.

综上所述，通过本实施例基于截断误差估计的EKF组合导航保护级确定方法，针对EKF线性化过程中产生的截断误差使保护级计算不够准确的问题，通过估计出的截断误差的大小和其协方差，减小了截断误差的不确定性，修正定位误差及其分布，从而得到更准确的保护级，使其满足给定的完好性风险；通过采用极值理论包络实际采样的伪距、伪距率等误差、截断残余误差分布，定位出更准确的误差尾部概率分布。结合更准确的保护级和更准确的误差尾部概率分布，保证了保护级的可靠性。进而保障了非线性GNSS/INS组合导航系统的完好性。To sum up, through the method for determining the protection level of EKF integrated navigation based on truncation error estimation in this embodiment, in view of the problem that the truncation error generated during the EKF linearization process makes the calculation of the protection level not accurate enough, through the estimated size of the truncation error and Its covariance reduces the uncertainty of the truncation error, corrects the positioning error and its distribution, and obtains a more accurate protection level to meet the given integrity risk; Errors such as range and pseudorange rate, and truncated residual error distribution, locate a more accurate error tail probability distribution. Combined with a more accurate protection level and a more accurate error tail probability distribution, the reliability of the protection level is guaranteed. Thus, the integrity of the nonlinear GNSS/INS integrated navigation system is guaranteed.

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。The above is only a preferred embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any person skilled in the art within the technical scope disclosed in the present invention can easily think of changes or Replacement should be covered within the protection scope of the present invention.

Claims (10)

1. An EKF integrated navigation protection level calculation method based on truncation error estimation is characterized by comprising the following steps:
establishing a GNSS/INS tight combination navigation system state space model considering truncation errors, and estimating truncation error values and covariance thereof;
modifying the model with the estimated truncation error and its covariance; calculating the statistical characteristic of the corrected positioning error by considering the residual amount of the truncation error to obtain a protection level expression comprising the truncation error;
acquiring actual pseudo-range, pseudo-range rate observed quantity samples and corresponding residual truncation error sample data;
respectively enveloping pseudo range, pseudo range rate observation error and residual truncation error sample data based on an extreme value theory to obtain corresponding amplification factors;
and substituting the amplification factor into the protection level expression to obtain the protection level of the GNSS/INS tight combination navigation system.
2. The EKF integrated navigation protection level calculation method based on truncation error estimation according to claim 1,
the GNSS/INS tightly combined navigation system state space model considering truncation errors comprises the following steps:
Figure 297192DEST_PATH_IMAGE001
；
wherein,kis the time of day or the like,
Figure 620857DEST_PATH_IMAGE002
is thatnA state vector of dimensions;
Figure 666173DEST_PATH_IMAGE003
is thatmA measurement vector of dimensions;
Figure 326962DEST_PATH_IMAGE004
is composed ofnThe state of the stage is transferred to the matrix in one step,
Figure 298591DEST_PATH_IMAGE005
is composed ofm×nA measurement matrix of orders;
Figure 386633DEST_PATH_IMAGE006
to truncate the error vector;
Figure 919245DEST_PATH_IMAGE007
is composed ofnThe system noise vector of the dimension(s),
Figure 524670DEST_PATH_IMAGE008
is composed ofmThe measurement noise vector of the dimension.
3. The method of claim 2, wherein the state vector of the GNSS/INS tightly-integrated navigation system state space model is a 17-dimensional error state vector:
Figure 459128DEST_PATH_IMAGE009
wherein,
Figure 718071DEST_PATH_IMAGE010
is the three-dimensional velocity error;
Figure 737980DEST_PATH_IMAGE011
is the three-dimensional attitude angle error;
Figure 130784DEST_PATH_IMAGE012
is the three-dimensional position error;
Figure 919748DEST_PATH_IMAGE013
is the three-dimensional accelerometer error;
Figure 615172DEST_PATH_IMAGE014
is the three-dimensional gyroscope error;
Figure 263322DEST_PATH_IMAGE015
is a one-dimensional receiver constant offset;
Figure 804025DEST_PATH_IMAGE016
is a one-dimensional receiver clock drift.
4. The EKF combined navigation protection level calculation method based on truncation error estimation as claimed in claim 3, wherein,
measurement of GNSS/INS tightly combined navigation system state space model quantity as GNSS observation pseudo range
Figure 713075DEST_PATH_IMAGE017
Pseudorange rates
Figure 939919DEST_PATH_IMAGE018
Pseudoranges to INS predictions
Figure 934420DEST_PATH_IMAGE019
Pseudorange rate
Figure 544393DEST_PATH_IMAGE020
The difference between them;
measuring noise vectors
Figure 183316DEST_PATH_IMAGE021
Including GNSS observed pseudoranges and pseudorange rate noise vectors.
5. The EKF integrated navigation protection level calculation method based on truncation error estimation according to claim 4,
estimating a truncation error value and a covariance thereof by using an NPF method, specifically comprising:
1) Obtaining a residual error of the current-time observed quantity and the estimation of the last-time observed quantity;
residual error
Figure 220542DEST_PATH_IMAGE022
；
Wherein,
Figure 436760DEST_PATH_IMAGE023
is composed ofkAt the first momentiThe pseudorange observations of the satellites,
Figure 975057DEST_PATH_IMAGE024
is composed ofkAt the first momentiPseudo range rate observed quantity of the particle;
Figure 593120DEST_PATH_IMAGE025
and
Figure 801248DEST_PATH_IMAGE026
is composed ofk-1 time of dayiEstimated value of observed quantity of particles, equal to function
Figure 645707DEST_PATH_IMAGE027
In thatk-values of the high order terms omitted after taylor expansion at the estimated localization solution at time 1;i=1,…,p，pnumber of visible stars;
2) Estimating truncation error using residual errorIs estimated by
Figure 331903DEST_PATH_IMAGE028
And its covariance
Figure 70052DEST_PATH_IMAGE029
：
Estimation quantity
Figure 329040DEST_PATH_IMAGE030
；
Covariance
Figure 519850DEST_PATH_IMAGE031
；
Wherein,
Figure 275316DEST_PATH_IMAGE032
is a matrix of weights that is semi-positively determined,
Figure 477759DEST_PATH_IMAGE033
to measure the noise matrix.
6. The EKF integrated navigation protection level calculation method based on truncation error estimation according to claim 5,
the GNSS/INS tightly combined navigation system state space model corrected by the estimated value of the truncation error and the covariance thereof is as follows:
Figure 27689DEST_PATH_IMAGE034
；
in the formula,
Figure 971374DEST_PATH_IMAGE035
in order to estimate the truncation error residual quantity by an NPF method,
Figure 389586DEST_PATH_IMAGE036
。
7. the EKF integrated navigation protection level calculation method based on truncation error estimation according to claim 6, wherein,
the determination of the protection level expression including the truncation error comprises:
1) Determining a positioning error according to error propagation of Kalman filtering;
2) Performing integrity hypothesis distribution to obtain a vertical positioning error expression of the GNSS/INS tight combination navigation system under the H0 hypothesis;
3) Determining the distribution characteristic of the vertical positioning error according to the vertical positioning error expression;
4) And obtaining a protection level expression comprising the truncation error under the assumption of H0 according to the error distribution characteristics.
8. The EKF integrated navigation protection level calculation method based on truncation error estimation according to claim 7,
under the assumption of H0, the vertical positioning error expression is as follows:
Figure 836748DEST_PATH_IMAGE037
；
wherein,
Figure 557579DEST_PATH_IMAGE038
；
Figure 598347DEST_PATH_IMAGE039
is the Kalman gain;vas error state vectorsXThe line number corresponding to the height information is obtained;
Figure 695616DEST_PATH_IMAGE040
is composed of
Figure 262864DEST_PATH_IMAGE041
Corresponding to height information in a matrixThe elements in the row of (a) and (b),
Figure 780695DEST_PATH_IMAGE042
is composed of
Figure 433393DEST_PATH_IMAGE043
The elements in the row of the matrix corresponding to the height information,
Figure 68774DEST_PATH_IMAGE044
is composed of
Figure 100315DEST_PATH_IMAGE045
Elements in rows of the matrix corresponding to the height information;
vertical positioning error
Figure 428528DEST_PATH_IMAGE046
Obeying the characteristics of Gaussian distribution;
Figure 302943DEST_PATH_IMAGE047
；
h0 assumes that the guard level expression including truncation error:
Figure 866649DEST_PATH_IMAGE048
；
wherein
Figure 408488DEST_PATH_IMAGE049
The failure-free undetected coefficient.
9. The EKF combined navigation protection level calculation method based on truncation error estimation as claimed in claim 8,
the method for obtaining the corresponding amplification factor based on the extreme value theory envelope pseudo range sample data comprises the following steps:
1) Arranging the obtained actual pseudo-range samples from small to large after taking absolute values, and calculating the standard deviation of the samples;
2) Selecting a pseudo range sample as a threshold value by using an average excess function graph method;
3) Re-sampling the pseudo-range samples for B times in a Bootstrap method to obtain a group B of samples;
3) According to an extreme value theory, performing parameter estimation on each group of samples to determine a parameter confidence limit value;
4) And substituting the parameter confidence limit value into the sample distribution, and determining a pseudo-range error amplification factor according to the integrity of the probability distribution of the combined navigation pseudo-range error.
10. The EKF combined navigation protection level calculation method based on truncation error estimation as claimed in claim 9, wherein,
substituting the amplification factor into the protection level expression to obtain the protection level of the GNSS/INS tightly combined navigation system
Figure 642024DEST_PATH_IMAGE050
(ii) a Wherein,
Figure 879101DEST_PATH_IMAGE051
；
wherein,
Figure 856284DEST_PATH_IMAGE052
is as followsiThe pseudorange error amplification factor for the visible satellites,
Figure 987051DEST_PATH_IMAGE053
is as followsiThe pseudo-range rate error magnification factor for the visible satellites,
Figure 283166DEST_PATH_IMAGE054
is as followsiThe pseudorange residue of the visible satellites truncates the error amplification factor,
Figure 132173DEST_PATH_IMAGE055
is as followsiThe pseudorange rate residual truncation error magnification factor for the visible satellites.
CN202211341296.4A 2022-10-31 2022-10-31 EKF (extended Kalman Filter) integrated navigation protection level determination method based on truncation error estimation Active CN115408652B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211341296.4A CN115408652B (en) 2022-10-31 2022-10-31 EKF (extended Kalman Filter) integrated navigation protection level determination method based on truncation error estimation

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211341296.4A CN115408652B (en) 2022-10-31 2022-10-31 EKF (extended Kalman Filter) integrated navigation protection level determination method based on truncation error estimation

Publications (2)

Publication Number Publication Date
CN115408652A true CN115408652A (en) 2022-11-29
CN115408652B CN115408652B (en) 2023-02-03

Family

ID=84167738

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211341296.4A Active CN115408652B (en) 2022-10-31 2022-10-31 EKF (extended Kalman Filter) integrated navigation protection level determination method based on truncation error estimation

Country Status (1)

Country Link
CN (1) CN115408652B (en)

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104898140A (en) * 2015-05-28 2015-09-09 北京航空航天大学 Extreme value theory-based error envelope method of satellite navigation ground-based augmentation system
EP3009860A1 (en) * 2014-10-16 2016-04-20 GMV Aerospace and Defence S.A. Method for computing an error bound of a Kalman filter based GNSS position solution
CN105806339A (en) * 2016-05-14 2016-07-27 中卫物联成都科技有限公司 Integrated navigation method and device based on GNSS, INS and time keeping systems
CN108646277A (en) * 2018-05-03 2018-10-12 山东省计算中心（国家超级计算济南中心） The Beidou navigation method adaptively merged with Extended Kalman filter based on robust
CN112034491A (en) * 2020-08-31 2020-12-04 中国电子科技集团公司第五十四研究所 Integrity protection level calculation method based on error core distribution
US20210072407A1 (en) * 2019-09-10 2021-03-11 Trimble Inc. Protection level generation methods and systems for applications using navigation satellite system (nss) observations
CN113483759A (en) * 2021-06-29 2021-10-08 北京航空航天大学 Integrity protection level calculation method for GNSS, INS and Vision integrated navigation system
CN114721017A (en) * 2022-03-04 2022-07-08 北京理工大学 GNSS/INS integrated navigation autonomous integrity monitoring method

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP3009860A1 (en) * 2014-10-16 2016-04-20 GMV Aerospace and Defence S.A. Method for computing an error bound of a Kalman filter based GNSS position solution
CN104898140A (en) * 2015-05-28 2015-09-09 北京航空航天大学 Extreme value theory-based error envelope method of satellite navigation ground-based augmentation system
CN105806339A (en) * 2016-05-14 2016-07-27 中卫物联成都科技有限公司 Integrated navigation method and device based on GNSS, INS and time keeping systems
CN108646277A (en) * 2018-05-03 2018-10-12 山东省计算中心（国家超级计算济南中心） The Beidou navigation method adaptively merged with Extended Kalman filter based on robust
US20210072407A1 (en) * 2019-09-10 2021-03-11 Trimble Inc. Protection level generation methods and systems for applications using navigation satellite system (nss) observations
CN112034491A (en) * 2020-08-31 2020-12-04 中国电子科技集团公司第五十四研究所 Integrity protection level calculation method based on error core distribution
CN113483759A (en) * 2021-06-29 2021-10-08 北京航空航天大学 Integrity protection level calculation method for GNSS, INS and Vision integrated navigation system
CN114721017A (en) * 2022-03-04 2022-07-08 北京理工大学 GNSS/INS integrated navigation autonomous integrity monitoring method

Non-Patent Citations (5)

* Cited by examiner, † Cited by third party
Title
HAITAO JIANG ETC.: "An Effective Integrity Monitoring Scheme for GNSS/INS/Vision Integration Based on Error State EKF Model", 《IEEE SENSORS JOURNAL》 *
NI ZHU ETC.: "Extended Kalman Filter (EKF) Innovation-based Integrity Monitoring Scheme with C/N0 Weighting", 《IEEE》 *
丛丽等: "基于INS辅助CLAMBDA与AFM的GPS/INS组合导航测姿方法", 《系统工程与电子技术》 *
于耕等: "基于卡尔曼滤波的卫星地基增强系统位置域完好性监测分析", 《科学技术与工程》 *
周召发等: "考虑里程计截断误差的SINS/OD组合导航算法", 《中国惯性技术学报》 *

Also Published As

Publication number Publication date
CN115408652B (en) 2023-02-03

Similar Documents

Publication Publication Date Title
CN110196443B (en) Fault-tolerant integrated navigation method and system of aircraft
CN108519615B (en) Mobile robot autonomous navigation method based on combined navigation and feature point matching
US7711482B2 (en) Hybrid INS/GNSS system with integrity monitoring and method for integrity monitoring
CN110715659A (en) Zero-speed detection method, pedestrian inertial navigation method, device and storage medium
US20240402361A1 (en) System and method for computing positioning protection levels
US8082099B2 (en) Aircraft navigation using the global positioning system and an attitude and heading reference system
CN110140065A (en) GNSS receiver protection class
CN113670337A (en) Method for detecting slow-changing fault of GNSS/INS combined navigation satellite
CN114384569A (en) Terminal positioning method, apparatus, device and medium
CN115291253B (en) A vehicle positioning integrity monitoring method and system based on residual detection
CN114179825A (en) Method for obtaining confidence of measurement value through multi-sensor fusion and automatic driving vehicle
Anbu et al. Integration of inertial navigation system with global positioning system using extended Kalman filter
CN116358566B (en) Coarse detection combined navigation method based on robust adaptive factor
US20240094415A1 (en) System and method for detecting outliers in gnss observations
WO2021236454A1 (en) Multi-inertial measurement unit fusion for fine-resolution position estimation
CN116399351A (en) Vehicle position estimation method
CN115408652B (en) EKF (extended Kalman Filter) integrated navigation protection level determination method based on truncation error estimation
US11788842B2 (en) Hybrid AHRS system comprising a device for measuring the integrity of the calculated attitude
CN114199236A (en) Positioning data processing method and device, electronic equipment and automatic driving vehicle
CN118031951A (en) Multisource fusion positioning method, multisource fusion positioning device, electronic equipment and storage medium
JP5058594B2 (en) POSITIONING POINT ESTIMATION DEVICE, POSITIONING POINT ESTIMATION METHOD, PROGRAM THEREOF, AND RECORDING MEDIUM
CN112985456B (en) Navigation positioning correction method integrating airborne height measurement data and topographic data
CN115597594A (en) Mileage correction factor-based satellite-inertial navigation method and system
CN115327587A (en) Method and system for orbit error correction of low-orbit satellites based on GNSS positioning information
Chu et al. Performance comparison of tight and loose INS-Camera integration

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