Part1 Tm(代码获取方式在文章最后)
Tm 是 GNSS 反演 PWV 的关键性因素,由 Tm 可以求得转换因素 Π,Π*ZWD(天顶湿延迟)可以的得到 PWV。
Part 2Tm 的计算方法
Tm 的计算方法有两种,下面进行分别介绍
1根据加权平均温度的定义
使用水汽压和温度计算 Tm
和 表示第 i 层大气的水汽压和温度, 是第 i 层的厚度。这种方法使用探空站数据计算得到的 Tm 一般为真值与其他数据进行比较。
2利用线性关系
使用最多一般是 Bevis 和 Yao 的线性公式,不一一列举了,这个精度较差。
3还有其他方法,GPT2w 等,目前我没用到,后期用到我会写出来教程
Part3 ERA5 数据集解算 Tm
我一般使用 ERA5 数据集根据第一种方法得到 Tm。以探空站求得 Tm 为真值,进行对比。ERA5 数据集覆盖率高,时间和空间分辨率也高,很适合用来计算 Tm。只需要求得水汽压,带入第一种方法的公式就可以得到 Tm。
水汽压计算公式
根据饱和水汽压和相对湿度求得,公式如下
es 为饱和水汽压(hPa),用下式求得
式中,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。
Part4 Matlab 获取代码可关注公众号WZZHHH,或者咸鱼关注:WZZHHH123
Part5部分代码展示:
% 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');