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

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

Info

Publication number
CN107728171A
CN107728171A CN201710790184.XA CN201710790184A CN107728171A CN 107728171 A CN107728171 A CN 107728171A CN 201710790184 A CN201710790184 A CN 201710790184A CN 107728171 A CN107728171 A CN 107728171A
Authority
CN
China
Prior art keywords
mrow
msub
msubsup
particle
msup
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
CN201710790184.XA
Other languages
English (en)
Other versions
CN107728171B (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

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)
  • Mathematical Physics (AREA)
  • Data Mining & Analysis (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Computational Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Software Systems (AREA)
  • Operations Research (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • General Engineering & Computer Science (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Power Engineering (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Evolutionary Biology (AREA)
  • Probability & Statistics with Applications (AREA)
  • Position Fixing By Use Of Radio Waves (AREA)
  • Other Investigation Or Analysis Of Materials By Electrical Means (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，其定义如下：
其中：xk为k时刻粒子值，wk为k时刻粒子对应权值，Nh为各组粒子个数，h为组号，h＝1，2；
S3：计算粒子组重心距离d：
d＝|g1-g2|；
S4：判断距离d与载波波长λ之间的差值是否小于eps，若|d-λ|＜eps成立，则将两个粒子组进行合并。
进一步的，所述步骤3具体过程如下：
(1)采样生成初始粒子集对于第k时刻粒子集由上一时刻滤波结果生成，其中：x为粒子值，w为粒子对应权值，N为粒子个数，i＝1，2…N为粒子序号；
(2)计算粒子均方根，判断均方根是否大于给定阈值，若大于阈值则对粒子进行聚簇分析，通过聚簇分析整合分簇粒子；
(3)对于每一个粒子值，改正GNSS观测值法方程系统间偏差；解算法方程，获得未知量的负电解和相应的协方差阵；使用LAMBDA方法进行模糊度固定并输出对应粒子的RATIO值；
(4)建立关于RATIO值的似然函数，用函数值进行粒子滤波权更新，并标准化粒子权值，作为新的粒子权值；
(5)计算粒子的期望值作为相位偏差的估计值
计算粒子的方差
(6)判断粒子滤波是否收敛，实时追踪时，判断均方根是否小于设定阈值stdthd，若是则输出相位偏差的估值和粒子的方差作为估计结果，进入步骤(8)；若否则直接进入步骤(8)；
(7)判断粒子滤波是否收敛，精确估计时，判断均方根是否小于设定阈值std2thd，若是则输出相位偏差的估计值和粒子的方差作为估值结果，退出滤波，进入步骤4；若否则进入步骤(8)；
(8)若满足重采样的条件，根据更新的权值重新采样，转入步骤(9)；
(9)在精确估计时，判断均方根是否小于设定阈值std1thd，若是则对重采样的粒子实施正则化，转入步骤(11)；若否则转入步骤(10)；实时追踪时，直接转入步骤(10)；
(10)预计下一时刻粒子，对重采样的粒子实时离散化：
其中：e为离散化时所加的随机噪声；
(11)推算下一历元的粒子值。
进一步的，所述步骤(4)中粒子滤波权更新过程如下：
S11：建立似然函数与RATIO之间的函数关系：
式中：f(RATIO)是关于RATIO值的函数；
S12：根据建立的函数关系和第i个粒子对应的RATIO值RATIOi计算对应粒子的似然函数值
S13：将似然函数值与对应粒子的权值相乘，获得更新后的粒子权
S14：标准化粒子权值，即将每个粒子的权与所有粒子权之和的比值，作为新的粒子权值
进一步的，所述步骤(8)中重新采样过程如下：
S21：根据序号累加粒子权值，获得各粒子的累积分布函数值集：
S22：计算所需粒子数Nk+1：
式中：n为单元方差对应的粒子个数，为最小粒子个数；
S23：生成均匀的或随机的累积分布函数值：
S24：依次将粒子序号对应的累积分布函数值，和均匀或随机的累积分布函数值进行对比；对于m＝1，i＝1，如则删除第i个粒子，i＝i+1，否则复制第i个粒子到新的粒子集，m＝m+1；直到m＝Nk+1，得到新的粒子集为
S25：设置新的粒子集为等权：
得到新的粒子集及权值。
进一步的，所述步骤(8)中正则化的过程如下：
S31：确定粒子正则化的核函数：
式中：nx为未知向量的维数，当x是标量时值为1；为nx空间单元球体的体积；
S32：根据粒子滤波中的维数计算最优带宽hopt：
S33：根据粒子方差计算均方根
S34：对核函数不为零部分采样，产生集合根据采样值，最优带宽和均方根，根据下式获得新的粒子值：
进一步的，所述步骤2中单历元法方程或与之前历元累加法方程建立过程如下：
GLONASS系统伪距非差观测方程为：
GLONASS系统相位非差观测方程为：
式中：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系统内双差伪距观测方程：
GNSS系统内双差载波相位观测方程：
式中：s1为任一卫星系统，b为双差观测值的另一测站的测站号，j为组成双差观测值的另一GNSS卫星的卫星号；
GNSS系统间双差伪距观测方程：
GNSS系统间双差载波相位观测方程：
式中：s2为s1之外的另一卫星系统，d为伪距系统间偏差，μ为所求的相位系统间偏差；
GNSS系统双差伪距观测方程和GNSS系统双差载波相位观测方程线性化后可转化为：
v＝Ax+Db+Cγ+l
式中：x为除模糊度和频间偏差外其他未知量包括测站坐标分量组成的矢量，b为接收机间单差模糊度未知数矢量，γ为频间偏差未知量，A、D和C分别为未知量对应系数矩阵，l为常数项矢量，P为权矩阵，v为观测值残差矢量；
根据线性化方程可可单历元法方程或与之前历元累加法方程：
本发明的有益效果是：
(1)本发明在GNSS系统间偏差发生变化时仍然适用，可高效的实现系统间偏差的实时追踪和系统无差别的GNSS多系统组合精密定位；
(2)本发明在权更新之前加入聚簇分析，可在当前历元解决粒子分周期收敛的问题；
(3)本发明滤波所用粒子个数依据粒子分布状况实时调节，降低计算量；
(4)本发明精确估计时，采用二次收敛方法，第一次收敛后的粒子随着滤波的推进，可以拥有足够小的均方根，从而获得均方根小于第二次收敛阈值的粒子，达到精确估计的目的。
附图说明
图1为本发明步骤3中实时追踪的流程示意图。
图2为本发明步骤3中精确估计的流程示意图。
图3为本发明实施例中GNSS相位系统间偏差实时追踪的频间偏差率追踪结果；
图4为本发明实施例中GNSS相位系统间偏差精确估计的频间偏差率估计结果；
图5为本发明实施例中GNSS相位系统间偏差实时追踪和精密估计方法改正相位偏差后的基线解算结果。
具体实施方式
下面结合附图和具体实施例对本发明做进一步说明。
如图1-2所示，一种基于粒子滤波的GNSS相位系统间偏差实时追踪和精密估计方法，包括以下步骤：
步骤1：对卫星导航数据进行预处理，导入卫星星历、当前历元伪距观测值和当前历元相位观测值；预处理包括数据格式转换、粗差探测与剔除，周跳探测与修复等处理；
GNSS的伪距非差观测方程表示为：
GNSS系统相位非差观测方程为：
式中：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系统内双差伪距观测方程：
GNSS系统内双差载波相位观测方程：
式中：s1为任一卫星系统，b为双差观测值的另一测站的测站号，j为组成双差观测值的另一GNSS卫星的卫星号；
对于不同卫星系统的观测值，得到GNSS系统间双差伪距观测方程：
GNSS系统间双差载波相位观测方程，
式中：s2为s1之外的另一卫星系统，d为伪距系统间偏差，又伪距观测值直接求得，μ为所求的相位系统间偏差，需要精确估计。
对上述GNSS系统双差伪距观测方程和GNSS系统双差载波相位观测方程进行精确的电离层延迟误差和对流程延迟误差改正，在短基线的情况下，可忽略此改正；将方程线性化后可表示为：
v＝Ax+Db+Cγ+l
式中：x为除模糊度和频间偏差外其他未知量包括测站坐标分量组成的矢量，b为接收机间单差模糊度未知数矢量，γ为频间偏差未知量，A、D和C分别为未知量对应系数矩阵，l为常数项矢量，P为权矩阵，v为观测值残差矢量；
根据线性化方程可可单历元法方程或与之前历元累加法方程：
对于已知频间偏差可通过下式进行更正：
步骤3：依次使用粒子值改正法方程中的相位偏差小数部分(Fractional ISB，F-ISB)；解算法方程；使用LAMBDA方法进行模糊度固定并输出RATIO值；建立关于RATIO值的函数，用函数值更新粒子滤波权；根据带权粒子，计算相位系统间偏差小数部分的数值及粒子均方根；在精密估计时采取二次收敛的方法，设定两级收敛阈值；其中实时跟踪和精确估计的流程如图1和图2所示。
具体步骤如下：
(1)在粒子滤波准备阶段，依据偏差的取值范围和分布特征，采样生成初始粒子集对于第k时刻粒子集由上一时刻滤波结果生成，其中：x为粒子值代表GLONASS频间偏差率γ的数值，w为粒子对应权值，N为粒子个数，i＝1，2…N为粒子序号；
(2)计算粒子均方根，判断均方根是否大于给定阈值stdgroup，若大于阈值则对粒子进行聚簇分析，通过聚簇分析整合分簇粒子；即把粒子分成两个群，如果粒子群重心之间的距离近似于一个载波波长，合并两个粒子群。
具体步骤如下：
S1：选择粒子值最大和最小的两个粒子作为起始粒子，计算其它粒子到起始粒子的距离，按照距离大小将粒子分为两个组。
S2：计算两个粒子组的重心g1、g2，其定义如下：
其中：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之间的函数关系：
式中：f(RATIO)是关于RATIO值的函数；
S12：根据建立的函数关系和第i个粒子对应的RATIO值RATIOi计算对应粒子的似然函数值可取：
S13：将似然函数值与对应粒子的权值相乘，获得更新后的粒子权
S14：标准化粒子权值，即将每个粒子的权与所有粒子权之和的比值，作为新的粒子权值
(5)计算粒子的期望值作为相位偏差的估计值
计算粒子的方差
(6)判断粒子滤波是否收敛，实时追踪时，判断均方根是否小于设定阈值stdthd，若是则输出相位偏差的估值和粒子的方差作为估计结果，进入步骤(8)；若否则直接进入步骤(8)；
(7)判断粒子滤波是否收敛，精确估计时，判断均方根是否小于设定阈值std2thd，若是则输出相位偏差的估计值和粒子的方差作为估值结果，退出滤波，进入步骤4；若否则进入步骤(8)；
(8)若满足重采样条件，则对得到的粒子，设置总的粒子数Nk+1，根据更新的权值重新采样；转入步骤(9)；
重采样过程如下：
S21：根据序号累加粒子权值，获得各粒子的累积分布函数值集：
S22：计算所需粒子数Nk+1：
式中：n为单元方差对应的粒子个数，为最小粒子个数，即粒子个数的下限；
S23：生成均匀的或随机的累积分布函数值：
S24：依次将粒子序号对应的累积分布函数值，和均匀或随机的累积分布函数值进行对比；对于m＝1，i＝1，如则删除第i个粒子，i＝i+1，否则复制第i个粒子到新的粒子集，m＝m+1；直到m＝Nk+1，得到新的粒子集为
S25：设置新的粒子集为等权：
得到新的粒子集及权值，将生成的新的粒子集及其权值输出到下一步骤。
(9)在精确估计时，判断均方根是否小于设定阈值std1thd，若是则对重采样的粒子实施正则化，转入步骤(11)；若否则转入步骤(10)；实时追踪时，直接转入步骤(10)；
正则化过程如下：
S31：确定粒子正则化的核函数：
式中：nx为未知向量的维数，当x是标量时值为1；为nx空间单元球体的体积；
S32：根据粒子滤波中的维数计算最优带宽hopt：
S33：根据粒子方差计算均方根
S34：对核函数不为零部分采样，产生集合根据采样值，最优带宽和均方根，根据下式获得新的粒子值：
将粒子值组成的新的集合输出到下一步骤。
(10)预计下一时刻粒子，对重采样的粒子实时离散化：
其中：e为离散化时所加的随机噪声；
(11)预计下一时刻粒子过程如下：
S41：建立偏差的系统方程，设定系统方程噪声：
xk+1＝h(xk，σ)
式中：h为系统方程，ζ为噪声水平；
S42：根据系统方程及系统噪声，推算下一历元的粒子值：
这里可直接取因为ζ很小甚至为零，所以在实施追踪时，(10)和(11)可简化为
S33：将推算的新的粒子集输出到下一步骤(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 (7)

1.一种基于粒子滤波的GNSS相位系统间偏差实时追踪和精密估计方法，其特征在于，包括以下步骤：
步骤1：对卫星导航数据进行预处理，导入卫星星历、当前历元伪距观测值和当前历元相位观测值；
步骤2：建立系统内和系统间的观测值方程，并线性化，得到单历元法方程或与之前历元累加法方程；
步骤3：计算粒子均方根，判断均方根是否大于给定阈值，若大于阈值则对粒子进行聚簇分析，通过聚簇分析整合分簇粒子；对每个粒子值，改正GNSS观测值法方程系统间偏差；解算法方程，使用LAMBDA方法进行模糊度固定并输出对应粒子的RATIO值；建立关于RATIO值的函数，用函数值更新粒子滤波权；根据带权粒子，计算相位系统间偏差小数部分的数值及粒子均方根；
步骤4：重复步骤1-3，实时追踪GNSS系统间相位偏差小数部分的值，实现实时追踪；重复步骤1-3，待滤波收敛后，通过估计的GNSS系统间相位偏差小数部分，改正观测值或法方程中相位系统间偏差值，固定系统内和系统间的双差整周模糊度，实现精密定位。
2.根据权利要求1所述的一种基于粒子滤波的GNSS相位系统间偏差实时追踪和精密估计方法，其特征在于，所述步骤3中通过聚簇分析整合分簇粒子过程如下：
S1：选择粒子值最大和最小的两个粒子作为起始粒子，计算其它粒子到起始粒子的距离，按照距离大小将粒子分为两组；
S2：分别计算两个粒子组的加权重心g1、g2，其定义如下：
<mrow> <msub> <mi>g</mi> <mi>h</mi> </msub> <mo>=</mo> <mfrac> <mrow> <msubsup> <mi>&amp;Sigma;</mi> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>N</mi> <mi>h</mi> </msub> </msubsup> <msubsup> <mi>x</mi> <mi>k</mi> <mi>i</mi> </msubsup> <msubsup> <mi>w</mi> <mi>k</mi> <mi>i</mi> </msubsup> </mrow> <mrow> <msubsup> <mi>&amp;Sigma;</mi> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>N</mi> <mi>h</mi> </msub> </msubsup> <msubsup> <mi>w</mi> <mi>k</mi> <mi>i</mi> </msubsup> </mrow> </mfrac> </mrow>
其中：xk为k时刻粒子值，wk为k时刻粒子对应权值，Nh为各组粒子个数，h为组号，h＝1，2；
S3：计算粒子组重心距离d：
d＝|g1-g2|；
S4：判断距离d与载波波长λ之间的差值是否小于eps，若|d-λ|＜eps成立，则将两个粒子组进行合并。
3.根据权利要求1所述的一种基于粒子滤波的GNSS相位系统间偏差实时追踪和精密估计方法，其特征在于，所述步骤3具体过程如下：
(1)采样生成初始粒子集对于第k时刻粒子集由上一时刻滤波结果生成，其中：x为粒子值，w为粒子对应权值，N为粒子个数，i＝1，2…N为粒子序号；
(2)计算粒子均方根，判断均方根是否大于给定阈值，若大于阈值则对粒子进行聚簇分析，通过聚簇分析整合分簇粒子；
(3)对于每一个粒子值，改正GNSS观测值法方程系统间偏差；解算法方程，获得未知量的负电解和相应的协方差阵；使用LAMBDA方法进行模糊度固定并输出对应粒子的RATIO值；
(4)建立关于RATIO值的似然函数，用函数值进行粒子滤波权更新，并标准化粒子权值，作为新的粒子权值；
(5)计算粒子的期望值作为相位偏差的估计值
<mrow> <msub> <mover> <mi>x</mi> <mo>^</mo> </mover> <mi>k</mi> </msub> <mo>&amp;ap;</mo> <msubsup> <mi>&amp;Sigma;</mi> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>N</mi> </msubsup> <msubsup> <mi>w</mi> <mi>k</mi> <mi>i</mi> </msubsup> <msubsup> <mi>x</mi> <mi>k</mi> <mi>i</mi> </msubsup> </mrow>
计算粒子的方差
<mrow> <msub> <mover> <mi>P</mi> <mo>^</mo> </mover> <mi>x</mi> </msub> <mo>&amp;ap;</mo> <munderover> <mi>&amp;Sigma;</mi> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>N</mi> </munderover> <msup> <mrow> <mo>(</mo> <mrow> <msubsup> <mi>x</mi> <mi>k</mi> <mi>i</mi> </msubsup> <mo>-</mo> <msub> <mover> <mi>x</mi> <mo>^</mo> </mover> <mi>k</mi> </msub> </mrow> <mo>)</mo> </mrow> <mn>2</mn> </msup> <msubsup> <mi>w</mi> <mi>k</mi> <mi>i</mi> </msubsup> <mo>;</mo> </mrow>
(6)判断粒子滤波是否收敛，实时追踪时，判断均方根是否小于设定阈值stdthd，若是则输出相位偏差的估值和粒子的方差作为估计结果，进入步骤(8)；若否则直接进入步骤(8)；
(7)判断粒子滤波是否收敛，精确估计时，判断均方根是否小于设定阈值std2thd，若是则输出相位偏差的估计值和粒子的方差作为估值结果，退出滤波，进入步骤4；若否则进入步骤(8)；
(8)如果满足重采样条件，根据更新的权值重新采样，转入步骤(9)；
(9)在精确估计时，判断均方根是否小于设定阈值std1thd，若是则对重采样的粒子实施正则化，转入步骤(11)；若否则转入步骤(10)；实时追踪时，直接转入步骤(10)；
(10)预计下一时刻粒子，对重采样的粒子实施离散化：
<mrow> <msubsup> <mi>x</mi> <mi>k</mi> <mi>i</mi> </msubsup> <mo>=</mo> <msubsup> <mi>x</mi> <mi>k</mi> <mi>i</mi> </msubsup> <mo>+</mo> <mi>e</mi> </mrow>
其中：e为离散化时所加的随机噪声；
(11)推算下一历元的粒子值。
4.根据权利要求3所述的一种基于粒子滤波的GNSS相位系统间偏差实时追踪和精密估计方法，其特征在于，所述步骤(4)中粒子滤波权更新过程如下：
S11：建立似然函数与RATIO之间的函数关系：
<mrow> <mi>p</mi> <mrow> <mo>(</mo> <mrow> <msub> <mover> <mi>b</mi> <mo>~</mo> </mover> <mi>k</mi> </msub> <mo>|</mo> <mi>x</mi> </mrow> <mo>)</mo> </mrow> <mo>=</mo> <mi>f</mi> <mrow> <mo>(</mo> <mrow> <mi>R</mi> <mi>A</mi> <mi>T</mi> <mi>I</mi> <mi>O</mi> </mrow> <mo>)</mo> </mrow> </mrow>
式中：f(RATIO)是关于RATIO值的函数；
S12：根据建立的函数关系和第i个粒子对应的RATIO值RATIOi计算对应粒子的似然函数值
<mrow> <mi>p</mi> <mrow> <mo>(</mo> <mrow> <msub> <mover> <mi>b</mi> <mo>~</mo> </mover> <mi>k</mi> </msub> <mo>|</mo> <msup> <mi>x</mi> <mi>i</mi> </msup> </mrow> <mo>)</mo> </mrow> <mo>=</mo> <mi>f</mi> <mrow> <mo>(</mo> <mrow> <msub> <mi>RATIO</mi> <mi>i</mi> </msub> </mrow> <mo>)</mo> </mrow> <mo>;</mo> </mrow>
S13：将似然函数值与对应粒子的权值相乘，获得更新后的粒子权
<mrow> <msubsup> <mover> <mi>w</mi> <mo>&amp;OverBar;</mo> </mover> <mi>k</mi> <mi>i</mi> </msubsup> <mo>=</mo> <msubsup> <mi>w</mi> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> </mrow> <mi>i</mi> </msubsup> <mi>p</mi> <mrow> <mo>(</mo> <mrow> <msub> <mover> <mi>b</mi> <mo>~</mo> </mover> <mi>k</mi> </msub> <mo>|</mo> <msubsup> <mi>x</mi> <mi>k</mi> <mi>i</mi> </msubsup> </mrow> <mo>)</mo> </mrow> <mo>;</mo> </mrow>
S14：标准化粒子权值，即将每个粒子的权与所有粒子权之和的比值，作为新的粒子权值
<mrow> <msubsup> <mi>w</mi> <mi>k</mi> <mi>i</mi> </msubsup> <mo>=</mo> <mfrac> <msubsup> <mover> <mi>w</mi> <mo>&amp;OverBar;</mo> </mover> <mi>k</mi> <mi>i</mi> </msubsup> <mrow> <msubsup> <mi>&amp;Sigma;</mi> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>N</mi> </msubsup> <msubsup> <mover> <mi>w</mi> <mo>&amp;OverBar;</mo> </mover> <mi>k</mi> <mi>j</mi> </msubsup> </mrow> </mfrac> <mo>.</mo> </mrow>
5.根据权利要求3所述的一种基于粒子滤波的GNSS相位系统间偏差实时追踪和精密估计方法，其特征在于，所述步骤(7)中重新采样过程如下：
S21：根据序号累加粒子权值，获得各粒子的累积分布函数值集：
<mrow> <msubsup> <mrow> <mo>{</mo> <mrow> <msubsup> <mi>x</mi> <mi>k</mi> <mi>i</mi> </msubsup> <mo>,</mo> <msubsup> <mi>W</mi> <mi>k</mi> <mi>i</mi> </msubsup> </mrow> <mo>}</mo> </mrow> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>N</mi> </msubsup> <mo>;</mo> </mrow>
S22：计算所需粒子数Nk+1：
式中：n为单元方差对应的粒子个数，为最小粒子个数；
S23：生成均匀的或随机的累积分布函数值：
<mrow> <msubsup> <mrow> <mo>{</mo> <msubsup> <mi>U</mi> <mi>k</mi> <mi>i</mi> </msubsup> <mo>}</mo> </mrow> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>N</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> </msubsup> <mo>;</mo> </mrow>
S24：依次将粒子序号对应的累积分布函数值，和均匀或随机的累积分布函数值进行对比；对于m＝1，i＝1，如则删除第i个粒子，i＝i+1，否则复制第i个粒子到新的粒子集，m＝m+1；直到m＝Nk+1，得到新的粒子集为
S25：设置新的粒子集为等权：
<msubsup> <mrow> <mo>{</mo> <mrow> <msubsup> <mover> <mi>x</mi> <mo>&amp;OverBar;</mo> </mover> <mi>k</mi> <mi>i</mi> </msubsup> <mo>,</mo> <mfrac> <mn>1</mn> <msub> <mi>N</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> </mfrac> </mrow> <mo>}</mo> </mrow> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>N</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> </msubsup>
得到新的粒子集及权值。
6.根据权利要求3所述的一种基于粒子滤波的GNSS相位系统间偏差实时追踪和精密估计方法，其特征在于，所述步骤(8)中正则化的过程如下：
S31：确定粒子正则化的核函数：
式中：nx为未知向量的维数，当x是标量时值为1；为nx空间单元球体的体积；
S32：根据粒子滤波中的维数计算最优带宽hopt：
<mrow> <msub> <mi>h</mi> <mrow> <mi>o</mi> <mi>p</mi> <mi>t</mi> </mrow> </msub> <mo>=</mo> <msup> <mrow> <mo>&amp;lsqb;</mo> <mn>8</mn> <msup> <msub> <mi>c</mi> <msub> <mi>n</mi> <mi>x</mi> </msub> </msub> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> <mrow> <mo>(</mo> <msub> <mi>n</mi> <mi>x</mi> </msub> <mo>+</mo> <mn>4</mn> <mo>)</mo> </mrow> <msup> <mrow> <mo>(</mo> <mn>2</mn> <msqrt> <mi>&amp;pi;</mi> </msqrt> <mo>)</mo> </mrow> <msub> <mi>n</mi> <mi>x</mi> </msub> </msup> <mo>&amp;rsqb;</mo> </mrow> <mfrac> <mn>1</mn> <mrow> <msub> <mi>n</mi> <mi>x</mi> </msub> <mo>+</mo> <mn>4</mn> </mrow> </mfrac> </msup> <msup> <mi>N</mi> <mrow> <mo>-</mo> <mfrac> <mn>1</mn> <mrow> <msub> <mi>n</mi> <mi>x</mi> </msub> <mo>+</mo> <mn>4</mn> </mrow> </mfrac> </mrow> </msup> <mo>;</mo> </mrow>
S33：根据粒子方差计算均方根
S34：对核函数不为零部分采样，产生集合根据采样值，最优带宽和均方根，根据下式获得新的粒子值：
<mrow> <msubsup> <mi>x</mi> <mi>k</mi> <mi>i</mi> </msubsup> <mo>=</mo> <msubsup> <mi>x</mi> <mi>k</mi> <mi>i</mi> </msubsup> <mo>+</mo> <msub> <mi>h</mi> <mrow> <mi>o</mi> <mi>p</mi> <mi>t</mi> </mrow> </msub> <msub> <mi>A</mi> <mi>k</mi> </msub> <msup> <mo>&amp;Element;</mo> <mi>i</mi> </msup> <mo>.</mo> </mrow>
7.根据权利要求1所述的一种基于粒子滤波的GNSS相位系统间偏差实时追踪和精密估计方法，其特征在于，所述步骤2中单历元法方程或与之前历元累加法方程建立过程如下：
GLONASS系统伪距非差观测方程为：
<mrow> <msubsup> <mi>P</mi> <mi>a</mi> <mi>i</mi> </msubsup> <mo>=</mo> <msubsup> <mi>&amp;rho;</mi> <mi>a</mi> <mi>i</mi> </msubsup> <mo>-</mo> <mi>c</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;delta;t</mi> <mi>a</mi> </msub> <mo>-</mo> <msup> <mi>&amp;delta;t</mi> <mi>i</mi> </msup> <mo>)</mo> </mrow> <mo>+</mo> <msubsup> <mi>d</mi> <mi>a</mi> <mi>i</mi> </msubsup> <mo>-</mo> <msup> <mi>d</mi> <mi>i</mi> </msup> <mo>+</mo> <msubsup> <mi>I</mi> <mi>a</mi> <mi>i</mi> </msubsup> <mo>+</mo> <msubsup> <mi>T</mi> <mi>a</mi> <mi>i</mi> </msubsup> <mo>+</mo> <msubsup> <mi>&amp;epsiv;</mi> <mi>a</mi> <mi>i</mi> </msubsup> </mrow>
GLONASS系统相位非差观测方程为：
<mrow> <msup> <mi>&amp;lambda;</mi> <mi>i</mi> </msup> <msubsup> <mi>&amp;Phi;</mi> <mi>a</mi> <mi>i</mi> </msubsup> <mo>=</mo> <msubsup> <mi>&amp;rho;</mi> <mi>a</mi> <mi>i</mi> </msubsup> <mo>-</mo> <mi>c</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;delta;t</mi> <mi>a</mi> </msub> <mo>-</mo> <msup> <mi>&amp;delta;t</mi> <mi>i</mi> </msup> <mo>)</mo> </mrow> <mo>+</mo> <msubsup> <mi>&amp;mu;</mi> <mi>a</mi> <mi>i</mi> </msubsup> <mo>-</mo> <msup> <mi>&amp;mu;</mi> <mi>i</mi> </msup> <mo>+</mo> <msup> <mi>&amp;lambda;</mi> <mi>i</mi> </msup> <msubsup> <mi>N</mi> <mi>a</mi> <mi>i</mi> </msubsup> <mo>-</mo> <msubsup> <mi>I</mi> <mi>a</mi> <mi>i</mi> </msubsup> <mo>+</mo> <msubsup> <mi>T</mi> <mi>a</mi> <mi>i</mi> </msubsup> <mo>+</mo> <msubsup> <mi>&amp;xi;</mi> <mi>a</mi> <mi>i</mi> </msubsup> </mrow>
式中：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系统内双差伪距观测方程：
<mrow> <msubsup> <mi>P</mi> <mrow> <mi>a</mi> <mi>b</mi> </mrow> <mrow> <msub> <mi>s</mi> <mn>1</mn> </msub> <msub> <mi>s</mi> <mn>1</mn> </msub> <mo>,</mo> <mi>i</mi> <mi>j</mi> </mrow> </msubsup> <mo>=</mo> <msubsup> <mi>&amp;rho;</mi> <mrow> <mi>a</mi> <mi>b</mi> </mrow> <mrow> <msub> <mi>s</mi> <mn>1</mn> </msub> <msub> <mi>s</mi> <mn>1</mn> </msub> <mo>,</mo> <mi>i</mi> <mi>j</mi> </mrow> </msubsup> <mo>+</mo> <msubsup> <mi>I</mi> <mrow> <mi>a</mi> <mi>b</mi> </mrow> <mrow> <msub> <mi>s</mi> <mn>1</mn> </msub> <msub> <mi>s</mi> <mn>1</mn> </msub> <mo>,</mo> <mi>i</mi> <mi>j</mi> </mrow> </msubsup> <mo>+</mo> <msubsup> <mi>T</mi> <mrow> <mi>a</mi> <mi>b</mi> </mrow> <mrow> <msub> <mi>s</mi> <mn>1</mn> </msub> <msub> <mi>s</mi> <mn>1</mn> </msub> <mo>,</mo> <mi>i</mi> <mi>j</mi> </mrow> </msubsup> <mo>+</mo> <msubsup> <mi>&amp;epsiv;</mi> <mrow> <mi>a</mi> <mi>b</mi> </mrow> <mrow> <msub> <mi>s</mi> <mn>1</mn> </msub> <msub> <mi>s</mi> <mn>1</mn> </msub> <mo>,</mo> <mi>i</mi> <mi>j</mi> </mrow> </msubsup> </mrow>
GNSS系统内双差载波相位观测方程：
<mrow> <msup> <mi>&amp;lambda;</mi> <mrow> <msub> <mi>s</mi> <mn>1</mn> </msub> <mo>,</mo> <mi>j</mi> </mrow> </msup> <msubsup> <mi>&amp;Phi;</mi> <mrow> <mi>a</mi> <mi>b</mi> </mrow> <mrow> <msub> <mi>s</mi> <mn>1</mn> </msub> <mo>,</mo> <mi>j</mi> </mrow> </msubsup> <mo>-</mo> <msup> <mi>&amp;lambda;</mi> <mrow> <msub> <mi>s</mi> <mn>1</mn> </msub> <mo>,</mo> <mi>i</mi> </mrow> </msup> <msubsup> <mi>&amp;Phi;</mi> <mrow> <mi>a</mi> <mi>b</mi> </mrow> <mrow> <msub> <mi>s</mi> <mn>1</mn> </msub> <mo>,</mo> <mi>i</mi> </mrow> </msubsup> <mo>=</mo> <msubsup> <mi>&amp;rho;</mi> <mrow> <mi>a</mi> <mi>b</mi> </mrow> <mrow> <msub> <mi>s</mi> <mn>1</mn> </msub> <msub> <mi>s</mi> <mn>1</mn> </msub> <mo>,</mo> <mi>i</mi> <mi>j</mi> </mrow> </msubsup> <mo>+</mo> <msup> <mi>&amp;lambda;</mi> <mrow> <msub> <mi>s</mi> <mn>1</mn> </msub> <mo>,</mo> <mi>j</mi> </mrow> </msup> <msubsup> <mi>N</mi> <mrow> <mi>a</mi> <mi>b</mi> </mrow> <mrow> <msub> <mi>s</mi> <mn>1</mn> </msub> <mo>,</mo> <mi>j</mi> </mrow> </msubsup> <mo>-</mo> <msup> <mi>&amp;lambda;</mi> <mrow> <msub> <mi>s</mi> <mn>1</mn> </msub> <mo>,</mo> <mi>i</mi> </mrow> </msup> <msubsup> <mi>N</mi> <mrow> <mi>a</mi> <mi>b</mi> </mrow> <mrow> <msub> <mi>s</mi> <mn>1</mn> </msub> <mo>,</mo> <mi>i</mi> </mrow> </msubsup> <mo>-</mo> <msubsup> <mi>I</mi> <mrow> <mi>a</mi> <mi>b</mi> </mrow> <mrow> <msub> <mi>s</mi> <mn>1</mn> </msub> <msub> <mi>s</mi> <mn>1</mn> </msub> <mo>,</mo> <mi>i</mi> <mi>j</mi> </mrow> </msubsup> <mo>+</mo> <msubsup> <mi>T</mi> <mrow> <mi>a</mi> <mi>b</mi> </mrow> <mrow> <msub> <mi>s</mi> <mn>1</mn> </msub> <msub> <mi>s</mi> <mn>1</mn> </msub> <mo>,</mo> <mi>i</mi> <mi>j</mi> </mrow> </msubsup> <mo>+</mo> <msubsup> <mi>&amp;xi;</mi> <mrow> <mi>a</mi> <mi>b</mi> </mrow> <mrow> <msub> <mi>s</mi> <mn>1</mn> </msub> <msub> <mi>s</mi> <mn>1</mn> </msub> <mo>,</mo> <mi>i</mi> <mi>j</mi> </mrow> </msubsup> </mrow>
式中：s1为任一卫星系统，b为双差观测值的另一测站的测站号，j为组成双差观测值的另一GNSS卫星的卫星号；
GNSS系统间双差伪距观测方程：
<mrow> <msubsup> <mi>P</mi> <mrow> <mi>a</mi> <mi>b</mi> </mrow> <mrow> <msub> <mi>s</mi> <mn>1</mn> </msub> <msub> <mi>s</mi> <mn>2</mn> </msub> <mo>,</mo> <mi>i</mi> <mi>j</mi> </mrow> </msubsup> <mo>=</mo> <msubsup> <mi>&amp;rho;</mi> <mrow> <mi>a</mi> <mi>b</mi> </mrow> <mrow> <msub> <mi>s</mi> <mn>1</mn> </msub> <msub> <mi>s</mi> <mn>2</mn> </msub> <mo>,</mo> <mi>i</mi> <mi>j</mi> </mrow> </msubsup> <mo>+</mo> <msubsup> <mi>I</mi> <mrow> <mi>a</mi> <mi>b</mi> </mrow> <mrow> <msub> <mi>s</mi> <mn>1</mn> </msub> <msub> <mi>s</mi> <mn>2</mn> </msub> <mo>,</mo> <mi>i</mi> <mi>j</mi> </mrow> </msubsup> <mo>+</mo> <msubsup> <mi>T</mi> <mrow> <mi>a</mi> <mi>b</mi> </mrow> <mrow> <msub> <mi>s</mi> <mn>1</mn> </msub> <msub> <mi>s</mi> <mn>2</mn> </msub> <mo>,</mo> <mi>i</mi> <mi>j</mi> </mrow> </msubsup> <mo>+</mo> <msubsup> <mi>&amp;epsiv;</mi> <mrow> <mi>a</mi> <mi>b</mi> </mrow> <mrow> <msub> <mi>s</mi> <mn>1</mn> </msub> <msub> <mi>s</mi> <mn>2</mn> </msub> <mo>,</mo> <mi>i</mi> <mi>j</mi> </mrow> </msubsup> </mrow>
GNSS系统间双差载波相位观测方程：
<mrow> <msup> <mi>&amp;lambda;</mi> <mrow> <msub> <mi>s</mi> <mn>2</mn> </msub> <mo>,</mo> <mi>j</mi> </mrow> </msup> <msubsup> <mi>&amp;Phi;</mi> <mrow> <mi>a</mi> <mi>b</mi> </mrow> <mrow> <msub> <mi>s</mi> <mn>2</mn> </msub> <mo>,</mo> <mi>j</mi> </mrow> </msubsup> <mo>-</mo> <msup> <mi>&amp;lambda;</mi> <mrow> <msub> <mi>s</mi> <mn>1</mn> </msub> <mo>,</mo> <mi>i</mi> </mrow> </msup> <msubsup> <mi>&amp;Phi;</mi> <mrow> <mi>a</mi> <mi>b</mi> </mrow> <mrow> <msub> <mi>s</mi> <mn>1</mn> </msub> <mo>,</mo> <mi>i</mi> </mrow> </msubsup> <mo>=</mo> <msubsup> <mi>&amp;rho;</mi> <mrow> <mi>a</mi> <mi>b</mi> </mrow> <mrow> <msub> <mi>s</mi> <mn>1</mn> </msub> <msub> <mi>s</mi> <mn>2</mn> </msub> <mo>,</mo> <mi>i</mi> <mi>j</mi> </mrow> </msubsup> <mo>+</mo> <msubsup> <mi>&amp;mu;</mi> <mrow> <mi>a</mi> <mi>b</mi> </mrow> <mrow> <msub> <mi>s</mi> <mn>1</mn> </msub> <msub> <mi>s</mi> <mn>2</mn> </msub> <mo>,</mo> <mi>i</mi> <mi>j</mi> </mrow> </msubsup> <mo>+</mo> <msup> <mi>&amp;lambda;</mi> <mrow> <msub> <mi>s</mi> <mn>2</mn> </msub> <mo>,</mo> <mi>j</mi> </mrow> </msup> <msubsup> <mi>N</mi> <mrow> <mi>a</mi> <mi>b</mi> </mrow> <mrow> <msub> <mi>s</mi> <mn>2</mn> </msub> <mo>,</mo> <mi>j</mi> </mrow> </msubsup> <mo>-</mo> <msup> <mi>&amp;lambda;</mi> <mrow> <msub> <mi>s</mi> <mn>1</mn> </msub> <mo>,</mo> <mi>i</mi> </mrow> </msup> <msubsup> <mi>N</mi> <mrow> <mi>a</mi> <mi>b</mi> </mrow> <mrow> <msub> <mi>s</mi> <mn>1</mn> </msub> <mo>,</mo> <mi>i</mi> </mrow> </msubsup> <mo>-</mo> <msubsup> <mi>I</mi> <mrow> <mi>a</mi> <mi>b</mi> </mrow> <mrow> <msub> <mi>s</mi> <mn>1</mn> </msub> <msub> <mi>s</mi> <mn>2</mn> </msub> <mo>,</mo> <mi>i</mi> <mi>j</mi> </mrow> </msubsup> <mo>+</mo> <msubsup> <mi>T</mi> <mrow> <mi>a</mi> <mi>b</mi> </mrow> <mrow> <msub> <mi>s</mi> <mn>1</mn> </msub> <msub> <mi>s</mi> <mn>2</mn> </msub> <mo>,</mo> <mi>i</mi> <mi>j</mi> </mrow> </msubsup> <mo>+</mo> <msubsup> <mi>&amp;xi;</mi> <mrow> <mi>a</mi> <mi>b</mi> </mrow> <mrow> <msub> <mi>s</mi> <mn>1</mn> </msub> <msub> <mi>s</mi> <mn>2</mn> </msub> <mo>,</mo> <mi>i</mi> <mi>j</mi> </mrow> </msubsup> </mrow>
式中：s2为s1之外的另一卫星系统，d为伪距系统间偏差，μ为所求的相位系统间偏差；
GNSS系统双差伪距观测方程和GNSS系统双差载波相位观测方程线性化后可转化为：
v＝Ax+Db+Cγ+l
式中：x为除模糊度和频间偏差外其他未知量包括测站坐标分量组成的矢量，b为接收机间单差模糊度未知数矢量，γ为系统间偏差未知量，A、D和C分别为未知量对应系数矩阵，l为常数项矢量，P为权矩阵，v为观测值残差矢量；
根据线性化方程可可单历元法方程或与之前历元累加法方程：
<mrow> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mrow> <msup> <mi>A</mi> <mi>T</mi> </msup> <mi>P</mi> <mi>A</mi> </mrow> </mtd> <mtd> <mrow> <msup> <mi>A</mi> <mi>T</mi> </msup> <mi>P</mi> <mi>D</mi> </mrow> </mtd> <mtd> <mrow> <msup> <mi>A</mi> <mi>T</mi> </msup> <mi>P</mi> <mi>C</mi> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow></mrow> </mtd> <mtd> <mrow> <msup> <mi>D</mi> <mi>T</mi> </msup> <mi>P</mi> <mi>D</mi> </mrow> </mtd> <mtd> <mrow> <msup> <mi>D</mi> <mi>T</mi> </msup> <mi>P</mi> <mi>C</mi> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mi>s</mi> <mi>y</mi> <mi>m</mi> </mrow> </mtd> <mtd> <mrow></mrow> </mtd> <mtd> <mrow> <msup> <mi>C</mi> <mi>T</mi> </msup> <mi>P</mi> <mi>C</mi> </mrow> </mtd> </mtr> </mtable> </mfenced> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mi>x</mi> </mtd> </mtr> <mtr> <mtd> <mi>b</mi> </mtd> </mtr> <mtr> <mtd> <mi>&amp;gamma;</mi> </mtd> </mtr> </mtable> </mfenced> <mo>=</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mrow> <msup> <mi>A</mi> <mi>T</mi> </msup> <mi>P</mi> <mi>l</mi> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msup> <mi>D</mi> <mi>T</mi> </msup> <mi>P</mi> <mi>l</mi> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msup> <mi>C</mi> <mi>T</mi> </msup> <mi>P</mi> <mi>l</mi> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>.</mo> </mrow>
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 true CN107728171A (zh) 2018-02-23
CN107728171B 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)

Cited By (10)

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

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102353969A (zh) * 2011-09-02 2012-02-15 东南大学 精密单点定位技术中相位偏差的估计方法
US20130120187A1 (en) * 2010-02-26 2013-05-16 NavCorn Technology, Inc. c/o Deere & Company Method and system for estimating position with bias compensation
US20160109579A1 (en) * 2014-10-16 2016-04-21 Gmv Aerospace And Defence, S.A. Device and method for computing an error bound of a kalman filter based gnss position solution
CN105738926A (zh) * 2016-03-30 2016-07-06 武汉大学 一种glonass系统接收机间相位频间偏差标定方法
CN106249256A (zh) * 2016-07-08 2016-12-21 辽宁工程技术大学 基于粒子群优化算法的实时glonass相位偏差估计方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20130120187A1 (en) * 2010-02-26 2013-05-16 NavCorn Technology, Inc. c/o Deere & Company Method and system for estimating position with bias compensation
CN102353969A (zh) * 2011-09-02 2012-02-15 东南大学 精密单点定位技术中相位偏差的估计方法
US20160109579A1 (en) * 2014-10-16 2016-04-21 Gmv Aerospace And Defence, S.A. Device and method for computing an error bound of a kalman filter based gnss position solution
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
YUMIAO TIAN ET AL.: ""Particle filter-based estimation of inter-system phase bias for real-time integer ambiguity resolution"", 《SPRINGER》 *
YUMIAO TIAN: "" Online Estimation of Inter-Frequency/System Phase Biases in Precise Positioning "", 《SPRINGER》 *

Cited By (14)

* 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 北京邮电大学 一种相位差确定方法、装置、电子设备及存储介质
CN108594277A (zh) * 2018-04-27 2018-09-28 北京邮电大学 一种相位差确定方法、装置、电子设备及存储介质
CN112014862A (zh) * 2019-05-30 2020-12-01 上海海积信息科技股份有限公司 一种载波相位观测数据生成方法及装置
CN112014862B (zh) * 2019-05-30 2024-03-29 上海海积信息科技股份有限公司 一种载波相位观测数据生成方法及装置
CN110988935A (zh) * 2019-11-25 2020-04-10 重庆市地理信息和遥感应用中心(重庆市测绘产品质量检验测试中心) 基于接收机端偏差聚类优化的多系统组合精密定位方法
CN113556304B (zh) * 2021-06-02 2022-11-04 北京大学 基于粒子滤波器的时变频偏估计方法、系统及介质
CN113556304A (zh) * 2021-06-02 2021-10-26 北京大学 基于粒子滤波器的时变频偏估计方法、系统及介质
CN113534210A (zh) * 2021-06-07 2021-10-22 湖南北斗微芯数据科技有限公司 一种基于混合卡尔曼滤波的模糊度固定方法
CN113504557A (zh) * 2021-06-22 2021-10-15 北京建筑大学 面向实时应用的gps频间钟差新预报方法
CN113589349A (zh) * 2021-09-30 2021-11-02 中国地质大学（武汉） 基于粒子滤波的gnss实时高精度单差测姿方法
CN114488227A (zh) * 2022-01-26 2022-05-13 西南交通大学 一种基于空间相关性的多路径误差改正方法
CN114488227B (zh) * 2022-01-26 2023-10-20 西南交通大学 一种基于空间相关性的多路径误差改正方法
CN117991303A (zh) * 2024-04-03 2024-05-07 武汉大学 一种天线环境变化情况下的多路径误差修正方法及装置

Also Published As

Publication number Publication date
CN107728171B (zh) 2020-10-23

Similar Documents

Publication Publication Date Title
CN107728171A (zh) 基于粒子滤波的gnss相位系统间偏差实时追踪和精密估计方法
CN107678050A (zh) 基于粒子滤波的glonass相位频间偏差实时追踪和精密估计方法
CN104714244B (zh) 一种基于抗差自适应Kalman滤波的多系统动态PPP解算方法
CN109581452B (zh) 一种gnss参考站载波相位整周模糊度解算方法
CN106168672B (zh) 一种gnss多模单频rtk周跳探测方法及装置
CN107728180B (zh) 一种基于多维粒子滤波偏差估计的gnss精密定位方法
CN107356947A (zh) 基于单频导航卫星数据确定卫星差分伪距偏差的方法
CN101403790B (zh) 单频gps接收机的精密单点定位方法
CN109061696A (zh) 一种确定导航卫星轨道和钟差的方法
CN107193029A (zh) 北斗三频信号的网络rtk基准站间模糊度快速确定方法
CN108919321B (zh) 一种基于尝试法的gnss定位粗差探测方法
CN104614741B (zh) 一种不受glonass码频间偏差影响的实时精密卫星钟差估计方法
CN111596321B (zh) 利用非差改正数的多gnss多路径误差恒星日滤波方法及系统
CN102998681A (zh) 一种卫星导航系统的高频钟差估计方法
CN102033236A (zh) 一种卫星导航位置速度联合估计方法
CN107966722A (zh) 一种gnss钟差解算方法
CN104459722B (zh) 一种基于多余观测分量的整周模糊度可靠性检验方法
CN112285745B (zh) 基于北斗三号卫星导航系统的三频模糊度固定方法及系统
CN106125113A (zh) 一种利用多系统gnss观测值的高精度基线解算方法
CN107422342A (zh) Gnss卫星钟差实时估计质量控制方法
CN109212563A (zh) 北斗/gps三频周跳探测与修复方法
CN107121689A (zh) Glonass频间偏差单历元快速估计方法
Banville et al. Defining the basis of an integer-levelling procedure for estimating slant total electron content
CN105974441A (zh) 一种接收机观测噪声的获取方法和装置
CN105954772A (zh) 一种稳健无偏的导航信号矢量跟踪方法

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