Tm 是 GNSS 反演 PWV 的关键性因素,由 Tm 可以求得转换因素 Π,Π*ZWD(天顶湿延迟)可以的得到 PWV。
Tm 的计算方法有两种,下面进行分别介绍
使用水汽压和温度计算 Tm
和 表示第 i 层大气的水汽压和温度, 是第 i 层的厚度。这种方法使用探空站数据计算得到的 Tm 一般为真值与其他数据进行比较。
使用最多一般是 Bevis 和 Yao 的线性公式,不一一列举了,这个精度较差。
我一般使用 ERA5 数据集根据第一种方法得到 Tm。以探空站求得 Tm 为真值,进行对比。ERA5 数据集覆盖率高,时间和空间分辨率也高,很适合用来计算 Tm。只需要求得水汽压,带入第一种方法的公式就可以得到 Tm。
根据饱和水汽压和相对湿度求得,公式如下
式中,T 为温度(K)。es 采用 ECMWF IFS 报告(IFS Documentation CY31R1 Part II )给出的模型,对水的不同状态做了区别 (ECMWF, 2007) :
(1)温度大于 0℃, R2 = 611.21 hPa,R3 = 17.502 K 和 R4 = 32.19 K;
(2)温度小于-23℃,R2 = 611.21 hPa,R3 = 22.587 K 和 R4 = -0.7 K;
(3)温度介于-23℃ 和 0℃ 之间,则用下式计算:
式中,T0 = 273.16 K,Ti = 250.16 K。
%?GNSS上方的ERA5数据集的求得Tm
%?根据ERA5最底层高于或者低于GNSS站点,使用补偿方法得到
%?前置数据需要下载ERA5的位势、温度、相对湿度。
%?更多代码获取微信公众号WZZHHH
%?-----------------------------------------------------------------------
%?你也可以把数据下载好打包发我,我来负责处理,不过费用贵
%?需要ERA5、站点经纬度名称海拔的excel表等,具体需要商量
%%?--------------------------需要修改的路径-----------------------------
clc;clear;
%?计算的年份以2020年为例
nc_path????=?'D:\paper_write\paper_code\2\ERA5\2020\';?%?原始ERA5的存放地址
Edata_path?=?'D:\paper_write\paper_code\2\mat\ERA5\';???%?ERA5读取后数据存放
%?Excel第一列站点名字、第二列纬度、第三列经度、第四列高程
%?station_x(:,1)=纬度、station_x(:,2)=经度、station_x(:,3)=高程
[station_x,~]?=?xlsread('D:\paper_write\paper_code\2\xls\使用站点的经纬度高程.xlsx');?%?GNSS站点经纬度坐标
resolution?=?0.25;????????%?ERA5分辨率0.25;ERAinterm分辨率0.75
%%?----------------------------NC数据读取-------------------------------
%?读取原始ERA5文件夹下所有nc数据
List?=?dir(fullfile(nc_path,'*.nc'));
Tm=[];
for?I?=?1:size(List,1)
????%?nc的具体路径
????filen?=?[nc_path?List(I).name];
????%?用ncinfo读取nc数据里面的元素,找到我们下载数据的缩写
????%?aa?=?ncinfo(filen);
????%?例如:z:位势、t:温度、r:相对湿度、q:比湿度、levels:气压
????z?=?ncread(filen,'z');
????t?=?ncread(filen,'t');
????r?=?ncread(filen,'r');
????levels?=?ncread(filen,'level');
????Time_num?=?ncread(filen,'time');
????latitude?=?ncread(filen,'latitude');
????longitude?=?ncread(filen,'longitude');