CN111126466A - 一种多源pwv数据融合方法 - Google Patents

一种多源pwv数据融合方法 Download PDF

Info

Publication number
CN111126466A
CN111126466A CN201911291390.1A CN201911291390A CN111126466A CN 111126466 A CN111126466 A CN 111126466A CN 201911291390 A CN201911291390 A CN 201911291390A CN 111126466 A CN111126466 A CN 111126466A
Authority
CN
China
Prior art keywords
pwv
model
data
different data
difference
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
CN201911291390.1A
Other languages
English (en)
Other versions
CN111126466B (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.)
Xian University of Science and Technology
Original Assignee
Xian University of Science and Technology
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 Xian University of Science and Technology filed Critical Xian University of Science and Technology
Priority to CN201911291390.1A priority Critical patent/CN111126466B/zh
Publication of CN111126466A publication Critical patent/CN111126466A/zh
Application granted granted Critical
Publication of CN111126466B publication Critical patent/CN111126466B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/25—Fusion techniques
    • G06F18/251—Fusion techniques of input or preprocessed data

Landscapes

  • Engineering & Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Artificial Intelligence (AREA)
  • Evolutionary Biology (AREA)
  • Evolutionary Computation (AREA)
  • Physics & Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种多源PWV数据融合方法，包括步骤：一、多源PWV数据的高程基准统一；二、GPT2w模型的PWV初值计算；三、获取高程基准统一的多源PWV数据；四、多项式PWV差值拟合；五、球谐函数模型PWV残差模拟；六、建立PWV混合模型。本发明引入GPT2w模型计算PWV初值，结合多项式PWV差值拟合和球谐函数模型PWV残差模拟，获取PWV混合模型，结合多项式PWV差值拟合和球谐函数模型PWV残差模拟过程中提出不同数据源的方程权重确定方法，并使用岭估计方法解决病态矩阵的最小二乘系数求解，建立的PWV混合模型，可得到任意位置以及高精度的事后PWV数据集。

Description

一种多源PWV数据融合方法
技术领域
本发明属于数据融合技术领域，具体涉及一种多源PWV数据融合方法。
背景技术
传统获取水汽数据主要是利用水汽辐射计、无线电探空仪和卫星观测等。但是该类方法易受观测成本和环境条件等限制，而且时空分辨率较差。全球定位系统(GlobalPositioning System，GPS)反演大气可降水量(Precipitable Water Vapor，PWV)后便因其低成本及高时空分辨率等优点使该方法得到了广泛应用。实验表明由全球导航卫星系统(Global Navigation Satellite System，GNSS)反演的PWV可达到与传统方法相同的精度。
但上述基于站点的测量手段由于稀疏且不均匀的分布以及观测的不连续性，很难得到准确的PWV趋势等重要信息。此外另一种水汽数据是基于大气环流模型将不同类型的观测数据同化的再分析数据集，具有记录均匀，全局覆盖和空间完整等优点。但是再分析数据在有限甚至没有观测数据同化的地区并不可靠。不同类型的水汽数据各有其优缺点，但是目前无论是模型联合还是数据融合同化的相关内容并不多。现急需建立好的的PWV产品融合模型。
发明内容
本发明所要解决的技术问题在于针对上述现有技术中的不足，提供一种多源PWV数据融合方法，引入GPT2w模型计算PWV初值，结合多项式PWV差值拟合和球谐函数模型PWV残差模拟，获取PWV混合模型，结合多项式PWV差值拟合和球谐函数模型PWV残差模拟过程中提出不同数据源的方程权重确定方法，并使用岭估计方法解决病态矩阵的最小二乘系数求解，建立的PWV混合模型，可得到任意位置以及高精度的事后PWV数据集，便于推广使用。
为解决上述技术问题，本发明采用的技术方案是：一种多源PWV数据融合方法，其特征在于，该方法包括以下步骤：
步骤一、多源PWV数据的高程基准统一：获取多源PWV数据，所述多源PWV数据包括GNSSPWV数据、ECMWF再分析数据和探空水汽数据，所述GNSSPWV数据的高程基准为大地高，ECMWF再分析数据的高程基准为位势，探空水汽数据的高程基准为位势高；
将ECMWF再分析数据的位势和探空水汽数据的位势高均统一到大地高；
步骤二、GPT2w模型的PWV初值计算：根据公式
Figure BDA0002319231770000021
计算GPT2w模型的天顶湿延迟ZWDGPT2w，其中，k3和k'2均为气体常数，k'2＝16.529K·mb-1，k3＝373900K2·mb-1，Tm为加权平均温度，λ为水汽递减因子，gm为当地的实际平均重力加速度，单位为m/s2，es为水汽压；
根据公式PWV0＝Π·ZWDGPT2w，计算GPT2w模型的PWV初值PWV0，其中，Π为转换因子且
Figure BDA0002319231770000022
Rv为水汽气体常数且Rv＝461.495J·kg-1·K-1；
步骤三、获取高程基准统一的多源PWV数据：分别获取高程基准统一后的GNSSPWV数据、ECMWF再分析数据和探空水汽数据；
步骤四、多项式PWV差值拟合，过程如下：
步骤401、根据公式
Figure BDA0002319231770000031
对GNSS数据中PWV值和GPT2w模型的PWV初值之间的差值方程
Figure BDA0002319231770000032
探空水汽数据的PWV值和GPT2w模型的PWV初值之间的差值方程
Figure BDA0002319231770000033
以及ECMWF再分析数据的PWV值和GPT2w模型的PWV初值之间的差值方程
Figure BDA0002319231770000034
进行多项式差值拟合，其中，
Figure BDA0002319231770000035
为纬度，φ为经度，h为高度，i表示GNSS数据源，j表示探空水汽数据源，q表示ECMWF再分析数据源，a0-a9表示多项式函数模型的拟合系数；
步骤402、初始化不同数据源的差值方程权重，令不同数据源的差值方程权重均为1；
步骤403、根据公式
Figure BDA0002319231770000036
计算不同数据源差值方程的后验单位权方差
Figure BDA0002319231770000037
其中，PMODEL,D为不同数据源的差值方程权重矩阵且MODEL＝i、j、q，VMODEL,D为不同数据源差值方程的后验误差矩阵，m为不同数据源的个数，
Figure BDA0002319231770000038
Figure BDA0002319231770000039
PD为包含Pi,D、Pj,D和Pq,D的对角权阵，
Figure BDA0002319231770000041
tr(·)为矩阵的秩函数；
步骤404、根据公式
Figure BDA0002319231770000042
计算不同数据源差值方程的检验统计量
Figure BDA0002319231770000043
其中，nMODEL为不同数据源的样本数，
Figure BDA0002319231770000044
为各数据源差值方程的初始后验单位权方差均值且
Figure BDA0002319231770000045
步骤405、利用Bartlett检验判断不同数据源差值方程的检验统计量
Figure BDA0002319231770000046
是否小于自由度为2时查找表检索的临界值
Figure BDA0002319231770000047
当不同数据源差值方程的检验统计量
Figure BDA0002319231770000048
小于自由度为2时查找表检索的临界值
Figure BDA0002319231770000049
时，则不同数据源差值方程的后验单位权方差在统计学上相等，进而确定不同数据源的差值方程权重，然后执行步骤407；否则执行步骤406；
步骤406、根据公式
Figure BDA00023192317700000410
更新不同数据源的差值方程权重，获取更新后的不同数据源的差值方程权重P'MODEL,D，并将更新后的不同数据源的差值方程权重P'MODEL,D视为新的不同数据源的差值方程权重矩阵PMODEL,D循环步骤403，其中，c为数据更新常数；
步骤407、根据公式
Figure BDA0002319231770000051
估计多项式函数模型的拟合系数xD，其中，
Figure BDA0002319231770000052
ηD为第一正则化参数，I为单位矩阵，
Figure BDA0002319231770000053
步骤五、球谐函数模型PWV残差模拟，过程如下：
步骤501、将步骤407获取的多项式函数模型的拟合系数xD代入步骤401，获取虚拟差值，再通过虚拟差值和真值差值之间的差值获取初始PWV残差
Figure BDA0002319231770000054
然后对初始PWV残差
Figure BDA0002319231770000055
进行高程面统一，获取PWV残差
Figure BDA0002319231770000056
其中，
Figure BDA0002319231770000057
κ和
Figure BDA0002319231770000058
均为球谐函数的阶次且
Figure BDA0002319231770000059
Figure BDA00023192317700000517
为勒让德函数且
Figure BDA00023192317700000510
Θ为κ的最大阶次取值，Ω为
Figure BDA00023192317700000511
的最大阶次取值，
Figure BDA00023192317700000512
θ为中间变量，
Figure BDA00023192317700000513
和
Figure BDA00023192317700000514
均为球谐函数模型的拟合系数；
步骤502、根据公式
Figure BDA00023192317700000515
对各数据源的样本的PWV残差方程进行球谐函数模型拟合，其中，nmodel为各数据源的样本总数且nmodel＝ni+nj+nq，ζ为样本编号且ζ＝1，2，...，nmodel，
Figure BDA00023192317700000516
为第ζ个样本数据的PWV残差，
Figure BDA0002319231770000061
Aζ为第ζ个样本数据的球谐函数模型的中间参数矩阵且
Figure BDA0002319231770000062
Figure BDA0002319231770000063
和
Figure BDA0002319231770000064
均为第ζ个样本数据的球谐函数模型的中间参数且
Figure BDA0002319231770000065
κ'和
Figure BDA0002319231770000066
均为球谐函数的阶次编号且且κ'＝0,1,...,Θ，
Figure BDA0002319231770000067
步骤503、初始化不同数据源的残差方程权重，令不同数据源的残差方程权重均为1；
步骤504、根据公式
Figure BDA0002319231770000068
计算不同数据源残差方程的后验单位权方差
Figure BDA0002319231770000069
其中，PMODEL,R为不同数据源的残差方程权重矩阵，VMODEL,R为不同数据源残差方程的后验误差矩阵，
Figure BDA00023192317700000610
Figure BDA00023192317700000611
PR为包含Pi,R、Pj,R和Pq,R的对角权阵，
Figure BDA00023192317700000612
AMODEL,R为AR中不同数据源的样本数据的球谐函数模型的中间参数矩阵；
步骤505、根据公式
Figure BDA00023192317700000613
计算不同数据源残差方程的检验统计量
Figure BDA00023192317700000614
其中，
Figure BDA00023192317700000615
为各数据源残差方程的初始后验单位权方差均值且
Figure BDA00023192317700000616
步骤506、利用Bartlett检验判断不同数据源残差方程的检验统计量
Figure BDA00023192317700000617
是否小于自由度为2时查找表检索的临界值
Figure BDA00023192317700000618
不同数据源残差方程的检验统计量
Figure BDA0002319231770000071
小于自由度为2时查找表检索的临界值
Figure BDA0002319231770000072
时，则不同数据源残差方程的后验单位权方差在统计学上相等，进而确定不同数据源的残差方程权重，然后执行步骤508；否则执行步骤507；
步骤507、根据公式
Figure BDA0002319231770000073
更新不同数据源的残差方程权重，获取更新后的不同数据源的残差方程权重P'MODEL,R，并将更新后的不同数据源的残差方程权重P'MODEL,R视为新的不同数据源的残差方程权重矩阵PMODEL,R循环步骤504；
步骤508、根据公式
Figure BDA0002319231770000074
估计球谐函数模型残差拟合系数矩阵xR，ηR为第二正则化参数，
Figure BDA0002319231770000075
步骤六、建立PWV混合模型：根据公式
Figure BDA0002319231770000076
建立PWV混合模型PWVMODEL，其中，
Figure BDA0002319231770000077
为不同数据源的GPT2w模型计算的PWV初值，DPWVMODEL为不同数据源的多项式拟合的PWV差值，RPWVMODEL为不同数据源的球谐函数拟合的PWV残差。
上述的一种多源PWV数据融合方法，其特征在于：步骤一中将ECMWF再分析数据的位势和探空水汽数据的位势高均统一到大地高的方法如下：
步骤101、根据公式
Figure BDA0002319231770000078
将位势GP转换为位势高GPH，其中，g为当地的实际重力加速度，单位为m/s2；
步骤102、根据公式
Figure BDA0002319231770000079
将位势高GPH转换为正高
Figure BDA00023192317700000710
其中，
Figure BDA00023192317700000711
为地球椭球化表面的正常重力值，
Figure BDA00023192317700000712
为纬度为
Figure BDA0002319231770000081
时的地球有效半径，Y45为地球椭球化表面上纬度为45°的重力值；
步骤103、根据公式
Figure BDA0002319231770000082
利用EGM2008将正高
Figure BDA0002319231770000083
转化为大地高h'，其中，N为大地水准面差距。
上述的一种多源PWV数据融合方法，其特征在于：步骤三中GNSSPWV数据中的PWV数据获取过程如下：
利用精密单点定位原理对接收到的GNSS数据进行解算得到天顶对流层总延迟ZTD，根据公式
Figure BDA0002319231770000084
计算天顶静力学延迟ZHD和GNSS数据中天顶湿延迟量ZWDGNSS，其中，P为地表气压，H为测站大地高；
根据公式PWVGNSS＝Π·ZWDGNSS，计算GNSSPWV数据中的PWV值PWVGNSS；
步骤三中ECMWF再分析数据中的PWV数据从European entre for Medium-RangeWeather Forecasts官网上下载；
步骤三中探空水汽数据中的PWV数据从Integrated Global Radiosonde Archive官网上下载。
上述的一种多源PWV数据融合方法，其特征在于：步骤501中根据公式
Figure BDA0002319231770000085
对高度为hδ的初始PWV残差进行高程面转换，实现高度hδ到高度为hβ的初始PWV残差的高程面统一。
上述的一种多源PWV数据融合方法，其特征在于：步骤407中第一正则化参数ηD和步骤508中的第二正则化参数ηR均通过L曲线法获取。
本发明与现有技术相比具有以下优点：
1、本发明引入GPT2w模型计算PWV初值，结合多项式PWV差值拟合和球谐函数模型PWV残差模拟，获取PWV混合模型，充分利用现有数据，极大提高了数据水平分辨率，满足PWV应用于气象临近预报的精度要求，建立的PWV混合模型，可得到任意位置以及高精度的事后PWV数据集，对监测和评价天气系统的变化具有重要的研究意义和应用价值，便于推广使用。
2、本发明多项式PWV差值拟合过程中提出不同数据源的差值方程权重确定方法，通过Bartlett检验评估验后单位权方差，确定不同数据源的最优差值方程权重，并使用岭估计方法解决病态矩阵的最小二乘系数求解，获取可靠的多项式函数模型的拟合系数，使用效果好。
3、本发明球谐函数模型PWV残差模拟过程中提出不同数据源的残差方程权重确定方法，通过Bartlett检验评估验后单位权方差，确定不同数据源的最优残差方程权重，并使用岭估计方法解决大区域球谐函数模型病态矩阵的系数求解问题，获取可靠的球谐函数模型的拟合系数，确保融合后PWV数据集的精度和可靠性，另外，球谐函数模型PWV残差模拟之前通过虚拟差值和真值差值之间的差值获取初始PWV残差，然后对初始PWV残差进行高程面统一，获取PWV残差，保证PWV残差数据均具有相同的高程面，便于推广使用。
4、本发明方法步骤简单，将获取的GNSSPWV数据、ECMWF再分析数据和探空水汽数据的高程基准统一，为后续计算提供了计算的基础，利用GPT2w模型的PWV初值计算为PWV混合模型提供了可靠的初值，结合多项式PWV差值拟合和球谐函数模型PWV残差模拟，建立PWV混合模型，可得到任意位置以及高精度的事后PWV数据集，弥补了现缺少PWV产品融合模型的缺陷。
综上所述，本发明引入GPT2w模型计算PWV初值，结合多项式PWV差值拟合和球谐函数模型PWV残差模拟，获取PWV混合模型，结合多项式PWV差值拟合和球谐函数模型PWV残差模拟过程中提出不同数据源的方程权重确定方法，并使用岭估计方法解决病态矩阵的最小二乘系数求解，建立的PWV混合模型，可得到任意位置以及高精度的事后PWV数据集，便于推广使用。
下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。
附图说明
图1为本发明的方法流程框图。
具体实施方式
如图1所示，本发明的一种多源PWV数据融合方法，包括以下步骤：
步骤一、多源PWV数据的高程基准统一：获取多源PWV数据，所述多源PWV数据包括GNSSPWV数据、ECMWF再分析数据和探空水汽数据，所述GNSSPWV数据的高程基准为大地高，ECMWF再分析数据的高程基准为位势，探空水汽数据的高程基准为位势高；
将ECMWF再分析数据的位势和探空水汽数据的位势高均统一到大地高；
本实施例中，步骤一中将ECMWF再分析数据的位势和探空水汽数据的位势高均统一到大地高的方法如下：
步骤101、根据公式
Figure BDA0002319231770000101
将位势GP转换为位势高GPH，其中，g为当地的实际重力加速度，单位为m/s2；
步骤102、根据公式
Figure BDA0002319231770000102
将位势高GPH转换为正高
Figure BDA0002319231770000103
其中，
Figure BDA0002319231770000104
为地球椭球化表面的正常重力值，
Figure BDA0002319231770000105
为纬度为
Figure BDA0002319231770000106
时的地球有效半径，Y45为地球椭球化表面上纬度为45°的重力值；
步骤103、根据公式
Figure BDA0002319231770000107
利用EGM2008将正高
Figure BDA0002319231770000108
转化为大地高h'，其中，N为大地水准面差距。
步骤二、GPT2w模型的PWV初值计算：根据公式
Figure BDA0002319231770000109
计算GPT2w模型的天顶湿延迟ZWDGPT2w，其中，k3和k'2均为气体常数，k'2＝16.529K·mb-1，k3＝373900K2·mb-1，Tm为加权平均温度，λ为水汽递减因子，gm为当地的实际平均重力加速度，单位为m/s2，es为水汽压；
根据公式PWV0＝Π·ZWDGPT2w，计算GPT2w模型的PWV初值PWV0，其中，Π为转换因子且
Figure BDA0002319231770000111
Rv为水汽气体常数且Rv＝461.495J·kg-1·K-1；
需要说明的是，GPT2w(Global pressure and temperature model 2wet)为全球温度气压湿度模型，引入GPT2w模型计算PWV初值，结合多项式PWV差值拟合和球谐函数模型PWV残差模拟，获取PWV混合模型，充分利用现有数据，极大提高了数据水平分辨率，满足PWV应用于气象临近预报的精度要求，建立的PWV混合模型，可得到任意位置以及高精度的事后PWV数据集，对监测和评价天气系统的变化具有重要的研究意义和应用价值。
步骤三、获取高程基准统一的多源PWV数据：分别获取高程基准统一后的GNSSPWV数据、ECMWF再分析数据和探空水汽数据；
本实施例中，步骤三中GNSSPWV数据中的PWV数据获取过程如下：
利用精密单点定位原理对接收到的GNSS数据进行解算得到天顶对流层总延迟ZTD，根据公式
Figure BDA0002319231770000112
计算天顶静力学延迟ZHD和GNSS数据中天顶湿延迟量ZWDGNSS，其中，P为地表气压，H为测站大地高；
根据公式PWVGNSS＝Π·ZWDGNSS，计算GNSSPWV数据中的PWV值PWVGNSS；
步骤三中ECMWF再分析数据中的PWV数据从European entre for Medium-RangeWeather Forecasts官网上下载；
步骤三中探空水汽数据中的PWV数据从Integrated Global Radiosonde Archive官网上下载。
步骤四、多项式PWV差值拟合，过程如下：
步骤401、根据公式
Figure BDA0002319231770000121
对GNSS数据中PWV值和GPT2w模型的PWV初值之间的差值方程
Figure BDA0002319231770000122
探空水汽数据的PWV值和GPT2w模型的PWV初值之间的差值方程
Figure BDA0002319231770000123
以及ECMWF再分析数据的PWV值和GPT2w模型的PWV初值之间的差值方程
Figure BDA0002319231770000124
进行多项式差值拟合，其中，
Figure BDA0002319231770000125
为纬度，φ为经度，h为高度，i表示GNSS数据源，j表示探空水汽数据源，q表示ECMWF再分析数据源，a0-a9表示多项式函数模型的拟合系数；
步骤402、初始化不同数据源的差值方程权重，令不同数据源的差值方程权重均为1；
步骤403、根据公式
Figure BDA0002319231770000126
计算不同数据源差值方程的后验单位权方差
Figure BDA0002319231770000127
其中，PMODEL,D为不同数据源的差值方程权重矩阵且MODEL＝i、j、q，VMODEL,D为不同数据源差值方程的后验误差矩阵，m为不同数据源的个数，
Figure BDA0002319231770000128
Figure BDA0002319231770000129
PD为包含Pi,D、Pj,D和Pq,D的对角权阵，
Figure BDA0002319231770000131
tr(·)为矩阵的秩函数；
步骤404、根据公式
Figure BDA0002319231770000132
计算不同数据源差值方程的检验统计量
Figure BDA0002319231770000133
其中，nMODEL为不同数据源的样本数，
Figure BDA0002319231770000134
为各数据源差值方程的初始后验单位权方差均值且
Figure BDA0002319231770000135
步骤405、利用Bartlett检验判断不同数据源差值方程的检验统计量
Figure BDA0002319231770000136
是否小于自由度为2时查找表检索的临界值
Figure BDA0002319231770000137
当不同数据源差值方程的检验统计量
Figure BDA0002319231770000138
小于自由度为2时查找表检索的临界值
Figure BDA0002319231770000139
时，则不同数据源差值方程的后验单位权方差在统计学上相等，进而确定不同数据源的差值方程权重，然后执行步骤407；否则执行步骤406；
步骤406、根据公式
Figure BDA00023192317700001310
更新不同数据源的差值方程权重，获取更新后的不同数据源的差值方程权重P'MODEL,D，并将更新后的不同数据源的差值方程权重P'MODEL,D视为新的不同数据源的差值方程权重矩阵PMODEL,D循环步骤403，其中，c为数据更新常数；
步骤407、根据公式
Figure BDA0002319231770000141
估计多项式函数模型的拟合系数xD，其中，
Figure BDA0002319231770000142
ηD为第一正则化参数，I为单位矩阵，
Figure BDA0002319231770000143
需要说明的是，多项式PWV差值拟合过程中提出不同数据源的差值方程权重确定方法，通过Bartlett检验评估验后单位权方差，确定不同数据源的最优差值方程权重，并使用岭估计方法解决病态矩阵的最小二乘系数求解，获取可靠的多项式函数模型的拟合系数，使用效果好。
步骤五、球谐函数模型PWV残差模拟，过程如下：
步骤501、将步骤407获取的多项式函数模型的拟合系数xD代入步骤401，获取虚拟差值，再通过虚拟差值和真值差值之间的差值获取初始PWV残差
Figure BDA0002319231770000144
然后对初始PWV残差
Figure BDA0002319231770000145
进行高程面统一，获取PWV残差
Figure BDA0002319231770000146
其中，
Figure BDA0002319231770000147
κ和
Figure BDA0002319231770000148
均为球谐函数的阶次且
Figure BDA0002319231770000149
Figure BDA00023192317700001414
为勒让德函数且
Figure BDA00023192317700001410
Θ为κ的最大阶次取值，Ω为
Figure BDA00023192317700001411
的最大阶次取值，
Figure BDA00023192317700001412
θ为中间变量，和
Figure BDA00023192317700001416
均为球谐函数模型的拟合系数；
本实施例中，步骤501中根据公式
Figure BDA00023192317700001413
对高度为hδ的初始PWV残差进行高程面转换，实现高度hδ到高度为hβ的初始PWV残差的高程面统一。
步骤502、根据公式
Figure BDA0002319231770000151
对各数据源的样本的PWV残差方程进行球谐函数模型拟合，其中，nmodel为各数据源的样本总数且nmodel＝ni+nj+nq，ζ为样本编号且ζ＝1，2，...，nmodel，
Figure BDA0002319231770000152
为第ζ个样本数据的PWV残差，
Figure BDA0002319231770000153
Aζ为第ζ个样本数据的球谐函数模型的中间参数矩阵且
Figure BDA0002319231770000154
Figure BDA0002319231770000155
和
Figure BDA0002319231770000156
均为第ζ个样本数据的球谐函数模型的中间参数且
Figure BDA0002319231770000157
κ'和
Figure BDA0002319231770000158
均为球谐函数的阶次编号且且κ'＝0,1,...,Θ，
Figure BDA0002319231770000159
步骤503、初始化不同数据源的残差方程权重，令不同数据源的残差方程权重均为1；
步骤504、根据公式
Figure BDA00023192317700001510
计算不同数据源残差方程的后验单位权方差
Figure BDA00023192317700001511
其中，PMODEL,R为不同数据源的残差方程权重矩阵，VMODEL,R为不同数据源残差方程的后验误差矩阵，
Figure BDA00023192317700001512
Figure BDA00023192317700001513
PR为包含Pi,R、Pj,R和Pq,R的对角权阵，
Figure BDA00023192317700001514
AMODEL,R为AR中不同数据源的样本数据的球谐函数模型的中间参数矩阵；
步骤505、根据公式
Figure BDA0002319231770000161
计算不同数据源残差方程的检验统计量
Figure BDA0002319231770000162
其中，
Figure BDA0002319231770000163
为各数据源残差方程的初始后验单位权方差均值且
Figure BDA0002319231770000164
步骤506、利用Bartlett检验判断不同数据源残差方程的检验统计量
Figure BDA0002319231770000165
是否小于自由度为2时查找表检索的临界值
Figure BDA0002319231770000166
不同数据源残差方程的检验统计量
Figure BDA0002319231770000167
小于自由度为2时查找表检索的临界值
Figure BDA0002319231770000168
时，则不同数据源残差方程的后验单位权方差在统计学上相等，进而确定不同数据源的残差方程权重，然后执行步骤508；否则执行步骤507；
步骤507、根据公式
Figure BDA0002319231770000169
更新不同数据源的残差方程权重，获取更新后的不同数据源的残差方程权重P'MODEL,R，并将更新后的不同数据源的残差方程权重P'MODEL,R视为新的不同数据源的残差方程权重矩阵PMODEL,R循环步骤504；
步骤508、根据公式
Figure BDA00023192317700001610
估计球谐函数模型残差拟合系数矩阵xR，ηR为第二正则化参数，
Figure BDA00023192317700001611
需要说明的是，球谐函数模型PWV残差模拟过程中提出不同数据源的残差方程权重确定方法，通过Bartlett检验评估验后单位权方差，确定不同数据源的最优残差方程权重，并使用岭估计方法解决大区域球谐函数模型病态矩阵的系数求解问题，获取可靠的球谐函数模型的拟合系数，确保融合后PWV数据集的精度和可靠性，另外，球谐函数模型PWV残差模拟之前通过虚拟差值和真值差值之间的差值获取初始PWV残差，然后对初始PWV残差进行高程面统一，获取PWV残差，保证PWV残差数据均具有相同的高程面，便于推广使用。
步骤六、建立PWV混合模型：根据公式
Figure BDA0002319231770000172
建立PWV混合模型PWVMODEL，其中，
Figure BDA0002319231770000171
为不同数据源的GPT2w模型计算的PWV初值，DPWVMODEL为不同数据源的多项式拟合的PWV差值，RPWVMODEL为不同数据源的球谐函数拟合的PWV残差。
本实施例中，步骤407中第一正则化参数ηD和步骤508中的第二正则化参数ηR均通过L曲线法获取。
本发明使用时，将获取的GNSSPWV数据、ECMWF再分析数据和探空水汽数据的高程基准统一，为后续计算提供了计算的基础，利用GPT2w模型的PWV初值计算为PWV混合模型提供了可靠的初值，结合多项式PWV差值拟合和球谐函数模型PWV残差模拟，建立PWV混合模型，可得到任意位置以及高精度的事后PWV数据集，弥补了现缺少PWV产品融合模型的缺陷。
以上所述，仅是本发明的较佳实施例，并非对本发明作任何限制，凡是根据本发明技术实质对以上实施例所作的任何简单修改、变更以及等效结构变化，均仍属于本发明技术方案的保护范围内。

Claims (5)

1.一种多源PWV数据融合方法，其特征在于，该方法包括以下步骤：
步骤一、多源PWV数据的高程基准统一：获取多源PWV数据，所述多源PWV数据包括GNSSPWV数据、ECMWF再分析数据和探空水汽数据，所述GNSS PWV数据的高程基准为大地高，ECMWF再分析数据的高程基准为位势，探空水汽数据的高程基准为位势高；
将ECMWF再分析数据的位势和探空水汽数据的位势高均统一到大地高；
步骤二、GPT2w模型的PWV初值计算：根据公式
Figure FDA0002319231760000011
计算GPT2w模型的天顶湿延迟ZWDGPT2w，其中，k3和k'2均为气体常数，k'2＝16.529K·mb-1，k3＝373900K2·mb-1，Tm为加权平均温度，λ为水汽递减因子，gm为当地的实际平均重力加速度，单位为m/s2，es为水汽压；
根据公式PWV0＝Π·ZWDGPT2w，计算GPT2w模型的PWV初值PWV0，其中，Π为转换因子且
Figure FDA0002319231760000012
Rv为水汽气体常数且Rv＝461.495J·kg-1·K-1；
步骤三、获取高程基准统一的多源PWV数据：分别获取高程基准统一后的GNSS PWV数据、ECMWF再分析数据和探空水汽数据；
步骤四、多项式PWV差值拟合，过程如下：
步骤401、根据公式
Figure FDA0002319231760000021
对GNSS数据中PWV值和GPT2w模型的PWV初值之间的差值方程
Figure FDA0002319231760000022
探空水汽数据的PWV值和GPT2w模型的PWV初值之间的差值方程
Figure FDA0002319231760000023
以及ECMWF再分析数据的PWV值和GPT2w模型的PWV初值之间的差值方程
Figure FDA0002319231760000024
进行多项式差值拟合，其中，
Figure FDA0002319231760000025
为纬度，φ为经度，h为高度，i表示GNSS数据源，j表示探空水汽数据源，q表示ECMWF再分析数据源，a0-a9表示多项式函数模型的拟合系数；
步骤402、初始化不同数据源的差值方程权重，令不同数据源的差值方程权重均为1；
步骤403、根据公式
Figure FDA0002319231760000026
计算不同数据源差值方程的后验单位权方差
Figure FDA0002319231760000027
其中，PMODEL,D为不同数据源的差值方程权重矩阵且MODEL＝i、j、q，VMODEL,D为不同数据源差值方程的后验误差矩阵，m为不同数据源的个数，
Figure FDA0002319231760000028
Figure FDA0002319231760000029
PD为包含Pi,D、Pj,D和Pq,D的对角权阵，
Figure FDA0002319231760000031
tr(·)为矩阵的秩函数；
步骤404、根据公式
Figure FDA0002319231760000032
计算不同数据源差值方程的检验统计量
Figure FDA0002319231760000033
其中，nMODEL为不同数据源的样本数，
Figure FDA0002319231760000034
为各数据源差值方程的初始后验单位权方差均值且
Figure FDA0002319231760000035
步骤405、利用Bartlett检验判断不同数据源差值方程的检验统计量
Figure FDA0002319231760000036
是否小于自由度为2时查找表检索的临界值
Figure FDA0002319231760000037
当不同数据源差值方程的检验统计量
Figure FDA0002319231760000038
小于自由度为2时查找表检索的临界值
Figure FDA0002319231760000039
时，则不同数据源差值方程的后验单位权方差在统计学上相等，进而确定不同数据源的差值方程权重，然后执行步骤407；否则执行步骤406；
步骤406、根据公式
Figure FDA00023192317600000310
更新不同数据源的差值方程权重，获取更新后的不同数据源的差值方程权重P'MODEL,D，并将更新后的不同数据源的差值方程权重P'MODEL,D视为新的不同数据源的差值方程权重矩阵PMODEL,D循环步骤403，其中，c为数据更新常数；
步骤407、根据公式
Figure FDA0002319231760000041
估计多项式函数模型的拟合系数xD，其中，
Figure FDA0002319231760000042
ηD为第一正则化参数，I为单位矩阵，
Figure FDA0002319231760000043
步骤五、球谐函数模型PWV残差模拟，过程如下：
步骤501、将步骤407获取的多项式函数模型的拟合系数xD代入步骤401，获取虚拟差值，再通过虚拟差值和真值差值之间的差值获取初始PWV残差
Figure FDA0002319231760000044
然后对初始PWV残差
Figure FDA0002319231760000045
进行高程面统一，获取PWV残差
Figure FDA0002319231760000046
其中，
Figure FDA0002319231760000047
κ和
Figure FDA0002319231760000048
均为球谐函数的阶次且
Figure FDA0002319231760000049
Figure FDA00023192317600000415
为勒让德函数且
Figure FDA00023192317600000410
Θ为κ的最大阶次取值，Ω为
Figure FDA00023192317600000411
的最大阶次取值，
Figure FDA00023192317600000412
θ为中间变量，
Figure FDA00023192317600000416
和
Figure FDA00023192317600000417
均为球谐函数模型的拟合系数；
步骤502、根据公式
Figure FDA00023192317600000413
对各数据源的样本的PWV残差方程进行球谐函数模型拟合，其中，nmodel为各数据源的样本总数且nmodel＝ni+nj+nq，ζ为样本编号且ζ＝1，2，...，nmodel，
Figure FDA00023192317600000414
为第ζ个样本数据的PWV残差，
Figure FDA0002319231760000051
Aζ为第ζ个样本数据的球谐函数模型的中间参数矩阵且
Figure FDA0002319231760000052
Figure FDA0002319231760000053
和
Figure FDA0002319231760000054
均为第ζ个样本数据的球谐函数模型的中间参数且
Figure FDA0002319231760000055
κ'和
Figure FDA0002319231760000056
均为球谐函数的阶次编号且κ'＝0,1,...,Θ，
Figure FDA0002319231760000057
步骤503、初始化不同数据源的残差方程权重，令不同数据源的残差方程权重均为1；
步骤504、根据公式
Figure FDA0002319231760000058
计算不同数据源残差方程的后验单位权方差
Figure FDA0002319231760000059
其中，PMODEL,R为不同数据源的残差方程权重矩阵，VMODEL,R为不同数据源残差方程的后验误差矩阵，
Figure FDA00023192317600000510
Figure FDA00023192317600000511
PR为包含Pi,R、Pj,R和Pq,R的对角权阵，
Figure FDA00023192317600000512
AMODEL,R为AR中不同数据源的样本数据的球谐函数模型的中间参数矩阵；
步骤505、根据公式
Figure FDA00023192317600000513
计算不同数据源残差方程的检验统计量
Figure FDA00023192317600000514
其中，
Figure FDA00023192317600000515
为各数据源残差方程的初始后验单位权方差均值且
Figure FDA00023192317600000516
步骤506、利用Bartlett检验判断不同数据源残差方程的检验统计量
Figure FDA00023192317600000517
是否小于自由度为2时查找表检索的临界值
Figure FDA00023192317600000518
不同数据源残差方程的检验统计量
Figure FDA0002319231760000061
小于自由度为2时查找表检索的临界值
Figure FDA0002319231760000062
时，则不同数据源残差方程的后验单位权方差在统计学上相等，进而确定不同数据源的残差方程权重，然后执行步骤508；否则执行步骤507；
步骤507、根据公式
Figure FDA0002319231760000063
更新不同数据源的残差方程权重，获取更新后的不同数据源的残差方程权重P'MODEL,R，并将更新后的不同数据源的残差方程权重P'MODEL,R视为新的不同数据源的残差方程权重矩阵PMODEL,R循环步骤504；
步骤508、根据公式
Figure FDA0002319231760000064
估计球谐函数模型残差拟合系数矩阵xR，ηR为第二正则化参数，
Figure FDA0002319231760000065
步骤六、建立PWV混合模型：根据公式
Figure FDA0002319231760000066
建立PWV混合模型PWVMODEL，其中，
Figure FDA0002319231760000067
为不同数据源的GPT2w模型计算的PWV初值，DPWVMODEL为不同数据源的多项式拟合的PWV差值，RPWVMODEL为不同数据源的球谐函数拟合的PWV残差。
2.按照权利要求1所述的一种多源PWV数据融合方法，其特征在于：步骤一中将ECMWF再分析数据的位势和探空水汽数据的位势高均统一到大地高的方法如下：
步骤101、根据公式
Figure FDA0002319231760000068
将位势GP转换为位势高GPH，其中，g为当地的实际重力加速度，单位为m/s2；
步骤102、根据公式
Figure FDA0002319231760000071
将位势高GPH转换为正高
Figure FDA0002319231760000072
其中，
Figure FDA0002319231760000073
为地球椭球化表面的正常重力值，
Figure FDA0002319231760000074
为纬度为
Figure FDA0002319231760000075
时的地球有效半径，Y45为地球椭球化表面上纬度为45°的重力值；
步骤103、根据公式
Figure FDA0002319231760000076
利用EGM 2008将正高
Figure FDA0002319231760000077
转化为大地高h'，其中，N为大地水准面差距。
3.按照权利要求1所述的一种多源PWV数据融合方法，其特征在于：步骤三中GNSS PWV数据中的PWV数据获取过程如下：
利用精密单点定位原理对接收到的GNSS数据进行解算得到天顶对流层总延迟ZTD，根据公式
Figure FDA0002319231760000078
计算天顶静力学延迟ZHD和GNSS数据中天顶湿延迟量ZWDGNSS，其中，P为地表气压，H为测站大地高；
根据公式PWVGNSS＝Π·ZWDGNSS，计算GNSS PWV数据中的PWV值PWVGNSS；
步骤三中ECMWF再分析数据中的PWV数据从European entre for Medium-RangeWeather Forecasts官网上下载；
步骤三中探空水汽数据中的PWV数据从Integrated Global Radiosonde Archive官网上下载。
4.按照权利要求1所述的一种多源PWV数据融合方法，其特征在于：步骤501中根据公式
Figure FDA0002319231760000079
对高度为hδ的初始PWV残差进行高程面转换，实现高度hδ到高度为hβ的初始PWV残差的高程面统一。
5.按照权利要求1所述的一种多源PWV数据融合方法，其特征在于：步骤407中第一正则化参数ηD和步骤508中的第二正则化参数ηR均通过L曲线法获取。
CN201911291390.1A 2019-12-16 2019-12-16 一种多源pwv数据融合方法 Expired - Fee Related CN111126466B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911291390.1A CN111126466B (zh) 2019-12-16 2019-12-16 一种多源pwv数据融合方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911291390.1A CN111126466B (zh) 2019-12-16 2019-12-16 一种多源pwv数据融合方法

Publications (2)

Publication Number Publication Date
CN111126466A true CN111126466A (zh) 2020-05-08
CN111126466B CN111126466B (zh) 2020-09-22

Family

ID=70499085

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911291390.1A Expired - Fee Related CN111126466B (zh) 2019-12-16 2019-12-16 一种多源pwv数据融合方法

Country Status (1)

Country Link
CN (1) CN111126466B (zh)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114019585A (zh) * 2021-10-11 2022-02-08 武汉大学 一种大高差地区高精度定位cors网fkp解算方法
CN115993668A (zh) * 2023-03-22 2023-04-21 成都云智北斗科技有限公司 一种基于多项式改正和神经网络的pwv重建方法及系统
CN117405968A (zh) * 2023-10-30 2024-01-16 广东海洋大学 基于大数据的设备能耗检测方法和系统

Citations (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20090189802A1 (en) * 2008-01-25 2009-07-30 Tillotson Brian J System and method for using iridium satellite signals for meteorological measurements
JP2010060443A (ja) * 2008-09-04 2010-03-18 Japan Weather Association 気象予測装置、方法及びプログラム
CN103323888A (zh) * 2013-04-24 2013-09-25 东南大学 Gnss大气探测数据中对流层延迟误差的消除方法
CN106324620A (zh) * 2016-08-02 2017-01-11 中国人民解放军空军工程大学 一种不依赖地表气象数据实时测量的对流层天顶延迟方法
US20170045624A1 (en) * 2015-08-14 2017-02-16 Trimble Navigation Limited Navigation satellite system positioning involving the generation of advanced correction information
CN107843943A (zh) * 2017-10-30 2018-03-27 西安科技大学 一种基于函数基的三维水汽探测方法
CN108229537A (zh) * 2017-12-07 2018-06-29 深圳市宏电技术股份有限公司 一种降水量的奇异值检测方法、装置及设备
CN108387169A (zh) * 2018-02-11 2018-08-10 羲和时空（武汉）网络科技有限公司 一种基于实时大气产品的gnss形变监测系统
CN108761574A (zh) * 2018-05-07 2018-11-06 中国电建集团北京勘测设计研究院有限公司 基于多源信息融合的降雨量估算方法
CN109344865A (zh) * 2018-08-24 2019-02-15 山东省环境规划研究院 一种多数据源的数据融合方法
CN109902346A (zh) * 2019-01-24 2019-06-18 东南大学 基于神经网络的区域加权平均温度信息获取方法

Patent Citations (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20090189802A1 (en) * 2008-01-25 2009-07-30 Tillotson Brian J System and method for using iridium satellite signals for meteorological measurements
JP2010060443A (ja) * 2008-09-04 2010-03-18 Japan Weather Association 気象予測装置、方法及びプログラム
CN103323888A (zh) * 2013-04-24 2013-09-25 东南大学 Gnss大气探测数据中对流层延迟误差的消除方法
US20170045624A1 (en) * 2015-08-14 2017-02-16 Trimble Navigation Limited Navigation satellite system positioning involving the generation of advanced correction information
CN106324620A (zh) * 2016-08-02 2017-01-11 中国人民解放军空军工程大学 一种不依赖地表气象数据实时测量的对流层天顶延迟方法
CN107843943A (zh) * 2017-10-30 2018-03-27 西安科技大学 一种基于函数基的三维水汽探测方法
CN108229537A (zh) * 2017-12-07 2018-06-29 深圳市宏电技术股份有限公司 一种降水量的奇异值检测方法、装置及设备
CN108387169A (zh) * 2018-02-11 2018-08-10 羲和时空（武汉）网络科技有限公司 一种基于实时大气产品的gnss形变监测系统
CN108761574A (zh) * 2018-05-07 2018-11-06 中国电建集团北京勘测设计研究院有限公司 基于多源信息融合的降雨量估算方法
CN109344865A (zh) * 2018-08-24 2019-02-15 山东省环境规划研究院 一种多数据源的数据融合方法
CN109902346A (zh) * 2019-01-24 2019-06-18 东南大学 基于神经网络的区域加权平均温度信息获取方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
FEI YANG ET AL；: "《Establishment and Assessment of a New GNSS PrecipitableWater Vapor Interpolation Scheme Based on the GPT2w Model 》", 《REMOTE SENSING》 *
李松青: "《基于 GPT2w 模型的区域大气可降水汽研究》", 《硕士学位论文库》 *

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114019585A (zh) * 2021-10-11 2022-02-08 武汉大学 一种大高差地区高精度定位cors网fkp解算方法
CN114019585B (zh) * 2021-10-11 2024-06-11 武汉大学 一种大高差地区高精度定位cors网fkp解算方法
CN115993668A (zh) * 2023-03-22 2023-04-21 成都云智北斗科技有限公司 一种基于多项式改正和神经网络的pwv重建方法及系统
CN117405968A (zh) * 2023-10-30 2024-01-16 广东海洋大学 基于大数据的设备能耗检测方法和系统
CN117405968B (zh) * 2023-10-30 2024-05-31 广东海洋大学 基于大数据的设备能耗检测方法和系统

Also Published As

Publication number Publication date
CN111126466B (zh) 2020-09-22

Similar Documents

Publication Publication Date Title
CN114518586B (zh) 一种基于球谐展开的gnss精密单点定位方法
CN111126466B (zh) 一种多源pwv数据融合方法
CN109543353B (zh) 三维水汽反演方法、装置、设备和计算机可读存储介质
CN103323888B (zh) Gnss大气探测数据中对流层延迟误差的消除方法
CN110031877B (zh) 一种基于grnn模型的区域nwp对流层延迟改正方法
Wilgan et al. Multi-observation meteorological and GNSS data comparison with Numerical Weather Prediction model
CN104965207B (zh) 一种区域对流层天顶延迟的获取方法
CN112034490B (zh) 一种nwp反演对流层延迟改进方法
CN107843943B (zh) 一种基于函数基的三维水汽探测方法
CN109917424B (zh) 多因子约束下的nwp反演对流层延迟的残差改正方法
CN109613582B (zh) 一种车载实时单频米级伪距定位方法
CN116182795B (zh) 普速铁路纵断面精密测量方法
Xia et al. Assessing water vapor tomography in Hong Kong with improved vertical and horizontal constraints
CN114910939B (zh) 短距离大高差rtk中对流层延迟实测气象改正方法
CN111123345B (zh) 一种基于gnss测量的经验电离层模型数据驱动方法
CN114357770B (zh) 一种对流层层析方法
CN116609859A (zh) 一种气象灾害高分辨率区域模式预报系统及方法
CN106908815B (zh) 一种基于探空数据的北半球对流层延迟改正方法
CN115097490A (zh) 耦合gnss及地面气象站监测数据的实时水汽场生成方法
CN112711022B (zh) 一种GNSS层析技术辅助的InSAR大气延迟改正方法
Raquet et al. Use of a Covariance Analysis Technique for Predicting Performance of Regional‐Area Differential Code and Carrier‐Phase Networks
CN115980317B (zh) 基于修正相位的地基gnss-r数据土壤水分估算方法
Seko et al. The meso-γ scale water vapor distribution associated with a thunderstorm calculated from a dense network of GPS receivers
Zhang et al. A new four-layer inverse scale height grid model of china for zenith tropospheric delay correction
CN115346128A (zh) 一种光学立体卫星dem高程改正和融合方法

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
CF01 Termination of patent right due to non-payment of annual fee
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20200922

Termination date: 20211216