根据ERA5数据解算站点高度的PWV

1.文中适用精度对比GNSS PWV，采用ERA5相对湿度、气压、温度以积分方法得到PWV。

2.也可用于ERA5插值温度、气压的精度对比。

2.由于GNSS与ERA5并不并址，需要进行水平和垂直方向的补偿。

3.水平方向补偿：在最接近GPS站经纬度坐标的4个网格点，利用反距离加权水平内插到每个水平压力级的GPS站位置  ，得到GPS测站点垂直剖面数据。

4.接下来，在垂直方向上将气象数据内插（外推）到每个GPS天线位置上。

5.下面是部分代码：

% 利用ERA5气压层数据积分计算PWV
%ERA5数据使用位势、温度、相对湿度
clc;clear;


% 加载站点的经纬度
[station_x,~] = xlsread('F:\predict_pm\xls\air_surface.xlsx');

List = dir('F:\predict_pm\ERA5\bj\allyear\leves\');%气压层数据文件夹
filename = 'F:\predict_pm\result\ERA5_data\21年数据.mat';%生成文件存放位置

% 经纬度边界（左上角点的位置）ERA5边界,ERA5分辨率0.25.
lon_t = 114.5;  lat_t = 41.25;
resolution = 0.25;


ERA_PR=[];ERA_T=[];pwvera5=[];esat=[];
for I = 3:size(List,1)

    % z：位势、t：温度、r：相对湿度
    filen = fullfile(List(I).folder,List(I).name);
    year = str2num(List(I).name(1:4));
    a = ncinfo(filen);

    %读取时间
    Time_num = ncread(filen,'time');
    times = length(Time_num);

    %读取经纬度
    longitude = ncread(filen,'longitude');
    latitude = ncread(filen,'latitude');

    %判断边界是否正确
    if longitude(1,1)~=lon_t||latitude(1,1)~=lat_t
        warning('边界出错！！！');
    end

    % 读取位势、温度、相对湿度、气压层
    z = ncread(filen,'z');
    t = ncread(filen,'t');
    r = ncread(filen,'r');
    levels = ncread(filen,'level');

    % 提前分配空间
    pwera = zeros(times,size(station_x,1));
    ERAPR = zeros(times,size(station_x,1));
    ERAT = zeros(times,size(station_x,1));
    rh = zeros(times,size(station_x,1));
    es = zeros(times,size(station_x,1));

    % 根据站点数量分别计算
    for P = 1:size(station_x,1)

        % 读取站点经纬度、海拔
        lon = station_x(P,1);
        lat = station_x(P,2);
        Height = station_x(P,3);