CN115755103A - 一种抗差自适应的gnss水汽层析方法 - Google Patents

一种抗差自适应的gnss水汽层析方法 Download PDF

Info

Publication number
CN115755103A
CN115755103A CN202211430336.2A CN202211430336A CN115755103A CN 115755103 A CN115755103 A CN 115755103A CN 202211430336 A CN202211430336 A CN 202211430336A CN 115755103 A CN115755103 A CN 115755103A
Authority
CN
China
Prior art keywords
gnss
water vapor
equation
chromatography
matrix
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
CN202211430336.2A
Other languages
English (en)
Other versions
CN115755103B (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.)
China University of Mining and Technology Beijing CUMTB
Original Assignee
China University of Mining and Technology Beijing CUMTB
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 China University of Mining and Technology Beijing CUMTB filed Critical China University of Mining and Technology Beijing CUMTB
Priority to CN202211430336.2A priority Critical patent/CN115755103B/zh
Publication of CN115755103A publication Critical patent/CN115755103A/zh
Application granted granted Critical
Publication of CN115755103B publication Critical patent/CN115755103B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • 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

  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种抗差自适应的GNSS水汽层析方法，其方法包括：S1、根据研究区构建GNSS水汽层析格网；S2、采集研究区内GNSS测站的观测数据、气象数据及精密星历数据，通过计算得到斜路径水汽含量SWV；S3、通过计算卫星射线在每个层析网格内的截距构建GNSS水汽层析斜路径方程，利用大气折射率参数随着高度的增加而指数递减的特性构建垂直约束方程，将GNSS水汽层析斜路径方程及垂直约束方程共同组成Kalman滤波的矩阵形式的观测方程，然后构建状态方程；S4、Kalman滤波解算模型进行依序解算更新得到大气水汽密度值。本发明能够改善GNSS水汽层析解算过程中观测值误差的影响及滤波解算的不稳定性，为GNSS气象学研究提供高精度稳定的大气水汽信息。

Description

一种抗差自适应的GNSS水汽层析方法
技术领域
本发明涉及GNSS数据及气象数据处理领域，尤其涉及一种抗差自适应的GNSS水汽层析方法。
背景技术
大气水汽的高精度、多维度、连续稳定监测对于全球温室效应、气候变化、改善数值天气预报初始场及卫星对地观测精度等研究都具有重要意义。随着GNSS测站卫星技术的发展成熟，GNSS数据与气象数据均具有全天候、连续运行、时空分辨率高等优点，如何有效利用两者数据反演得到水汽数据(包括水汽含量、水汽密度)，是GNSS气象学研究热点，也是现有技术需要解决的技术难题，目前现有技术均难以获取稳定的水汽结果，也难以进行基于水汽数据的科学研究，这严重制约水汽数据的获取、处理及研究。
发明内容
本发明的目的在于克服背景技术所指出的技术问题，提供一种抗差自适应的GNSS水汽层析方法，对研究区进行水平经纬度与垂直高程划分得到三维GNSS水汽层析格网，同时采集研究区内GNSS测站的观测数据、气象数据及精密星历数据计算得到斜路径水汽含量SWV，通过对GNSS斜路径水汽含量及卫星射线网格截距的获取，构建水汽层析观测方程，同时建立基于IGG3的抗差估计、基于预测残差的自适应因子，并开展GNSS水汽层析方程的卡尔曼滤波解算，能够改善GNSS水汽层析解算过程中观测值误差的影响及滤波解算的不稳定性，为GNSS气象学研究提供高精度稳定的大气水汽信息，能够更好地基于地基GNSS反演大气水汽信息。
本发明的目的通过下述技术方案实现：
一种抗差自适应的GNSS水汽层析方法，其方法包括：
S1、确定研究区，根据研究区构建GNSS水汽层析格网，GNSS水汽层析格网由若干个水平格网呈垂直方向分布，水平格网包括若干个层析网格；
S2、采集研究区内GNSS测站的观测数据、气象数据及精密星历数据，通过如下公式计算得到斜路径水汽含量SWV：
Figure BDA0003944190620000021
其中，Π为转换系数，ρw表示液态水密度，R为通用气体常数，mw是湿空气摩尔质量，k′2和k3是大气折射率常数，Tm为大气加权平均温度，SWD为GNSS斜路径湿延迟；
S3、以斜路径水汽含量SWV为观测量、各个层析网格内的水汽密度为待求参数，通过计算卫星射线在每个层析网格内的截距表示如下GNSS水汽层析斜路径方程：
Figure BDA0003944190620000022
SWVq表示第q条卫星射线的斜路径水汽含量截距，n表示GNSS水汽层析网格的总数，
Figure BDA0003944190620000023
表示第q条卫星射线穿过网格i的截距，xi是第i个层析网格内的水汽密度值；
利用大气折射率参数随着高度的增加而指数递减的特性构建如下垂直约束方程：
Figure BDA0003944190620000024
其中xi和xi+l分别表示第i个和第i+l个层析网格的大气水汽密度；hi和hi+l分别代表第i和第i+l个层析网格的高程；
将GNSS水汽层析斜路径方程及垂直约束方程共同组成Kalman滤波的矩阵形式的观测方程，观测方程公式如下：
L＝A·X+Δ，其中L为层析斜路径方程、垂直约束方程构成的观测值列向量，A为系数矩阵，X为所有层析网格的待求水汽密度值，Δ观测噪声矩阵；
构建如下GNSS水汽层析Kalman滤波的状态方程：
Xk＝Xk-1+wk，其中Xk表示所有水汽层析网格k时刻的大气水汽密度状态向量，wk为状态噪声；
S4、构建包含状态方程、观测方程的Kalman滤波解算模型；Kalman滤波解算模型分别计算抗差因子
Figure BDA0003944190620000031
自适应因子αk，基于k-1时刻状态初值Xk-1及其协方差矩阵
Figure BDA0003944190620000032
解算更新k时刻的状态矩阵
Figure BDA0003944190620000033
和协方差矩阵
Figure BDA0003944190620000034
其方法如下：
S41、通过如下公式计算k时刻的预测残差vk和标准化残差
Figure BDA0003944190620000035
Figure BDA0003944190620000036
Figure BDA0003944190620000037
其中Lk为k时刻层析斜路径方程、垂直约束方程构成的观测值列向量，
Figure BDA0003944190620000038
为利用k-1时刻的初值基于状态方程得到k时刻的中间预测值，σ为单位权重误差；Ak为k时刻的系数矩阵；
S42、按照如下方法得到观测方程的观测值的抗差权因子：
Figure BDA0003944190620000039
其中，
Figure BDA00039441906200000310
表示k时刻的观测方程中的观测值的抗差权因子，抗差权因子
Figure BDA00039441906200000311
与协方差矩阵
Figure BDA00039441906200000312
互为逆矩阵；pk表示k时刻的观测值原始权因子，c0和c1为常数；
S43、基于预测残差构造k时刻判别自适应因子的统计量，方法如下：
Figure BDA00039441906200000313
其中，tr表示求矩阵的迹，
Figure BDA00039441906200000314
表示预测残差vk的矩阵转置，
Figure BDA0003944190620000041
表示系数矩阵Ak的矩阵转置，
Figure BDA0003944190620000042
表示抗差权因子
Figure BDA0003944190620000043
的逆矩阵；Ak表示k时刻的观测方程系数矩阵；
Figure BDA0003944190620000044
表示k时刻的协方差矩阵；
然后构造出如下自适应因子函数模型：
Figure BDA0003944190620000045
由此得到k时刻的自适应因子αk；
S44、通过如下公式计算得到Kalman滤波增益矩阵
Figure BDA0003944190620000046
Figure BDA0003944190620000047
S45、按照如下方法进行层析k时刻的状态矩阵
Figure BDA0003944190620000048
和协方差矩阵
Figure BDA0003944190620000049
Figure BDA00039441906200000410
Figure BDA00039441906200000411
其中I表示单位阵；
以k时刻的状态矩阵
Figure BDA00039441906200000412
和协方差矩阵
Figure BDA00039441906200000413
作为初值按照方法S41～S45解算更新k+1时刻的状态矩阵
Figure BDA00039441906200000414
和协方差矩阵
Figure BDA00039441906200000415
以此类推，直至最后一个时刻完成上述循环计算，得到研究区每个层析网格的大气水汽密度值。
优选地，本发明在Kalman滤波解算模型中，状态方程中的状态噪声的协方差矩阵表示为Q，观测方程的观测噪声的协方差矩阵表示为R。
优选地，本发明GNSS斜路径湿延迟SWD通过如下公式得到：
Figure BDA00039441906200000416
其中ele表示GNSS卫星高度角；mw(ele)和mΔ(ele)分别为湿延迟和湿梯度的映射函数；
Figure BDA0003944190620000051
表示南北方向的湿大气水平梯度，
Figure BDA0003944190620000052
表示东西方向的湿大气水平梯度；Re表示GNSS观测值的非差验后残差，ZWD为测站对流层湿延迟。
进一步地，测站对流层湿延迟ZWD通过如下方法得到：
S21、采集研究区内GNSS测站的气象数据，通过如下公式计算得到天顶静力学延迟ZHD：
Figure BDA0003944190620000053
其中P为GNSS测站的地面气压，
Figure BDA0003944190620000054
为GNSS测站的地理纬度，h为GNSS测站的高程；
S22、采集研究区内GNSS测站的观测数据及精密星历数据，利用Gamit软件进行双差解算，得到GNSS天顶对流层延迟ZTD；测站对流层湿延迟ZWD通过如下公式得到：
ZWD＝ZTD-ZHD。
优选地，本发明水平格网的层析网格满足如下关系式：
s≤h1/tan(ele)，其中s表示水平格网任一边的边长，ele表示GNSS卫星高度角，hl表示研究区的层析高度。
优选地，本发明GNSS斜路径湿延迟SWD通过联合GAMIT软件利用测站对流层湿延迟ZWD处理得到。
本发明较现有技术相比，具有以下优点及有益效果：
(1)本发明对研究区进行水平经纬度与垂直高程划分得到三维GNSS水汽层析格网，同时采集研究区内GNSS测站的观测数据、气象数据及精密星历数据计算得到斜路径水汽含量SWV，通过对GNSS斜路径水汽含量及卫星射线网格截距的获取，构建水汽层析观测方程，同时建立基于IGG3的抗差估计、基于预测残差的自适应因子，并开展GNSS水汽层析方程的卡尔曼滤波解算，能够改善GNSS水汽层析解算过程中观测值误差的影响及滤波解算的不稳定性，为GNSS气象学研究提供高精度稳定的大气水汽信息，能够更好地基于地基GNSS反演大气水汽信息。
(2)本发明在GNSS水汽层析的卡尔曼滤波解算中，利用IGGIII计算抗差因子，来减弱GNSS水汽层析模型中观测值异常对水汽层析解算结果的影响，基于预测残差利用两段法构造了自适应因子，来调节GNSS水汽层析中预测量与观测量的权重；本发明能够提升水汽层析解算的精度和稳定性。
附图说明
图1为本发明的方法流程示意图；
图2为实施例中抗差自适应GNSS水汽层析方法的方法流程原理图；
图3为实施例中举例的研究区范围的水平网格划分示意图；
图4为图3举例研究区的GNSS水汽层析大气水汽密度空间分布图。
具体实施方式
下面结合实施例对本发明作进一步地详细说明：
实施例
如图1～图4所示，一种抗差自适应的GNSS水汽层析方法，其方法包括：
S1、确定研究区，根据研究区构建GNSS水汽层析格网，GNSS水汽层析格网由若干个水平格网呈垂直方向分布(也即垂直方向分辨率，垂直方向的范围选取根据实际需求确定，一般垂直方向研究区域探空历史数据进行范围确定，GNSS水汽层析格网按照所确定垂直方向的范围来确定垂直方向的划分，一般水平格网最小厚度不低于300m)，水平格网包括若干个层析网格(水平格网中的层析网格密度为水平分辨率)。优选地，水平格网的层析网格满足如下关系式：
s≤h1/tan(ele)其中s表示水平格网任一边的边长，ele表示GNSS卫星高度角，h1表示研究区的层析高度。
如图1所示，GNSS水汽层析格网为三维层析格网，包括垂直划分和水平划分，对研究区进行垂直划分得到呈垂直方向分布的若干个水平格网，对研究区的水平格网进行水平划分，水平格网上就会划分出若干个层析网格。如图3所示，以选取香港地区为研究区为例，研究区的范围是东西方向113.87°E～114.35°E，南北方向22.19N°～22.54°N，分布着14个GNSS观测站和，根据研究区测站分布特点及大小，东西方向划分8列的经度网格线(网格分辨率为0.06°)；南北方向划分为7行的纬度网格线(网格分辨率为0.05°)。根据探空历史资料，研究区高程方向的范围选定0～8千米，每层水平格网的高度为800米，参见图3，研究区被划分为8×7×10个层析网格。
S2、采集研究区内GNSS测站的观测数据、气象数据及精密星历数据，通过如下公式计算得到斜路径水汽含量SWV：
Figure BDA0003944190620000071
其中，Π为转换系数，ρw表示液态水密度，R为通用气体常数，mw是湿空气摩尔质量，k′2和k3是大气折射率常数，Tm为大气加权平均温度，SWD为GNSS斜路径湿延迟。
在一些实施例中，本发明GNSS斜路径湿延迟SWD通过如下公式得到：
Figure BDA0003944190620000072
其中ele表示GNSS卫星高度角；mw(ele)和mΔ(ele)分别为湿延迟和湿梯度的映射函数；
Figure BDA0003944190620000073
表示南北方向的湿大气水平梯度，
Figure BDA0003944190620000074
表示东西方向的湿大气水平梯度；Re表示GNSS观测值的非差验后残差，ZWD为测站对流层湿延迟。本发明GNSS斜路径湿延迟SWD通过联合GAMIT软件利用测站对流层湿延迟ZWD处理得到。
在一些实施例中，测站对流层湿延迟ZWD通过如下方法得到：
S21、采集研究区内GNSS测站的气象数据，通过如下公式计算得到天顶静力学延迟ZHD：
Figure BDA0003944190620000075
其中P为GNSS测站的地面气压，
Figure BDA0003944190620000076
为GNSS测站的地理纬度，h为GNSS测站的高程；
S22、采集研究区内GNSS测站的观测数据及精密星历数据，利用Gamit软件进行双差解算，得到GNSS天顶对流层延迟ZTD；测站对流层湿延迟ZWD通过如下公式得到：
ZWD＝ZTD-ZHD。
以图3中香港地区2021年5月1日的GNSS观测数据、气象数据及精密星历数据，首先利用Gamit软件计算得到14个测站的对流层参数(包括天顶对流层延迟ZTD和水平梯度信息)，然后根据测站气象文件记录的地表气压数据计算ZHD，表达式如下：
Figure BDA0003944190620000081
将从ZTD扣除ZHD得到的ZWD，然后与水平梯度一起来计算得到GNSS斜路径湿延迟SWD，表达式如下：
Figure BDA0003944190620000082
再将GNSS斜路径湿延迟SWD通过转化因子的转化得到GNSS斜路径水汽含量SWV：
Figure BDA0003944190620000083
S3、以斜路径水汽含量SWV为观测量、各个层析网格内的水汽密度为待求参数，通过计算卫星射线在每个层析网格内的截距表示如下GNSS水汽层析斜路径方程：
Figure BDA0003944190620000084
SWVq表示第q条卫星射线的斜路径水汽含量截距(上标q代表第q条用于GNSS水汽层析的卫星射线)，n表示GNSS水汽层析网格的总数，
Figure BDA0003944190620000085
表示第q条卫星射线穿过网格i的截距，xi是第i个层析网格内的水汽密度值。以图3的香港地区为例，在2021年5月1日，每个层析时段的斜路径方程平均数为1890个。
利用大气折射率参数随着高度的增加而指数递减的特性构建如下垂直约束方程(采取指数相关的方法构建)：
Figure BDA0003944190620000091
其中xi和xi+l分别表示第i个和第i+l个层析网格的大气水汽密度；hi和hi+l分别代表第i和第i+l个层析网格的高程。以图3的香港地区为例，选取的最大高程h为8km，每个层析时段560个层析网格建立了504个垂直约束方程。
将GNSS水汽层析斜路径方程及垂直约束方程共同组成Kalman滤波的矩阵形式的观测方程，观测方程公式如下：
L＝A·X+Δ，其中L为层析斜路径方程、垂直约束方程构成的观测值列向量，A为系数矩阵，X为所有层析网格的待求水汽密度值，Δ观测噪声矩阵；
构建如下GNSS水汽层析Kalman滤波的状态方程(基于一定层析时间内所有层析网格内水汽密度参数符合高斯-马尔科夫随机游走平稳过程，构建GNSS水汽层析Kalman滤波的状态方程)：
Xk＝Xk-1+wk，其中Xk表示所有水汽层析网格k时刻的大气水汽密度状态向量，wk为状态噪声；
S4、构建包含状态方程、观测方程的Kalman滤波解算模型；在Kalman滤波解算模型中，状态方程中的状态噪声的协方差矩阵表示为Q，观测方程的观测噪声的协方差矩阵表示为R。Kalman滤波解算模型分别计算抗差因子
Figure BDA0003944190620000092
自适应因子αk，基于k-1时刻状态初值Xk-1及其协方差矩阵
Figure BDA0003944190620000093
解算更新k时刻的状态矩阵
Figure BDA0003944190620000094
和协方差矩阵
Figure BDA0003944190620000095
其方法如下：
S41、通过如下公式计算k时刻的预测残差vk和标准化残差
Figure BDA0003944190620000096
Figure BDA0003944190620000097
Figure BDA0003944190620000101
其中Lk为k时刻层析斜路径方程、垂直约束方程构成的观测值列向量，
Figure BDA0003944190620000102
为利用k-1时刻的初值基于状态方程得到k时刻的中间预测值，σ为单位权重误差；Ak为k时刻的系数矩阵；
S42、按照如下方法得到观测方程的观测值的抗差权因子：
Figure BDA0003944190620000103
其中，
Figure BDA0003944190620000104
表示k时刻的观测方程中的观测值的抗差权因子(抗差权矩阵
Figure BDA0003944190620000105
的对角线元素)，抗差权因子
Figure BDA0003944190620000106
与协方差矩阵
Figure BDA0003944190620000107
互为逆矩阵；pk表示k时刻的观测值原始权因子，c0和c1为常数；
S43、基于预测残差构造k时刻判别自适应因子的统计量，方法如下：
Figure BDA0003944190620000108
其中，tr表示求矩阵的迹，
Figure BDA0003944190620000109
表示预测残差vk的矩阵转置，
Figure BDA00039441906200001010
表示系数矩阵Ak的矩阵转置，
Figure BDA00039441906200001011
表示抗差权因子
Figure BDA00039441906200001012
的逆矩阵；Ak表示k时刻的观测方程系数矩阵；
Figure BDA00039441906200001013
表示k时刻的协方差矩阵；上述公式中，上标T表示相应矩阵的转置，上标-1表示相应矩阵的逆矩阵。
然后构造出如下自适应因子函数模型：
Figure BDA00039441906200001014
由此得到k时刻的自适应因子αk；
S44、通过如下公式计算得到Kalman滤波增益矩阵
Figure BDA0003944190620000111
Figure BDA0003944190620000112
S45、按照如下方法进行层析k时刻的状态矩阵
Figure BDA0003944190620000113
和协方差矩阵
Figure BDA0003944190620000114
Figure BDA0003944190620000115
Figure BDA0003944190620000116
其中I表示单位阵；
以k时刻的状态矩阵
Figure BDA0003944190620000117
和协方差矩阵
Figure BDA0003944190620000118
作为初值按照方法S41～S45解算更新k+1时刻的状态矩阵
Figure BDA0003944190620000119
和协方差矩阵
Figure BDA00039441906200001110
以此类推，直至最后一个时刻完成上述循环计算(若最后一个时刻为k+p，则就会解算更新k+p时刻的状态矩阵
Figure BDA00039441906200001111
和协方差矩阵
Figure BDA00039441906200001112
)，得到研究区每个层析网格的大气水汽密度值(以k-1为起始时刻，解算更新得到k+p时刻的状态矩阵
Figure BDA00039441906200001113
研究区的时间范围为k-1～k+p，状态矩阵
Figure BDA00039441906200001114
就为结算更新得到k+p时刻每个层析网格的大气水汽密度值)。以图3的香港地区为例，选取k-1～k+p时间为2021年5月1日11点～12点，则按照步骤S1～S4的计算，就能得到研究区2021年5月1日12点所有层析网格的大气水汽密度值，其根据大气水汽密度值制作成大气水汽密度空间分布图(如图4所示)，参见图4，大气水汽密度按照研究区(香港地区)的经纬度、垂直方向按照层析网格进行叠加表示其分布情况，颜色越深，其大气水汽密度越高，大气水汽含量也越大。
以上所述仅为本发明的较佳实施例而L，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (6)

1.一种抗差自适应的GNSS水汽层析方法，其特征在于：其方法包括：
S1、确定研究区，根据研究区构建GNSS水汽层析格网，GNSS水汽层析格网由若干个水平格网呈垂直方向分布，水平格网包括若干个层析网格；
S2、采集研究区内GNSS测站的观测数据、气象数据及精密星历数据，通过如下公式计算得到斜路径水汽含量SWV：
Figure FDA0003944190610000011
其中，Π为转换系数，ρw表示液态水密度，R为通用气体常数，mw是湿空气摩尔质量，k′2和k3是大气折射率常数，Tm为大气加权平均温度，SWD为GNSS斜路径湿延迟；
S3、以斜路径水汽含量SWV为观测量、各个层析网格内的水汽密度为待求参数，通过计算卫星射线在每个层析网格内的截距表示如下GNSS水汽层析斜路径方程：
Figure FDA0003944190610000012
SWVq表示第q条卫星射线的斜路径水汽含量截距，n表示GNSS水汽层析网格的总数，
Figure FDA0003944190610000013
表示第q条卫星射线穿过网格i的截距，xi是第i个层析网格内的水汽密度值；
利用大气折射率参数随着高度的增加而指数递减的特性构建如下垂直约束方程：
Figure FDA0003944190610000014
其中xi和xi+l分别表示第i个和第i+l个层析网格的大气水汽密度；hi和hi+l分别代表第i和第i+l个层析网格的高程；
将GNSS水汽层析斜路径方程及垂直约束方程共同组成Kalman滤波的矩阵形式的观测方程，观测方程公式如下：
L＝A·X+Δ，其中L为层析斜路径方程、垂直约束方程构成的观测值列向量，A为系数矩阵，X为所有层析网格的待求水汽密度值，Δ观测噪声矩阵；
构建如下GNSS水汽层析Kalman滤波的状态方程：
Xk＝Xk-1+wk，其中Xk表示所有水汽层析网格k时刻的大气水汽密度状态向量，wk为状态噪声；
S4、构建包含状态方程、观测方程的Kalman滤波解算模型；Kalman滤波解算模型分别计算抗差因子
Figure FDA0003944190610000021
自适应因子αk，基于k-1时刻状态初值Xk-1及其协方差矩阵
Figure FDA0003944190610000022
解算更新k时刻的状态矩阵
Figure FDA0003944190610000023
和协方差矩阵
Figure FDA0003944190610000024
其方法如下：
S41、通过如下公式计算k时刻的预测残差vk和标准化残差
Figure FDA0003944190610000025
Figure FDA0003944190610000026
Figure FDA0003944190610000027
其中Lk为k时刻层析斜路径方程、垂直约束方程构成的观测值列向量，
Figure FDA0003944190610000028
为利用k-1时刻的初值基于状态方程得到k时刻的中间预测值，σ为单位权重误差；Ak为k时刻的系数矩阵；
S42、按照如下方法得到观测方程的观测值的抗差权因子：
Figure FDA0003944190610000029
其中，
Figure FDA00039441906100000210
表示k时刻的观测方程中的观测值的抗差权因子，抗差权因于
Figure FDA00039441906100000211
与协方差矩阵
Figure FDA00039441906100000212
互为逆矩阵；pk表示k时刻的观测值原始权因子，c0和c1为常数；
S43、基于预测残差构造k时刻判别自适应因子的统计量，方法如下：
Figure FDA0003944190610000031
其中，tr表示求矩阵的迹，
Figure FDA0003944190610000032
表示预测残差vk的矩阵转置，
Figure FDA0003944190610000033
表示系数矩阵Ak的矩阵转置，
Figure FDA0003944190610000034
表示抗差权因子
Figure FDA0003944190610000035
的逆矩阵；Ak表示k时刻的观测方程系数矩阵；
Figure FDA0003944190610000036
表示k时刻的协方差矩阵；
然后构造出如下自适应因子函数模型：
Figure FDA0003944190610000037
其中c为常数；
由此得到k时刻的自适应因子αk；
S44、通过如下公式计算得到Kalman滤波增益矩阵
Figure FDA0003944190610000038
Figure FDA0003944190610000039
S45、按照如下方法进行层析k时刻的状态矩阵
Figure FDA00039441906100000310
和协方差矩阵
Figure FDA00039441906100000311
Figure FDA00039441906100000312
Figure FDA00039441906100000313
其中I表示单位阵；
以k时刻的状态矩阵
Figure FDA00039441906100000314
和协方差矩阵
Figure FDA00039441906100000315
作为初值按照方法S41～S45解算更新k+1时刻的状态矩阵
Figure FDA00039441906100000316
和协方差矩阵
Figure FDA00039441906100000317
以此类推，直至最后一个时刻完成上述循环计算，得到研究区每个层析网格的大气水汽密度值。
2.按照权利要求1所述的一种抗差自适应的GNSS水汽层析方法，其特征在于：在Kalman滤波解算模型中，状态方程中的状态噪声的协方差矩阵表示为Q，观测方程的观测噪声的协方差矩阵表示为R。
3.按照权利要求1所述的一种抗差自适应的GNSS水汽层析方法，其特征在于：GNSS斜路径湿延迟SWD通过如下公式得到：
Figure FDA0003944190610000041
其中ele表示GNSS卫星高度角；mw(ele)和mΔ(ele)分别为湿延迟和湿梯度的映射函数；
Figure FDA0003944190610000042
表示南北方向的湿大气水平梯度，
Figure FDA0003944190610000043
表示东西方向的湿大气水平梯度；Re表示GNSS观测值的非差验后残差，ZWD为测站对流层湿延迟。
4.按照权利要求3所述的一种抗差自适应的GNSS水汽层析方法，其特征在于：测站对流层湿延迟ZWD通过如下方法得到：
S21、采集研究区内GNSS测站的气象数据，通过如下公式计算得到天顶静力学延迟ZHD：
Figure FDA0003944190610000044
其中P为GNSS测站的地面气压，
Figure FDA0003944190610000045
为GNSS测站的地理纬度，h为GNSS测站的高程；
S22、采集研究区内GNSS测站的观测数据及精密星历数据，利用Gamit软件进行双差解算，得到GNSS天顶对流层延迟ZTD；测站对流层湿延迟ZWD通过如下公式得到：
ZWD＝ZTD-ZHD。
5.按照权利要求1所述的一种抗差自适应的GNSS水汽层析方法，其特征在于：所述水平格网的层析网格满足如下关系式：
s≤h1/tan(ele)其中s表示水平格网任一边的边长，ele表示GNSS 卫星高度角，h1表示研究区的层析高度。
6.按照权利要求3所述的一种抗差自适应的GNSS水汽层析方法，其特征在于：GNSS斜路径湿延迟SWD通过联合GAMIT软件利用测站对流层湿延迟ZWD处理得到。
CN202211430336.2A 2022-11-15 2022-11-15 一种抗差自适应的gnss水汽层析方法 Active CN115755103B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211430336.2A CN115755103B (zh) 2022-11-15 2022-11-15 一种抗差自适应的gnss水汽层析方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211430336.2A CN115755103B (zh) 2022-11-15 2022-11-15 一种抗差自适应的gnss水汽层析方法

Publications (2)

Publication Number Publication Date
CN115755103A true CN115755103A (zh) 2023-03-07
CN115755103B CN115755103B (zh) 2023-07-07

Family

ID=85371841

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211430336.2A Active CN115755103B (zh) 2022-11-15 2022-11-15 一种抗差自适应的gnss水汽层析方法

Country Status (1)

Country Link
CN (1) CN115755103B (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116068670B (zh) * 2023-03-30 2023-06-06 中国科学院精密测量科学与技术创新研究院 一种无地面网络区适用的北斗水汽场实时重构方法及装置
CN118193905B (zh) * 2024-05-17 2024-07-09 武汉大学 北斗水汽层析垂直约束构建与自适应调节方法及装置

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114859389A (zh) * 2022-04-18 2022-08-05 华力智芯(成都)集成电路有限公司 一种gnss多系统抗差自适应融合rtk解算方法

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114859389A (zh) * 2022-04-18 2022-08-05 华力智芯(成都)集成电路有限公司 一种gnss多系统抗差自适应融合rtk解算方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
杨飞: "地基GNSS反演水汽空间分布关键技术研究", 中国博士学位论文全文数据库 基础科学辑 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116068670B (zh) * 2023-03-30 2023-06-06 中国科学院精密测量科学与技术创新研究院 一种无地面网络区适用的北斗水汽场实时重构方法及装置
CN118193905B (zh) * 2024-05-17 2024-07-09 武汉大学 北斗水汽层析垂直约束构建与自适应调节方法及装置

Also Published As

Publication number Publication date
CN115755103B (zh) 2023-07-07

Similar Documents

Publication Publication Date Title
Sun et al. An ERA5‐based model for estimating tropospheric delay and weighted mean temperature over China with improved spatiotemporal resolutions
Chen et al. Assessment of ZTD derived from ECMWF/NCEP data with GPS ZTD over China
Chen et al. Seasonal global water mass budget and mean sea level variations
Zhao et al. High-precision ZTD model of altitude-related correction
CN110335355B (zh) 一种大型浅水湖泊水面高自动计算方法
CN112034490B (zh) 一种nwp反演对流层延迟改进方法
Yeh et al. Determining the precipitable water vapor thresholds under different rainfall strengths in Taiwan
CN113009531A (zh) 一种小尺度高精度的低空对流层大气折射率模型
CN101866021A (zh) 并行化Abel变换大气参数数据处理方法
Ssessanga et al. On imaging South African regional ionosphere using 4D‐var technique
CN115755103A (zh) 一种抗差自适应的gnss水汽层析方法
CN114415208B (zh) 一种附加外部数据集信息的地基gnss对流层顶探测方法
CN113960635B (zh) 一种顾及日变化的对流层延迟改正方法
CN113639893B (zh) 一种基于多气象因子的近地加权平均温度信息获取方法
Xia et al. Establishing a high-precision real-time ZTD model of China with GPS and ERA5 historical data and its application in PPP
CN112711022B (zh) 一种GNSS层析技术辅助的InSAR大气延迟改正方法
Mateus et al. Using InSAR data to improve the water vapor distribution downstream of the core of the South American low‐level Jet
Tang et al. High-spatial-resolution mapping of precipitable water vapour using SAR interferograms, GPS observations and ERA-Interim reanalysis
Izanlou et al. Enhanced troposphere tomography: Integration of GNSS and remote sensing data with optimal vertical constraints
Huang et al. A global grid model for the estimation of zenith tropospheric delay considering the variations at different altitudes
CN115099159B (zh) 基于神经网络且顾及地表差异的modis水汽反演方法
CN116466376A (zh) 数值预报模式辅助的实时ppp改善方法
Bakış et al. Analysis and comparison of spatial rainfall distribution applying different interpolation methods in Porsuk river basin, Turkey
CN115980317A (zh) 基于修正相位的地基gnss-r数据土壤水分估算方法
CN111538943B (zh) 新的高时空分辨率全球ztd垂直剖面格网模型构建方法

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