CN113281793A - 单历元位置界限的方法和设备 - Google Patents
单历元位置界限的方法和设备 Download PDFInfo
- Publication number
- CN113281793A CN113281793A CN202110138235.7A CN202110138235A CN113281793A CN 113281793 A CN113281793 A CN 113281793A CN 202110138235 A CN202110138235 A CN 202110138235A CN 113281793 A CN113281793 A CN 113281793A
- Authority
- CN
- China
- Prior art keywords
- probability density
- model
- state
- parameter
- distribution
- 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.)
- Pending
Links
- 238000000034 method Methods 0.000 title claims abstract description 108
- 238000005259 measurement Methods 0.000 claims abstract description 100
- 238000013442 quality metrics Methods 0.000 claims abstract description 41
- 238000009826 distribution Methods 0.000 claims description 56
- 238000005070 sampling Methods 0.000 claims description 33
- 230000006870 function Effects 0.000 claims description 29
- 241001061260 Emmelichthys struhsakeri Species 0.000 claims description 13
- 230000008859 change Effects 0.000 claims description 13
- 230000008569 process Effects 0.000 claims description 9
- 238000005315 distribution function Methods 0.000 claims description 8
- 238000007476 Maximum Likelihood Methods 0.000 claims description 7
- 230000003993 interaction Effects 0.000 claims description 5
- 230000010354 integration Effects 0.000 claims description 4
- 230000001419 dependent effect Effects 0.000 claims description 3
- 230000000644 propagated effect Effects 0.000 claims description 3
- 230000001902 propagating effect Effects 0.000 claims description 2
- 238000012937 correction Methods 0.000 description 17
- 238000010586 diagram Methods 0.000 description 11
- 230000000694 effects Effects 0.000 description 10
- 238000004364 calculation method Methods 0.000 description 8
- 238000012545 processing Methods 0.000 description 8
- 238000013178 mathematical model Methods 0.000 description 5
- 239000000306 component Substances 0.000 description 4
- 238000010606 normalization Methods 0.000 description 4
- 238000005457 optimization Methods 0.000 description 4
- 238000012360 testing method Methods 0.000 description 4
- 238000013461 design Methods 0.000 description 3
- 239000011159 matrix material Substances 0.000 description 3
- 238000002156 mixing Methods 0.000 description 3
- 230000010287 polarization Effects 0.000 description 3
- 238000013459 approach Methods 0.000 description 2
- 230000001934 delay Effects 0.000 description 2
- 230000002708 enhancing effect Effects 0.000 description 2
- 239000000463 material Substances 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 230000002829 reductive effect Effects 0.000 description 2
- 230000002441 reversible effect Effects 0.000 description 2
- 230000002159 abnormal effect Effects 0.000 description 1
- 230000002411 adverse Effects 0.000 description 1
- 230000003321 amplification Effects 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 238000003491 array Methods 0.000 description 1
- 230000003416 augmentation Effects 0.000 description 1
- 238000013398 bayesian method Methods 0.000 description 1
- 230000002902 bimodal effect Effects 0.000 description 1
- 230000001413 cellular effect Effects 0.000 description 1
- 238000004891 communication Methods 0.000 description 1
- 239000008358 core component Substances 0.000 description 1
- 238000001914 filtration Methods 0.000 description 1
- 239000004973 liquid crystal related substance Substances 0.000 description 1
- 238000013507 mapping Methods 0.000 description 1
- 239000000203 mixture Substances 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 230000005404 monopole Effects 0.000 description 1
- 238000003199 nucleic acid amplification method Methods 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 230000036961 partial effect Effects 0.000 description 1
- 238000012805 post-processing Methods 0.000 description 1
- 239000004065 semiconductor Substances 0.000 description 1
- 230000003068 static effect Effects 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 230000001131 transforming effect Effects 0.000 description 1
- 230000001960 triggered effect Effects 0.000 description 1
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/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/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
-
- 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/20—Integrity monitoring, fault detection or fault isolation of space segment
-
- 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/24—Acquisition or tracking or demodulation of signals transmitted by the system
- G01S19/243—Demodulation of navigation message
Landscapes
- Engineering & Computer Science (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Computer Networks & Wireless Communication (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Computer Security & Cryptography (AREA)
- Signal Processing (AREA)
- Position Fixing By Use Of Radio Waves (AREA)
Abstract
本发明涉及单历元位置界限的方法和设备。一种使用GNSS测量的单个历元来确定位置估计的保护级别的方法,该方法包括:指定状态x的先验概率密度P(x);指定将状态x与测量的可观测量z联系起来的***模型h(x);将与测量关联的质量度量q量化;指定非高斯残留误差概率密度模型f(r|θ,q)并使用一组实验数据对模型参数θ进行拟合;定义后验概率密度P(x|z,q,θ);估计状态x;以及通过在状态x上对后验概率密度P(x|z,q,θ)进行积分来计算保护级别。
Description
技术领域
本公开总体上涉及基于卫星的导航,更具体地,涉及使用全球导航卫星***(GNSS)观测数据的单个历元和/或来自其它相关传感器的测量来确定位置估计的保护级别或者与确定有关的参数的方法及其操作设备。
背景技术
GNSS接收器接收从一个或多个GNSS卫星群发送的卫星信号并且使用包含在卫星信号中的信息来估计其位置等等。通常,由于例如卫星信号所遭受的多径干扰,此估计无法单独提供足够的准确性和可靠性。提供误差超过可容范围的位置估计的后果会是显著的。因此,需要确定与位置估计关联的误差上的概率界限。
可通过对一组GNSS可观测量的误差进行建模来确定与位置估计关联的误差上的概率界限。然而,当可观测量彼此不独立时,可观测量误差的建模是困难的,导致难以确定与位置估计关联的误差上的概率界限。
发明内容
根据本公开的一些实施方式,提供了一种使用全球导航卫星***(GNSS)测量的单个历元来确定位置估计的保护级别的方法。该方法包括:预先指定状态x的先验概率密度P(x);预先指定将状态x与测量的可观测量z联系起来的***模型h(x);在操作期间,将与测量关联的质量度量q量化;预先指定非高斯残留误差概率密度模型f(r|θ,q)并对先验离线确定的模型参数θ进行拟合;在操作期间,定义后验状态概率密度P(x|z,q,θ);在操作期间,估计状态x;以及在操作期间,通过在状态x上对后验概率密度P(x|z,q,θ)进行积分来计算保护级别。
可选地,***模型可被指定为z=h(x)+r,其中r是残留误差。
在一些布置方式中,将与测量关联的质量度量q量化可包括使用随机样本共识(RANSAC)技术来排除测量的一个或更多个异常值。
非高斯残留误差概率密度模型可以是student-t分布函数。指定残留误差概率密度模型f(r|θ,q)还可包括使用student-t分布函数将伪距的模型识别为:
其中r是残留误差,v是自由度参数,σ是缩放参数,并且Γ是伽马函数。
可选地,非高斯残留误差概率密度模型可以是student-t分布函数。指定残留误差概率密度模型f(r|θ,q)还可包括使用student-t分布函数将载波相位的模型识别为:
其中r是残留误差,ν是自由度参数,σ是缩放参数,w是权重,并且Γ是伽马函数。
使用一组实验数据对模型参数θ进行拟合可包括:使用无约束student-t参数定义自由度参数和缩放参数;将无约束student-t参数定义为质量度量的函数;以及寻找无约束student-t参数的最大似然值。
在示例性方法中,使用一组实验数据对模型参数θ进行拟合可包括:使用无约束参数来定义自由度参数和缩放参数和均匀权重参数;将无约束参数定义为质量度量的函数;以及寻找无约束参数的最大似然值。
定义后验概率密度可基于下式:
其中i是自然数。
可选地,在状态上对后验概率密度进行积分可包括:使用具有修改的密度函数的多个交互链来应用马尔科夫链蒙特卡洛(MCMC)数值积分。
在一些布置方式中,在状态上对后验概率密度进行积分可包括:通过从无偏分布提取第一样本集合并在第一样本集合上估计分位数的第一位置来执行第一轮采样;通过从基于分位数的第一位置确定的第一约束分布提取第二样本集合来执行第二轮采样;估计第二样本集合上的分位数的第二位置,其中,分位数的第二位置用于为第三轮采样设定第二约束分布;使用在先前采样中设定的约束分布来重复采样,直至第n轮采样完成,其中n是大于2的自然数;确定n轮采样的样本是否足以计算保护级别;以及当n轮样本足以计算保护级别时,将n轮采样的样本与适当的权重组合。
该方法还可包括:使用由一个或更多个传感器获得的δ相误差分布或数据使最近的单历元界限向前传播到当前时间。
传播还可包括:使用非高斯分布函数对δ相上的误差分布进行建模并且使用马尔科夫链蒙特卡洛技术数值上评估δ位置上的界限;或者使用高斯过界方法对δ位置上的误差分布进行建模。
可选地,状态x可被定义为:
其中p是流动站的位置,c是真空中的光速,Δdtr是流动站与基站之间的接收器时钟偏差的差,ΔkP,r是伪距上的接收器代码仪器延迟的差,ΔkL,r是载波相位上的接收器代码仪器延迟的差,并且ΔTzenith是天顶对流层延迟的差。
根据本公开的一些实施方式,提供了一种量化GNSS测量的质量的方法。该方法包括以下步骤:在围绕感兴趣的历元的窗口中执行GNSS测量;确定GNSS测量的改变;识别位置和时钟偏差的改变的共识解;以及通过与共识的偏离来表征GNSS测量的质量。GNSS测量可包括相位测量、伪距测量和多普勒测量中的至少一个。
可选地,识别位置和时钟偏差的改变的共识解还可包括:使用表示相对于窗口的第一历元的位置和时钟偏差偏移的***来对GNSS测量数据进行建模;以及将该***作为加权线性最小二乘问题求解。
将***作为加权线性最小二乘问题求解还可包括引入残留误差上的损失函数作为:
其中ρ(s)是损失函数,wi是各个测量的权重,ri是各个测量的残留误差,i是自然数。
可选地,损失函数可以是Huber损失函数:
其中δ表示阈值参数。
根据本公开的一些实施方式,还提供了一种被配置为使用GNSS测量的单个历元来确定位置估计的保护级别的设备。该设备包括:接收器,其被配置为接收GNSS信号并处理所接收的信号以生成测量;以及处理器,其被配置为:指定状态x的先验概率密度P(x);指定将状态x与测量的可观测量z联系起来的***模型h(x);将与测量关联的质量度量q量化;指定非高斯残留误差概率密度模型f(r|θ,q)并对先验离线确定的模型参数θ进行拟合;定义后验概率密度P(x|z,q,θ);例如使用后验概率密度P(x|z,q,θ)来估计状态x;以及通过在状态x上对后验概率密度P(x|z,q,θ)进行积分来计算保护级别。
可选地,处理器还可被配置为与外部传感器通信以获得传感器所跟踪的位置改变数据以确定界限传播。
根据本公开的一些实施方式,还提供了一种存储有指令的非暂时性计算机可读介质,所述指令在由处理器执行时执行一种使用GNSS测量的单个历元来确定位置估计的保护级别的方法,该方法包括以下步骤:指定状态x的先验概率密度P(x);指定将状态x与测量的可观测量z联系起来的***模型h(x);将与测量关联的质量度量q量化;指定非高斯残留误差概率密度模型f(r|θ,q)并对先验离线确定的模型参数θ进行拟合;定义后验概率密度P(x|z,q,θ);例如使用后验概率密度P(x|z,q,θ)来估计状态x;以及通过在状态x上对后验概率密度P(x|z,q,θ)进行积分来计算保护级别。
附图说明
图1是与本公开的一些实施方式一致的装置的框图。
图2是示出与本公开的一些实施方式一致的确定位置估计的保护级别的方法的流程图。
图3是示出与本公开的一些实施方式一致的识别GNSS测量的异常值的方法的流程图。
图4是示出与本公开的一些实施方式一致的量化与GNSS测量关联的质量度量的方法的流程图。
图5示出与本公开的一些实施方式一致的观测的残留误差分布、拟合到观测的残留误差分布的高斯分布和拟合到观测的残留误差分布的student-t分布的比较。
图6是示出与本公开的一些实施方式一致的使用并行交互链的多峰采样方法的示意图。
图7是示出与本公开的一些实施方式一致的执行伞形采样的方法的流程图。
图8A和图8B示出显示与本公开的一些实施方式一致的在顺轨方向上从实验数据获得的界限半宽和真实误差的斯坦福图(Stanford plot),图8C和图8D示出显示在跨轨维度上从实验数据获得的界限半宽和真实误差的斯坦福图,图8E和图8F示出显示在垂直维度上从实验数据获得的界限半宽和真实误差的斯坦福图。
具体实施方式
现在将详细参考示例性实施方式,其示例示出于附图中。以下描述参考了附图,其中除非另外表示,否则不同附图中的相同标号表示相同或相似的元件。示例性实施方式的以下描述中阐述的实现方式并非表示与本公开一致的所有实现方式。相反,它们仅是与所附权利要求中陈述的本公开相关的方面一致的***、设备和方法的示例。
GNSS接收器通过天线接收从一个或多个GNSS卫星群发送的卫星信号并且使用包含在卫星信号中的信息来估计其位置。许多GNSS应用需要所估计的位置的特定级别的准确性和可靠性。例如对于依赖于可信赖的全球位置信息的自主车辆,提供误差超过可容范围的位置估计(或者与估计有关的任何其它参数)的后果会是显著的。
与沿着特定方向的位置估计关联的误差可使用该方向上的“保护级别”的概念来量化。保护级别可被定义为位置估计的误差上的概率界限,指定约束区域并且建议真实位置在该区域之外的可能性低于给定值。该区域可按照间隔、圆圈或矩形或任何其它形状的形式提供。例如,可作为围绕GNSS接收器所估计的水平位置的圆提供保护级别。保护级别可用于确定何时警示用户接收器无法以足够的准确性提供位置估计。例如,如果保护级别超过预定警示极限,则可触发警报,以使得客户端***可进行适当的调节。当***位置估计不可靠时***向用户提供警告的能力可使用“完整性”的概念来描述。如果保护级别没有超过警示极限,则用户可对位置估计的准确性有信心。
从一组GNSS可观测量(例如,伪距或载波相位或多普勒)确定保护级别是主要取决于可观测量的误差概率分布的统计问题。给定可观测量的一组已知误差分布,贝耶斯方法可用于计算后验概率密度(即,位置上的概率密度)。然后,可通过对后验概率密度进行积分来确定保护级别。
然而,当可观测量统计上彼此不独立时,需要在联合分布中对可观测量误差进行建模。这使得误差分布更加难以建模,并且使得后验概率密度更加难以确定,从而导致难以计算保护级别。
本公开的实施方式提供了一种使用GNSS测量的单个历元来确定位置估计的保护级别的方法。历元可被表示为在单个时间点来自GNSS接收器的一组测量。在实施方式中,指定状态的先验概率密度;指定将状态与测量的可观测量联系起来的***模型;并且基于先验离线确定的拟合模型参数来指定非高斯测量误差概率密度模型。在实施方式中,使用基于窗口的方法或RANSAC技术,推导质量度量并且将测量的异常值排除或减权重。使用马尔科夫链蒙特卡洛(MCMC)或重要性采样方法,定义后验概率密度并且通过对后验概率密度进行积分来确定与位置估计关联的保护级别。在实施方式中,使用δ相误差分布或通过传感器获得的数据来实现界限传播。
本文所公开的实施方式具有一个或更多个技术效果。通过考虑仅来自GNSS测量的单个历元的可观测量,去除主要的非独立性(例如,时间相关性)来源,并且各个信号上的误差分布被单独地建模为一维概率密度函数,从而导致快速且严谨地确定保护级别。使用非高斯误差概率密度模型允许适当地考虑尾概率密度,从而确保处理的准确性。将测量的异常值排除或减权重允许确定的准确性和一致性增强。使用具有修改的密度函数的多个MCMC链对后验概率密度进行积分允许考虑多峰分布上的尾概率密度,从而进一步增强处理的准确性。执行界限传播允许***维持高可用性和低延迟。
图1是与本公开的一些实施方式一致的装置100的框图。参照图1,装置100可被安装在移动车辆中或固定位置中。装置100可采取任何形式,包括(但不限于)膝上型计算机、GNSS接收器、包括移动电话、无线手持装置或无线个人装置的无线终端或者任何其它形式。装置100包括天线102、联接到天线102的接收器104、处理器106、存储器108、本地时钟110和输入/输出装置112。
天线102被配置为接收GNSS信号。GNSS信号可从单个卫星接收,或者是分别在多个频带中从多个卫星发送的多个卫星信号。GNSS信号还可包括源自反射和/或散射卫星信号的一个或更多个虚拟源的信号。然而,天线102所接收的信号不限于卫星信号,可以是从任何源发送的任何电磁波(例如,无线蜂窝信号)。天线102可以是任何类型的天线,例如贴片天线、螺旋天线、十字弓形天线、正交放置的单极天线等。天线102可以是天线阵列。
接收器104联接到天线102并且被配置为经由天线102接收GNSS信号。由天线102接收的GNSS信号可经由同轴RF线缆或适合于发送RF信号的任何其它线缆发送到接收器104。在实施方式中,接收器104可以是收发器调制解调器的一部分,其包括被配置为向外部装置发送数据的发送器。
接收器104可对所接收的GNSS信号执行诸如放大、滤波、混合和数字化的操作。接收器104还可处理数字化的信号并生成诸如伪距测量和载波相位测量的测量。接收器104还可确定测量质量指标(也称为质量度量)。质量度量可提供关于执行测量的环境的质量的信息。质量度量的示例包括(但不限于)载波噪声密度比、卫星仰角、锁相时间、锁码时间和周跳。接收器104所进行的测量的准确性可与测量所关联的质量度量相关。接收器104还被配置为与时钟110通信。
处理器106可包括一个或更多个专用处理单元、专用集成电路(ASIC)、现场可编程门阵列(FPGA)或者各种其它类型的处理器或处理单元。处理器106可从接收器104接收伪距和载波相位测量的结果以及与测量关联的质量度量,并且还处理信息以估计装置100的当前位置和位置估计的保护级别。处理器106还可提供其自己的测量质量度量。在实施方式中,处理器106可包括导航滤波器(例如,卡尔曼滤波器或递归LS滤波器)以用于确定保护级别。处理器106可联接到诸如加速度计、陀螺仪或车轮速度传感器的外部运动传感器,以从运动传感器获得位置跟踪信息并且有助于确定保护级别。处理器106还被配置为与输入/输出装置112和存储器108通信。
在实施方式中,接收器104可包括执行处理器106的全部或部分功能的内置处理器(未示出)。在实施方式中,接收器104的内置处理器可以是控制接收器104中的信号处理的前端处理器,并且处理器106可以是基于接收器104中的信号处理执行进一步计算的后端处理器。在实施方式中,处理器106可向远程计算机(未示出)指派计算任务,以使得远程计算机执行一部分计算并将计算结果发送到处理器106。
存储器108可以是任何类型的计算机可读存储介质,包括易失性或非易失性存储器装置或其组合。存储器108可存储与装置100的标识和天线102所接收的GNSS信号有关的信息。存储器108还可存储后处理信号。存储器108还可存储伪距和载波相位测量以及与这些测量和来自其它传感器的测量关联的质量度量。存储器108还可存储用于接收器104中的信号处理和处理器106中的计算的计算机可读程序指令、数学模型和算法。存储器108还可存储计算机可读程序指令以供处理器106执行以操作装置100。
本地时钟110提供装置100所处的本地的时间。本地时钟110可用于确定用于伪距测量和位置估计的信号的到达时间。
输入/输出装置112可用于将信号处理和计算的结果通信给用户或另一装置。输入/输出装置112可包括用户接口,其包括显示器和输入装置以将用户命令发送到处理器106。显示器可被配置为显示装置100处的信号接收状态、存储器108处存储的数据、信号处理状态和计算结果等。例如,显示器可向用户显示确定保护级别的结果,以使得用户可更好地理解装置100的位置。显示器可包括(但不限于)阴极射线管(CRT)、液晶显示器(LCD)、发光二极管(LED)、气体等离子体显示器、触摸屏或用于将信息显示给用户的其它图像投影装置。输入装置可以是用于接收来自用户的数据和控制信号的任何类型的计算机硬件设备。输入装置可包括(但不限于)键盘、鼠标、扫描仪、数字相机、操纵杆、轨迹球、光标方向键、触摸屏监视器或音频/视频命令器等。输入/输出装置112还可包括诸如电总线连接或无线通信链路的机器接口。
在实施方式中,用于确定位置估计的保护级别的数学模型和算法可基于贝叶斯定理。
贝叶斯定理提供了一种基于与事件有关的已知先验信息来确定该事件的概率密度的方式。贝叶斯定理可利用下式来表示:
其中P(A|B)是假定B为真时发生事件A的条件概率,P(B|A)是假定A为真时发生事件B的条件概率,P(A)和P(B)分别是独立地考虑时A和B的概率。
当给定一组连续域测量应用于推理问题时,贝叶斯公式可被转换为下式:
其中P(state|data)是给定一组观测的特定状态的概率密度,P(data|state)是给定状态的特定一组观测的概率密度,P(state)是状态的先验概率密度,P(data)是一组观测的概率密度。关于状态的信息可从可观测量(例如,伪距或载波相位)推断,并且P(state|data)对应于为了计算位置估计的保护级别而需要确定的后验概率密度。P(data)可被当作未知归一化因子,并且可利用所有状态上的后验概率密度P(state|data)的积分等于1的事实来推断。P(data|state)与测量误差概率分布密切相关并且可如下所述由数学模型指定。
在实施方式中,对GNSS测量应用上式(2)。状态(x)可包括(但不限于)位置、时钟偏差、仪器偏差或大气参数。时钟偏差可包括接收器时钟偏差和卫星时钟偏差。GNSS测量可包括诸如伪距测量和载波相位测量的一组观测。测量还可包括一组质量度量q(元数据),其指示关于观测质量的信息,例如载波噪声密度比。
一组观测z可被分成两个分量:将状态x与一组观测z中的可观测量联系起来的***模型h(x)以及随机测量误差分量(残差)r:
z=h(x)+r; (3)
残差r的概率密度可由函数fjoint(r|θ,q)定义,其在给定误差模型参数θ和一组质量度量q下,指定一组指定的残差r的概率密度。返回式(2),后验概率密度P(x|z,q,θ)可被表示为:
以分母P(z|q,θ)为未知归一化常数,并且r=z-h(x),后验概率密度P(x|z,q,θ)可如下表示:
P(x|z,q,θ)∝P(z|x,q,θ)P(x) (5)
P(x|z,q,θ)∝fjoint(z-h(x)|θ,q)P(x) (6)
在实施方式中,在确定后验概率密度P(x|z,q,θ)时,考虑来自单个测量历元的可观测量。对于单个测量历元,由于测量之间不存在时间相关性,所以测量可被视为独立的。因此,联合概率密度fjoint(r|θ,q)可被计算为各个观测的独立概率密度f(r|θ,q)的简单乘积:
fjoint(z-h(x)|θ,q)=∏if(zi-hi(x)|θ,qi) (7)
因此,式(6)可如下简化:
P(x|z,q,θ)∝P(x)∏if(zi-hi(x)|θ,qi) (8)
这样,测量误差分布被简化为残差r的一维函数f(r|θ,q),该误差分布可使用实验数据来确定。还可例如通过在状态x上对后验概率密度P(x|z,q,θ)进行积分来从后验概率密度P(x|z,q,θ)进一步确定位置估计的保护级别。
图2是示出与本公开的一些实施方式一致的确定位置估计的保护级别的方法200的流程图。该方法可由诸如图1的装置100的装置执行。参照图2,方法200包括指定关于状态x的先验信息的步骤S210。例如,可在进行测量之前使用关于状态x的知识来指定状态x的先验概率密度P(x)。在实施方式中,指定先验信息可在测量之前执行。
方法200包括指定将状态x与测量的可观测量联系起来的***模型h(x)的步骤S220。可观测量可包括伪距R和载波相位Φ。
在实施方式中,伪距R的数学模型可表示如下:
R=ρ+c(dtr-dts)+T+αfSTEC+kP,r+kP,s+mP+∈P (9)
其中ρ是几何距离,c是真空中的光速,dtr和dts是与GNSS时间标度的接收器和卫星时钟偏差,T是对流层延迟,αfSTEC是电离层延迟(STEC是沿着信号路径对电子密度进行积分的倾斜总电子计数,αf是将STEC与延迟联系起来的频率相关转换因子),kP,r和kP,s是取决于代码和频率的接收器和卫星代码仪器延迟,mP是多径干扰对伪距的影响,∈P是接收器噪声对伪距的影响。
在实施方式中,载波相位Φ的数学模型可表示如下:
Φ=ρ+c(dtr-dts)+T-αfSTEC+kL,r-kL,s+λn+λw+mL+∈L (10)
除了如上定义的项之外,kL,r和kL,s分别是接收器和卫星代码仪器延迟,mL是多径干扰对载波相位的影响,λ是所发送的GNSS信号的波长,n是整周模糊度,w是圆偏振缠绕(wind-up),∈L是接收器噪声对载波相位的影响。载波相位测量可以是传统单频测量或者从不同频率的信号获得的载波相位测量的组合(例如,宽波组合)。
在实施方式中,可使用从一个或多个物理或虚拟基站或校正服务或增强***获得的数据来计算对伪距R的校正和对载波相位Φ的校正。具体地考虑来自基站的校正,基站可以是准确位置已知的参考装置。该校正可去除卫星偏差项并减小电离层和对流层延迟的影响,并且减轻时钟和轨道中的误差。基于基站数据的校正可表示如下:
Rcorrection=Rbase-ρpredicted (11)
Φcorrection=Φbase-ρpredicted/λ (12)其中Rbase和Φbase分别是从基站观测的伪距和载波相位,Rcorrection和Φcorrection分别是对伪距和载波相位的校正,ρpredicted是给定基站的已知位置而预测的几何距离。
基于流动站(例如,图1的装置100)所获得的数据的校正然后被计算为:
Rcorrected=Rrover-Rcorrection (13)
Φcorrected=Φrover-Φcorrection (14)
校正的伪距可观测量可表示如下:
Rcorrected=ρ+cΔdtr+ΔT+αfΔSTEC+ΔkP,r+mP,rover-mP,base+∈P (15)
其中,除了如上定义的项之外,Δdtr是流动站与基站之间的接收器时钟偏差的差(Δdtr=dtrover–dtbase),ΔT是对流层延迟的差,αfΔSTEC是电离层延迟的差,ΔkP,r是接收器代码仪器延迟的差,mP,rover是多径对流动站处测量的伪距的影响,mP,base是多径对基站处测量的伪距的影响。流动站和基站接收器噪声二者已被归入∈P中。
类似地,校正的相位可观测量可表示如下:
Φcorrected=ρ+cΔdtr+ΔT-αfΔSTEC+ΔkL,r+λ(nrover-nbase)+λΔw+mL,rover-mL,base+∈L (16)
其中,除了上面定义的项之外,ΔkL,r是接收器代码仪器延迟的差,mL,rover是多径对流动站处测量的载波相位的影响,mL,base是多径对基站处测量的载波相位的影响,n是整周模糊度,Δw是圆偏振缠绕的差。流动站和基站接收器噪声二者已被归入∈L中。
在实施方式中,由卫星旋转导致的圆偏振缠绕项的部分可被抵消,并且对于单历元分析,由于流动站旋转而引起的缠绕可被仪器偏差项ΔkL,r吸收,因此被忽略。如果仅使用短基线(<20km),则电离层项αfΔSTEC也可被忽略。对流层延迟使用仅取决于卫星仰角E的倾斜因子M(E)来建模。因此,最终可观测量模型可如下简化:
Rcorrected=ρ+cΔdtr+ΔTzenithM(E)+ΔkP,r+mP,rover-mP,base+∈P (17)
Φcorrected=ρ+cΔdtr+ΔTzenithM(E)+ΔkL,r+λ(nrover-nbase)+mL,rover-mL,base+∈L(18)
其中,除了如上定义的项之外,ΔTzenith是天顶对流层延迟的差。
这样,可观测量(伪距R和载波相位Φ)被表示为状态x的函数。在实施方式中,可观测量的校正可为可选的。在实施方式中,指定***模型可在测量之前执行。
在实施方式中,状态x可被定义为以下状态向量:
其中,除了上面定义的项之外,p是流动站的位置。仪器偏差ΔkP,r和ΔkL,r每GNSS信号类型包括一个偏差,然而根据定义,一个GNSS上的伪距偏差被取为零,因此在状态向量中省略。在实施方式中,GNSS包括GPS、Galileo(由欧盟创建的GNSS)和GLONASS(由俄罗斯创建的GNSS),因此,典型历元在x中具有10个状态(在p中三个,一个用于Δdtr,两个用于Galileo和GLONASS伪距偏差ΔkP,r,三个用于GNSS特定相位偏差ΔkL,r并且一个用于ΔTzenith)。
方法200包括将与测量关联的质量度量q量化的步骤S230。在实施方式中,一些质量度量(例如,载波噪声密度比、多径干扰的影响和卫星仰角)可由GNSS装置(例如,图1的装置100)提供;一些其它质量度量(例如,GNSS信号的更多信息)可通过检查不同卫星之间的一致性来获得。在实施方式中,量化质量度量可在操作期间(例如,在测量期间)执行。在检查不同卫星之间的一致性时,主界限计算未使用的量(例如,紧邻感兴趣的历元的伪距和载波相位的改变)可能是检查的重点。
在实施方式中,影响方法的准确性和一致性的一些测量被排除或减权重。例如,识别并排除受到多径干扰严重影响的一些信号。还排除可能造成大误差的一些非视线信号。
在实施方式中,基于随机样本共识(RANSAC)技术排除不利地影响方法的准确性和一致性的测量。例如,一秒间隔内的一组δ相测量(包括频带L1和L2二者中的信号)可以是RANSAC技术的输入。仅锁相达至少一秒的信号可参与此方法,因此所有其它信号被立即排除。在一秒短间隔内,仪器偏差ΔkL,r、对流层延迟ΔTzenithM(E)和基站多径mL,base的改变极其小,因此可被忽略。然后,相位的改变ΔΦcorrected由下式给出:
ΔΦcorrected=Δρ+cΔΔdtr+λ(Δnrover-Δnbase)+ΔmL,rover+∈L (20)
其中,Δρ是一秒间隔内的距离改变,ΔΔdtr是残留时钟偏差的改变,Δnrover和Δnbase是模糊度的改变(即,周跳),ΔmL,rover是多径干扰贡献的改变。
一些信号可受到非视线(NLOS)信号和多径干扰的严重影响,这里,多径干扰项ΔmL,rover可较大。一些信号可受到周跳的严重影响,因此λ(Δnrover–Δnbase)可较大。多径干扰项ΔmL,rover大于阈值多径干扰的信号或者周跳项λ(Δnrover–Δnbase)大于阈值周跳的信号可被识别为异常值并被排除。
图3是示出与本公开的一些实施方式一致的识别GNSS测量的异常值的方法300的流程图。参照图3,方法300包括选择四个随机卫星并为各个卫星随机地选择L1频带或L2频带δ相测量的步骤S310。
方法300包括使用四个δ相测量来求解位置的改变和残留时钟偏差的改变的步骤S320。这是可解的,因为***是具有四个测量和四个未知数的线性***。
方法300包括使用步骤S320中的解来预测所有卫星和频带的δ相的步骤S330。基于在步骤S320中获得的解,可预测所有卫星和频带的δ相。方法300包括确定预测的信号与实际信号之间的差并对差小于阈值差的信号的数量进行计数的步骤S340。预测的值与实际值之间的差小于阈值差的信号被视为测量的正常值。
方法300包括将步骤310-340重复特定迭代次数的步骤S350。迭代次数可基于确保有效收敛到解的置信度来确定。方法300包括识别具有数量最多的正常值的测量的步骤S360。在步骤S350中的所有迭代完成之后执行步骤S360。
方法300包括确定最高正常值数量是否大于阈值正常值数量(例如,预定最少数量的正常值)的步骤S370。如果确定结果为“是”,则方法300执行形成测量的正常值的最终输出的步骤S380。另一方面,如果确定结果为“否”,则方法300执行将所有信号作为测量的异常值排除的步骤S390。在实施方式中,步骤S350中的迭代次数可为100,步骤S340中的阈值差为0.02m,并且步骤S370中的正常值的阈值数量为15。
这样,使用RANSAC技术来识别并排除测量的异常值,从而增强测量的准确性和一致性。
返回参照图2,在方法200的步骤S230中,代替使用RANSAC技术,可使用基于窗口的技术来获得推导的质量度量。可通过在小时间窗口(例如,至多四秒长)内寻找δ相和δ伪距二者的一致性来指定基于窗口的质量度量,如下面参照图4详细描述的。
图4是示出与本公开的一些实施方式一致的将与GNSS测量关联的质量度量量化的方法400的流程图。参照图4,方法400包括围绕感兴趣的历元定义窗口的步骤S410。窗口的宽度是可配置的。在实施方式中,窗口宽度可以是二秒的半宽。窗口可以是对称的并且可具有四秒的全宽,并且可包含来自例如五个历元(即,样本数量nsamples=5)的数据。使用感兴趣的历元之后的历元意味着当在实时***中实现时,算法具有等于窗口半宽的延迟。
方法400包括确定存在于所有指定的历元(例如,五个历元)处并且对于所有历元具有锁相连续性的信号的步骤S420。信号总数nsig是L1频带中的信号数量nL1与L2频带中的信号数量nL2之和(nsig=nL1+nL2)。为了被包括,对于所有历元,信号需要存在并且具有锁相连续性。不满足这些条件的信号可被排除。
方法400包括获得相位和伪距测量数据的步骤S430。相位测量的总数由nsamples×(nL1+nL2)给出。类似地,伪距测量的总数为nsamples×(nL1+nL2),提供数据点的总数m=2nsamples×(nL1+nL2)。在实施方式中,代替相位和伪距测量,在围绕感兴趣的历元的窗口中执行多普勒测量,并且步骤S430中的获得是获得多普勒测量数据。
方法400包括通过表示相对于窗口的第一历元的位置和时钟偏差偏移的***连同各个频带中的相位和伪距的固定偏差项一起对测量数据进行建模的步骤S440。因此,自由参数的总数为n=4(nsamples–1)+4。
方法400包括例如使用加权线性最小二乘问题来求解***的步骤S450。例如,在列向量y中放置m个数据点,定义列向量x以保持n个自由参数,并且为线性模型y=Ax定义m×n大小的矩阵A。接下来,定义权重w的列向量,每个测量一个权重,伪距给予权重1/σPR,相位给予权重1/σphase。对于给定状态x,残差r可表示如下:
r=y-Ax (21)
可通过将该问题变换为标准形式:
A′=diag(w)A (23-1)
y′=diag(w)y (23-2)
来在数值上求解式(22)。
在实施方式中,残差上的损失函数可如下定义:
其中ρ(s)是损失函数。
例如,可如下引入Huber损失函数:
其中,阈值参数δ被设定为2。这具有对测量的异常值减权重的效果。然后,可使用鲁棒的高斯-牛顿算法来寻找最佳拟合解步骤S440和S450相当于识别共识解,例如,对于位置和时钟偏差的改变的共识(consensus)。可使用相对于共识的偏离来如下所述表征GNSS测量的质量。
方法400以计算各个信号的残差的均方根(RMS)以推导质量度量的步骤S460结束。这相当于确定相对于共识的偏离。一旦找到式(22)的解,可计算残差。对于各个信号,可计算残差的RMS,并且这形成推导的质量度量。相位和伪距度量可单独地计算。在实施方式中,这些度量可用于排除信号。作为示例,方法400的参数可如下设定:最大允许伪距RMS为3m;最大允许相位RMS为0.05循环;并且至少15个信号必须通过这两个检查,否则将排除整个历元。
这样,实现测量的异常值的减权重,并且使用δ伪距和δ相二者来获得推导的质量度量,而无需大量迭代以便于收敛,从而导致处理的效率和准确性增强。
返回参照图2,方法200包括指定非高斯误差概率密度模型f(r|θ,q)并使用实验数据对误差概率密度模型参数θ进行拟合的步骤S240。在实施方式中,指定非高斯误差概率密度模型并对模型参数进行拟合可在测量之前(例如,先验离线)执行。在实施方式中,非高斯模型是student-t分布模型,并且伪距误差被建模为student-t分布,其中概率密度被表示为:
其中,r是残差,v是自由度参数,σ是定义核心分布的宽度的缩放参数,Γ是伽马函数。
参数v和σ不是常数,并且可取决于质量度量q。在实施方式中,可使用三个质量度量:卫星仰角(el)、载波噪声密度比(cno)和锁相时间(t)。为了简化拟合处理,首先根据对任何实数均有效的量定义基本student-t参数σ和v(称为无约束中间参数σ′和v′):
σ=log(1+e(σ′)) (28)
ν=log(1+e(ν′)) (29)
使用这些无约束中间参数σ′和v′允许将模型拟合处理作为无约束优化问题处理,从而避免表示σ和v必须均为正数的约束的需要。
然后,无约束student-t参数被定义为质量度量的函数:
其中,X0是2×1矩阵,X1是2×2矩阵。此模型总共具有六个自由参数:两个在X0中,四个在X1中。将所有这些无约束参数捆绑为向量θ并将所有质量度量捆绑为q,具有质量度量qi的残差ri的概率密度被写为:
fpr(ri|θ,qi) (31)
通过给定一组测量观测数据,寻找九个参数(六个自由参数加残差ri、向量θ和质量度量qi)的最大似然值,来识别模型。使用贝叶斯公式:
最大似然模型是使P(model|data)最大化的模型。使用θ表示式(32)中描述的模型,使用r表示观测的残差,并使用q表示关联的质量度量,式(32)变为:
这是针对θ的最大化问题。P(r|q)是归一化因子并且可被忽略。然后,式(35)变为:
先验模型参数P(θ)被认为均匀,因此可被省略。因此,式(36)可被简化为:
整个数据集fpr(r|θ,q)的组合概率密度可被评估为各个残差的所有单独概率密度的乘积:
在实施方式中,代替直接评估整个数据集以进行数值优化,使用对数似然,并且式(39)变为:
在实施方式中,类似于上述伪距模型,载波相位残差被建模为student-t分布。在相位残差模型中,误差分布可被建模为student-t加均匀混合(以范围±0.5为界)。由于这里可观测量是相位,所以不会发生大于0.5循环的误差。
混合的基本概率密度如下:
如伪距一样,基本参数w、σ和ν不是常数,而是质量度量的函数。首先根据无约束值定义基本参数:
σ=log(1+e(σ′)) (42)
ν=log(1+e(ν′)) (43)
w′和w之间的映射允许w′取任何值,同时将w约束在0和1之间。然后,根据质量度量来定义无约束参数。这里仅使用卫星仰角(el)信息:
剩余无约束参数w′被认为与卫星仰角无关。总共存在五个自由参数,两个在X0中,两个在X1中,一个用于w′。然后,模型拟合处理与用于伪距的方法相似。因为车辆的运动不显著影响相位误差分布,所以单个模型覆盖静止历元和移动历元二者。
为了性能优化,可通过取fphase(r)的两个分量的最大值来近似均匀student混合:
类似于上述伪距,式(46)也可使用对数似然来评估。
图5示出与本公开的一些实施方式一致的观测的残留误差分布(410)、对观测的残留误差分布的高斯分布拟合(420)以及对观测的残留误差分布的student-t分布拟合(430)的比较。观测的残留误差分布是使用设置在移动车辆中的接收器获得的伪距残留误差分布,其中卫星(Galileo)仰角介于10-30度之间。如图5所示,高斯分布拟合(420)是对观测的残留误差分布(410)的较差拟合,因为观测的残留误差分布(410)的尾部远比高斯分布拟合严重。另一方面,student-t分布拟合(430)是对观测的残留误差分布(410)的改进拟合,指示与高斯分布相比,student-t分布是对残留误差分布进行建模的更好选择。
返回参照图2,方法200包括定义状态x上的未归一化后验概率密度的步骤S250。在实施方式中,定义未归一化后验概率密度可在操作期间执行。在实施方式中,未归一化后验概率密度使用上式(8)来定义。在定义误差分布的数学形式并且误差模型参数θ拟合到一组实验数据的情况下,可给定状态x计算特定一组测量的似然:
P(z|θ,x)=∏if(ri|θ,qi) (47)
=∏if(zi-hi(x)|θ,qi) (48)
这里,hi(x)是给定***的状态x,预测可观测量zi的***模型。以P(z)为未知归一化常数,最终后验概率密度由下式给出:
P(x|θ,z,q)∝P(x)∏if(zi-hi(x)|θ,qi) (49)其中,P(x)是先验状态概率。在实施方式中,除了被认为具有以零为中心的高斯分布并且标准偏差为5cm的ΔTzenith之外,在所有状态下P(x)被认为均匀。
方法200包括估计状态x并通过在状态x上对后验概率密度进行积分来计算保护级别的步骤S260。可例如使用后验概率密度P(x|z,q,θ)来估计状态x。积分可在数值上例如使用马尔科夫链蒙特卡洛(MCMC)或重要性采样方法来完成。在实施方式中,估计状态x和计算保护级别可在操作期间执行。
在实施方式中,后验概率密度在状态x上的积分基于MCMC来执行,但是使用具有修改的密度函数的多个链来给出多峰分布上的尾密度的准确估计。多峰可由存在多个不同的分布模式来表征。例如,在第一阶段,使用一组并行交互链来解决对多峰后验采样的问题。代替直接从后验概率密度提取样本,定义修改的概率密度函数,以使得相位概率分布(式(42))的σ参数可改变为:
σ=max(log(1+e(σ′)),σmin) (50)
这具有拓宽后验概率密度的模式的效果。如果σmin被设定为较大值(例如,0.5m),则模式合并并且分布变为单峰。例如,一组10个MCMC链并行执行,各个MCMC链具有不同的σmin值。这些链可被视为形成链层叠,其中大大拓宽的单峰链在顶部,未拓宽的链(σmin设定为0)在底部。相邻链可经由副本交换来交互,即,两个相邻链的状态有时使用常见的Metropolis标准来交换。在采样结束时,仅使用底部(未拓宽)链的状态。该底部链可准确地对多峰后验概率密度采样。
图6是示出与本公开的一些实施方式一致的使用并行交互链的多峰采样方法的示意图。10个MCMC链并行执行(为了清晰,图6中仅示出三个链),各个MCMC链具有不同的σmin值。历元具有在-1.5m和0m处具有模式的双峰后验。未修改的链(σmin=0)可通过与上面逐渐拓宽的链交互在-1.5m和0m处的两个模式之间移动。在最高级别(σmin=0.5),分布为单峰的。这样,使用具有修改的密度函数的多个链,可准确地估计多峰分布上的尾密度。
在实施方式中,代替直接从后验概率分布提取大量样本,使用重要性采样。例如,使用伞形采样来提取偏向感兴趣的区域(例如,分布的尾部)的样本。图7是示出与本公开的一些实施方式一致的执行伞形采样的方法700的流程图。如图7所示,方法700包括步骤S710:通过从无偏分布提取第一样本集合并在第一样本集合上估计分位数的第一位置来执行第一轮采样。例如,分位数可被设定为中等值(例如,0.2)。该分位数将用于在第二轮采样中约束采样池。
方法700包括步骤S720:通过从基于分位数的第一位置确定的第一约束分布提取第二样本集合来执行第二轮采样。第一约束分布可被设定为无偏分布的尾部。方法700包括步骤S730:在第二样本集合上估计分位数的第二位置。分位数的第二位置用于为第三轮采样设定第二约束分布。
方法700包括步骤S740:使用在先前采样中设定的约束分布来重复采样,直至第n轮采样完成。例如,可使用六轮采样。对于第n样本集合,在第(n-1)采样中设定的约束分布强制仅采样尾部的qn部分。例如,对于分位数值0.2,在第一轮中仅采样尾部的0.2部分,在第二轮中仅采样尾部的0.04部分,在第三轮中仅采样尾部的0.008部分,依此类推。
方法700包括步骤S750:确定对于所需保护级别计算,采样是否足够了。如果步骤S740中的确定结果为“否”,则该方法从步骤S740重复并提取另一轮样本。另一方面,如果步骤S740中的确定结果为“是”,则在步骤S760中,所有样本与适当的权重组合,并且确定最终保护级别。
在实施方式中,使用一小组试验数据测试了包括接收器(例如,图1中的装置100)和卫星的***。接收器设置在车辆中,并且这一组试验数据包括车辆在公路上行驶约25小时而获得的数据。使用本地(例如,20km内)参考站并固定L1和L2整周模糊度,使用实时运动学(RTK)方法获得测量残差。使用相同的数据集来对误差模型参数θ进行拟合并测试保护级别计算。数据包括来自GPS、Galileo和GLONASS的测量。相位观测是宽波的,在不发送民用L2信号的GPS卫星上不可用。与本公开的一些实施方式一致,图8A和图8B示出显示从顺轨方向上的实验数据获得的界限半宽和真实误差的斯坦福图,图8C和图8D示出显示从跨轨维度的实验数据获得的界限半宽和真实误差的斯坦福图,图8E和图8F示出显示从垂直维度的实验数据获得的界限半宽和真实误差的斯坦福图。图8A和图8B示出相同的数据,但是具有不同的轴极限。类似地,图8C和图8D示出相同的数据,但是具有不同的轴极限,图8E和图8F示出相同的数据,但是具有不同的轴极限。图8A至图8F示出测试数据的保护级别和实际误差,完整性风险被设定为每25小时的数据预期一个事件。下表1中示出跨轨界限的可用性。
表1
这里,可用性仅考虑足够的GNSS信号可用于计算界限的历元。
上述单历元位置界限算法形成完整性方法的核心组件。该算法还可用于确定全(多峰)后验分布,其可用于其它应用,而不仅用于界限或位置估计。在实施方式中,以低速率(例如,每10秒一次)评估界限传播(例如,单历元界限),并且并行地,使用δ相观测和传感器数据使最近的单历元界限和位置向前传播至当前时间。数学上,当前时间的界限是某一先前历元处的界限之和加上自该时间以来位置改变上的界限。为了保持完整性,在绝对界限计算与传播计算之间共享误导信息的风险。因此,以与主界限计算相似的完整性级别来完成传播。通过这种方式,***可维持高可用性、低延迟和减少的计算负荷。
在实施方式中,使用原型传播方案。原型传播方案基于两个不同的方法:(1)使用非高斯分布(例如,student-t分布)对δ相上的误差分布进行建模,并且使用MCMC在数值上评估δ位置上的界限;以及(2)使用高斯过界(Gaussian over-bound)方法对误差分布进行建模,这在航空领域中由接收器自主完整性监测(RAIM)算法使用。
在实施方式中,与GNSS信号一起或者当GNSS信号不可用时,可使用基于传感器的传播。例如,由运动传感器(例如,加速度计、陀螺仪、车轮速度传感器或者设置在车辆中的任何其它传感器)进行的测量可用于测量。在基于传感器的传播中,可使用扩展卡尔曼滤波器来跟踪位置的改变。在断电之前,GNSS信号可用于估计***的状态(例如,取向和传感器偏差)。
在实施方式中,如果多径干扰严重,则可添加一个或更多个模型参数。模型参数可与位置同时估计。在实施方式中,代替本地参考站,可使用状态空间表示(SSR)校正服务。与SSR校正服务关联的误差可被当作解中的坐标移位,因此可通过扩大极限以允许该移位的可能大小来恢复界限有效性。可选地,可以明确地对SSR误差进行建模。
可读存储介质可以是可存储指令以供指令执行装置执行的有形装置。计算机可读存储介质可以是例如(但不限于)电子存储装置、磁存储装置、光学存储装置、电磁存储装置、半导体存储装置或者前述的任何合适的组合。计算机可读存储介质的更具体示例的非穷尽列表包括以下:便携式计算机磁盘、硬盘、随机存取存储器(RAM)、只读存储器(ROM)、可擦除可编程只读存储器(EPROM)、静态随机存取存储器(SRAM)、便携式紧凑盘只读存储器(CD-ROM)、数字多功能盘(DVD)、记忆棒、软盘、机械编码装置(例如,打孔卡或记录有指令的凹槽中的凸起结构)以及前述的任何合适的组合。
本公开的计算机可读程序指令可以是汇编指令、指令集架构(ISA)指令、机器指令、机器相关指令、微代码、固件指令、状态设定数据或者以一种或更多种编程语言的任何组合编写的源代码或目标代码,包括面向对象的编程语言和传统过程编程语言。计算机可读程序指令可作为独立软件包完全在计算装置上执行,或者部分地在第一计算装置上,部分地在远离第一计算装置的第二计算装置上执行。在后一种情况下,第二远程计算装置可通过任何类型的网络连接到第一计算装置,包括局域网(LAN)或广域网(WAN)。
附图中的流程图和框图示出根据各种实施方式的***、方法和装置的可能实现方式的架构、功能和操作的示例。应该注意的是,在一些替代实现方式中,方框中指出的功能可不按图中指出的顺序发生。例如,根据所涉及的功能,连续示出的两个方框可事实上基本上同时执行,或者方框有时可按相反顺序执行。
将理解,所描述的实施方式并不互斥,结合一个示例实施方式描述的元件、组件、材料或步骤可按照合适的方式与其它实施方式组合或从其它实施方式去除,以实现期望的设计目的。
本文中提及“一些实施方式”或“一些示例性实施方式”意指结合实施方式描述的特定特征、结构或特性可包括在至少一个实施方式中。在本公开的各种地方出现短语“一个实施方式”、“一些实施方式”或“另一实施方式”未必全部指同一实施方式,也未必是与其它实施方式互斥的单独或替代实施方式。
应该理解,本文中阐述的示例方法的步骤未必需要按所描述的顺序执行,这些方法的步骤的顺序应该被理解为仅是示例。例如,根据所涉及的功能,连续示出的两个方框可事实上基本上同时执行,或者方框有时可按相反顺序执行。同样,这些方法中可包括附加步骤,并且在与各种实施方式一致的方法中,某些步骤可被省略或组合。
如本公开中所使用,词语“示例性”在本文中用于意指充当示例、实例或例示。本文中描述为“示例性”的任何方面或设计未必被解释为比其它方面或设计优选或有优势。相反,使用该词语旨在以具体方式呈现概念。
如本公开中所使用,除了具体地另外说明,否则术语“或”涵盖除了不可行之外的所有可能的组合。例如,如果称数据库可包括A或B,则除非具体地另外说明或者不可行,否则数据库可包括A、或B、或A和B。作为第二示例,如果称数据库可包括A、B或C,则除非具体地另外说明或者不可行,否则数据库可包括A、或B、或C、或A和B、或A和C、或B和C、或A和B和C。
另外,除非另外指明或者从上下文清楚是指单数形式,否则如本公开和所附权利要求中使用的冠词通常应该被解释为意指“一个或更多个”。
除非明确地另外说明,否则各个数值和范围应该被解释为近似的,如同该值或范围的值之前有词语“约”或“大约”一样。
尽管随附方法权利要求(如果有的话)中的元素以特定顺序叙述,但除非权利要求叙述另外暗指实现那些元素中的一些或全部的特定顺序,否则那些元素未必旨在限于以该特定顺序实现。
将理解,为了清晰而在单独实施方式的上下文中描述的本公开的某些特征也可在单个实施方式中组合提供。相反,为了清晰而在单个实施方式的上下文中描述的说明书的各种特征也可单独地或以任何合适的子组合或如说明书的任何其它描述的实施方式中那样合适地提供。除非如此指出,否则在各种实施方式的上下文中描述的某些特征不是那些实施方式的必要特征。
还将理解,在不脱离随附权利要求的范围的情况下,本领域技术人员可进行为了说明所描述的实施方式的性质而描述和示出的部件的细节、材料和布置方式的各种修改、替代和变化。因此,随附权利要求涵盖落在权利要求的条款内的所有这些替代、修改和变化。
Claims (15)
1.一种使用全球导航卫星***GNSS测量的单个历元来确定位置估计的保护级别的方法(200),该方法包括以下步骤:
预先指定(S210)状态x的先验概率密度P(x);
预先指定(S220)将所述状态x与所述测量的可观测量z联系起来的***模型h(x);
在操作期间,将与所述测量关联的质量度量q量化(S230);
预先指定(S240)非高斯残留误差概率密度模型f(r|θ,q)并对先验离线确定的模型参数θ进行拟合;
在操作期间,定义(S250)后验概率密度P(x|z,q,θ);
在操作期间,估计(S260)所述状态x;以及
在操作期间,通过在所述状态x上对所述后验概率密度P(x|z,q,θ)进行积分来计算(S260)所述保护级别。
2.根据权利要求1所述的方法,其中,所述***模型被指定为:
z=h(x)+r;
其中r是残留误差。
3.根据权利要求1所述的方法,其中,将与所述测量关联的所述质量度量q量化(S230)的步骤包括使用随机样本共识RANSAC技术来排除所述测量的一个或更多个异常值。
6.根据权利要求4或5所述的方法,其中,使用一组实验数据对所述模型参数θ进行拟合的步骤包括:
使用无约束student-t参数定义所述自由度参数和所述缩放参数,并且当从属于权利要求5时,定义所述权重参数;
将所述无约束student-t参数定义为所述质量度量的函数;以及
寻找所述无约束student-t参数的最大似然值。
8.根据前述权利要求中的任一项所述的方法,其中,在所述状态上对所述后验概率密度进行积分的步骤包括:
使用具有修改的密度函数的多个交互链来应用马尔科夫链蒙特卡洛MCMC数值积分。
9.根据前述权利要求中的任一项所述的方法,其中,在所述状态上对所述后验概率密度进行积分(S260)的步骤包括:
通过从无偏分布提取第一样本集合并在所述第一样本集合上估计分位数的第一位置来执行(S710)第一轮采样;
通过从基于所述分位数的所述第一位置确定的第一约束分布提取第二样本集合来执行(S720)第二轮采样;
估计(S730)所述第二样本集合上的所述分位数的第二位置,其中,所述分位数的所述第二位置用于为第三轮采样设定第二约束分布;
使用在先前采样中设定的约束分布来重复(S740)所述采样,直至第n轮采样完成,其中n是大于2的自然数;
确定(S750)n轮采样的样本是否足以计算所述保护级别;以及
当n轮的样本足以计算所述保护级别时,将所述n轮采样的样本与适当的权重组合(S760)。
10.根据前述权利要求中的任一项所述的方法,该方法还包括以下步骤:
使用由一个或更多个传感器获得的δ相误差分布或数据使最近的单历元界限向前传播到当前时间。
11.根据权利要求10所述的方法,其中,所述传播的步骤还包括:
使用非高斯分布函数对δ相上的误差分布进行建模并且使用马尔科夫链蒙特卡洛技术在数值上评估δ位置上的界限;或者
使用高斯过界方法对δ位置上的所述误差分布进行建模。
13.一种被配置为使用GNSS测量的单个历元来确定位置估计的保护级别的设备(100),该设备包括:
接收器(104),该接收器被配置为接收GNSS信号并处理所接收的信号以进行测量;以及
处理器(106),该处理器被配置为:
指定(S210)状态x的先验概率密度P(x);
指定(S220)将所述状态x与所述测量的可观测量z联系起来的***模型h(x);
将与所述测量关联的质量度量q量化(S230);
指定(S240)非高斯残留误差概率密度模型f(r|θ,q)并使用一组实验数据对模型参数θ进行拟合;
定义(S250)后验概率密度P(x|z,q,θ);
估计(S260)所述状态x;以及
通过在所述状态x上对所述后验概率密度P(x|z,q,θ)进行积分来计算(S260)所述保护级别。
14.根据权利要求13所述的设备,其中,所述处理器(106)还被配置为:
与外部传感器通信以获得由所述传感器跟踪的位置改变数据以确定界限传播。
15.一种存储有指令的非暂时性计算机可读介质,所述指令在由处理器(106)执行时执行一种使用GNSS测量的单个历元来控制GNSS装置的方法,该方法包括以下步骤:
指定(S210)状态x的先验概率密度P(x);
指定(S220)将所述状态x与所述测量的可观测量z联系起来的***模型h(x);
将与所述测量关联的质量度量q量化(S230);
指定(S240)非高斯残留误差概率密度模型f(r|θ,q),其中模型参数θ先验离线确定;
定义(S250)后验概率密度P(x|z,q,θ);以及
估计(S260)所述状态x并且通过在所述状态x上对所述后验概率密度P(x|z,q,θ)进行积分来计算(S260)所述保护级别。
Applications Claiming Priority (4)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
US16/779,464 US11885890B2 (en) | 2020-01-31 | 2020-01-31 | Method and apparatus of single epoch position bound |
US16/779,464 | 2020-01-31 | ||
EP20170246.1A EP3859397A1 (en) | 2020-01-31 | 2020-04-17 | Method and apparatus of single epoch position bound |
EP20170246.1 | 2020-04-17 |
Publications (1)
Publication Number | Publication Date |
---|---|
CN113281793A true CN113281793A (zh) | 2021-08-20 |
Family
ID=70333810
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202110138235.7A Pending CN113281793A (zh) | 2020-01-31 | 2021-02-01 | 单历元位置界限的方法和设备 |
Country Status (3)
Country | Link |
---|---|
US (2) | US11885890B2 (zh) |
EP (1) | EP3859397A1 (zh) |
CN (1) | CN113281793A (zh) |
Cited By (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN113792371A (zh) * | 2021-09-27 | 2021-12-14 | 江西科技学院 | 基于锁相值的轨道异常匹配的诊断方法 |
WO2023050760A1 (en) * | 2021-09-28 | 2023-04-06 | Qualcomm Incorporated | Precise positioning engine (ppe) base station swap handling |
Families Citing this family (11)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
WO2020247839A1 (en) * | 2019-06-06 | 2020-12-10 | Star Ally International Limited | Single-epoch pseudo-range positioning under varying ionosphere delays |
US20220080991A1 (en) * | 2020-09-11 | 2022-03-17 | Beijing Wodong Tianjun Information Technology Co., Ltd. | System and method for reducing uncertainty in estimating autonomous vehicle dynamics |
CN114329342B (zh) * | 2021-12-20 | 2024-02-06 | 中交第二公路勘察设计研究院有限公司 | 基于极值-高斯混合分布模型的高速公路通行能力计算方法 |
WO2023132987A1 (en) * | 2022-01-07 | 2023-07-13 | Qualcomm Incorporated | Integrity information for radio access technology dependent positioning assistance data |
EP4307012A1 (en) | 2022-07-11 | 2024-01-17 | u-blox AG | Gnss measurement processing and residual error model estimation |
EP4307008A1 (en) | 2022-07-11 | 2024-01-17 | u-blox AG | Gnss measurement processing with carrier phase residual error model |
EP4307009A1 (en) | 2022-07-11 | 2024-01-17 | u-blox AG | Gnss measurement processing to identify modes |
CN115144411B (zh) * | 2022-09-05 | 2022-12-06 | 国家卫星海洋应用中心 | 基于卫星散射计的海冰检测的方法、装置、设备和介质 |
US20240201392A1 (en) * | 2022-12-15 | 2024-06-20 | Swift Navigation Inc. | System and method for bounding a satellite positioning solution integrity |
CN116794695B (zh) * | 2023-08-22 | 2023-10-24 | 中国科学院空天信息创新研究院 | 基于双重优化的高程辅助低轨卫星机会信号定位方法 |
CN116819585B (zh) * | 2023-08-31 | 2023-12-29 | 长沙金维信息技术有限公司 | 基于非线性优化的gnss单点定位方法及导航方法 |
Family Cites Families (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN106569233B (zh) * | 2016-11-03 | 2019-03-19 | 中国铁路总公司 | 基于学生t分布的接收机自主完好性检测和故障排除方法 |
JP6855580B2 (ja) * | 2016-12-30 | 2021-04-07 | ユー−ブロックス、アクチエンゲゼルシャフトu−blox AG | Gnss受信機保護レベル |
-
2020
- 2020-01-31 US US16/779,464 patent/US11885890B2/en active Active
- 2020-04-17 EP EP20170246.1A patent/EP3859397A1/en active Pending
-
2021
- 2021-02-01 CN CN202110138235.7A patent/CN113281793A/zh active Pending
-
2023
- 2023-12-13 US US18/537,921 patent/US20240118433A1/en active Pending
Cited By (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN113792371A (zh) * | 2021-09-27 | 2021-12-14 | 江西科技学院 | 基于锁相值的轨道异常匹配的诊断方法 |
CN113792371B (zh) * | 2021-09-27 | 2024-01-26 | 江西科技学院 | 基于锁相值的轨道异常匹配的诊断方法 |
WO2023050760A1 (en) * | 2021-09-28 | 2023-04-06 | Qualcomm Incorporated | Precise positioning engine (ppe) base station swap handling |
Also Published As
Publication number | Publication date |
---|---|
US20240118433A1 (en) | 2024-04-11 |
EP3859397A1 (en) | 2021-08-04 |
US20210239845A1 (en) | 2021-08-05 |
US11885890B2 (en) | 2024-01-30 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
US11885890B2 (en) | Method and apparatus of single epoch position bound | |
US11237276B2 (en) | System and method for gaussian process enhanced GNSS corrections generation | |
US11860260B2 (en) | Systems and methods for high-integrity satellite positioning | |
US9766349B1 (en) | Localization and tracking using location, signal strength, and pseudorange data | |
JP4781313B2 (ja) | マルチパス検出装置、測位装置、姿勢方位標定装置、マルチパス検出方法およびマルチパス検出プログラム | |
US20220107427A1 (en) | System and method for gaussian process enhanced gnss corrections generation | |
JP2010528320A (ja) | リアルタイムキネマティック(rtk)測位における距離依存性誤差の軽減 | |
JP6234549B2 (ja) | 測位装置 | |
CN110662987B (zh) | 定位方法和定位终端 | |
CN112313542A (zh) | 定位方法及定位终端 | |
JP2008039691A (ja) | 搬送波位相式測位装置 | |
CN118339480A (zh) | 一种数据预处理方法、数据预处理装置和芯片 | |
US11112508B2 (en) | Positioning method and positioning terminal | |
JP2023508604A (ja) | 整数なしのgnss測位のためのシステムおよび方法 | |
JP2023508119A (ja) | Gnssアンビギュイティ決定のためのシステムおよび方法 | |
JP2010197353A (ja) | 計算機装置 | |
CN110678781B (zh) | 定位方法和定位终端 | |
CN114355410B (zh) | 基于并行计算的卫星导航实时精密单点定位***及方法 | |
Madhu Krishna et al. | Improving the position accuracy of rover receiver using differential positioning in indian regional navigation satellite system | |
KR101803652B1 (ko) | 광역 보강 시스템의 전리층 불균질 위협모델을 위한 기하학적 측정 메트릭 개발을 통한 sbas 가용성 향상 방법 및 장치 | |
Dekkiche et al. | Differential ionosphere modelling for single-reference long-baseline GPS kinematic positioning | |
Viandier et al. | Studies on dpm for the density estimation of pseudorange noises and evaluations on real data |
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 |