CN113791434B - 一种gnss对流层湿延迟估计方法、系统、设备及介质 - Google Patents

一种gnss对流层湿延迟估计方法、系统、设备及介质 Download PDF

Info

Publication number
CN113791434B
CN113791434B CN202111075538.5A CN202111075538A CN113791434B CN 113791434 B CN113791434 B CN 113791434B CN 202111075538 A CN202111075538 A CN 202111075538A CN 113791434 B CN113791434 B CN 113791434B
Authority
CN
China
Prior art keywords
zwd
files
observation
random walk
delay
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
CN202111075538.5A
Other languages
English (en)
Other versions
CN113791434A (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.)
National Time Service Center of CAS
Original Assignee
National Time Service Center 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 National Time Service Center of CAS filed Critical National Time Service Center of CAS
Priority to CN202111075538.5A priority Critical patent/CN113791434B/zh
Publication of CN113791434A publication Critical patent/CN113791434A/zh
Application granted granted Critical
Publication of CN113791434B publication Critical patent/CN113791434B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO 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/00Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/01Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/13Receivers
    • G01S19/35Constructional details or hardware or software details of the signal processing chain
    • G01S19/37Hardware or software details of the signal processing chain
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO 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/00Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/38Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system
    • G01S19/39Determining 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/396Determining accuracy or reliability of position or pseudorange measurements
    • YGENERAL 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
    • Y02TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02ATECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Signal Processing (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Position Fixing By Use Of Radio Waves (AREA)

Abstract

本发明公开了一种GNSS对流层湿延迟估计方法、系统、设备及介质,所述方法包括以下步骤:获取不同频点的观测值,在每个历元形成无电离层组合观测方程;根据先验ZWD随机游走参数进行Kalman滤波状态更新和观测更新,得到未知变量预报值和其方差协方差阵Pt|t‑1;计算给定先验ZWD随机游走噪声和观测噪声时的数学期望;以所述数学期望作为伪观测值,应用滑动窗口里的数据估计新的ZWD随机游走参数;用新的ZWD随机游走参数更新先验ZWD随机游走参数。本发明的方法可以实时估计随时间变化的随机模型参数,可改善对流层湿延迟估计值。

Description

一种GNSS对流层湿延迟估计方法、系统、设备及介质
技术领域
本发明属于卫星定位技术领域,涉及对流层湿延迟估计领域,特别涉及一种GNSS对流层湿延迟估计方法、系统、设备及介质。
背景技术
对流层是全球导航卫星系统(GNSS)、甚长基线干涉测量(VLBI)和多普勒定轨与无线电定位系统(DORIS)等空间大地测量技术的主要误差来源之一,目前对流层已经成为实现毫米级高精度地面参考系,可靠地测量和解释冰后回弹,监测海平面和全球变暖等气候变化影响的制约因素。
目前,通常将对流层延迟分为干延迟(Zenith Hydrostatic Delay,ZHD)和湿延迟(Zenith Wet Delay,ZWD)两部分。
(1)干延迟可以通过经验模型准确建模,经验模型可以分为两类:第一类是以地面气象数据作为输入,根据Saastamoinen(1972)和Hopfield(1969)等经验气象公式计算天顶方向的对流层延迟(Zenith Total Delay,ZTD);第二类需要时间和近似坐标,以使用数值天气模型(Numerical Weather Model,NWM)的周年项、半周年项和温度变化率等参数的平均值计算对流层延迟,这类模型包括GPT和GPT2等。
(2)相比于干延迟ZHD,由于对流层湿延迟ZWD依赖于大气水汽含量,随着时间和空间剧烈变化,导致ZWD难以准确地预测和模型化,只能作为未知参数进行估计。一般是将GNSS信号在倾斜方向的湿延迟通过映射函数映射到天顶方向,逐历元估计对天顶湿延迟ZWD、接收机坐标和钟差;将倾斜对流层延迟映射到天顶方向的映射函数包括Neill(1996)和VMF1等。另外,有时还需要考虑大气方位角不对称性的影响而使用三维映射函数,即在估计对流层天顶延迟的同时估计对流层梯度。
随着新的GNSS数据处理技术发展,与事后处理的ZWD产品相比,GNSS轨道和时钟的(近)实时高精度产品使人们能够在大约1厘米的精度实时估算ZWD,仍需要改善(近)实时ZWD精度。目前主流的方法包括:
(1)从单GPS星座演变为多GNSS星座,以便可以利用更多的观测值和更好的星座几何构型。
(2)改善离散的经验映射函数,使得不规则变化ZWD的一部分进入经验模型,并且该映射函数可以在特定仰角下表现更好。这两种方法已经发展得相对成熟。
(3)通过实时估计ZWD随机模型改进实时ZWD估计值。在处理GNSS数据时,通常将ZWD模型化随机游走噪声(Random Walk Noise,RWN),目前不同的软件采用不同的固定RWN值,取值范围从0.01到但是,这种固定RWN的方法有时并不能区分正常天气事件和极端天气事件,导致极端天气条件下的ZWD估计值不准确,故有必要根据GNSS台站位置和当地的天气状况确定最佳的RWN。
针对上述现有方法(3),根据RWN定义可知,可从数值天气模型NWM获得RWN;然而,通过NWM计算RWN尚存在以下缺点:
(1)由于NWM不准确,且ZWD的大小和变化趋势取决于所使用的映射函数和GPS观测值随机模型,因此由NWM得出的RWN可能与软件的背景模型不一致。
(2)NWM并不能保证RWN估计的准确度。
(3)NWM的产品不是公开可用的。
(4)通过NWM计算RWN需要采用光线追踪法计算,导致计算复杂且耗时。
(5)NWM通常每6个小时发布一次,时间分辨率低;例如,时长仅2小时的暴风雨等强对流天气突然来临时,无法满足要求。
发明内容
本发明的目的在于提供一种GNSS对流层湿延迟估计方法、系统、设备及介质,以解决上述存在的一个或多个技术问题。本发明的方法以贝叶斯原理和极大似然估计法为基础,充分考虑了GNSS测站的地理位置和对流层快速变化的时变特性,可以实时估计随时间变化的随机模型参数,可改善对流层湿延迟估计值。
为达到上述目的,本发明采用以下技术方案:
本发明提供的一种GNSS对流层湿延迟估计方法,包括以下步骤:
获取不同频点的观测值,在每个历元形成无电离层组合观测方程;
根据先验ZWD随机游走参数进行Kalman滤波状态更新和观测更新,得到未知变量预报值和其方差协方差阵Pt|t-1;所述未知变量包括测站坐标、接收机钟、模糊度参数和ZWD;基于所述无电离层组合观测方程,计算得到新的参数滤波值/>和其方差协方差阵Pt|t,并计算获得事后残差/>做一次向后平滑,得到参数平滑估计值/>和其方差协方差阵Pt-1|t,以及/>和/>的协方差阵Pt,t-1|t
根据Pt|t、/>Pt-1|t、Pt,t-1|t和事后残差/>计算给定先验ZWD随机游走噪声和观测噪声时的数学期望;以所述数学期望作为伪观测值,并插入滑动窗口队列尾部,进行最大化,应用滑动窗口里的数据估计新的ZWD随机游走参数;删掉滑动窗口头部数据,用新的ZWD随机游走参数更新先验ZWD随机游走参数,完成GNSS对流层湿延迟估计。
本发明方法的进一步改进在于,所述获取不同频点的观测值,在每个历元形成无电离层组合观测方程的步骤具体包括:
从码伪距和载波观测文件提取不同频点的观测值,在每个历元形成以下无电离层组合:
其中,表示卫星位置;c表示光速;/>表示测站坐标;dti表示卫星钟差;dtr表示接收机钟差;/>表示对流层映射函数;/>表示对流层延迟修正值, 是干延迟映射函数,ZHD是天顶延迟静力学延迟,/>是湿延迟映射函数,ZWD如前所述为湿天顶延迟;cδtrel是相对论效应修正值,其中/>表示地球引力常量,/>和/>为无电离层组合的波长和模糊度;/>和/>分别是伪距和载波的噪声;δPCO和δPCV分别为天线相位中心和相位变化修正值;brcv和bi分别是接收机和卫星码硬件延迟;δWindup为相位缠绕修正;/> 修正后的测站坐标,/>是接收机标石坐标,/>是测站固体潮修正值,/>是海潮修正值,是极潮位置修正值。
本发明方法的进一步改进在于,在获取不同频点的观测值,在每个历元形成无电离层组合观测方程之前,还包括:
数据准备,包括:指定先验ZWD随机游走参数;下载码伪距和载波观测文件、码间差分DCB文件、精密轨道文件、钟差文件、行星星历文件、海潮文件、地球定向参数文件、接收机天线文件、卫星和接收机硬件延迟文件。
本发明方法的进一步改进在于,用于监测和预报极端天气事件。
本发明的一种GNSS对流层湿延迟估计系统,包括:
方程获取模块,用于获取不同频点的观测值,在每个历元形成无电离层组合观测方程;
参数获取模块,用于根据先验ZWD随机游走参数进行Kalman滤波状态更新和观测更新,得到未知变量预报值和其方差协方差阵Pt|t-1;所述未知变量包括测站坐标、接收机钟、模糊度参数和ZWD;基于所述无电离层组合观测方程,计算得到新的参数滤波值/>和其方差协方差阵Pt|t,并计算获得事后残差/>做一次向后平滑,得到参数平滑估计值和其方差协方差阵Pt-1|t,以及/>和/>的协方差阵Pt,t-1|t
参数更新模块,用于根据Pt|t、/>Pt-1|t、Pt,t-1|t和事后残差/>计算给定先验ZWD随机游走噪声和观测噪声时的数学期望;以所述数学期望作为伪观测值,并插入滑动窗口队列尾部,进行最大化,应用滑动窗口里的数据估计新的ZWD随机游走参数;删掉滑动窗口头部数据,用新的ZWD随机游走参数更新先验ZWD随机游走参数,完成GNSS对流层湿延迟估计。
本发明系统的进一步改进在于,所述方程获取模块中,获取不同频点的观测值,在每个历元形成无电离层组合观测方程的步骤具体包括:
从码伪距和载波观测文件提取不同频点的观测值,在每个历元形成以下无电离层组合:
其中,表示卫星位置;c表示光速;/>表示测站坐标;dti表示卫星钟差;dtr表示接收机钟差;/>表示对流层映射函数;/>表示对流层延迟修正值, 是干延迟映射函数,ZHD是天顶延迟静力学延迟,/>是湿延迟映射函数,ZWD如前所述为湿天顶延迟;cδtrel是相对论效应修正值,和/>为无电离层组合的波长和模糊度;和/>分别是伪距和载波的噪声;δPCO和δPCV分别为天线相位中心和相位变化修正值;brcv和bi分别是接收机和卫星码硬件延迟;δWindup为相位缠绕修正;修正后的测站坐标,/>是接收机标石坐标,/>是测站固体潮修正值,/>是海潮修正值,/>是极潮位置修正值。
本发明系统的进一步改进在于,还包括:
数据准备模块,用于指定先验ZWD随机游走参数;下载码伪距和载波观测文件、码间差分DCB文件、精密轨道文件、钟差文件、行星星历文件、海潮文件、地球定向参数文件、接收机天线文件、卫星和接收机硬件延迟文件。
本发明的一种计算机设备,包括存储器、处理器以及存储在所述存储器中并可在所述处理器上运行的计算机程序,所述处理器执行所述计算机程序时实现本发明任一项上述GNSS对流层湿延迟估计方法的步骤。
本发明的一种计算机可读存储介质,所述计算机可读存储介质存储有计算机程序,所述计算机程序被处理器执行时实现本发明任一项上述GNSS对流层湿延迟估计方法的步骤。
与现有技术相比,本发明具有以下有益效果:
当极端天气发生时,不同的随机游走噪声约束和GNSS观测权会导致GNSS对流层估计值有很大差异,当前根据经验方法设置GNSS对流层随机游走参数已经不能满足极端天气下的对流层估计需求。本发明公开的方法独立于数值天气模型,采用了适合于实时估计GNSS对流层随机游走噪声的方差分量估计算法,能够为监测极端天气和气候变化提供更加成熟可靠的GNSS对流层产品。具体的,本发明的方法以贝叶斯原理和极大似然估计法为基础,充分考虑了GNSS测站的地理位置和对流层快速变化的时变特性,可以实时估计随时间变化的随机模型参数,可改善对流层天顶延迟估计值,有利于监测和追踪极端天气事件。
进一步解释的,本发明方法独立于数值天气模型,则不需要购买实时数值天气产品;相比于依赖于数值天气模型的传统射线追踪算法,本发明容易实现,计算高效。兼顾GNSS观测噪声和多路径依赖于高度角的特点,充分考虑了GNSS测站的地理位置差异和对流层快速变化的时变特性,可以实时估计随时间变化的最优随机模型参数,能够将ZWD提高(示例性的,能够将ZWD提高约11%左右),且使得ZWD局部峰值更加显著,可以更加可靠地监测和预报恶劣的天气事件。
进一步解释的,应用传统的方差分量估计技术直接估计随机游走噪声时,难以收敛,而且只能事后处理。本发明以期望最大化算法基础,进行合理的算法假设,对GNSS载波相位方差分量进行约束,保证了算法的收敛性,并可以实时估计。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案,下面对实施例或现有技术描述中所需要使用的附图做简单的介绍;显而易见地,下面描述中的附图是本发明的一些实施例,对于本领域普通技术人员来说,在不付出创造性劳动的前提下,还可以根据这些附图获得其他的附图。
图1是本发明实施例的一种GNSS对流层湿延迟估计方法的流程示意图;
图2是本发明实施例中,随机噪声的条件数学期望滑动窗口队列示意图;
图3是本发明实施例中,用于进行算法数值测试验证的IGS测站分布图;
图4是本发明实施例中,应用本发明方法估计的各个IGS测站随机游走噪声时间序列示意图;
图5是本发明实施例中,极端天气下的对流层ZTD时间序列示意图;其中为了和JPL事后对流层产品进行比较,将ZWD转换为ZTD。
具体实施方式
为使本发明实施例的目的、技术效果及技术方案更加清楚,下面结合本发明实施例中的附图,对本发明实施例中的技术方案进行清楚、完整地描述;显然,所描述的实施例是本发明一部分实施例。基于本发明公开的实施例,本领域普通技术人员在没有做出创造性劳动的前提下所获得的其它实施例,都应属于本发明保护的范围。
请参阅图1,本发明实施例的一种改进GNSS对流层湿延迟估计方法,每一个历元估计ZWD时可以分解为以下三个步骤,其中Step1和Step2是传统ZWD估计算法;Step3是本发明实施例新加入的算法模块,是本发明实施例方法的创新核心。
Step 1.数据准备,包括指定先验随机参数、下载码伪距和载波观测文件和码间差分DCB文件、精密轨道和钟差文件、行星星历文件和海潮文件、地球定向参数文件、接收机天线文件、卫星和接收机硬件延迟文件。在准备完毕这些文件后,进入下一步提取观测值并修正观测值。
Step 2.精密单点定位(PPP)估计ZWD。下面以GPSL1和L2频点的观测值为例予以简要介绍。
首先,从码伪距和载波观测文件提取不同频点的观测值,在每个历元形成以下无电离层组合:
其中,卫星位置,可由精密轨道文件插值计算得到;c表示光速;/>是测站坐标,可由SINEX直接提取;dti为卫星钟差,可由精密卫星钟差文件计算得到;dtr是接收机钟差;是对流层映射函数,/>是对流层延迟修正值,分为干延迟和湿延迟:
是干延迟映射函数,ZHD是天顶延迟静力学延迟,/>是湿延迟映射函数,ZWD如前所述为湿天顶延迟。
cδtrel是相对论效应修正值,可由接收机、卫星位置和速度计算得到:
其中,是地球引力常量。
和/>为无电离层组合的波长和模糊度;/>和/>分别是伪距和载波的噪声。
δPCO和δPCV分别为天线相位中心和相位变化修正值,可由天线文件计算修正;brcv和bi分别是接收机和卫星码硬件延迟,可以由硬件延迟文件提取修正;δWindup为相位缠绕修正,可由卫星位置和手机位置计算得到,具体计算方法见Wu et al(1993).
在形成观测方程(1)和(2)时,需要应用精密轨道和钟差文件,对数据进行预处理,包括周跳探测、接收机毫秒跳修复、周跳探测和接收机毫秒跳修复。同时需要应用星历文件、海潮文件、地球定向参数和天线文件对GPS位置予以修正:
其中,修正后的测站坐标,/>是接收机标石坐标,/>是测站固体潮修正值,/>是海潮修正值,/>是极潮位置修正值。其它修正如计算相位中心改正、相位解缠改正、相对论效应改正、先验对流层改正等可以参考IERS2020协议,这里不再赘述。
然后,根据先验的对流层随机参数信息进行Kalman滤波状态更新和观测更新,得到新的未知变量(包括测站坐标、接收机钟、模糊度参数和ZWD)预报值和预报值的方差协方差阵Pt|t-1,下标表示由历元t-1前面的观测值(包括t-1时刻)得到历元t时刻的估计量,后面的下标含义以此类推;计算得到新的参数滤波值/>和其方差Pt|t,下标表示由历元t前面的观测值(包括t时刻)计算得到历元t时刻的估计量;并计算事后残差/>下标表示由历元t前面的观测值(包括t时刻)计算得到历元t时刻的估计量;最后做一次向后平滑(Lag-one smoothing),得到参数平滑估计值/>和其方差Pt-1|t,以及/>和/>的协方差阵Pt,t-1|t。以/>Pt|t、/>Pt-1|t、Ptt-1|t和事后残差/>作为第Step3的输入值,用来在下一步计算随机游走噪声,进入第三步。
Step 3.估计新的ZWD随机游走参数。首先,根据Pt|t、/>Pt-1|t、Pt,t-1|t和事后残差/>计算给定当前先验ZWD随机游走噪声和观测噪声时的数学期望。
然后,以此数学期望作为伪观测值,并插入滑动窗口队列尾部,紧接着进行最大化,即应用滑动窗口里的数据估计新的ZWD随机游走参数。最后删掉滑动窗口头部数据,将新的估计参数代入Step 1以更新先验ZWD随机游走参数,重复步骤Step 1至Step 3。
如前所述,Step 3是本发明算法的主要创新点,本发明方法与现有技术相比具有以下优点:
1)在现有算法中,ZWD的随机游走噪声被设置为固定值,没有考虑不同区域不同时间的ZWD的时空差异,导致实时ZWD估计值估计不准确,使得ZWD时间序列相位滞后。而本发明方法主要根据贝叶斯原理实时估计ZWD随机游走噪声和观测噪声方差分量,将ZWD提高约11%左右。
2)本发明的算法具有高时间分辨率,特别地,当极端天气发生时,本发明的方法会使得ZWD局部峰值更加显著,可以更加可靠地监测和预报恶劣的天气事件。
3)应用现有的方差分量估计技术直接估计随机游走噪声时,难以收敛,而且一般只能事后处理。由于本发明的算法不但能实时估计随机游走噪声模型,而且相比于同类的极大似然估计方法更加稳健,结果更加可靠。
为了更好地理解本申请的技术方案,下面结合一个具体实例详细说明本技术的实现步骤。该示例以GPS星座L1和L1载波,P1和P1伪随机噪声码作为基本观测值,形成无电离层组合LC和PC来模拟实时估计ZWD。
请参阅图1至图5,本发明实施例的一种改进GNSS对流层湿延迟估计方法,具体包括以下步骤:
Step 1:预先下载rinex观测文件,准备GPS精密卫星轨道和钟差文件、与钟差和轨道配套的地球旋转参数ERP、GPS码差分文件DCB、IGS天线文件、JPL行星星历文件和海潮FES2004文件。
Step 2:为详细阐述该步骤,需要首先给出ZWD时变随机游走噪声的定义,然后介绍与精密单点定位相对应的状态方程和观测方程。
1)时变随机游走噪声的定义
考虑对流层随机噪声序列{ZWDt|t=0,1,…},满足
ZWDt=ZWDt-1+wt (6)
其中,wt独立的并且与平均值的零正态分布,其协方差线性地依赖于比例因子q和ZWDt和ZWDt-1采样间隔Δt,即
E{wt}=0,Cov{wt,wt}=qt·Δt (7)
那么{ZWDt|t=0,1,…}被定义为随机游走噪声,参数qt随时间变化,是描述时变随机游走噪声特征的谱密度(Spectral Power Density,SPD)。
2)状态方程定义
假定在精密单点定位PPP中,应用卡尔曼滤波进行数据处理,其状态方程为:
其中右箭头表示向量,tandt-1是表示前后两个历元时刻,t=1,…,n;是r×1维状态向量,X,Y和Z分别是接收机的地心坐标分量;cdt表示接收机钟差,/>表示无电离层组合的模糊度参数,mt表示可见卫星数。Φt表示状态转移矩阵;/>是过程噪声向量,服从高斯分布,/>的均值为零,方差矩阵满足条件:
其中,Qk,t(k=1,…,p)为对角矩阵或者标量;αk(k=1,…,p)表示过程噪声Qt的方差分量,为待求参数。可以看出对流湿天顶延迟ZWD的随机游走参数qt满足公式(4)的要求。
3)观测方程定义
读取一个历元观测数据,进行数据预处理,形成以下无电离层组合:
其中,记号与公式(1)和(2)同。公式(8)和(9)是非线性的,可应用扩展的卡尔曼滤波进行估计。为了方便描述,将公式(8)和(9)写为一般的矩阵形式:
其中是包含接收机坐标、接收机钟、模糊度和对流层ZWD的未知向量。/>是由无电离层组合/>构成的观测向量,/>是观测噪声向量。/>Rt是t时刻的观测噪声矩阵:
其中,Rk,t(k=1,…,q)为对角矩阵或者标量;βk(k=1,…,q)表示过程噪声Rt的方差分量,为待求参数。
根据公式(6-11)可以按照扩展的卡尔曼滤波理论得到状态向量的预测值及其预测方差Pt|t-1;进一步进行观测更新可得到状态向量的估计值预测值/>和协方差方差Pt|t;最后进行一步平(Lag-one smooth)得到平滑状态向量/>和协方差方差Pt-1|t,进一步转入Step 3进行随机噪声或方差分量估计。
Step 3:
Step 3.1计算数学期望
假定满足/>t=1,…,n。
其中是观测数据,/>是未知变量,由于/>无法直接观察到,被称为隐藏变量。
其视为完备数据,并假设其似然函数/>依赖于未知参数向量Θ={αk(k=1,…,p),βj(j=1,…,q)},假定Θ(t-1)表示应用历史数据估计得到的先验未知参,则可以得到数学期望:
其中,是事后残差,/>和Pt|t为参数滤波值和其方差;/>和Pt-1|t为参数平滑估计值和其方差Pt-1|t,Pt,t-1|t是/>和/>的协方差阵,/>
Step 3.2将Xt和Yt分别插入各自滑动窗口队列尾部。
Step 3.3最大化M-step步骤。
根据滑动窗口中的队列元素,根据以下公式计算方差分量:
Step 3.3.1滑动窗口方差分量估计
公式(16)和公式(17)使用一定时间内的历史数据递推更新当前历元参数,当前历元的参数估计值会用下一迭代的先验值,滑动窗口不仅能够追踪迅速变化的时变参数也能快速过滤野值。
变量L和N和用于强调对过程噪声和观测噪声滑动窗口大小不相同,因为他们可能有不同的变化率。例如,由于用于观测值方差分量几乎恒定的,N应该取较大的值,如大于2000个历元,随机游走噪声L应该相对较小。根据实践经验,L应该超过10分钟,调整滑动窗口的大小将改变未知数的时间特性。
Step3.3.2边界约束
为了加强数值稳定性和保证算法的收敛性,有必要对待估参数设置上限和下限约束。由于GPS载波相位测量的精度约为2至3厘米,远高于随机行走噪声导致随机行走噪声很容易会吸收到GPS载波测量中,导致随机游走噪声估计值不算变小,而观测噪声方差分量不断增加。如果没有边界约束,则随机步行噪声会被低估,进而会导致对流层天顶总延迟被低估。通常,可以通过分析参数的物理特性来获得这样的参数的上下边界。然而,由于地理位置差异和天气变化复杂,导致难以具体量化模型误差。但是,观察GNSS相位残差给出一个提示以间接约束未知参数。当发生强对流天气时,GPS相位残差的15分钟RMS平均值会在11毫米至25毫米之间波动。因此,如果使用相位残差小于某个值(例如15毫米)的载波相位来更新参数以达到(间接)约束参数变化的目的。
本发明实施例的一种GNSS对流层时变随机游走噪声估计方法,涉及应用GNSS进行极端天气预测领域;极端天气下的GNSS对流层实时随机游走参数估计技术以期望最大化算法和卡尔曼滤波基础,可以独立于数值天气模型实时估计GNSS对流层随机游走参数。本发明独立于数值天气模型,仅仅应用GNSS精密单点定位技术PPP就可以实时估计对流层湿天顶延迟(Zenith Wet Delay,ZWD)随机游走参数;以递推方式对卡尔曼滤波过程噪声和观测噪声分量采用不同长度的滑动窗口长度提取随机游走噪声参数;并指定GNSS载波方差分量的边界约束条件,以保证整个算法的收敛性和可靠性。
Step 3.4删除滑动窗口队列头部元素为下一次估计做准备,如图2所示。
本发明使用精密单点定位(PPP)来计算天顶湿延迟。PPP映射函数是全局映射函数(Boehm et al,2006a),并且PPP是基于Kalman滤波器实现的,公式(16,17)用于估计随机步行噪声的参数。为确保收敛,(17)中用于观察的滑动窗口的大小应大于2000个历元,公式(16)随机行走噪音的滑动窗口长度约为10分钟。只用残差小于15毫米的载波相位观测量估计随机游走噪声参数qt
本发明实施例提供的一种测试方案,具体包括:选取13个IGS台站,被选择的台站随机且均匀地分布在地球上的不同位置具有不同天气(见图3)。其中测站SCIP上空有两次大气河经过,导致SCIP经历过两次强降雨过程。PPP观测方程是标准的GPS双频无电离层组合,用PL和LC分别表示为代码伪距和载波相位观测的无电离层组合(单位:米),它们已经校正了卫星时钟,相对论效应,固体地球潮汐,极地潮汐,海洋潮汐,相位缠绕,先验对流层延迟和相位结束。在(13)中,PC和LC观测量的随机特性被认为依赖于高度角,但是具有不同的方差分量,即:
其中,是历元t时刻卫星s的观测噪声,/>是从接收器到卫星的仰角,β1andβ2分别表示方差分量PC和LC方差分量,mt是视野中卫星数目。
对于除SCIP以外的所有GPS台站,均使用卫星轨道和时钟的最终GFZ多GNSS(GBM)产品(https://cddis.nasa.gov/archive/gnss/products/mgex/)。对测站SCIP,精密星历的轨道和卫星时钟从该中心CODE下载(http://ftp.aiub.unibe.ch/CODE/2017/)。精确的轨道校正的采样间隔为15分钟,而精确的时钟的采样间隔为30分钟,均在GFZ处生成。本发明使用CODE差分码偏差(DCB)乘积将C1对齐为P1,使用的频段是GPS的L1和L2波段。
图4显示了一个星期内随机行走参数的估计值(不包括测站SCIP)。横轴是历元数,谱密度qt平方根结果表明不同地理位置的IGS台站的估计值明显不同,大概范围是/>最小且最稳定值是HRAG站,而且它的时间序列比其他时间序列平滑得多,这是因为HRAG使用氢主时钟,并且HRAG位于沙漠地区,其上空的对流层相对稳定。与HRAG相比,其他站点显示出较大的干扰,KOUR的幅度大约是HRAG的三倍。另外,在ZIM2站,ZWD改进可以达到11%。
图5显示测站SCIP天顶对流层延迟ZTD时间序列,两个峰值表示极端天气-大气河流两次经过测站SCIP。○表示应用普通PPP算法估计的ZTD时间序列,◇表示应用本发明的算法得到的结果,黑色+曲线JPL发布的事后ZTD产品。由于大气河流2017年1月两次在北美西海岸登陆,导致SCIP经历了两次极端降雨天气。可以看出,本发明的算法估计的ZTD和JPL时间序列呈现现出相似的行为,两者在竖轴方向上平均偏移为3.9mm,RMS约为7.8mm,这两种解都明显地捕获了对流层的突然变化。从2017年1月18日至1月19日,ZTD首次急剧上升了约10cm,然后由于大气河流离去ZTD值随后立即下降。同样,由于第二条大气河道在2017年1月22日至2017年1月23日之间到达,ZTD再一次急剧增加。图中可以清晰地看出,本发明的ZTD时间序列比JPL ZTD强降水时段峰值更加明显,他们在峰值差异可以高达约2厘米,因此,应用本发明的算法在监测和预报恶劣天气事件时更加优越。此外,还可以看出,当极端天气发生时,普通的PPP的ZTD值会产生时间滞后,ZTD的也小于本发明的ZTD估计值。在三种解中,本发明实施例的方案明显更好。
综上,本发明的核心创新点在于,1)根据期望-最大化算法和Kalman滤波理论,本发明开发出了一种能够实时估计GNSS对流层随机游走噪声的算法框架。2)本发明算法独立于数值天气模型,计算高效。3)该算法采用滑动窗口和间接约束载波相位方差分量来估计GNSS对流层随机游走噪声,从而保证了算法的收敛性和数值可靠性。4)相比于传统的GNSS精密单点定位技术估计的ZTD值,本发明的算法可以将ZTD提高约11%左右。5)特别地,本发明是算法可以显著提高极端天气下的ZTD估计值,可以更加有效地利用GNSS进行极端天气的监测和预报。
本发明公开了一种应用GNSS进行极端天气预测时,能够实时估计其随机游走噪声(Random Walk Noise,RWN)的新算法框架。在当前的GNSS数据处理软件中,RWN参数通常被视作为常数按照经验给出。但是不同的RWN参数和相对应的观测权会导致极端天气下的ZWD估计值存在巨大差异,实时估计ZWD时变最优RWN参数,是准确估计极端天气下的ZWD值,监测和预报极端天气事件的关键。本发明了一种以卡尔曼滤波和期望最大化算法为基础,兼顾GNSS观测噪声和多路径依赖于高度角的特点,充分考虑了GNSS测站的地理位置差异和对流层快速变化的时变特性,可以实时估计随时间变化的最优随机模型参数,将ZWD提高约11%左右,并且使得ZWD局部峰值更加显著,可以更加可靠地监测和预报恶劣的天气事件。
本领域内的技术人员应明白,本申请的实施例可提供为方法、系统、或计算机程序产品。因此,本申请可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且,本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。
本申请是参照根据本申请实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器,使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中,使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品,该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上,使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理,从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
以上实施例仅用以说明本发明的技术方案而非对其限制,尽管参照上述实施例对本发明进行了详细的说明,所属领域的普通技术人员依然可以对本发明的具体实施方式进行修改或者等同替换,这些未脱离本发明精神和范围的任何修改或者等同替换,均在申请待批的本发明的权利要求保护范围之内。

Claims (9)

1.一种GNSS对流层湿延迟估计方法,其特征在于,包括以下步骤:
获取不同频点的观测值,在每个历元形成无电离层组合观测方程;
根据先验ZWD随机游走参数进行Kalman滤波状态更新和观测更新,得到未知变量预报值和其方差协方差阵Pt|t-1;所述未知变量包括测站坐标、接收机钟、模糊度参数和ZWD;基于所述无电离层组合观测方程,计算得到新的参数滤波值/>和其方差协方差阵Pt|t,并计算获得事后残差/>做一次向后平滑,得到参数平滑估计值/>和其方差协方差阵Pt-1|t,以及/>和/>的协方差阵Pt,t-1|t
根据Pt|t、/>Pt-1|t、Pt,t-1|t和事后残差/>计算给定先验ZWD随机游走噪声和观测噪声时的数学期望;以所述数学期望作为伪观测值,并插入滑动窗口队列尾部,进行最大化,应用滑动窗口里的数据估计新的ZWD随机游走参数;删掉滑动窗口头部数据,用新的ZWD随机游走参数更新先验ZWD随机游走参数,完成GNSS对流层湿延迟估计。
2.根据权利要求1所述的一种GNSS对流层湿延迟估计方法,其特征在于,所述获取不同频点的观测值,在每个历元形成无电离层组合观测方程的步骤具体包括:
从码伪距和载波观测文件提取不同频点的观测值,在每个历元形成以下无电离层组合:
其中,表示卫星位置;c表示光速;/>表示测站坐标;dti表示卫星钟差;dtr表示接收机钟差;/>表示对流层映射函数;/>表示对流层延迟修正值,/> 是干延迟映射函数,ZHD是天顶延迟静力学延迟,/>是湿延迟映射函数,ZWD如前所述为湿天顶延迟;cδtrel是相对论效应修正值,/>其中/>表示地球引力常量,/>和/>为无电离层组合的波长和模糊度;/>和/>分别是伪距和载波的噪声;δPCO和δPCV分别为天线相位中心和相位变化修正值;brcv和bi分别是接收机和卫星码硬件延迟;δWindup为相位缠绕修正; 修正后的测站坐标,/>是接收机标石坐标,/>是测站固体潮修正值,/>是海潮修正值,/>是极潮位置修正值。
3.根据权利要求1所述的一种GNSS对流层湿延迟估计方法,其特征在于,在获取不同频点的观测值,在每个历元形成无电离层组合观测方程之前,还包括:
数据准备,包括:指定先验ZWD随机游走参数;下载码伪距和载波观测文件、码间差分DCB文件、精密轨道文件、钟差文件、行星星历文件、海潮文件、地球定向参数文件、接收机天线文件、卫星和接收机硬件延迟文件。
4.根据权利要求1所述的一种GNSS对流层湿延迟估计方法,其特征在于,用于监测和预报极端天气事件。
5.一种GNSS对流层湿延迟估计系统,其特征在于,包括:
方程获取模块,用于获取不同频点的观测值,在每个历元形成无电离层组合观测方程;
参数获取模块,用于根据先验ZWD随机游走参数进行Kalman滤波状态更新和观测更新,得到未知变量预报值和其方差协方差阵Pt|t-1;所述未知变量包括测站坐标、接收机钟、模糊度参数和ZWD;基于所述无电离层组合观测方程,计算得到新的参数滤波值/>和其方差协方差阵Pt|t,并计算获得事后残差/>做一次向后平滑,得到参数平滑估计值/>和其方差协方差阵Pt-1|t,以及/>和/>的协方差阵Pt,t-1|t
参数更新模块,用于根据Pt|t、/>Pt-1|t、Pt,t-1|t和事后残差/>计算给定先验ZWD随机游走噪声和观测噪声时的数学期望;以所述数学期望作为伪观测值,并插入滑动窗口队列尾部,进行最大化,应用滑动窗口里的数据估计新的ZWD随机游走参数;删掉滑动窗口头部数据,用新的ZWD随机游走参数更新先验ZWD随机游走参数,完成GNSS对流层湿延迟估计。
6.根据权利要求5所述的一种GNSS对流层湿延迟估计系统,其特征在于,所述方程获取模块中,获取不同频点的观测值,在每个历元形成无电离层组合观测方程的步骤具体包括:
从码伪距和载波观测文件提取不同频点的观测值,在每个历元形成以下无电离层组合:
其中,表示卫星位置;c表示光速;/>表示测站坐标;dti表示卫星钟差;dtr表示接收机钟差;/>表示对流层映射函数;/>表示对流层延迟修正值,/> 是干延迟映射函数,ZHD是天顶延迟静力学延迟,/>是湿延迟映射函数,ZWD如前所述为湿天顶延迟;cδtrel是相对论效应修正值,/> 表示地球引力常量,/>和/>为无电离层组合的波长和模糊度;/>和/>分别是伪距和载波的噪声;δPCO和δPCV分别为天线相位中心和相位变化修正值;brcv和bi分别是接收机和卫星码硬件延迟;δWindup为相位缠绕修正;/> 修正后的测站坐标,/>是接收机标石坐标,/>是测站固体潮修正值,/>是海潮修正值,/>是极潮位置修正值。
7.根据权利要求5所述的一种GNSS对流层湿延迟估计系统,其特征在于,还包括:
数据准备模块,用于指定先验ZWD随机游走参数;下载码伪距和载波观测文件、码间差分DCB文件、精密轨道文件、钟差文件、行星星历文件、海潮文件、地球定向参数文件、接收机天线文件、卫星和接收机硬件延迟文件。
8.一种计算机设备,包括存储器、处理器以及存储在所述存储器中并可在所述处理器上运行的计算机程序,其特征在于,所述处理器执行所述计算机程序时实现如权利要求1至3任一项所述GNSS对流层湿延迟估计方法的步骤。
9.一种计算机可读存储介质,所述计算机可读存储介质存储有计算机程序,其特征在于,所述计算机程序被处理器执行时实现如权利要求1至3任一项所述GNSS对流层湿延迟估计方法的步骤。
CN202111075538.5A 2021-09-14 2021-09-14 一种gnss对流层湿延迟估计方法、系统、设备及介质 Active CN113791434B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111075538.5A CN113791434B (zh) 2021-09-14 2021-09-14 一种gnss对流层湿延迟估计方法、系统、设备及介质

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111075538.5A CN113791434B (zh) 2021-09-14 2021-09-14 一种gnss对流层湿延迟估计方法、系统、设备及介质

Publications (2)

Publication Number Publication Date
CN113791434A CN113791434A (zh) 2021-12-14
CN113791434B true CN113791434B (zh) 2023-08-01

Family

ID=78880234

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111075538.5A Active CN113791434B (zh) 2021-09-14 2021-09-14 一种gnss对流层湿延迟估计方法、系统、设备及介质

Country Status (1)

Country Link
CN (1) CN113791434B (zh)

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102331582A (zh) * 2010-05-30 2012-01-25 天宝导航有限公司 利用联合离子层滤波器的gnss大气估计
WO2012128980A1 (en) * 2011-03-22 2012-09-27 Trimble Navigation Limited Gnss signal processing with known position for reconvergence
CN105182388A (zh) * 2015-10-10 2015-12-23 安徽理工大学 一种快速收敛的精密单点定位方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102331582A (zh) * 2010-05-30 2012-01-25 天宝导航有限公司 利用联合离子层滤波器的gnss大气估计
CN102331583A (zh) * 2010-05-30 2012-01-25 天宝导航有限公司 利用模糊度固定的gnss大气估计
WO2012128980A1 (en) * 2011-03-22 2012-09-27 Trimble Navigation Limited Gnss signal processing with known position for reconvergence
CN105182388A (zh) * 2015-10-10 2015-12-23 安徽理工大学 一种快速收敛的精密单点定位方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
Dynamic Modelling of Zenith wet Delay in GNSS Measurements;EI-Mowafy et al.;Proceedings of the 2011 International Technical meeting of the institute of navigation;全文 *
估计对流层延迟的单频RTK卡尔曼滤波算法;徐彦田;程鹏飞;蔡艳辉;甄杰;徐宗秋;;测绘通报(08);全文 *
对流层湿延迟估计方法对PPP数据处理的影响;李盼 等;武汉大学学报(信息科学版);第35卷(第07期);全文 *

Also Published As

Publication number Publication date
CN113791434A (zh) 2021-12-14

Similar Documents

Publication Publication Date Title
Jakowski et al. Total electron content models and their use in ionosphere monitoring
US10078140B2 (en) Navigation satellite system positioning involving the generation of advanced correction information
Herman et al. Mountain glacier velocity variation during a retreat/advance cycle quantified using sub-pixel analysis of ASTER images
Rohm et al. Ground-based GNSS ZTD/IWV estimation system for numerical weather prediction in challenging weather conditions
CN103728643B (zh) 附有宽巷约束的北斗三频网络rtk模糊度单历元固定方法
CN107942346B (zh) 一种高精度gnss电离层tec观测值提取方法
Paziewski Study on desirable ionospheric corrections accuracy for network-RTK positioning and its impact on time-to-fix and probability of successful single-epoch ambiguity resolution
CN103760572A (zh) 一种基于区域cors的单频ppp电离层加权方法
CN103529461A (zh) 一种基于强跟踪滤波和埃尔米特插值法的接收机快速定位方法
Choi et al. The influence of grounding on GPS receiver differential code biases
CN115544785B (zh) 一种无资料梯级水库流域水文模拟方法和系统
CN107607971A (zh) 基于gnss共视时间比对算法的时间频率传递方法及接收机
CN108181633A (zh) 一种gnss时间频率传递接收机及接收方法
Andersen et al. The DTU21 global mean sea surface and first evaluation
Tolman et al. Absolute precise kinematic positioning with GPS and GLONASS
Saracoglu et al. Effect of meteorological seasons on the accuracy of GPS positioning
Chen et al. Near real-time global ionospheric modeling based on an adaptive Kalman filter state error covariance matrix determination method
CN117272812A (zh) 一种低纬小区域电离层模型构建方法
CN113791434B (zh) 一种gnss对流层湿延迟估计方法、系统、设备及介质
Wan et al. Toward terrain effects on GNSS interferometric reflectometry snow depth retrievals: geometries, modeling, and applications
Wagner et al. Assimilation of GNSS and synoptic data in a convection permitting limited area model: improvement of simulated tropospheric water vapor content
Chen et al. Undifferenced zenith tropospheric modeling and its application in fast ambiguity recovery for long-range network RTK reference stations
CN115435674A (zh) 一种现场观测和卫星遥感联合反演北极海冰积雪深度方法
Wang et al. Long-term time-varying characteristics of UPD products generated by a global and regional network and their interoperable application in PPP
Schaer GNSS ionosphere analysis at CODE

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