ERA5和GNSS站点不并址的处理方法二:垂直补偿(获得GNSS站点高度的PWV、温度和气压,基于matlab)

0. 码字不易,点赞加关注(Matlab 获取代码可关注公众号WZZHHH回复:ERA5合集,或者咸鱼关注:WZZHHH123

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));     % 插值得到的站点水汽压
  • 12
    点赞
  • 5
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

WZZHHH

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值