ERA5和GNSS站点不并址的处理方法二：垂直补偿（获得GNSS站点高度的PWV、温度和气压，基于matlab）

0. 码字不易，点赞加关注（公众号：WZZHHH，部分资料在公众号可以下载）

Part1 垂直补偿的原因

我们在获取GNSS站点位置的PWV有很多方法，其中一个是利用ERA5数据进行积分得到。由于ERA5格网位置和GNSS不并址，需要进行水平方向和垂直方向的补偿的补偿获得。这里讲的是垂直方向补偿，水平方向补偿已经写过反距离加权水平内插，附matlab代码(ERA5和GNSS站点不并址的处理方法之水平补偿)-CSDN博客，这里介绍垂直方向的补偿。

水平方向补偿方法

Part2 ERA5下载方法

ERA5多层温度、位势、相对湿度的下载方法在这面

Part3 具体方法

其中相对湿度取最接近层数据即可，温度和气压用下面方法

1. GNSS 站点高于最底层

使用内插估计站点参数即可。

2. GNSS 站点低于最底层

假定温度递减率为-6.5 K/km，外推获得温度。

ERA5数据中的气压所在的高度和测站高度通常不一致，这时候就需要对气压进行高程改正，用下面方法得到测站高度上的气压：

式中， Ps表示测站气压；P为最接近层的ERA5气压；g是重力加速度；Rd为干燥气体的气体常数；对于HM如下：

式中，hs表示测站高度； h为最接近层的 ERA5高度；T是最接近层的温度；Г为温度递减率-6.5 K/km。

Part4 ERA5 PWV的计算

插值得到测站处的气象参数温度、气压相对湿度后，利用公式求得PWV即可：

g是重力加速度，我考虑了纬度和高程变化。

a1=5.2885×10^-3；a2=-5.9×10^-6；

Re=6378.1km；Rp=Re-21.5km

Part5 Matlab部分代码如下（代码输出GNSS站点位置的ERA5 PWV）

这个代码包括了水平方向和垂直方向的补偿，获取方式在

详细代码在公众号WZZHHH获取

%?利用ERA5解算出GNSS站点位置出温度、气压、Tm以及PWV
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读取后数据存放
[station_x,~]?=?xlsread('D:\paper_write\paper_code\2\xls\使用站点的经纬度高程.xlsx');???????%?GNSS站点经纬度坐标

%%?----------------------------NC数据读取-------------------------------
%?读取原始ERA5文件夹下所有nc数据
List?=?dir(fullfile(nc_path,'*.nc'));

%?pwvera5：gnss站点积分得到的PWV
ERA_P=[];ERA_T=[];tmm_all=[];pwvera5=[];esat=[];
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');
????q?=?ncread(filen,'q');
????levels?=?ncread(filen,'level');
????Time_num?=?ncread(filen,'time');
????latitude?=?ncread(filen,'latitude');
????longitude?=?ncread(filen,'longitude');

????%?提前根据时间长度分配空间
????times?=?length(Time_num);
????Tm?=?zeros(times,size(station_x,1));?????%?插值得到的站点Tm
????pwera?=?zeros(times,size(station_x,1));??%?插值得到的站点pwv
????ERAT?=?zeros(times,size(station_x,1));???%?插值得到的站点温度
????ERAPR?=?zeros(times,size(station_x,1));??%?插值得到的站点气压
????es?=?zeros(times,size(station_x,1));?????%?插值得到的站点水汽压