CN115390114A - 一种基于idr粗差探测的gnss定位方法、装置及存储介质 - Google Patents

一种基于idr粗差探测的gnss定位方法、装置及存储介质 Download PDF

Info

Publication number
CN115390114A
CN115390114A CN202211039658.4A CN202211039658A CN115390114A CN 115390114 A CN115390114 A CN 115390114A CN 202211039658 A CN202211039658 A CN 202211039658A CN 115390114 A CN115390114 A CN 115390114A
Authority
CN
China
Prior art keywords
normal
abnormal
error
variance
unit weight
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
CN202211039658.4A
Other languages
English (en)
Other versions
CN115390114B (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.)
Tongji University
Original Assignee
Tongji 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 Tongji University filed Critical Tongji University
Priority to CN202211039658.4A priority Critical patent/CN115390114B/zh
Publication of CN115390114A publication Critical patent/CN115390114A/zh
Application granted granted Critical
Publication of CN115390114B publication Critical patent/CN115390114B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • 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
    • 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/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/13—Receivers
    • G01S19/35—Constructional details or hardware or software details of the signal processing chain
    • G01S19/37—Hardware or software details of the signal processing chain

Landscapes

  • Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Signal Processing (AREA)
  • Complex Calculations (AREA)

Abstract

本发明涉及一种基于IDR粗差探测的GNSS定位方法、装置及存储介质，其中方法包括：步骤1)获取GNSS观测值；步骤2)求解待估参数和观测值残差方差‑协方差矩阵的最小二乘估值；步骤3)初始化正常集合和异常集合；步骤4)确定待检集合，确定w统计量，筛选最小w统计量；步骤5)对最小w统计量进行显著性检验，若显著性检验通过，则将待检集合中的观测值并入正常集合，更新异常集合和正常集合，若未通过，则转步骤7)；步骤6)判断是否满足迭代终止条件，若是，则转步骤7)，若否，则返回步骤4)；步骤7)剔除异常集合中的值，并基于正常集合求解待估参数，得到定位结果。与现有技术相比，本发明具有提高了粗差探测的准确度和GNSS定位精度等优点。

Description

一种基于IDR粗差探测的GNSS定位方法、装置及存储介质
技术领域
本发明涉及大地测量数据处理技术，尤其是涉及一种基于IDR粗差探测的 GNSS定位方法、装置及存储介质。
背景技术
最小二乘(Least Square,LS)估计，是目前最常用的GNSS定位解算方法，能实现无偏性无偏性的前提是使用现实的函数模型和随机模型。也就是说，如果观测中只有包含随机误差，待估参数的最小二乘解就具有最优无偏性。然而，当观测值中存在粗差时，最小二乘解的精度会受到严重的影响。因此，在最小二乘平差之前，对于粗差的探测和识别是至关重要的。
经过半个多世纪的大地测量学发展，数据探测法是粗差检测和识别的最成熟方法之一，已被用作大地测量学数据质量控制的标准方法。通常，数据探测法是基于统计假设检验原理，用于检测和识别线性模型中的单个粗差。它起源于Baarda在 1967和1968年的开创性工作，后来由Pope将其扩展到观测值精度未知的情况。数据探测法通过对每个单独的观测值的筛选来确定其中是否存在粗差，也就是众所周知的w检验。w检验中使用的检验统计量通常是标准化或学生化的最小二乘残差。
数据探测法首先对原假设和一组备择假设模型进行假设检验。在检测过程中，通常遵循正态分布、τ分布、χ^2分布、F分布等检验统计量来决定观测系统中是否存在粗差。识别是一种模型选择，目的是选择最值得信赖的、能够正确排除未建模粗差的备选假设模型。最后，在被识别模型下进行自适应或参数估计，以消除或减轻对未知参数的潜在偏差。在应用数据探测法时，一个重要前提是需确定可疑粗差点的数量，这在实际应用中似乎是不现实的。
在现代大地测量学中，当数据量很大时平差系统中极有可能存在多个粗差。在这种情况下，理论上为单个粗差开发的数据探测法通常通过一个接一个地处理一个粗差来连续实施以处理多个粗差，这被称为迭代数据探测法(Iterative data snooping,IDS)。在不事先确定粗差个数的情况下，在单个粗差假设下进行连续检测、识别和自适应。一旦确定了一个含粗差的观测值，通过排除已确定粗差的观测值重新进行最小二乘平差。然而，这种方法在理论上并不严谨，因为它假设每次迭代中只有一个粗差，但这个假设在下一次迭代中立即被推翻。实际上，在多个粗差的情况下，多个粗差的假设会受到掩蔽和淹没效应的阻碍。一方面，多个粗差很容易相互掩盖，增加了检测和识别的难度。另一方面，如果要检测的粗差的数量很大时，则统计检验往往会声明比实际更多的粗差。随着迭代次数的增加，漏检和错误识别的概率将呈指数级增加，这将在参数估计中引入严重偏差。
发明内容
本发明的目的就是为了提供一种基于迭代数据精检法(Iterative datarefining， IDR)粗差探测的GNSS定位方法，将数据探测法中的w检验应用于数据精炼和检验过程，确定一种迭代数据精检法，并基于迭代数据精检法利用正常集合中的观测值对异常集合中的观测值进行逐个地筛选和检验，将异常集合中被检验为正常的观测值放入正常集合中，直到异常集合中的观测值都被检验为粗差或异常集合中没有观测值，得到可靠性更高的正常集合，并基于正常集合进行粗差剔除和参数估计，有效缓解掩蔽和淹没效应，提高粗差探测与识别的精度，得到高精度GNSS定位结果。
本发明的目的可以通过以下技术方案来实现：
一种基于IDR粗差探测的GNSS定位方法，包括以下步骤：
步骤1)获取GNSS的观测值；
步骤2)基于Gauss-Markov模型求解待估参数和观测值残差的方差-协方差矩阵的最小二乘估值；
步骤3)基于预配置的最大可疑粗差个数初始化观测值的正常集合和异常集合；
步骤4)确定待检集合，并基于异常集合与待估参数、观测值残差的方差-协方差矩阵的最小二乘估值确定w统计量，筛选最小w统计量；
步骤5)对最小w统计量进行显著性检验，若显著性检验通过，则将待检集合中的观测值并入正常集合，更新异常集合和正常集合，若显著性检验未通过，则结束迭代，转步骤7)；
步骤6)判断是否满足迭代终止条件，若是，则转步骤7)，若否，则返回步骤 4)进行迭代；
步骤7)剔除异常集合中的值，基于正常集合与Gauss-Markov模型进行参数估计，求解待估参数，得到无粗差影响的GNSS高精度定位结果。
所述步骤2)包括以下步骤：
步骤2-1)建立基于Gauss-Markov模型的GNSS线性观测方程
y＝Ax+e,e～N(0,∑)
式中，A为m×n阶非随机设计矩阵；x为n×1阶待估参数向量；y为m×1阶观测值向量；e为m×1阶偶然误差；∑为m×m阶表征观测值精度的方差-协方差矩阵：
∑＝σ2Q
式中，Q为m×m阶协因数矩阵；σ2为协因数矩阵对应的单位权方差；
步骤2-2)确定待估参数的最小二乘解：
Figure BDA0003819656650000031
步骤2-3)确定观测值残差：
Figure BDA0003819656650000032
步骤2-4)确定待估参数和观测值残差的方差-协方差矩阵；
当单位权方差σ2已知时，
待估参数的方差-协方差矩阵为：
Figure BDA0003819656650000033
观测值残差的方差-协方差矩阵为：
Figure BDA0003819656650000034
当单位权方差σ2未知时，
待估参数的方差-协方差矩阵为
Figure BDA0003819656650000035
观测值残差的方差-协方差矩阵为
Figure BDA0003819656650000036
其中，
Figure BDA0003819656650000037
Figure BDA0003819656650000038
Figure BDA0003819656650000039
为单位权方差σ2的最小二乘估值：
Figure BDA00038196566500000310
。
所述步骤3)包括以下步骤：
步骤3-1)假定平差系统中所有的观测值形成一个全集，则每一个观测值被唯一地划分为到正常集合或异常集合中，基于最大可疑粗差个数qmax构造原始异常集合：
Figure BDA0003819656650000041
与原始异常集合对应的原始正常集合记为：
Figure BDA0003819656650000042
步骤3-2)从原始异常-正常集合的组合中筛选满足参数可解条件的组合，所述参数可解条件为矩阵
Figure BDA0003819656650000043
或
Figure BDA0003819656650000044
可逆，其中，
Figure BDA0003819656650000045
是异常集合
Figure BDA0003819656650000046
对应的m×qmax阶且只包含0和1两种元素的设计矩阵，矩阵
Figure BDA0003819656650000047
的每一列代表异常集合
Figure BDA0003819656650000048
中的一个粗差；
步骤3-3)从满足参数可解条件的组合中筛选出唯一一对正常-异常集合作为初始正常-异常集合，完成初始化。
筛选得到的初始正常-异常集合为第j个异常集合
Figure BDA0003819656650000049
当单位权中误差σ已知时，
Figure BDA00038196566500000410
当单位权中误差σ未知时，
Figure BDA00038196566500000411
所述步骤4)包括以下步骤：
步骤4-1)在第k次迭代中，将所有的观测值分为互无交集的三个集合：异常集合
Figure BDA00038196566500000412
正常集合
Figure BDA00038196566500000413
和待检集合{yi}，其中，yi∈O(k-1)，O(k-1)为在第k-1次迭代中保留的异常集合，异常集合
Figure BDA00038196566500000414
的元素个数为qmax-k，正常集合
Figure BDA00038196566500000415
的元素个数为m-qmax+k-1，待检集合{yi}的元素个数为1，三个集合同时满足以下两个关系条件：
Figure BDA00038196566500000416
Figure BDA00038196566500000417
步骤4-2)对第k次迭代的异常集合
Figure BDA0003819656650000051
中的每个观测值构造w统计量
Figure BDA0003819656650000052
当单位权中误差σ已知时：
Figure BDA0003819656650000053
当单位权中误差σ未知时：
Figure BDA0003819656650000054
其中，
Figure BDA0003819656650000055
Figure BDA0003819656650000056
是一个m×(qmax-k)阶设计矩阵，由异常集合
Figure BDA0003819656650000057
所生成，矩阵
Figure BDA0003819656650000058
的每一列代表集合
Figure BDA0003819656650000059
中的一个异常观测值，
Figure BDA00038196566500000510
是一个由待检集合{yi}所生成的m×1 阶的设计向量，
Figure BDA00038196566500000511
是一个由集合
Figure BDA00038196566500000512
生成m×(qmax-k+1)阶设计矩阵，即Hi＝[Ci ci]；
步骤4-3)筛选最小w统计量：
当单位权中误差σ已知时：
Figure BDA00038196566500000513
当单位权中误差σ未知时：
Figure BDA00038196566500000514
所述步骤5)包括以下步骤：
步骤5-1)基于预配置的的显著水平α确定显著性检验阈值；
步骤5-2)比较最小w统计量的值与显著性检验阈值，
若最小w统计量的值小于等于显著性检验阈值，则显著性检验通过，将待检集合{yi}中的观测值yj视为正常观测值，将其放入正常集合，更新正常集合和异常集合，即：
I(k)＝I(k-1)∪{yj}
O(k)＝O(k-1)-{yj}
若最小w统计量的值大于显著性检验阈值，则显著性检验未通过，异常集合中的观测值都被视为异常观测值，结束迭代。
单位权中误差σ已知时，显著性检验阈值kα为：
kα＝Nα/2(0,1)
其中，Nα/2(0,1)表示标准正态分布分布在α/2处的分位值；
单位权中误差σ未知时，显著性检验阈值k′α为：
k′α＝tα/2(m-n-qmax+k-1,0)
其中，tα/2(m-n-qmax+k-1,0)表示自由度为m-n-qmax+k-1的学生分布在α/2处的分位值。
所述迭代终止条件为迭代次数等于最大可疑粗差个数。
一种基于IDR粗差探测的GNSS定位装置，包括存储器、处理器，以及存储于所述存储器中的程序，所述处理器执行所述程序时实现如上述所述的方法。
一种存储介质，其上存储有程序，所述程序被执行时实现如上述所述的方法。
与现有技术相比，本发明具有以下有益效果：
本发明通过初始化异常集合和正常集合，保证在迭代开始时所有的异常值均已被归入异常集合，再通过迭代的方法从异常值中找出最有可能正常的值并进行检验，从而判定其是否能归入正常集合，从而有效地缓解了掩蔽和淹没效应，保证了正常集合中观测值的数量和质量，有利于提高后续基于正常集合进行的参数估计的精度，从而提高了GNSS系统的定位精度，并且在粗差的探测和识别方面，具有更强的鲁棒性。
附图说明
图1为本发明的方法流程图；
图2为σ已知的条件下IDS和IDR的正确决策概率；
图3为IDS、IDR和LS在σ已知条件下的RMSE比较结果图；
图4为σ未知的条件下IDS和IDR的正确决策概率；
图5为IDS、IDR和LS在σ未知条件下的RMSE比较结果图。
具体实施方式
下面结合附图和具体实施例对本发明进行详细说明。本实施例以本发明技术方案为前提进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。
本发明将数据探测法中的w检验应用于数据精炼和检验过程，得到了一种基于Gauss-Markov模型的粗差探测方法，称为数据精检法(data refining)，为了其迭代形式也可以用来进行多粗差的探测和识别，又称为迭代数据精检法(Iterative data refining，IDR)，并在IDR粗差探测的基础上提出了一种高精度的GNSS定位方法，如图1所示，包括以下步骤：
步骤1)获取GNSS的观测值；
步骤2)基于Gauss-Markov模型求解待估参数和观测值残差的方差-协方差矩阵的最小二乘估值；
步骤2-1)建立基于Gauss-Markov模型的GNSS线性观测方程
y＝Ax+e,e～N(0,∑)
式中，A为m×n阶非随机设计矩阵；x为n×1阶待估参数向量；y为m×1阶观测值向量；e为m×1阶偶然误差；∑为m×m阶表征观测值精度的方差-协方差矩阵：
∑＝σ2Q
式中，Q为m×m阶协因数矩阵；σ2为协因数矩阵对应的单位权方差；
步骤2-2)确定待估参数的最小二乘解：
Figure BDA0003819656650000071
步骤2-3)确定观测值残差：
Figure BDA0003819656650000072
步骤2-4)确定待估参数和观测值残差的方差-协方差矩阵；
当单位权方差σ2已知时，
待估参数的方差-协方差矩阵为：
Figure BDA0003819656650000073
观测值残差的方差-协方差矩阵为：
Figure BDA0003819656650000081
当单位权方差σ2未知时，
待估参数的方差-协方差矩阵为
Figure BDA0003819656650000082
观测值残差的方差-协方差矩阵为
Figure BDA0003819656650000083
其中，
Figure BDA0003819656650000084
Figure BDA0003819656650000085
Figure BDA0003819656650000086
为单位权方差σ2的最小二乘估值：
Figure BDA0003819656650000087
步骤3)基于预配置的最大可疑粗差个数初始化观测值的正常集合和异常集合；
步骤3-1)假定平差系统中所有的观测值形成一个全集，则每一个观测值被唯一地划分为到正常集合或异常集合中，基于最大可疑粗差个数qmax构造原始异常集合：
Figure BDA0003819656650000088
与原始异常集合对应的原始正常集合记为：
Figure BDA0003819656650000089
步骤3-2)由于观测值和参数之间特殊的对应关系，比如某些观测值并不包含所有参数的信息，导致在一些异常-正常集合的组合条件下，参数不是完全可解的，此时只需要从原始异常-正常集合的组合中筛选出满足参数可解条件的组合即可，所述参数可解条件为矩阵
Figure BDA00038196566500000810
或
Figure BDA00038196566500000811
可逆，其中，
Figure BDA00038196566500000812
是异常集合
Figure BDA00038196566500000813
对应的m×qmax阶且只包含0和1两种元素的设计矩阵，矩阵
Figure BDA00038196566500000814
的每一列代表异常集合
Figure BDA00038196566500000815
中的一个粗差；例如，如果异常集合
Figure BDA00038196566500000816
中包含第k 个观测值yk，那么矩阵
Figure BDA00038196566500000817
就有某一列中的第k个元素为1，该列的其他元素都为 0。
步骤3-3)从满足参数可解条件的组合中筛选出唯一一对最有可能的正常-异常集合作为初始正常-异常集合，完成初始化；
可以选择第j个异常集合
Figure BDA00038196566500000818
作为初始异常集合，其中，
当单位权中误差σ已知时，
Figure BDA0003819656650000091
当单位权中误差σ未知时，
Figure BDA0003819656650000092
步骤4)确定待检集合，并基于异常集合与待估参数、观测值残差的方差-协方差矩阵的最小二乘估值确定w统计量，筛选最小w统计量；
假设
Figure BDA0003819656650000093
作为初始的异常集合O(0)，它的补集即为正常集合I(0)。在接下来的迭代中，粗差的探测识别和检测可以看成是是一个从异常集合挑选正常观测值到正常集合中的过程。
步骤4-1)在第k次迭代中，将所有的观测值分为互无交集的三个集合：异常集合
Figure BDA0003819656650000094
正常集合
Figure BDA0003819656650000095
和待检集合{yi}，其中，yi∈O(k-1)，O(k-1)为在第k-1次迭代中保留的异常集合，异常集合
Figure BDA0003819656650000096
的元素个数为qmax-k，正常集合
Figure BDA0003819656650000097
的元素个数为m-qmax+k-1，待检集合{yi}的元素个数为1，三个集合同时满足以下两个关系条件：
Figure BDA0003819656650000098
Figure BDA0003819656650000099
步骤4-2)对第k次迭代的异常集合
Figure BDA00038196566500000910
中的每个观测值构造w统计量
Figure BDA00038196566500000911
当单位权中误差σ已知时：
Figure BDA00038196566500000912
当单位权中误差σ未知时：
Figure BDA0003819656650000101
其中，
Figure BDA0003819656650000102
Figure BDA0003819656650000103
是一个m×(qmax-k)阶设计矩阵，由异常集合
Figure BDA0003819656650000104
所生成，矩阵
Figure BDA0003819656650000105
的每一列代表集合
Figure BDA0003819656650000106
中的一个异常观测值，
Figure BDA0003819656650000107
是一个由待检集合{yi}所生成的m×1 阶的设计向量，
Figure BDA0003819656650000108
是一个由集合
Figure BDA0003819656650000109
生成m×(qmax-k+1)阶设计矩阵，即Hi＝[Cici]；
步骤4-3)通过筛选最小w统计量，确定在异常集合O(k-1)中最有可能正常的观察值；
当单位权中误差σ已知时：
Figure BDA00038196566500001010
当单位权中误差σ未知时：
Figure BDA00038196566500001011
步骤5)对最小w统计量进行显著性检验，若显著性检验通过，则将待检集合中的观测值并入正常集合，更新异常集合和正常集合，若显著性检验未通过，则结束迭代，转步骤7)；
步骤5-1)基于预配置的的显著水平α确定显著性检验阈值；
单位权中误差σ已知时，显著性检验阈值kα为：
kα＝Nα/2(0,1)
其中，Nα/2(0,1)表示标准正态分布分布在α/2处的分位值；
单位权中误差σ未知时，显著性检验阈值k′α为：
k′α＝tα/2(m-n-qmax+k-1,0)
其中，tα/2(m-n-qmax+k-1,0)表示自由度为m-n-qmax+k-1的学生分布在α/2处的分位值。
步骤5-2)比较最小w统计量的值与显著性检验阈值，
若最小w统计量的值小于等于显著性检验阈值，即
Figure BDA0003819656650000111
或者
Figure BDA0003819656650000112
则显著性检验通过，将待检集合{yi}中的观测值yj视为正常观测值，将其放入正常集合，更新正常集合和异常集合，即：
I(k)＝I(k-1)∪{yj}
O(k)＝O(k-1)-{yj}
若最小w统计量的值大于显著性检验阈值，则显著性检验未通过，异常集合中的观测值都被视为异常观测值，结束迭代。
步骤6)判断迭代次数是否等于最大可疑粗差个数，若是，此时异常集合中没有观测值，结束迭代，转步骤7)，若否，则返回步骤4)进行迭代；
步骤7)剔除异常集合中的值，基于正常集合与Gauss-Markov模型进行参数估计，求解待估参数，得到无粗差影响的GNSS高精度定位结果。
上述功能如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分或者该技术方案的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例所述方法的全部或部分步骤。而前述的存储介质包括：U盘、移动硬盘、只读存储器(ROM，Read-Only Memory)、随机存取存储器(RAM，Random Access Memory)、磁碟或者光盘等各种可以存储程序代码的介质。
本实施例通过一个GNSS单点定位的数值实验来验证本发明有效性。其设计矩阵A和方差阵∑为：
Figure BDA0003819656650000121
∑
＝diag([8.123 7.790 8.004 14.571 8.722 6.690 8.867 6.955 6.89311.026])
数值实验通过蒙特卡罗模拟(MCS)的方式进行。在每次模拟中，除了正态分布的随机误差外，还将在观测值中加入不同数量和大小的粗差
Figure BDA0003819656650000122
其中，加入粗差的观测值(异常值)数量q*为1到3个，即分别在)y1)，(y1,y2)和)y1,y2,y3)的观测值组合中加入粗差，粗差大小的范围为对应观测精度的1到20倍。为了公平起见，IDR和IDS使用同样的参数，其中最大可疑异常值数量qmax均设置为3，显著水平均取常用的0.1％。针对不同的应用场景，考虑了观测值单位权中误差σ已知和未知的两种情况。为对比和评估IDR和IDS两种方法的异常检测和识别的能力，下文统计分析了多重假设检验的各类决策概率水平和定位解算的均方根误差 (RMSE)。
图2给出了σ已知时IDR和IDS的正确决策概率。当仅有单个异常值时，IDR
和IDS的正确决策概率相当。当存在多个异常值时，IDR的正确决策概率高于IDS，且两种方法下的正确决策概率变化趋势差异显著。随着粗差的增大，IDS的正确决策概率先增后降，但始终保持在较低水平，这意味着IDS在面对多个异常值的时候难以做出完全正确的决策。而IDR的正确决策概率将随粗差的增大而显著提高，由此可见当σ已知时，对于多个异常值，IDR较IDS具有更强的检测和识别能力。
图3给出了一种GNSS案例中IDR、IDS和LS的定位RMSE。可以看出，仅有单个异常值时，LS的RMSE也会随着粗差的增大而呈指数增长，这意味着LS
对于少量的异常值也不具有抗差性。此时，IDR、IDS这两种异常值检测和识别方法都可以一定程度地降低定位误差。然而，当存在多个异常值时，IDS的RMSE
也会严重增大，甚至高于LS。这说明了虽然IDS已在各种实际应用中得到广泛应用，但实际上面对多个异常值时可能会失去抗差性。相比之下，使用IDR的RMSE 要小得多，并且始终可以控制在一定范围。由此可见，当σ已知时，相较于LS和 IDS，IDR对于多个异常值具有更强的抗差性。
接下来评估σ未知时IDR和IDS对异常值的检测和识别的能力。图4给出了σ未知时IDR和IDS的正确决策概率。可以发现，在单个异常值的情况下，IDR 和IDS的正确决策概率相当，且都随着粗差增大而先增加然后保持稳定。对于多个异常值的情况，由于严重的掩蔽和淹没效应，导致IDS的正确决策率极低，几乎接近于0。此时，IDR的正确决策概率仍然随着粗差的增大而保持上升趋势。这证明了在处理多个异常值时IDR比IDS具有更强的检测和识别能力。
图5进一步比较了σ未知时IDR、IDS和LS的定位RMSE。在单个异常值的情况下，相较于LS，两种方法都在很大程度上有效地降低了定位误差。同时，IDS 与IDR的精度相当。然而，对于多个异常值的情况，由于受到掩蔽和淹没效应的影响，IDS难以正确地检测和识别异常值，导致其抗差性严重下降，RMSE几乎与 LS一致。相反，IDR的RMSE显著较小，并且总是保持在一定范围内。这说明了当σ未知时，相较于LS和IDS，IDR对于多个异常值具有更强的抗差性。
以上详细描述了本发明的较佳具体实施例。应当理解，本领域的普通技术人员无需创造性劳动就可以根据本发明的构思做出诸多修改和变化。因此，凡本技术领域中技术人员依据本发明的构思在现有技术的基础上通过逻辑分析、推理、或者有限的实验可以得到的技术方案，皆应在权利要求书所确定的保护范围内。

Claims (10)

1.一种基于IDR粗差探测的GNSS定位方法，其特征在于，包括以下步骤：
步骤1)获取GNSS的观测值；
步骤2)基于Gauss-Markov模型求解待估参数和观测值残差的方差-协方差矩阵的最小二乘估值；
步骤3)基于预配置的最大可疑粗差个数初始化观测值的正常集合和异常集合；
步骤4)确定待检集合，并基于异常集合与待估参数、观测值残差的方差-协方差矩阵的最小二乘估值确定w统计量，筛选最小w统计量；
步骤5)对最小w统计量进行显著性检验，若显著性检验通过，则将待检集合中的观测值并入正常集合，更新异常集合和正常集合，若显著性检验未通过，则结束迭代，转步骤7)；
步骤6)判断是否满足迭代终止条件，若是，则转步骤7)，若否，则返回步骤4)进行迭代；
步骤7)剔除异常集合中的值，基于正常集合与Gauss-Markov模型进行参数估计，求解待估参数，得到无粗差影响的GNSS高精度定位结果。
2.根据权利要求1所述的一种基于IDR粗差探测的GNSS定位方法，其特征在于，所述步骤2)包括以下步骤：
步骤2-1)建立基于Gauss-Markov模型的GNSS线性观测方程：
y＝Ax+e，e～N(O，Σ)
式中，A为m×n阶非随机设计矩阵；x为n×1阶待估参数向量；y为m×1阶观测值向量；e为m×1阶偶然误差；Σ为m×m阶表征观测值精度的方差-协方差矩阵：
∑＝σ2Q
式中，Q为m×m阶协因数矩阵；σ2为协因数矩阵对应的单位权方差；
步骤2-2)确定待估参数的最小二乘解：
Figure FDA0003819656640000011
步骤2-3)确定观测值残差：
Figure FDA0003819656640000012
步骤2-4)确定待估参数和观测值残差的方差-协方差矩阵；
当单位权方差σ2已知时，
待估参数的方差-协方差矩阵为：
Figure FDA0003819656640000021
观测值残差的方差-协方差矩阵为：
Figure FDA0003819656640000022
当单位权方差σ2未知时，
待估参数的方差-协方差矩阵为
Figure FDA0003819656640000023
观测值残差的方差-协方差矩阵为
Figure FDA0003819656640000024
其中，
Figure FDA0003819656640000025
Figure FDA0003819656640000026
Figure FDA0003819656640000027
为单位权方差σ2的最小二乘估值：
Figure FDA0003819656640000028
3.根据权利要求2所述的一种基于IDR粗差探测的GNSS定位方法，其特征在于，所述步骤3)包括以下步骤：
步骤3-1)假定平差系统中所有的观测值形成一个全集，则每一个观测值被唯一地划分为到正常集合或异常集合中，基于最大可疑粗差个数qmax构造原始异常集合：
Figure FDA0003819656640000029
与原始异常集合对应的原始正常集合记为：
Figure FDA00038196566400000210
步骤3-2)从原始异常-正常集合的组合中筛选满足参数可解条件的组合，所述参数可解条件为矩阵
Figure FDA00038196566400000211
或
Figure FDA00038196566400000212
可逆，其中，
Figure FDA00038196566400000213
是异常集合
Figure FDA00038196566400000214
对应的m×qmax阶且只包含0和1两种元素的设计矩阵，矩阵
Figure FDA00038196566400000215
的每一列代表异常集合
Figure FDA00038196566400000216
中的一个粗差；
步骤3-3)从满足参数可解条件的组合中筛选出唯一一对正常-异常集合作为初始正常-异常集合，完成初始化。
4.根据权利要求3所述的一种基于IDR粗差探测的GNSS定位方法，其特征在于，步骤3-3)筛选得到的初始正常-异常集合为第j个异常集合
Figure FDA0003819656640000031
当单位权中误差σ已知时，
Figure FDA0003819656640000032
当单位权中误差σ未知时，
Figure FDA0003819656640000033
5.根据权利要求4所述的一种基于IDR粗差探测的GNSS定位方法，其特征在于，所述步骤4)包括以下步骤：
步骤4-1)在第k次迭代中，将所有的观测值分为互无交集的三个集合：异常集合
Figure FDA0003819656640000034
正常集合
Figure FDA0003819656640000035
和待检集合{yi}，其中，yi∈O(k-1)，O(k-1)为在第k-1次迭代中保留的异常集合，异常集合
Figure FDA0003819656640000036
的元素个数为qmax-k，正常集合
Figure FDA0003819656640000037
的元素个数为m-qmax+k-1，待检集合{yi}的元素个数为1，三个集合同时满足以下两个关系条件：
Figure FDA0003819656640000038
Figure FDA0003819656640000039
步骤4-2)对第k次迭代的异常集合
Figure FDA00038196566400000310
中的每个观测值构造w统计量
Figure FDA00038196566400000311
当单位权中误差σ已知时：
Figure FDA00038196566400000312
当单位权中误差σ未知时：
Figure FDA0003819656640000041
其中，
Figure FDA0003819656640000042
Figure FDA0003819656640000043
是一个m×(qmax-k)阶设计矩阵，由异常集合
Figure FDA0003819656640000044
所生成，矩阵
Figure FDA0003819656640000045
的每一列代表集合
Figure FDA0003819656640000046
中的一个异常观测值，
Figure FDA0003819656640000047
是一个由待检集合{yi}所生成的m×1阶的设计向量，
Figure FDA0003819656640000048
是一个由集合
Figure FDA0003819656640000049
生成m×(qmax-k+1)阶设计矩阵，即Hi＝[Ci ci]；
步骤4-3)筛选最小w统计量：
当单位权中误差σ已知时：
Figure FDA00038196566400000410
当单位权中误差σ未知时：
Figure FDA00038196566400000411
6.根据权利要求5所述的一种基于IDR粗差探测的GNSS定位方法，其特征在于，所述步骤5)包括以下步骤：
步骤5-1)基于预配置的的显著水平α确定显著性检验阈值；
步骤5-2)比较最小w统计量的值与显著性检验阈值，
若最小w统计量的值小于等于显著性检验阈值，则显著性检验通过，将待检集合{yi}中的观测值yj视为正常观测值，将其放入正常集合，更新正常集合和异常集合，即：
I(k)＝I(k-1)∪{yj}
O(k)＝O(k-1)-{yj}
若最小w统计量的值大于显著性检验阈值，则显著性检验未通过，异常集合中的观测值都被视为异常观测值，结束迭代。
7.根据权利要求6所述的一种基于IDR粗差探测的GNSS定位方法，其特征在于，单位权中误差σ已知时，显著性检验阈值kα为：
kα＝Nα/2(0,1)
其中，Nα/2(0,1)表示标准正态分布分布在α/2处的分位值；
单位权中误差σ未知时，显著性检验阈值k′α为：
k′α＝tα/2(m-n-qmax+k-1,0)
其中，tα/2(m-n-qmax+k-1,0)表示自由度为m-n-qmax+k-1的学生分布在α/2处的分位值。
8.根据权利要求1所述的一种基于IDR粗差探测的GNSS定位方法，其特征在于，所述迭代终止条件为迭代次数等于最大可疑粗差个数。
9.一种基于IDR粗差探测的GNSS定位装置，包括存储器、处理器，以及存储于所述存储器中的程序，其特征在于，所述处理器执行所述程序时实现如权利要求1-8中任一所述的方法。
10.一种存储介质，其上存储有程序，其特征在于，所述程序被执行时实现如权利要求1-8中任一所述的方法。
CN202211039658.4A 2022-08-29 2022-08-29 一种基于idr粗差探测的gnss定位方法、装置及存储介质 Active CN115390114B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211039658.4A CN115390114B (zh) 2022-08-29 2022-08-29 一种基于idr粗差探测的gnss定位方法、装置及存储介质

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211039658.4A CN115390114B (zh) 2022-08-29 2022-08-29 一种基于idr粗差探测的gnss定位方法、装置及存储介质

Publications (2)

Publication Number Publication Date
CN115390114A true CN115390114A (zh) 2022-11-25
CN115390114B CN115390114B (zh) 2024-11-26

Family

ID=84123484

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211039658.4A Active CN115390114B (zh) 2022-08-29 2022-08-29 一种基于idr粗差探测的gnss定位方法、装置及存储介质

Country Status (1)

Country Link
CN (1) CN115390114B (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117991308A (zh) * 2024-04-03 2024-05-07 长江三峡集团实业发展(北京)有限公司 一种面向复杂环境gnss数据质量的高精度数据处理方法
CN119471742A (zh) * 2025-01-13 2025-02-18 深圳华大北斗科技股份有限公司 基于w检验的M估计GNSS质量控制的定位方法

Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2012079616A1 (en) * 2010-12-13 2012-06-21 Telecom Italia S.P.A. Method and system for localizing mobile communications terminals
WO2014026074A2 (en) * 2012-08-09 2014-02-13 Bae Systems Information And Electronic Systems Integration Inc. Integrated data registration
CN103809161A (zh) * 2014-01-09 2014-05-21 中国人民解放军海军航空工程学院 雷达网抗距离欺骗+soj复合干扰方法
CN107180450A (zh) * 2017-06-06 2017-09-19 广西师范学院 一种基于dem的河谷横断面形态的算法
CN108919321A (zh) * 2018-05-18 2018-11-30 长安大学 一种基于尝试法的gnss定位粗差探测方法
CN109298389A (zh) * 2018-08-29 2019-02-01 东南大学 基于多粒子群优化的室内行人组合位姿估计方法
WO2019228439A1 (zh) * 2018-06-01 2019-12-05 浙江亚特电器有限公司 基于gnss-rtk的定位方法
CN110632625A (zh) * 2019-08-19 2019-12-31 中国矿业大学 一种gnss时间序列阶跃探测与修复方法
CN114117802A (zh) * 2021-11-29 2022-03-01 同济大学 一种基于极大后验估计的多粗差探测方法、装置及介质

Patent Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2012079616A1 (en) * 2010-12-13 2012-06-21 Telecom Italia S.P.A. Method and system for localizing mobile communications terminals
WO2014026074A2 (en) * 2012-08-09 2014-02-13 Bae Systems Information And Electronic Systems Integration Inc. Integrated data registration
CN103809161A (zh) * 2014-01-09 2014-05-21 中国人民解放军海军航空工程学院 雷达网抗距离欺骗+soj复合干扰方法
CN107180450A (zh) * 2017-06-06 2017-09-19 广西师范学院 一种基于dem的河谷横断面形态的算法
CN108919321A (zh) * 2018-05-18 2018-11-30 长安大学 一种基于尝试法的gnss定位粗差探测方法
WO2019228439A1 (zh) * 2018-06-01 2019-12-05 浙江亚特电器有限公司 基于gnss-rtk的定位方法
CN109298389A (zh) * 2018-08-29 2019-02-01 东南大学 基于多粒子群优化的室内行人组合位姿估计方法
CN110632625A (zh) * 2019-08-19 2019-12-31 中国矿业大学 一种gnss时间序列阶跃探测与修复方法
CN114117802A (zh) * 2021-11-29 2022-03-01 同济大学 一种基于极大后验估计的多粗差探测方法、装置及介质

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
喻杨康 等: "Bernoulli–Gaussian Model with Model Parameter Estimation", JOURNAL OF SURVEYING ENGINEERING, vol. 150, no. 4, 24 August 2024 (2024-08-24) *
孙志鹏;: "基于验后方差原理的选择权迭代法定位粗差", 北京测绘, no. 04, 25 August 2016 (2016-08-25) *

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117991308A (zh) * 2024-04-03 2024-05-07 长江三峡集团实业发展(北京)有限公司 一种面向复杂环境gnss数据质量的高精度数据处理方法
CN117991308B (zh) * 2024-04-03 2024-06-07 长江三峡集团实业发展(北京)有限公司 一种面向复杂环境gnss数据质量的高精度数据处理方法
CN119471742A (zh) * 2025-01-13 2025-02-18 深圳华大北斗科技股份有限公司 基于w检验的M估计GNSS质量控制的定位方法

Also Published As

Publication number Publication date
CN115390114B (zh) 2024-11-26

Similar Documents

Publication Publication Date Title
CN107103171B (zh) 机器学习模型的建模方法及装置
US20190166024A1 (en) Network anomaly analysis apparatus, method, and non-transitory computer readable storage medium thereof
CN111723865B (zh) 评估图像识别模型、攻击方法性能的方法、装置和介质
CN110717824A (zh) 基于知识图谱的银行对公客群风险传导测算的方法及装置
Dudczyk Radar emission sources identification based on hierarchical agglomerative clustering for large data sets
CN104869126B (zh) 一种网络入侵异常检测方法
CN115390114A (zh) 一种基于idr粗差探测的gnss定位方法、装置及存储介质
Sun et al. Quantifying variable interactions in continuous optimization problems
CN110852755A (zh) 针对交易场景的用户身份识别方法和装置
CN110570312B (zh) 样本数据获取方法、装置、计算机设备和可读存储介质
CN112437053B (zh) 入侵检测方法及装置
CN108919059A (zh) 一种电网故障诊断方法、装置、设备及可读存储介质
CN114662602B (zh) 一种离群点检测方法、装置、电子设备及存储介质
Evans et al. Outlier identification in model-based cluster analysis
CN117376228B (zh) 一种网络安全测试工具确定方法及装置
CN112632609A (zh) 异常检测方法、装置、电子设备及存储介质
CN110008972A (zh) 用于数据增强的方法和装置
Mameche et al. Learning causal models under independent changes
Lim et al. More powerful selective kernel tests for feature selection
CN115600105A (zh) 基于mic-lstm的水体缺失数据插补方法及装置
CN105989095B (zh) 顾及数据不确定性的关联规则显著性检验方法及装置
CN118152567A (zh) 一种员工画像异常检测方法及终端
CN111209567A (zh) 提高检测模型鲁棒性的可知性判断方法及装置
CN117235434B (zh) 林业碳汇项目基线构建方法、系统、终端及介质
CN115267704B (zh) 运动目标rcs模拟数据的生成方法、装置、设备及介质

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