CN110059361A - 一种基于抗差卡尔曼滤波算法的实时区域对流层建模方法 - Google Patents

一种基于抗差卡尔曼滤波算法的实时区域对流层建模方法 Download PDF

Info

Publication number
CN110059361A
CN110059361A CN201910221535.4A CN201910221535A CN110059361A CN 110059361 A CN110059361 A CN 110059361A CN 201910221535 A CN201910221535 A CN 201910221535A CN 110059361 A CN110059361 A CN 110059361A
Authority
CN
China
Prior art keywords
real
formula
time
ztd
epoch
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
CN201910221535.4A
Other languages
English (en)
Other versions
CN110059361B (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.)
Institute of Geodesy and Geophysics of CAS
Original Assignee
Institute of Geodesy and Geophysics of CAS
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 Institute of Geodesy and Geophysics of CAS filed Critical Institute of Geodesy and Geophysics of CAS
Priority to CN201910221535.4A priority Critical patent/CN110059361B/zh
Publication of CN110059361A publication Critical patent/CN110059361A/zh
Application granted granted Critical
Publication of CN110059361B publication Critical patent/CN110059361B/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/14—Receivers specially adapted for specific applications
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/08—Probabilistic or stochastic CAD
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • Position Fixing By Use Of Radio Waves (AREA)

Abstract

一种基于抗差卡尔曼滤波算法的实时区域对流层建模方法，包括：利用实时精密单点定位技术，实时估计GNSS参考站网内各参考站的天顶对流层延迟参数ZTD；选定数学函数模型、随机模型及动态模型的过程噪声；基于所述参考站网内各参考站的位置信息及实时ZTD信息，采用抗差卡尔曼滤波算法实时逐历元递归估计模型参数；对求解得到的模型参数进行编码，基于NTRIP协议通过互联网实时播发。不仅利用各GNSS参考站当前历元的对流层延迟信息，还可顾及区域对流层延迟的短时平稳变化特性；此外，引入抗差机制，有效避免了因个别GNSS参考站实时对流层延迟解算异常而引起的区域对流层延迟拟合的整体偏差，提升了建模的精度及稳健性。

Description

一种基于抗差卡尔曼滤波算法的实时区域对流层建模方法
技术领域
本发明涉及GNSS区域对流层建模技术领域，尤其涉及一种基于抗差卡尔曼滤波算法的实时区域对流层建模方法，主要适用于提升区域对流层建模的精度和稳健性。
背景技术
利用区域GNSS参考站网数据实时解算天顶对流层延迟，并将其进行空间拟合，生成区域对流层延迟模型，可用于增强区域内的GNSS定位、导航、授时等应用。目前实时区域对流层建模方法一般是基于各GNSS参考站当前历元解算的对流层延迟信息，采用合适的空间拟合函数，进行单历元最小二乘估计(本发明中简称其为LSETROP)。现有方法主要不足表现为以下两个方面：①建模过程将区域对流层模型参数随时间的变化视为白噪声过程，未顾及区域对流层延迟的短时平稳变化特性，时间域信息利用不充分；②建模过程对GNSS实时观测数据、卫星轨道与钟差产品的稳定度具有强依赖性，实时数据流传输的不稳定，会严重影响实时建模的精度。
发明内容
本发明的目的是克服现有技术中存在的建模精度低、建模稳健性差的缺陷与问题，提供一种建模精度高、建模稳健性好的基于抗差卡尔曼滤波算法的实时区域对流层建模方法。
为实现以上目的，本发明的技术解决方案是：一种基于抗差卡尔曼滤波算法的实时区域对流层建模方法，该方法包括以下步骤：
A、利用实时精密单点定位技术，实时估计GNSS参考站网内各参考站的天顶对流层延迟参数ZTD；
实时精密单点定位技术解算ZTD的公式如下：
式(1)中，s和r分别代表卫星和接收机号，f代表频率号，和为经精密卫星轨道与钟差产品改正后的伪距和相位观测值，为测站至卫星方向的单位向量，i为历元标识，Δx为待估测站坐标，dtr为待估接收机钟差，为待估电离层延迟，为待估模糊度参数，εP和εΦ分别表示伪距与相位观测噪声，Tr s为站星方向对流层总延迟；
Tr s进一步参数化为天顶方向对流层干延迟ZHD和湿延迟ZWD与相应干投影函数Mh、湿投影函数Mw乘积的形式，如下：
式(2)中，ZHD由Saastamoinen模型计算得出，ZWD作为未知参数进行求解；
B、选定对流层建模的数学函数模型、随机模型及动态模型的过程噪声；
数学函数模型为：
ZTD＝a00+a01·y+a10·x+a11·xy+b1·h (3)
式(3)中，ZTD表示参考站处天顶对流层延迟参数，x、y和h表示参考站的位置信息参数，a00为区域中心位置处天顶对流层延迟，a01和a10分别为ZTD沿x和y方向变化率，a11为ZTD沿xy综合变化率，b1为ZTD高度递减率；
随机模型采用等权方案，设定建模所用实时ZTD观测量的标准差为10mm；
模型参数随时间变化均模型化为随机游走过程，其过程噪声设置如下表：
C、基于所述参考站网内各参考站的位置信息及实时ZTD信息，采用抗差卡尔曼滤波算法实时逐历元递归估计区域对流层模型参数a00、a01、a10、a11和b1；
卡尔曼滤波的状态方程和观测方程表示为：
式(4)中，Xk表示k历元时刻的状态向量，对应待估模型参数a00、a01、a10、a11和b1，其协方差阵记为Yk表示k历元时刻的观测向量，对应各GNSS参考站实时解算的ZTD信息，Ak为设计矩阵，Yk＝AkXk即为选定的空间拟合数学函数模型；Xk-1表示k-1历元时刻的状态向量，其协方差阵记为Φk，k-1表示由k-1历元时刻到k历元时刻的状态转移矩阵；Γk-1表示过程噪声，其协方差阵记为Δk为k历元时刻观测噪声，其协方差阵记为
滤波值按下列步骤进行求解：
首先，一步预测：
式(5)中，Xk，k-1表示k历元时刻的状态向量预报值，其协方差阵记为ΦT k，k-1表示Φk，k-1的转置矩阵；
然后，计算滤波值：
式(6)中，为滤波增益矩阵，其中，表示Ak的转置；I表示单位矩阵；为k历元时刻的滤波值，其协方差阵记为
借鉴拟稳估计思想，采用如下IGGIII权函数调整上述卡尔曼滤波估计中的增益矩阵Kk，实现抗差估计；
式(7)中，表示为经IGGIII权函数调整过后Kk矩阵中第p行、q列对应元素；Kpq为Kk矩阵中第p行、q列对应元素；k0、k1为经验值，取值分别为1.5和2.5；rq为抗差因子，基于卡尔曼滤波残差及其协方差确定；
最后，给定***状态向量初值X0及其协方差阵根据不断更新的区域参考站实时ZTD观测信息，利用抗差卡尔曼滤波算法递归估计得到当前观测历元的最优滤波值及其协方差阵
D、对求解得到的进行编码，基于NTRIP协议通过互联网实时播发。
步骤C中，所述抗差因子rq的确定方法如下：
Vk，k-1＝Yk-AkXk，k-1 (9)
式(8)、式(9)中，为k历元时刻的滤波残差；Vk，k-1为新息向量，其协方差阵为：
将式(10)代入滤波增益矩阵Kk计算公式，则有：
根据式(6)、式(7)和式(9)，卡尔曼滤波残差计算公式如下：
根据误差传播定律，卡尔曼滤波残差的协方差阵为：
式(13)中，为的转置；
故抗差因子rq为：
式(14)中，vq为第q个观测量的滤波残差，其相应的协方差声
与现有技术相比，本发明的有益效果为：
本发明一种基于抗差卡尔曼滤波算法的实时区域对流层建模方法采用抗差卡尔曼滤波算法实现区域对流层模型参数的实时逐历元递归估计，使得建模过程不仅利用各GNSS参考站当前历元解算的对流层延迟信息，还可充分顾及区域对流层延迟的短时平稳变化特性，充分利用了时间域信息；此外，在建模过程引入抗差机制，有效避免了因个别GNSS参考站实时对流层延迟解算异常而引起的区域对流层延迟拟合的整体偏差，提升了建模的精度及稳健性。因此，本发明提升了区域对流层建模的精度及稳健性。
附图说明
图1是一种基于抗差卡尔曼滤波算法的实时区域对流层建模方法的流程图。
图2是本发明的实施例中所用的香港SatRef跟踪站网的测站分布图。
图3是本发明的实施例中的分别利用图1所示EKFTROP建模方法与传统LSETROP建模方法实时计算的连续72小时的模型参数时间序列对比图。
图4是本发明的实施例中的每个建模历元所用的有效参考站数量序列图。
图5是本发明的实施例中的分别利用图1所示EKFTROP建模方法和传统LSETROP建模方法计算的两个测试站处连续72小时的实时ZTD序列与IGS最终ZTD产品序列比较图。
图3中，1号线表示传统LSETROP建模方法实时计算的连续72小时的模型参数时间序列，2号线表示EKFTROP建模方法实时计算的连续72小时的模型参数时间序列。
图5中，3号线表示IGS最终ZTD产品序列，4号线表示传统LSETROP建模方法计算的两个测试站处连续72小时的实时ZTD序列，5号线表示EKFTROP建模方法计算的两个测试站处连续72小时的实时ZTD序列。
具体实施方式
以下结合附图说明和具体实施方式对本发明作进一步详细的说明。
参见图1，一种基于抗差卡尔曼滤波算法的实时区域对流层建模方法，该方法包括以下步骤：
A、利用实时精密单点定位技术，实时估计GNSS参考站网内各参考站的天顶对流层延迟参数ZTD；
实时精密单点定位技术解算ZTD的公式如下：
式(1)中，s和r分别代表卫星和接收机号，f代表频率号，和为经精密卫星轨道与钟差产品改正后的伪距和相位观测值，为测站至卫星方向的单位向量，i为历元标识，Δx为待估测站坐标，dtr为待估接收机钟差，为待估电离层延迟，为待估模糊度参数，εP和εΦ分别表示伪距与相位观测噪声，Tr s为站星方向对流层总延迟；
Tr s进一步参数化为天顶方向对流层干延迟ZHD和湿延迟ZWD与相应干投影函数Mh、湿投影函数Mw乘积的形式，如下：
式(2)中，ZHD由Saastamoinen模型计算得出，ZWD作为未知参数进行求解；
B、选定数学函数模型、随机模型及动态模型的过程噪声；
数学函数模型为：
ZTD＝a00+a01·y+a10·x+a11·xy+b1·h (3)
式(3)中，ZTD表示参考站处天顶对流层延迟参数，x、y和h表示参考站的位置信息参数，a00为区域中心位置处天顶对流层延迟，a01和a10分别为ZTD沿x和y方向变化率，a11为ZTD沿xy综合变化率，b1为ZTD高度递减率；
随机模型采用等权方案，设定建模所用ZTD观测量的标准差为10mm；
模型参数随时间变化均模型化为随机游走过程，其过程噪声设置如下表：
C、基于所述参考站网内各参考站的位置信息及实时ZTD信息，采用抗差卡尔曼滤波算法实时逐历元递归估计区域对流层模型参数a00、a01、a10、a11和b1；
卡尔曼滤波的状态方程和观测方程表示为：
式(4)中，Xk表示k历元时刻的状态向量，对应待估模型参数a00、a01、a10、a11和b1，其协方差阵记为Yk表示k历元时刻的观测向量，对应各GNSS参考站实时解算的ZTD信息，Ak为设计矩阵，Yk＝AkXk即为选定的空间拟合数学函数模型；Xk-1表示k-1历元时刻的状态向量，其协方差阵记为Φk，k-1表示由k-1历元时刻到k历元时刻的状态转移矩阵；Γk-1表示过程噪声，其协方差阵记为Δk为k历元时刻观测噪声，其协方差阵记为
滤波值按下列步骤进行求解：
首先，一步预测：
式(5)中，Xk，k-1表示k历元时刻的状态向量预报值，其协方差阵记为ΦT k，k-1表示Φk，k-1的转置矩阵；
然后，计算滤波值：
式(6)中，为滤波增益矩阵，其中，表示Ak的转置；I表示单位矩阵；为k历元时刻的滤波值，其协方差阵记为
借鉴拟稳估计思想，采用如下IGGIII权函数调整上述卡尔曼滤波估计中的增益矩阵Kk，实现抗差估计；
式(7)中，表示为经IGGIII权函数调整过后Kk矩阵中第p行、q列对应元素；Kpq为Kk矩阵中第p行、q列对应元素；k0、k1为经验值，取值分别为1.5和2.5；rq为抗差因子，基于卡尔曼滤波残差及其协方差确定；
最后，给定***状态向量初值X0及其协方差阵根据不断更新的区域参考站ZTD观测信息，利用抗差卡尔曼滤波算法递归估计得到当前观测历元的最优滤波值及其协方差阵
D、对求解得到的进行编码，基于NTRIP协议通过互联网实时播发。
步骤C中，所述抗差因子rq的确定方法如下：
Vk，k-1＝Yk-AkXk，k-1 (9)
式(8)、式(9)中，为k历元时刻的滤波残差；Vk，k-1为新息向量，其协方差阵为：
将式(10)代入滤波增益矩阵Kk计算公式，则有：
根据式(6)、式(7)和式(9)，卡尔曼滤波残差计算公式如下：
根据误差传播定律，卡尔曼滤波残差的协方差阵为：
式(13)中，为的转置；
故抗差因子rq为：
式(14)中，vq为第q个观测量的滤波残差，其相应的协方差为
本发明的原理说明如下：
ZHD通过Saastamoinen模型计算得出，其输入参数为标准气象参数。
等权方案是指各观测量(本设计中具体指GNSS参考站网内各参考站解算的ZTD)在平差过程中赋予的权重相同。
拟稳估计思想指的是测量数据处理中相关平差理论的公知。
卡尔曼滤波基于观测信息序列及动力学模型信息，通过不断预测与修正过程，实现函数模型中未知参数X的递归估计。
实施例：
参见图1，一种基于抗差卡尔曼滤波算法的实时区域对流层建模方法，该方法包括以下步骤：
A、利用实时精密单点定位技术(PPP)，实时估计GNSS参考站网内各参考站的天顶对流层延迟参数ZTD；
实时精密单点定位技术解算ZTD的公式如下：
式(1)中，s和r分别代表卫星和接收机号，f代表频率号，和为经精密卫星轨道与钟差产品改正后的伪距和相位观测值，为测站至卫星方向的单位向量，i为历元标识，Δx为待估测站坐标，dt r为待估接收机钟差，为待估电离层延迟，为待估模糊度参数，εP和εΦ分别表示伪距与相位观测噪声，Tr s为站星方向对流层总延迟；
Tr s进一步参数化为天顶方向对流层干延迟ZHD和湿延迟ZWD与相应干投影函数Mh、湿投影函数Mw乘积的形式，如下：
式(2)中，由于ZHD较稳定，ZHD由Saastamoinen模型计算得出，ZWD作为未知参数进行求解；
利用式(1)和式(2)可实现ZTD参数的精确解算，ZTD参数解算过程中所用数据及主要参数处理策略如表1所示，数据时间为2018年年积日第195天至197天；
表1.ZTD参数估计过程中所用数据及主要参数处理策略
B、选定对流层建模的数学函数模型、随机模型及动态模型的过程噪声，为后续抗差卡尔曼滤波估计做准备；
数学函数模型为：
ZTD＝a00+a01·y+a10·x+a11·xy+b1·h (3)
式(3)中，ZTD表示参考站处天顶对流层延迟参数，x、y和h表示参考站的位置信息参数，a00为区域中心位置处天顶对流层延迟，a01和a10分别为ZTD沿x和y方向变化率，a11为ZTD沿xy综合变化率，b1为ZTD高度递减率；可以看出该公式中待估的模型参数共有5个，即a00、a01、a10、a11和b1，因此，应用时要求区域参考站网内有效GNSS参考站数量最少为5个；
随机模型采用等权方案，设定建模所用ZTD观测量的标准差为10mm；
模型参数随时间变化均模型化为随机游走过程，其过程噪声设置如表2：
表2.模型参数过程噪声
C、基于所述参考站网内各参考站的位置信息及实时ZTD信息，采用抗差卡尔曼滤波算法实时逐历元递归估计区域对流层模型参数a00、a01、a10、a11和b1；
与单历元最小二乘估计不同，卡尔曼滤波算法可实现区域对流层模型参数的实时逐历元递归估计，卡尔曼滤波的状态方程和观测方程表示为：
式(4)中，Xk表示k历元时刻的状态向量，对应待估模型参数a00、a01、a10、a11和b1，其协方差阵记为Yk表示k历元时刻的观测向量，对应各GNSS参考站实时解算的ZTD信息，Ak为设计矩阵，Yk＝AkXk即为选定的空间拟合数学函数模型；Xk-1表示k-1历元时刻的状态向量，其协方差阵记为Φk，k-1表示由k-1历元时刻到k历元时刻的状态转移矩阵；Γk-1表示过程噪声，其协方差阵记为Δk为k历元时刻观测噪声，其协方差阵记为
滤波值安下列步骤进行求解：
首先，一步预测：
式(5)中，Xk，k-1表示k历元时刻的状态向量预报值，其协方差阵记为ΦT k，k-1表示Φk，k-1的转置矩阵；
然后，计算滤波值：
式(6)中，为滤波增益矩阵，其中，表示Ak的转置；I表示单位矩阵；为k历元时刻的滤波值，其协方差阵记为
借鉴拟稳估计思想，采用如下IGGIII权函数调整上述卡尔曼滤波估计中的增益矩阵Kk，实现抗差估计；
式(7)中，表示为经IGGIII权函数调整过后Kk矩阵中第p行、q列对应元素；Kpq为Kk矩阵中第p行、q列对应元素；k0、k1为经验值，取值分别为1.5和2.5；rq为抗差因子，基于卡尔曼滤波残差及其协方差确定；
最后，给定***状态向量初值X0及其协方差阵根据不断更新的区域参考站ZTD观测信息，利用抗差卡尔曼滤波算法递归估计得到当前观测历元的最优滤波值及其协方差阵
所述抗差因子rq的确定方法如下：
Vk，k-1＝Yk-AkXk，k-1 (9)
式(8)、式(9)中，为k历元时刻的滤波残差；Vk，k-1为新息向量，其协方差阵为：
将式(10)代入滤波增益矩阵Kk计算公式，则有：
根据式(6)、式(7)和式(9)，卡尔曼滤波残差计算公式如下：
根据误差传播定律，卡尔曼滤波残差的协方差阵为：
式(13)中，为的转置；
故抗差因子rq为：
式(14)中，vq为第q个观测量的滤波残差，其相应的协方差为
结合式(12)和式(14)可以看出，当第q个观测量存在粗差时，则其滤波残差vq变大，对应的抗差因子rq增大，Kpq则相应减小，从而减弱了观测粗差对滤波估值的影响，实现了抗差估计；
D、对求解得到的(对应对流层模型参数)进行编码，基于NTRIP协议通过互联网实时播发。
本实施例选取中国香港SatRef GNSS观测网络部分站的实时GPS观测数据，测站分布如图2所示(图片源于http：//www.geodetic.gov.hk/en/satref/satref.htm)，测站总数为15个(HKTK、HKFN T430、HKKT、HKLT、HKCL、HKNP、HKMW、HKPC、HKSC、HKST、HKSS、HKKS、HKQT、HKOH、HKLM)，其中参与建模的测站数量为13个(HKTK、HKFN T430、HKKT、HKLT、HKCL、HKNP、HKMW、HKPC、HKSC、HKST、HKSS、HKOH、HKLM)，另外2个选作测试站(HKSL和HKWS，均为IGS站)，测区面积约为60km×30km，测站间最大高差约为345m。
图3显示了在上述测区分别利用图1所示对流层建模方法(EKFTROP)和传统LSETROP建模方法实时计算的连续72小时的模型参数时间序列对比图。从图中可明显看出，传统的LSETROP模型系数具有较大噪声，这主要是由于LSETROP建模方法仅利用当前历元的实时ZTD信息进行空间建模，而实时GNSS轨道与钟差产品估计的实时ZTD通常具有较大噪声，加之，实时数据流传输的不稳定性，更增加了LSETROP方法实时建模结果的不确定性。
图4给出了整个测试段内每个建模历元所用的有效参考站数量序列图。从图中可看出，在第8个小时附近，建模可用有效测站数由13个急剧减少为5个。由于参与建模的参考站数量及其位置分布在此历元的急剧变化，因此，造成了图3中LSETROP模型系数(b1：ZTD的高度递减)在此历元出现了不合理的急剧减小。而图1所示对流层建模方法不仅利用各参考站当前历元的ZTD信息，还充分顾及了对流层延迟的短时平稳变化特性，因此，图1所示对流层建模方法系数随时间变化较为平稳，还可有效避免个别参考站实时数据传输的不稳定对建模过程的不利影响，使得实时区域对流层建模更加稳健。
为进一步比较分析图1所示EKFTROP建模方法与传统LSETROP建模方法的应用效果，本设计选取香港地区两个未参与建模的IGS站(HKSL和HKWS)作为测试站，并以IGS提供的最终ZTD产品(IGS-ZTD)为参考，评估分析了应用上述两种方法建立的实时区域对流层模型的应用精度。图5显示了分别利用图1所示EKFTROP建模方法和传统LSETROP建模方法计算的两个测试站处连续72小时的实时ZTD序列(时间分辨率5min)与IGS-ZTD序列比较图。从图中可明显看出，利用传统LSETROP建模方法计算的ZTD序列具有较大噪声，且个别历元存在不合理的短时急剧变化，主要原因同上述的分析内容。而图1所示EKFTROP建模方法充分顾及了对流层延迟的短时平稳变化特性，在时间域对建模过程增加了合理的约束，因此，利用图1所示EKFTROP建模方法计算的ZTD序列随时间变化较为平稳，有效避免了不合理的短时急剧变化。此外，以IGS-ZTD产品为参考，统计上述两种建模方法的应用精度，结果见表3。
表3.EKFTROP与LSETROP方法在HKWS和HWSL两个测试站处的应用精度统计结果(mm)
结合图5和表3可以看出，在模型稳健性及精度方面，图1所示EKFTROP建模方法(5号线)显示出了一定的优越性，最大瞬时精度提升约达40mm。
综上，图1所示EKFTROP建模方法相对于现有方法可有效提升现有实时区域对流层建模的精度及稳定度。

Claims (2)

1.一种基于抗差卡尔曼滤波算法的实时区域对流层建模方法，其特征在于，该方法包括以下步骤：
A、利用实时精密单点定位技术，实时估计GNSS参考站网内各参考站的天顶对流层延迟参数ZTD；
实时精密单点定位技术解算ZTD的公式如下：
式(1)中，s和r分别代表卫星和接收机号，f代表频率号，和为经精密卫星轨道与钟差产品改正后的伪距和相位观测值，为测站至卫星方向的单位向量，i为历元标识，Δx为待估测站坐标，dtr为待估接收机钟差，为待估电离层延迟，为待估模糊度参数，εP和εΦ分别表示伪距与相位观测噪声，Tr s为站星方向对流层总延迟；
Tr s进一步参数化为天顶方向对流层干延迟ZHD和湿延迟ZWD与相应干投影函数Mh、湿投影函数Mw乘积的形式，如下：
式(2)中，ZHD由Saastamoinen模型计算得出，ZWD作为未知参数进行求解；
B、选定对流层建模的数学函数模型、随机模型及动态模型的过程噪声；
数学函数模型为：
ZTD＝a00+a01·y+a10·x+a11·xy+b1·h (3)
式(3)中，ZTD表示参考站处天顶对流层延迟参数，x、y和h表示参考站的位置信息参数，a00为区域中心位置处天顶对流层延迟，a01和a10分别为ZTD沿x和y方向变化率，a11为ZTD沿xy综合变化率，b1为ZTD高度递减率；
随机模型采用等权方案，设定建模所用实时ZTD观测量的标准差为10mm；
模型参数随时间变化均模型化为随机游走过程，其过程噪声设置如下表：
C、基于所述参考站网内各参考站的位置信息及实时ZTD信息，采用抗差卡尔曼滤波算法实时逐历元递归估计区域对流层模型参数a00、a01、a10、a11和b1；
卡尔曼滤波的状态方程和观测方程表示为：
式(4)中，Xk表示k历元时刻的状态向量，对应待估模型参数a00、a01、a10、a11和b1，其协方差阵记为Yk表示k历元时刻的观测向量，对应各GNSS参考站实时解算的ZTD信息，Ak为设计矩阵，Yk＝AkXk即为选定的空间拟合数学函数模型；Xk-1表示k-1历元时刻的状态向量，其协方差阵记为Φk，k-1表示由k-1历元时刻到k历元时刻的状态转移矩阵；Γk-1表示过程噪声，其协方差阵记为Δk为k历元时刻观测噪声，其协方差阵记为
滤波值按下列步骤进行求解：
首先，一步预测：
式(5)中，Xk，k-1表示k历元时刻的状态向量预报值，其协方差阵记为ΦT k，k-1表示Φk，k-1的转置矩阵；
然后，计算滤波值：
式(6)中，为滤波增益矩阵，其中，表示Ak的转置；I表示单位矩阵；为k历元时刻的滤波值，其协方差阵记为
借鉴拟稳估计思想，采用如下IGGIII权函数调整上述卡尔曼滤波估计中的增益矩阵Kk，实现抗差估计；
式(7)中，表示为经IGGIII权函数调整过后Kk矩阵中第p行、q列对应元素；Kpq为Kk矩阵中第p行、q列对应元素；k0、k1为经验值，取值分别为1.5和2.5；rq为抗差因子，基于卡尔曼滤波残差及其协方差确定；
最后，给定***状态向量初值X0及其协方差阵根据不断更新的区域参考站实时ZTD观测信息，利用抗差卡尔曼滤波算法递归估计得到当前观测历元的最优滤波值及其协方差阵
D、对求解得到的进行编码，基于NTRIP协议通过互联网实时播发。
2.根据权利要求1所述的一种基于抗差卡尔曼滤波算法的实时区域对流层建模方法，其特征在于：步骤C中，所述抗差因子rq的确定方法如下：
Vk，k-1＝Yk-AkXk，k-1 (9)
式(8)、式(9)中，为k历元时刻的滤波残差；Vk，k-1为新息向量，其协方差阵为：
将式(10)代入滤波增益矩阵Kk计算公式，则有：
根据式(6)、式(7)和式(9)，卡尔曼滤波残差计算公式如下：
根据误差传播定律，卡尔曼滤波残差的协方差阵为：
式(13)中，为的转置；
故抗差因子rq为：
式(14)中，vq为第q个观测量的滤波残差，其相应的协方差为
CN201910221535.4A 2019-03-22 2019-03-22 一种基于抗差卡尔曼滤波算法的实时区域对流层建模方法 Active CN110059361B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910221535.4A CN110059361B (zh) 2019-03-22 2019-03-22 一种基于抗差卡尔曼滤波算法的实时区域对流层建模方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910221535.4A CN110059361B (zh) 2019-03-22 2019-03-22 一种基于抗差卡尔曼滤波算法的实时区域对流层建模方法

Publications (2)

Publication Number Publication Date
CN110059361A true CN110059361A (zh) 2019-07-26
CN110059361B CN110059361B (zh) 2021-01-15

Family

ID=67316265

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910221535.4A Active CN110059361B (zh) 2019-03-22 2019-03-22 一种基于抗差卡尔曼滤波算法的实时区域对流层建模方法

Country Status (1)

Country Link
CN (1) CN110059361B (zh)

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111942618A (zh) * 2020-07-08 2020-11-17 北京控制工程研究所 一种适用于动中成像的基于gnss数据的轨道获取方法
CN114527497A (zh) * 2022-02-08 2022-05-24 国汽大有时空科技(安庆)有限公司 一种基于ppp-rtk的定位增强信息传输方法
CN114912551A (zh) * 2022-07-18 2022-08-16 中国铁路设计集团有限公司 面向桥梁变形监测的gnss和加速度计实时融合算法
CN115047505A (zh) * 2022-08-17 2022-09-13 长沙金维信息技术有限公司 基于载波相位差分辅助的gnss定位方法及导航方法
CN117992706A (zh) * 2024-04-07 2024-05-07 武汉大学 面向实时对流层天顶延迟的点面转换方法及***

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1361409A (zh) * 2000-12-23 2002-07-31 林清芳 增强型导航定位之方法及其***
CN101403790A (zh) * 2008-11-13 2009-04-08 浙江师范大学 单频gps接收机的精密单点定位方法
CN104180781A (zh) * 2014-09-10 2014-12-03 中南大学 一种单双频gps混合网变形监测数据处理方法
CN105866807A (zh) * 2016-04-05 2016-08-17 南信大影像技术工程（苏州）有限公司 一种提高gnss实时监测数据精度的算法
CN108254773A (zh) * 2017-11-24 2018-07-06 中国测绘科学研究院 一种多gnss的实时钟差解算方法
CN108920414A (zh) * 2018-05-18 2018-11-30 中国人民解放军61540部队 一种利用气象资料计算局部天顶对流层湿延迟的新方法

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1361409A (zh) * 2000-12-23 2002-07-31 林清芳 增强型导航定位之方法及其***
CN101403790A (zh) * 2008-11-13 2009-04-08 浙江师范大学 单频gps接收机的精密单点定位方法
CN104180781A (zh) * 2014-09-10 2014-12-03 中南大学 一种单双频gps混合网变形监测数据处理方法
CN105866807A (zh) * 2016-04-05 2016-08-17 南信大影像技术工程（苏州）有限公司 一种提高gnss实时监测数据精度的算法
CN108254773A (zh) * 2017-11-24 2018-07-06 中国测绘科学研究院 一种多gnss的实时钟差解算方法
CN108920414A (zh) * 2018-05-18 2018-11-30 中国人民解放军61540部队 一种利用气象资料计算局部天顶对流层湿延迟的新方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
熊永良 等: "虚拟参考站技术中对流层误差建模方法研究", 《测绘学报》 *
罗天宇 等: "虚拟参考站对流层延迟算法", 《测绘科学》 *
许长辉 等: "精密单点定位的抗差卡尔曼滤波研究", 《中国矿业大学学报》 *

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111942618A (zh) * 2020-07-08 2020-11-17 北京控制工程研究所 一种适用于动中成像的基于gnss数据的轨道获取方法
CN114527497A (zh) * 2022-02-08 2022-05-24 国汽大有时空科技(安庆)有限公司 一种基于ppp-rtk的定位增强信息传输方法
CN114912551A (zh) * 2022-07-18 2022-08-16 中国铁路设计集团有限公司 面向桥梁变形监测的gnss和加速度计实时融合算法
CN115047505A (zh) * 2022-08-17 2022-09-13 长沙金维信息技术有限公司 基于载波相位差分辅助的gnss定位方法及导航方法
CN115047505B (zh) * 2022-08-17 2022-11-22 长沙金维信息技术有限公司 基于载波相位差分辅助的gnss定位方法及导航方法
CN117992706A (zh) * 2024-04-07 2024-05-07 武汉大学 面向实时对流层天顶延迟的点面转换方法及***
CN117992706B (zh) * 2024-04-07 2024-06-11 武汉大学 面向实时对流层天顶延迟的点面转换方法及***

Also Published As

Publication number Publication date
CN110059361B (zh) 2021-01-15

Similar Documents

Publication Publication Date Title
CN110059361A (zh) 一种基于抗差卡尔曼滤波算法的实时区域对流层建模方法
CN105629263B (zh) 一种对流层大气延迟误差估计改正方法和改正***
CN104714244B (zh) 一种基于抗差自适应Kalman滤波的多***动态PPP解算方法
CN106468774B (zh) 一种应用于星基增强***的星历星钟改正参数及空间信号完好性参数方法
CN104614741B (zh) 一种不受glonass码频间偏差影响的实时精密卫星钟差估计方法
CN107678050B (zh) 基于粒子滤波的glonass相位频间偏差实时追踪和精密估计方法
CN110275186A (zh) Leo卫星增强的gnss电离层归一化与融合建模方法
CN109001771A (zh) 导航卫星和低轨卫星实时钟差确定及预报方法和***
CN107422351A (zh) 一种基于虚拟网格的gnss分米级差分定位方法
CN109683186B (zh) 一种消除多卫星导航***载波相位时间传递天跳变的方法
CN111596321B (zh) 利用非差改正数的多gnss多路径误差恒星日滤波方法及***
CN104777498B (zh) 一种基于卡尔曼滤波的gnss单点定位的方法及装置
CN106125110A (zh) 基于分区改正的星基增强***定位精度提高方法
CN108120994A (zh) 一种基于星载gnss的geo卫星实时定轨方法
CN104965207A (zh) 一种区域对流层天顶延迟的获取方法
CN109116394A (zh) 一种适用于不同长度基线的实时动态定位方法
CN108254773A (zh) 一种多gnss的实时钟差解算方法
TWI459016B (zh) 移動資訊確定裝置、方法以及接收器
CN113253314B (zh) 一种低轨卫星间时间同步方法及***
CN107728171A (zh) 基于粒子滤波的gnss相位***间偏差实时追踪和精密估计方法
CN113568020A (zh) 一种顾及硬件频间差的卫星导航定位误差修正方法和装置
Pan et al. GPS inter-frequency clock bias modeling and prediction for real-time precise point positioning
CN107121689B (zh) Glonass频间偏差单历元快速估计方法
CN107607971A (zh) 基于gnss共视时间比对算法的时间频率传递方法及接收机
CN111381264A (zh) 网络rtk中长基线模糊度固定方法和平台

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