CN107728171B - 基于粒子滤波的gnss相位系统间偏差实时追踪和精密估计方法 - Google Patents

基于粒子滤波的gnss相位系统间偏差实时追踪和精密估计方法 Download PDF

Info

Publication number
CN107728171B
CN107728171B CN201710790184.XA CN201710790184A CN107728171B CN 107728171 B CN107728171 B CN 107728171B CN 201710790184 A CN201710790184 A CN 201710790184A CN 107728171 B CN107728171 B CN 107728171B
Authority
CN
China
Prior art keywords
value
particle
particles
equation
gnss
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201710790184.XA
Other languages
English (en)
Other versions
CN107728171A (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.)
Southwest Jiaotong University
Original Assignee
Southwest Jiaotong 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 Southwest Jiaotong University filed Critical Southwest Jiaotong University
Priority to CN201710790184.XA priority Critical patent/CN107728171B/zh
Publication of CN107728171A publication Critical patent/CN107728171A/zh
Application granted granted Critical
Publication of CN107728171B publication Critical patent/CN107728171B/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/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/27—Acquisition or tracking or demodulation of signals transmitted by the system creating, predicting or correcting ephemeris or almanac data within the receiver
    • 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
    • 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
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis

Landscapes

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

Abstract

本发明公开了一种基于粒子滤波的GNSS相位系统间偏差实时追踪和精密估计方法，包括以下步骤：步骤1：数据预处理，导入卫星星历、当前历元伪距观测值和当前历元相位观测值；步骤2：建立系统内和系统间的观测值方程，并线性化，得到单历元法方程或与之前历元累加法方程；步骤3：改正GNSS观测值法方程系统间偏差；解法方程，更新粒子滤波权；根据带权粒子，计算相位系统间偏差小数部分的数值及粒子均方根；步骤4：重复步骤1‑3，实时追踪GNSS系统间相位偏差小数部分的值，实现实时追踪；或重复步骤1‑3，精确估计相位系统间偏差值，固定系统内和系统间的双差整周模糊度，实现精密定位；本发明可高效的实现GNSS相位系统间偏差的实时追踪和精确估计，实现系统无差别的多模GNSS精密定位。

Description

基于粒子滤波的GNSS相位系统间偏差实时追踪和精密估计 方法
技术领域
本发明涉及卫星定位系统和定位测量技术领域，具体涉及一种基于粒子滤波的GNSS相位系统间偏差实时追踪和精密估计方法。
背景技术
全球导航卫星系统GNSS发展迅速，不久的将来，能够使用的系统至少包括美国的GPS、俄罗斯的GLONASS、中国的BDS和欧盟的Galileo等；以及局域系统日本的QZSS和印度NAVIC；可用卫星总数超过100颗，载波频率近20个；这些频率有很多重叠或者相近，在GNSS多系统组合定位中，对应的观测值可以放在一起组成系统间的双差观测值，从而充分利用观测信息；在接收机端，来自不同系统的卫星信号在硬件路径、数据处理以及初始相位等方面存在差异；因此双差之后会产生一个接收机端的系统间相位偏差(Inter-System Bias，ISB)；由于这个偏差的存在，使系统间的双差模糊度不是整数，也就不能在模糊度固定中被固定为整数；但是如果ISB已知，则可以在双差观测值或相应法方程中进行改正，恢复模糊度的整周特性，最终实现系统间双差模糊度的固定；系统间偏差ISB可分为两个部分，为波长整数倍的部分和剩余不足一周的部分(Fractional ISB，F-ISB)；前者可以忽略，后者需要精确估计或改正；传统的处理方法包括后处理法、假设偏差值和方差法；前者先精确解算整周模糊度再反推误差值，显然在偏差所在方程对整周模糊度解算很重要时，会因为不能成功固定整周模糊度而失败；后者假设的偏差值离真值较远时会起相反作用，影响偏差估计中的模糊度固定。
发明内容
本发明提供一种基于粒子滤波的GNSS相位系统间偏差实时追踪和精密估计方法。
一种基于粒子滤波的GNSS相位系统间偏差实时追踪和精密估计方法，包括以下步骤：
步骤1：对卫星导航数据进行预处理，导入卫星星历、当前历元伪距观测值和当前历元相位观测值；
步骤2：建立系统内和系统间的观测值方程，并线性化，得到单历元法方程或与之前历元累加法方程；
步骤3：计算粒子均方根，判断均方根是否大于给定阈值，若大于阈值则对粒子进行聚簇分析，通过聚簇分析整合分簇粒子；对每个粒子值，改正GNSS观测值法方程系统间偏差；解算法方程，使用LAMBDA方法进行模糊度固定并输出对应粒子的RATIO值；建立关于RATIO值的函数，用函数值更新粒子滤波权；根据带权粒子，计算相位系统间偏差小数部分的数值及粒子均方根；
步骤4：重复步骤1-3，实时追踪GNSS系统间相位偏差小数部分的值，实现实时追踪；重复步骤1-3，待滤波收敛后，通过估计的GNSS系统间相位偏差小数部分，改正观测值或法方程中相位系统间偏差值，固定系统内和系统间的双差整周模糊度，实现精密定位。
进一步的，所述步骤3中通过聚簇分析整合分簇粒子过程如下：
S1：选择粒子值最大和最小的两个粒子作为起始粒子，计算其它粒子到起始粒子的距离，按照距离大小将粒子分为两组；
S2：分别计算两个粒子组的加权重心g1、g2，其定义如下：
Figure GDA0002663941780000021
其中：xk为k时刻粒子值，wk为k时刻粒子对应权值，Nh为各组粒子个数，h为组号，h＝1，2；
S3：计算粒子组重心距离d：
d＝|g1-g2|；
S4：判断距离d与载波波长λ之间的差值是否小于eps，eps为较小的正数，是一个预设的阈值，若|d-λ|＜eps成立，则将两个粒子组进行合并。
进一步的，所述步骤3具体过程如下：
(1)采样生成初始粒子集
Figure GDA0002663941780000022
对于第k时刻粒子集由上一时刻滤波结果生成，其中：x为粒子值，w为粒子对应权值，N为粒子个数，i＝1，2…N为粒子序号；
(2)计算粒子均方根，判断均方根是否大于给定阈值，若大于阈值则对粒子进行聚簇分析，通过聚簇分析整合分簇粒子；
(3)对于每一个粒子值，改正GNSS观测值法方程系统间偏差；解算法方程，获得未知量的负电解和相应的协方差阵；使用LAMBDA方法进行模糊度固定并输出对应粒子的RATIO值；
(4)建立关于RATIO值的似然函数，用函数值进行粒子滤波权更新，并标准化粒子权值，作为新的粒子权值；
(5)计算粒子的期望值作为相位偏差的估计值
Figure GDA0002663941780000023
Figure GDA0002663941780000024
计算粒子的方差
Figure GDA0002663941780000025
Figure GDA0002663941780000031
(6)判断粒子滤波是否收敛，实时追踪时，判断均方根是否小于设定阈值stdthd，若是则输出相位偏差的估值和粒子的方差作为估计结果，进入步骤(8)；
(7)判断粒子滤波是否收敛，精确估计时，判断均方根是否小于设定阈值std2thd，若是则输出相位偏差的估计值和粒子的方差作为估值结果，退出滤波，进入步骤4；若否则进入步骤(8)；
(8)若满足重采样的条件，根据更新的权值重新采样，转入步骤(9)；
(9)在精确估计时，判断均方根是否小于设定阈值std1thd，若是则对重采样的粒子实施正则化，转入步骤(11)；若否则转入步骤(10)；实时追踪时，直接转入步骤(10)；
(10)预计下一时刻粒子，对重采样的粒子实时离散化：
Figure GDA0002663941780000032
其中：e为离散化时所加的随机噪声；
(11)推算下一历元的粒子值。
进一步的，所述步骤(4)中粒子滤波权更新过程如下：
S11：建立似然函数与RATIO之间的函数关系：
Figure GDA0002663941780000033
式中：f(RATIO)是关于RATIO值的函数；
S12：根据建立的函数关系和第i个粒子对应的RATIO值RATIOi计算对应粒子的似然函数值
Figure GDA0002663941780000034
Figure GDA0002663941780000035
S13：将似然函数值与对应粒子的权值相乘，获得更新后的粒子权
Figure GDA0002663941780000036
Figure GDA0002663941780000037
S14：标准化粒子权值，即将每个粒子的权与所有粒子权之和的比值，作为新的粒子权值
Figure GDA0002663941780000038
Figure GDA0002663941780000041
进一步的，所述步骤(8)中重新采样过程如下：
S21：根据序号累加粒子权值，获得各粒子的累积分布函数值集：
Figure GDA0002663941780000042
S22：计算所需粒子数Nk+1：
Figure GDA0002663941780000043
式中：
Figure GDA0002663941780000044
n为单元方差对应的粒子个数，
Figure GDA0002663941780000045
为最小粒子个数；
S23：生成均匀的或随机的累积分布函数值：
Figure GDA0002663941780000046
S24：依次将粒子序号对应的累积分布函数值，和均匀或随机的累积分布函数值进行对比；对于m＝1，i＝1，如
Figure GDA0002663941780000047
则删除第i个粒子，i＝i+1，否则复制第i个粒子到新的粒子集，m＝m+1；直到m＝Nk+1，得到新的粒子集为
Figure GDA0002663941780000048
S25：设置新的粒子集为等权：
Figure GDA0002663941780000049
得到新的粒子集及权值。
进一步的，所述步骤(9)中正则化的过程如下：
S31：确定粒子正则化的核函数：
Figure GDA00026639417800000410
式中：nx为未知向量的维数，当x是标量时值为1；
Figure GDA00026639417800000411
为nx空间单元球体的体积；
S32：根据粒子滤波中的维数计算最优带宽hopt：
Figure GDA0002663941780000051
S33：根据粒子方差计算均方根
Figure GDA0002663941780000052
S34：对核函数不为零部分采样，产生集合
Figure GDA0002663941780000053
根据采样值，最优带宽和均方根，根据下式获得新的粒子值：
Figure GDA0002663941780000054
进一步的，所述步骤2中单历元法方程或与之前历元累加法方程建立过程如下：
GLONASS系统伪距非差观测方程为：
Figure GDA0002663941780000055
GLONASS系统相位非差观测方程为：
Figure GDA0002663941780000056
式中：i为卫星序号，a为观测站序号，P为GNSS卫星的非差伪距观测值，Φ为GNSS卫星的非差相位观测值，c为光速，δta为GNSS观测站接收机钟差，ρ为观测站到GNSS卫星之间的距离，δti为GNSS卫星钟差，di a为接收机端伪距硬件延迟，di为GNSS卫星端伪距硬件延迟，I为电离层延迟误差，T为对流层延迟误差，ε为伪距观测值的观测噪声，μi a为接收机端相位硬件延迟，μi为GNSS卫星端相位硬件延迟，λi为第i颗卫星的载波波长，Ni a为整周模糊度，ζ为相位观测值的观测噪声。
对上述观测值进行双差组合，消除卫星钟差，接收机钟差，改正电离层延迟误差和对流层延迟误差；得到GNSS系统内双差伪距观测方程：
Figure GDA0002663941780000057
GNSS系统内双差载波相位观测方程：
Figure GDA0002663941780000058
式中：s1为任一卫星系统，b为双差观测值的另一测站的测站号，j为组成双差观测值的另一GNSS卫星的卫星号；
GNSS系统间双差伪距观测方程：
Figure GDA0002663941780000059
GNSS系统间双差载波相位观测方程：
Figure GDA0002663941780000061
式中：s2为s1之外的另一卫星系统，d为伪距系统间偏差，μ为所求的相位系统间偏差；
GNSS系统双差伪距观测方程和GNSS系统双差载波相位观测方程线性化后可转化为：
v＝Ax+Db+Cγ+l
式中：x为除模糊度和频间偏差外其他未知量包括测站坐标分量组成的矢量，b为接收机间单差模糊度未知数矢量，γ为频间偏差未知量，A、D和C分别为未知量对应系数矩阵，l为常数项矢量，P为权矩阵，v为观测值残差矢量；
根据线性化方程可可单历元法方程或与之前历元累加法方程：
Figure GDA0002663941780000062
本发明的有益效果是：
(1)本发明在GNSS系统间偏差发生变化时仍然适用，可高效的实现系统间偏差的实时追踪和系统无差别的GNSS多系统组合精密定位；
(2)本发明在权更新之前加入聚簇分析，可在当前历元解决粒子分周期收敛的问题；
(3)本发明滤波所用粒子个数依据粒子分布状况实时调节，降低计算量；
(4)本发明精确估计时，采用二次收敛方法，第一次收敛后的粒子随着滤波的推进，可以拥有足够小的均方根，从而获得均方根小于第二次收敛阈值的粒子，达到精确估计的目的。
附图说明
图1为本发明步骤3中实时追踪的流程示意图。
图2为本发明步骤3中精确估计的流程示意图。
图3为本发明实施例中GNSS相位系统间偏差实时追踪的频间偏差率追踪结果；
图4为本发明实施例中GNSS相位系统间偏差精确估计的频间偏差率估计结果；
图5为本发明实施例中GNSS相位系统间偏差实时追踪和精密估计方法改正相位偏差后的基线解算结果。
具体实施方式
下面结合附图和具体实施例对本发明做进一步说明。
如图1-2所示，一种基于粒子滤波的GNSS相位系统间偏差实时追踪和精密估计方法，包括以下步骤：
步骤1：对卫星导航数据进行预处理，导入卫星星历、当前历元伪距观测值和当前历元相位观测值；预处理包括数据格式转换、粗差探测与剔除，周跳探测与修复等处理；
GNSS的伪距非差观测方程表示为：
Figure GDA0002663941780000071
GNSS系统相位非差观测方程为：
Figure GDA0002663941780000072
式中：i为卫星序号，a为观测站序号，P为GNSS卫星的非差伪距观测值，Φ为GNSS卫星的非差相位观测值，c为光速，δta为GNSS观测站接收机钟差，ρ为观测站到GNSS卫星之间的距离，δti为GNSS卫星钟差，di a为接收机端伪距硬件延迟，di为GNSS卫星端伪距硬件延迟，I为电离层延迟误差，T为对流层延迟误差，ε为伪距观测值的观测噪声，μi a为接收机端相位硬件延迟，μi为GNSS卫星端相位硬件延迟，λi为第i颗卫星的载波波长，Ni a为整周模糊度，ζ为相位观测值的观测噪声；上述观测方程包括载波L1和L2上的观测值，此处并没有区分GNSSL1和L2频带的观测值；硬件延迟di a和μi a实际包含信号硬件路径延迟、数字信号处理过程延迟和接收机初始相位延迟等。
步骤2：建立系统内(Intra-System)和系统间(Inter-System)的观测值方程并线性化；得到单历元法方程或与之前历元累加法方程；
首先，对上述观测值进行双差组合，消除卫星钟差，接收机钟差，改正电离层延迟误差和对流层延迟误差；
对于同一卫星系统的观测值，得到GNSS系统内双差伪距观测方程：
Figure GDA0002663941780000073
GNSS系统内双差载波相位观测方程：
Figure GDA0002663941780000074
式中：s1为任一卫星系统，b为双差观测值的另一测站的测站号，j为组成双差观测值的另一GNSS卫星的卫星号；
对于不同卫星系统的观测值，得到GNSS系统间双差伪距观测方程：
Figure GDA0002663941780000075
GNSS系统间双差载波相位观测方程，
Figure GDA0002663941780000081
式中：s2为s1之外的另一卫星系统，d为伪距系统间偏差，又伪距观测值直接求得，μ为所求的相位系统间偏差，需要精确估计。
对上述GNSS系统双差伪距观测方程和GNSS系统双差载波相位观测方程进行精确的电离层延迟误差和对流程延迟误差改正，在短基线的情况下，可忽略此改正；将方程线性化后可表示为：
v＝Ax+Db+Cγ+l
式中：x为除模糊度和频间偏差外其他未知量包括测站坐标分量组成的矢量，b为接收机间单差模糊度未知数矢量，γ为频间偏差未知量，A、D和C分别为未知量对应系数矩阵，l为常数项矢量，P为权矩阵，v为观测值残差矢量；
根据线性化方程可可单历元法方程或与之前历元累加法方程：
Figure GDA0002663941780000082
对于已知频间偏差
Figure GDA0002663941780000083
可通过下式进行更正：
Figure GDA0002663941780000084
步骤3：依次使用粒子值改正法方程中的相位偏差小数部分(Fractional ISB，F-ISB)；解算法方程；使用LAMBDA方法进行模糊度固定并输出RATIO值；建立关于RATIO值的函数，用函数值更新粒子滤波权；根据带权粒子，计算相位系统间偏差小数部分的数值及粒子均方根；在精密估计时采取二次收敛的方法，设定两级收敛阈值；其中实时跟踪和精确估计的流程如图1和图2所示。
具体步骤如下：
(1)在粒子滤波准备阶段，依据偏差的取值范围和分布特征，采样生成初始粒子集
Figure GDA0002663941780000085
对于第k时刻粒子集由上一时刻滤波结果生成，其中：x为粒子值代表GLONASS频间偏差率γ的数值，w为粒子对应权值，N为粒子个数，i＝1，2…N为粒子序号；
(2)计算粒子均方根，判断均方根是否大于给定阈值stdgroup，若大于阈值则对粒子进行聚簇分析，通过聚簇分析整合分簇粒子；即把粒子分成两个群，如果粒子群重心之间的距离近似于一个载波波长，合并两个粒子群。
具体步骤如下：
S1：选择粒子值最大和最小的两个粒子作为起始粒子，计算其它粒子到起始粒子的距离，按照距离大小将粒子分为两个组。
S2：计算两个粒子组的重心g1、g2，其定义如下：
Figure GDA0002663941780000091
其中：xk为k时刻粒子值，wk为k时刻粒子对应权值，h＝1，2为组号，Nh为各组的粒子个数；
S3：计算粒子组重心距离：
d＝|g1-g2|；
S3：判断是否近似等于载波波长λ，即判断|d-λ|＜eps是否成立，其中eps是一个相对于载波波长较小的值，如果成立则将其中的一组粒子通过增加或减去一个波长，转移到另一组，实现两个粒子组的合并。
(3)对于每一个粒子值，改正GNSS观测值法方程系统间偏差；解算单历元法方程或与之前历元累加法方程，获得未知量的负电解和相应的协方差阵(variance-covariancematrix)；使用LAMBDA方法进行模糊度固定并输出对应粒子的RATIO值；
(4)建立关于RATIO值的似然函数，根据RATIO的值，计算似然函数的值；用函数值进行粒子滤波权更新，并标准化粒子权值，作为新的粒子权值；
粒子滤波权更新过程如下：
S11：建立似然函数与RATIO之间的函数关系：
Figure GDA0002663941780000092
式中：f(RATIO)是关于RATIO值的函数；
S12：根据建立的函数关系和第i个粒子对应的RATIO值RATIOi计算对应粒子的似然函数值
Figure GDA0002663941780000093
可取：
Figure GDA0002663941780000094
S13：将似然函数值与对应粒子的权值相乘，获得更新后的粒子权
Figure GDA0002663941780000095
Figure GDA0002663941780000096
S14：标准化粒子权值，即将每个粒子的权与所有粒子权之和的比值，作为新的粒子权值
Figure GDA0002663941780000101
Figure GDA0002663941780000102
(5)计算粒子的期望值作为相位偏差的估计值
Figure GDA0002663941780000103
Figure GDA0002663941780000104
计算粒子的方差
Figure GDA0002663941780000105
Figure GDA0002663941780000106
(6)判断粒子滤波是否收敛，实时追踪时，判断均方根是否小于设定阈值stdthd，若是则输出相位偏差的估值和粒子的方差作为估计结果，进入步骤(8)；
(7)判断粒子滤波是否收敛，精确估计时，判断均方根是否小于设定阈值std2thd，若是则输出相位偏差的估计值和粒子的方差作为估值结果，退出滤波，进入步骤4；若否则进入步骤(8)；
(8)若满足重采样条件，则对得到的粒子，设置总的粒子数Nk+1，根据更新的权值重新采样；转入步骤(9)；
重采样过程如下：
S21：根据序号累加粒子权值，获得各粒子的累积分布函数值集：
Figure GDA0002663941780000107
S22：计算所需粒子数Nk+1：
Figure GDA0002663941780000108
式中：
Figure GDA0002663941780000109
n为单元方差对应的粒子个数，
Figure GDA00026639417800001010
为最小粒子个数，即粒子个数的下限；
S23：生成均匀的或随机的累积分布函数值：
Figure GDA0002663941780000111
S24：依次将粒子序号对应的累积分布函数值，和均匀或随机的累积分布函数值进行对比；对于m＝1，i＝1，如
Figure GDA0002663941780000112
则删除第i个粒子，i＝i+1，否则复制第i个粒子到新的粒子集，m＝m+1；直到m＝Nk+1，得到新的粒子集为
Figure GDA0002663941780000113
S25：设置新的粒子集为等权：
Figure GDA0002663941780000114
得到新的粒子集及权值，将生成的新的粒子集及其权值输出到下一步骤。
(9)在精确估计时，判断均方根是否小于设定阈值std1thd，若是则对重采样的粒子实施正则化，转入步骤(11)；若否则转入步骤(10)；实时追踪时，直接转入步骤(10)；
正则化过程如下：
S31：确定粒子正则化的核函数：
Figure GDA0002663941780000115
式中：nx为未知向量的维数，当x是标量时值为1；
Figure GDA0002663941780000116
为nx空间单元球体的体积；
S32：根据粒子滤波中的维数计算最优带宽hopt：
Figure GDA0002663941780000117
S33：根据粒子方差计算均方根
Figure GDA0002663941780000118
S34：对核函数不为零部分采样，产生集合
Figure GDA0002663941780000119
根据采样值，最优带宽和均方根，根据下式获得新的粒子值：
Figure GDA00026639417800001110
将粒子值组成的新的集合输出到下一步骤。
(10)预计下一时刻粒子，对重采样的粒子实时离散化：
Figure GDA0002663941780000121
其中：e为离散化时所加的随机噪声；
(11)预计下一时刻粒子过程如下：
S41：建立偏差的系统方程，设定系统方程噪声：
xk+1＝h(xk，σ）
式中：h为系统方程，σ为噪声水平；
S42：根据系统方程及系统噪声，推算下一历元的粒子值：
Figure GDA0002663941780000122
这里可直接取
Figure GDA0002663941780000123
因为σ很小甚至为零，所以在实施追踪时，(10)和(11)可简化为
Figure GDA0002663941780000124
S33：将推算的新的粒子集
Figure GDA0002663941780000125
输出到下一步骤(11)推算下一历元的粒子值。
步骤4：实时追踪相位系统间偏差F-ISB的值；或精确估计偏差值，待滤波收敛后，用计算出的F-ISB值改正系统间观测方程或法方程中相位系统间偏差；固定系统内和系统间的双差整周模糊度，实现系统无差别的GNSS实时或事后精密定位。
利用本发明的方法对一条短基线进行处理，利用实时追踪程序，获得的某条短基线的GPS L1频率和Galileo E1频率的系统间偏差结果序列如图3所示，均方根均值为0.0032m；利用精确估计程序，获得的该基线的GPS L1频率和Galileo E1频率的系统间偏差结果序列如图4所示，均方根均值为0.0007m；利用估计的GPS L1和Galileo E1之间的系统间偏差值，获得的基线解算结果坐标序列如图5所示，单历元模糊度固定率由77.1％提高到81.2％。
本发明提出一种基于粒子滤波的实时追踪和精确估计GNSS相位系统间偏差的方法，即使在GNSS相位系统间偏差发生数值变化时仍然适用；可以高效率地实现相位系统间偏差的实时追踪和系统无差别的GNSS多系统精密定位；该方法在权更新之前加入聚簇分析方法，可以在当前历元解决粒子分周期收敛的问题；该方法中滤波所用粒子个数依据粒子分布状况实时调节，在需要较多粒子时能够增加粒子数，起到探索的目的；在需要较少粒子时，减少粒子数，降低计算量；该方法在精确估计时采用二次收敛操作，第一次收敛后的粒子随着滤波的推进，可以拥有足够小的均方根，从而获得均方根小于第二次收敛阈值的粒子集，达到精确估计的目的；这种基于粒子滤波的实时追踪和精确估计GNSS系统间偏差的方法是建立在完备的理论基础之上，具有普遍适用性。
文中：LAMBDA是Least-squares AMBiguity Decorrelation Adjustment的缩写，指最小二乘模糊度去相关平差；RATIO在GNSS领域单指整周模糊度固定时的一个可靠性检验指标。

Claims (5)

1.一种基于粒子滤波的GNSS相位系统间偏差实时追踪和精密估计方法，其特征在于，包括以下步骤：
步骤1：对卫星导航数据进行预处理，导入卫星星历、当前历元伪距观测值和当前历元相位观测值；
步骤2：建立系统内和系统间的观测值方程，并线性化，得到单历元法方程或与之前历元累加法方程；
步骤3：计算粒子均方根，判断均方根是否大于给定阈值，若大于阈值则对粒子进行聚簇分析，通过聚簇分析整合分簇粒子；对每个粒子值，改正GNSS观测值法方程系统间偏差；解算法方程，使用LAMBDA方法进行模糊度固定并输出对应粒子的RATIO值；建立关于RATIO值的函数，用函数值更新粒子滤波权；根据带权粒子，计算相位系统间偏差小数部分的数值及粒子均方根；
步骤4：重复步骤1-3，实时追踪GNSS系统间相位偏差小数部分的值，实现实时追踪；重复步骤1-3，待滤波收敛后，通过估计的GNSS系统间相位偏差小数部分，改正观测值或法方程中相位系统间偏差值，固定系统内和系统间的双差整周模糊度，实现精密定位；
所述步骤3中通过聚簇分析整合分簇粒子过程如下：
S1：选择粒子值最大和最小的两个粒子作为起始粒子，计算其它粒子到起始粒子的距离，按照距离大小将粒子分为两组；
S2：分别计算两个粒子组的加权重心g1、g2，其定义如下：
Figure FDA0002663941770000011
其中：xk为k时刻粒子值，wk为k时刻粒子对应权值，Nh为各组粒子个数，h为组号，h＝1，2；
S3：计算粒子组重心距离d：
d＝|g1-g2|；
S4：判断距离d与载波波长λ之间的差值是否小于eps，eps为较小的正数，是一个预设的阈值，若|d-λ|＜eps成立，则将两个粒子组进行合并；
所述步骤3具体过程如下：
(1)采样生成初始粒子集
Figure FDA0002663941770000012
对于第k时刻粒子集由上一时刻滤波结果生成，其中：x为粒子值，w为粒子对应权值，N为粒子个数，i＝1，2…N为粒子序号；
(2)计算粒子均方根，判断均方根是否大于给定阈值，若大于阈值则对粒子进行聚簇分析，通过聚簇分析整合分簇粒子；
(3)对于每一个粒子值，改正GNSS观测值法方程系统间偏差；解算法方程，获得未知量的负电解和相应的协方差阵；使用LAMBDA方法进行模糊度固定并输出对应粒子的RATIO值；
(4)建立关于RATIO值的似然函数，用函数值进行粒子滤波权更新，并标准化粒子权值，作为新的粒子权值；
(5)计算粒子的期望值作为相位偏差的估计值
Figure FDA0002663941770000021
Figure FDA0002663941770000022
计算粒子的方差
Figure FDA0002663941770000023
Figure FDA0002663941770000024
(6)判断粒子滤波是否收敛，实时追踪时，判断均方根是否小于设定阈值stdthd，若是则输出相位偏差的估值和粒子的方差作为估计结果，进入步骤(8)；
(7)判断粒子滤波是否收敛，精确估计时，判断均方根是否小于设定阈值std2thd，若是则输出相位偏差的估计值和粒子的方差作为估值结果，退出滤波，进入步骤4；若否则进入步骤(8)；
(8)如果满足重采样条件，根据更新的权值重新采样，转入步骤(9)；
(9)在精确估计时，判断均方根是否小于设定阈值std1thd，若是则对重采样的粒子实施正则化，转入步骤(11)；若否则转入步骤(10)；实时追踪时，直接转入步骤(10)；
(10)预计下一时刻粒子，对重采样的粒子实施离散化：
Figure FDA0002663941770000025
其中：e为离散化时所加的随机噪声；
(11)推算下一历元的粒子值。
2.根据权利要求1所述的一种基于粒子滤波的GNSS相位系统间偏差实时追踪和精密估计方法，其特征在于，所述步骤(4)中粒子滤波权更新过程如下：
S11：建立似然函数与RATIO之间的函数关系：
Figure FDA0002663941770000026
式中：f(RATIO)是关于RATIO值的函数；
S12：根据建立的函数关系和第i个粒子对应的RATIO值RATIOi计算对应粒子的似然函数值
Figure FDA0002663941770000031
Figure FDA0002663941770000032
S13：将似然函数值与对应粒子的权值相乘，获得更新后的粒子权
Figure FDA0002663941770000033
Figure FDA0002663941770000034
S14：标准化粒子权值，即将每个粒子的权与所有粒子权之和的比值，作为新的粒子权值
Figure FDA0002663941770000035
Figure FDA0002663941770000036
3.根据权利要求1所述的一种基于粒子滤波的GNSS相位系统间偏差实时追踪和精密估计方法，其特征在于，所述步骤(8)中重新采样过程如下：
S21：根据序号累加粒子权值，获得各粒子的累积分布函数值集：
Figure FDA0002663941770000037
S22：计算所需粒子数Nk+1：
Figure FDA0002663941770000038
式中：
Figure FDA0002663941770000039
n为单元方差对应的粒子个数，
Figure FDA00026639417700000310
为最小粒子个数；
S23：生成均匀的或随机的累积分布函数值：
Figure FDA00026639417700000311
S24：依次将粒子序号对应的累积分布函数值，和均匀或随机的累积分布函数值进行对比；对于m＝1，i＝1，如
Figure FDA00026639417700000312
则删除第i个粒子，i＝i+1，否则复制第i个粒子到新的粒子集，m＝m+1；直到m＝Nk+1，得到新的粒子集为
Figure FDA00026639417700000313
S25：设置新的粒子集为等权：
Figure FDA0002663941770000041
得到新的粒子集及权值。
4.根据权利要求1所述的一种基于粒子滤波的GNSS相位系统间偏差实时追踪和精密估计方法，其特征在于，所述步骤(9)中正则化的过程如下：
S31：确定粒子正则化的核函数：
Figure FDA0002663941770000042
式中：nx为未知向量的维数，当x是标量时值为1；
Figure FDA0002663941770000043
为nx空间单元球体的体积；
S32：根据粒子滤波中的维数计算最优带宽hopt：
Figure FDA0002663941770000044
S33：根据粒子方差计算均方根
Figure FDA0002663941770000045
S34：对核函数不为零部分采样，产生集合
Figure FDA0002663941770000046
根据采样值，最优带宽和均方根，根据下式获得新的粒子值：
Figure FDA0002663941770000047
5.根据权利要求1所述的一种基于粒子滤波的GNSS相位系统间偏差实时追踪和精密估计方法，其特征在于，所述步骤2中单历元法方程或与之前历元累加法方程建立过程如下：
GLONASS系统伪距非差观测方程为：
Figure FDA0002663941770000048
GLONASS系统相位非差观测方程为：
Figure FDA0002663941770000049
式中：i为卫星序号，a为观测站序号，P为GNSS卫星的非差伪距观测值，Φ为GNSS卫星的非差相位观测值，c为光速，δta为GNSS观测站接收机钟差，ρ为观测站到GNSS卫星之间的距离，δti为GNSS卫星钟差，di a为接收机端伪距硬件延迟，di为GNSS卫星端伪距硬件延迟，I为电离层延迟误差，T为对流层延迟误差，ε为伪距观测值的观测噪声，μi a为接收机端相位硬件延迟，μi为GNSS卫星端相位硬件延迟，λi为第i颗卫星的载波波长，Ni a为整周模糊度，ζ为相位观测值的观测噪声；
对上述观测值进行双差组合，消除卫星钟差，接收机钟差，改正电离层延迟误差和对流层延迟误差；得到GNSS系统内双差伪距观测方程：
Figure FDA0002663941770000051
GNSS系统内双差载波相位观测方程：
Figure FDA0002663941770000052
式中：s1为任一卫星系统，b为双差观测值的另一测站的测站号，j为组成双差观测值的另一GNSS卫星的卫星号；
GNSS系统间双差伪距观测方程：
Figure FDA0002663941770000053
GNSS系统间双差载波相位观测方程：
Figure FDA0002663941770000054
式中：s2为s1之外的另一卫星系统，d为伪距系统间偏差，μ为所求的相位系统间偏差；
GNSS系统双差伪距观测方程和GNSS系统双差载波相位观测方程线性化后可转化为：
v＝Ax+Db+Cγ+l
式中：x为除模糊度和频间偏差外其他未知量包括测站坐标分量组成的矢量，b为接收机间单差模糊度未知数矢量，γ为系统间偏差未知量，A、D和C分别为未知量对应系数矩阵，l为常数项矢量，P为权矩阵，v为观测值残差矢量；
根据线性化方程可可单历元法方程或与之前历元累加法方程：
Figure FDA0002663941770000055
CN201710790184.XA 2017-09-05 2017-09-05 基于粒子滤波的gnss相位系统间偏差实时追踪和精密估计方法 Active CN107728171B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710790184.XA CN107728171B (zh) 2017-09-05 2017-09-05 基于粒子滤波的gnss相位系统间偏差实时追踪和精密估计方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710790184.XA CN107728171B (zh) 2017-09-05 2017-09-05 基于粒子滤波的gnss相位系统间偏差实时追踪和精密估计方法

Publications (2)

Publication Number Publication Date
CN107728171A CN107728171A (zh) 2018-02-23
CN107728171B true CN107728171B (zh) 2020-10-23

Family

ID=61205759

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710790184.XA Active CN107728171B (zh) 2017-09-05 2017-09-05 基于粒子滤波的gnss相位系统间偏差实时追踪和精密估计方法

Country Status (1)

Country Link
CN (1) CN107728171B (zh)

Families Citing this family (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108802780A (zh) * 2018-03-09 2018-11-13 东南大学 一种gps/bds差分系统间偏差特性分析方法
CN108594277B (zh) * 2018-04-27 2020-06-12 北京邮电大学 一种相位差确定方法、装置、电子设备及存储介质
CN112014862B (zh) * 2019-05-30 2024-03-29 上海海积信息科技股份有限公司 一种载波相位观测数据生成方法及装置
CN110988935B (zh) * 2019-11-25 2021-11-12 重庆市地理信息和遥感应用中心(重庆市测绘产品质量检验测试中心) 基于接收机端偏差聚类优化的多系统组合精密定位方法
CN113556304B (zh) * 2021-06-02 2022-11-04 北京大学 基于粒子滤波器的时变频偏估计方法、系统及介质
CN113534210B (zh) * 2021-06-07 2022-05-31 湖南北斗微芯产业发展有限公司 一种基于混合卡尔曼滤波的模糊度固定方法
CN113504557B (zh) * 2021-06-22 2023-05-23 北京建筑大学 面向实时应用的gps频间钟差新预报方法
CN113589349A (zh) * 2021-09-30 2021-11-02 中国地质大学（武汉） 基于粒子滤波的gnss实时高精度单差测姿方法
CN114488227B (zh) * 2022-01-26 2023-10-20 西南交通大学 一种基于空间相关性的多路径误差改正方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102353969A (zh) * 2011-09-02 2012-02-15 东南大学 精密单点定位技术中相位偏差的估计方法
CN105738926A (zh) * 2016-03-30 2016-07-06 武汉大学 一种glonass系统接收机间相位频间偏差标定方法
CN106249256A (zh) * 2016-07-08 2016-12-21 辽宁工程技术大学 基于粒子群优化算法的实时glonass相位偏差估计方法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8803736B2 (en) * 2010-02-26 2014-08-12 Navcom Technology, Inc. Method and system for estimating position with bias compensation
ES2776723T3 (es) * 2014-10-16 2020-07-31 Gmv Aerospace And Defence S A Método para computar una cota de error de una solución de posición de GNSS basada en filtro de Kalman

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102353969A (zh) * 2011-09-02 2012-02-15 东南大学 精密单点定位技术中相位偏差的估计方法
CN105738926A (zh) * 2016-03-30 2016-07-06 武汉大学 一种glonass系统接收机间相位频间偏差标定方法
CN106249256A (zh) * 2016-07-08 2016-12-21 辽宁工程技术大学 基于粒子群优化算法的实时glonass相位偏差估计方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
"Particle filter-based estimation of inter-system phase bias for real-time integer ambiguity resolution";Yumiao Tian et al.;《Springer》;20161119;第950页右栏－第951页左栏，第956页 *
" Online Estimation of Inter-Frequency/System Phase Biases in Precise Positioning ";Yumiao Tian;《Springer》;20161024;第1.2节，第2.2.1节，第2.3.2节，第2.3.3节，第2.4节，第5.3.4节，第5.4节，第6.3节 *

Also Published As

Publication number Publication date
CN107728171A (zh) 2018-02-23

Similar Documents

Publication Publication Date Title
CN107728171B (zh) 基于粒子滤波的gnss相位系统间偏差实时追踪和精密估计方法
CN107678050B (zh) 基于粒子滤波的glonass相位频间偏差实时追踪和精密估计方法
CN107728180B (zh) 一种基于多维粒子滤波偏差估计的gnss精密定位方法
CN109581452B (zh) 一种gnss参考站载波相位整周模糊度解算方法
CN108931915B (zh) 利用导航卫星的授时方法和装置、计算机可读存储介质
US10895646B2 (en) Outlier-tolerant navigation satellite system positioning method and system
AU2015227414B2 (en) Precise GNSS positioning system with improved ambiguity estimation
CN108508461B (zh) 基于gnss载波相位高精度定位完好性监测方法
CN111045034B (zh) 基于广播星历的gnss多系统实时精密时间传递方法及系统
CN108254773B (zh) 一种多gnss的实时钟差解算方法
CN109001781B (zh) 一种顾及电离层约束的bds三频模糊度解算方法
CN108828640B (zh) 一种卫星导航定位观测值定权方法及装置
CN109581455B (zh) 一种bds和gps融合的三频宽巷紧组合定位方法
CN107966722B (zh) 一种gnss钟差解算方法
CN111694030A (zh) 一种基于格网虚拟观测值的bds局域差分方法及系统
CN110346823B (zh) 可用于北斗精密单点定位的三频模糊度解算方法
CN109061694B (zh) 一种基于gnss钟差固定的低轨导航增强定位方法及系统
CN113325446B (zh) 一种多模共频gnss载波相位时间传递方法及系统
CN110824505B (zh) Gnss卫星接收机的偏差估计方法及系统、定位方法及终端
CN114935770B (zh) 一种多历元加快精密单点定位收敛速度的方法及装置
CN112146557A (zh) 一种基于gnss的实时桥梁变形监测系统及方法
CN110646823A (zh) 一种基于Helmet后验定权法的GPS\BDS紧组合精密单点定位方法
CN110988935B (zh) 基于接收机端偏差聚类优化的多系统组合精密定位方法
CN116359968A (zh) 一种联合北斗二号和北斗三号的三频差分定位方法
CN115480279A (zh) Gnss导航方法和终端、组合导航系统、存储介质

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