学习教程

学习教程 | HGPT2模型计算大气参数

由 张 良承 没有评论

HGPT2（Hourly Global Pressure and Temperature 2）是一个基于 ERA5 再分析数据的全球模型，用于估算大气参数，包括表面压力（P）、表面空气温度（T）、相对湿度（RH）、加权平均温度（Tm）、天顶静力延迟（ZHD）、天顶湿延迟（ZWD）和大气可降水量（PWV）。该模型是 HGPT（Hourly Global Pressure and Temperature）模型的改进版，HGPT 模型最初于 2020 年提出，用于提供表面压力、温度、ZHD 和 Tm 等参数（参考文献 [1]）。HGPT2 在此基础上扩展了功能，特别引入了相对湿度的估算，使用 ERA5 的全空间分辨率（0.25° × 0.25°）和时间分辨率（1 小时），并考虑了年周期、半年周期和季度周期的振幅和相位变化，以更好地捕捉水汽动态（参考文献 [2]）。

HGPT2 模型的主要优势在于：

  • 高分辨率：利用 ERA5 的精细网格，优于早期模型如 GPT 系列（分辨率通常为 1° × 1° 或更粗）。
  • 实时/后处理适用：适用于 GNSS 气象学、InSAR 气象校正、水文模型、气候研究等领域，尤其在实时应用中，因为 ERA5 再分析数据无法用于即时预测（延迟 5 天发布）。
  • 气候变化考虑：引入线性趋势，以模拟全球气候变化场景。
  • 验证：模型使用 282 个湿度传感器和 35 个 IGS GNSS 站点进行验证，显示出比 GPT 系列更低的 RMSE 和偏差，尤其在高度组件上（参考文献 [2]）。

本教程以 MATLAB 为例，介绍如何使用 HGPT2 模型进行计算。其他语言（Python、Fortran）类似，可参考 GitHub 仓库。

参考文献：

  1. Mateus, P.; Catalão, J.; Mendes, V.B.; Nico, G. An ERA5-Based Hourly Global Pressure and Temperature (HGPT) Model. Remote Sens.2020, 12, 1098. https://doi.org/10.3390/rs12071098
  2. Mateus, P.; Mendes, V.B.; Plecha, S.M. HGPT2: An ERA5-Based Global Model to Estimate Relative Humidity. Remote Sens.2021, 13, 2179. https://doi.org/10.3390/rs13112179

准备工作

1. 下载必要文件

HGPT2 模型依赖四个二进制网格文件（.bin），这些文件包含模型系数，用于空间双线性插值计算。

  • 下载地址：https://github.com/pjmateus/hgpt_model/releases
    • 文件列表：
      • press_grid.bin：压力相关系数。
      • temp_grid.bin：温度相关系数。
      • rh_grid.bin：相对湿度相关系数（HGPT2 新增）。
      • tm_grid.bin：加权平均温度相关系数。

这些文件基于 20 年 ERA5 数据计算，年周期、半年周期和季度周期的振幅和相位存储在网格中（参考文献 [2] 中的 Section 2）。

2. 下载计算代码

  • GitHub 仓库：https://github.com/pjmateus/hgpt_model
    • 下载 MATLAB 版本的文件：hgpt2.m（主函数）和 es_wexler.m（用于计算饱和水汽压的 Wexler 方程辅助函数，参考文献 [2] 中的相对湿度计算）。
    • 其他版本：Python (hgpt2.py) 和 Fortran (hgpt2.f90) 可用类似方式使用。

3. 设置环境

  • 安装 MATLAB（推荐 R2020a 或更高版本）。
  • 将下载的四个 .bin 文件、hgpt2.m 和 es_wexler.m 放置在同一目录下。这是必需的，因为 hgpt2.m 会直接读取这些文件。
  • 如果使用其他语言，确保相应文件在工作目录中。

HGPT2 函数说明

输入参数

hgpt2.m 函数的调用格式：

[P, T, RH, Tm, ZHD, ZWD, PWV] = hgpt2(year, doy, hour, lat, lon, height, height_type)
  • year：整数，年份（e.g., 2024）。模型考虑线性趋势，因此年份影响气候变化校正。
  • doy：整数，年积日（Day of Year，1-366）。
  • hour：整数，小时（0-23）。
  • lat：浮点数，纬度（-90 到 90 度）。
  • lon：浮点数，经度（-180 到 180 度，或 0 到 360 度）。
  • height：浮点数，高度（米）。
  • height_type：字符串，高度类型：
    • ‘orth’：正高（orthometric height，相对于大地水准面的高度）。
    • ‘elli’：椭球高（ellipsoidal height，相对于椭球面的高度）。 模型会根据类型进行高度校正（参考文献 [1] 中的高度插值，文献 [2] 继承）。

输出参数

  • P：表面压力（hPa）。
  • T：表面空气温度（K）。
  • RH：相对湿度（%）。
  • Tm：加权平均温度（K，用于 ZWD 到 PWV 转换）。
  • ZHD：天顶静力延迟（m，使用 Saastamoinen 模型计算，参考文献 [1] Section 2.2）。
  • ZWD：天顶湿延迟（m，使用 Saastamoinen 湿公式，参考文献 [2] Section 2）。
  • PWV：大气可降水量（mm，使用 Bevis 公式转换，参考文献 [2] 中的 PWV 计算）。

函数内部使用双线性插值从 .bin 文件中获取系数，然后应用傅里叶分析计算周期性（年、半年、季度）和线性趋势（参考文献 [1] Section 2.1，文献 [2] 扩展到 RH）。

注意：

  • 输入可以是标量（单个点）或向量（批量计算，纬度/经度/高度等必须同维）。
  • 如果位置超出网格，函数会报错或返回 NaN。
  • RH 计算依赖 dew point temperature 的 ERA5 数据和 Wexler 方程（参考文献 [2] Section 2）。

示例：单个位置计算

假设计算 2024 年第 1 天 0 时的北京（纬度 39.9°，经度 116.4°，椭球高 50 m）参数。

% 示例：单个计算
year = 2024;
doy = 1;
hour = 0;
lat = 39.9;
lon = 116.4;
height = 50;
height_type = 'elli';

[P, T, RH, Tm, ZHD, ZWD, PWV] = hgpt2(year, doy, hour, lat, lon, height, height_type);

% 显示结果
fprintf('压力: %.2f hPa\n', P);
fprintf('温度: %.2f K\n', T);
fprintf('相对湿度: %.2f %%\n', RH);
fprintf('加权平均温度: %.2f K\n', Tm);
fprintf('ZHD: %.4f m\n', ZHD);
fprintf('ZWD: %.4f m\n', ZWD);
fprintf('PWV: %.2f mm\n', PWV);

输出示例（实际值取决于模型）：

压力: 1013.25 hPa
温度: 273.15 K
相对湿度: 50.00 %
加权平均温度: 280.00 K
ZHD: 2.3000 m
ZWD: 0.1000 m
PWV: 10.00 mm

示例：批量处理多个站点和时间

HGPT2 支持批量计算，常用于处理多个 GNSS 站点或时间序列。假设您有一个 Excel 文件 stations.xlsx，包含列：Station（站点名）、Lon（经度）、Lat（纬度）、Height（高度）、HeightType（’orth’ 或 ‘elli’）。

脚本示例：计算 2024 年全年的逐小时数据，并保存到 CSV。

% ========== 批量处理脚本 ==========

% 加载站点信息
stations = readtable('stations.xlsx');  % 假设列: Station, Lon, Lat, Height, HeightType

% 参数设置
year = 2024;  % 年份
start_doy = 1;  % 起始年积日
start_hour = 0;  % 起始小时 (0-23)
end_doy = 366;  % 结束年积日 (闰年)
end_hour = 23;  % 结束小时 (0-23)

% 初始化结果表
results = table();

% 循环每个站点
for i = 1:height(stations)
    station_name = stations.Station{i};
    lon = stations.Lon(i);
    lat = stations.Lat(i);
    hgt = stations.Height(i);
    hgt_type = stations.HeightType{i};
    
    % 循环时间
    for doy = start_doy:end_doy
        for hour = start_hour:end_hour
            % 调用 hgpt2
            [P, T, RH, Tm, ZHD, ZWD, PWV] = hgpt2(year, doy, hour, lat, lon, hgt, hgt_type);
            
            % 添加到结果表
            new_row = table({station_name}, year, doy, hour, lat, lon, hgt, {hgt_type}, P, T, RH, Tm, ZHD, ZWD, PWV, ...
                'VariableNames', {'Station', 'Year', 'DOY', 'Hour', 'Lat', 'Lon', 'Height', 'HeightType', ...
                'P', 'T', 'RH', 'Tm', 'ZHD', 'ZWD', 'PWV'});
            results = [results; new_row];
        end
    end
end

% 保存结果
writetable(results, 'hgpt2_results.csv');
disp('计算完成，结果保存到 hgpt2_results.csv');

解释：

  • 读取 Excel 文件，使用 readtable。
  • 双重循环：外层站点，内层时间（年积日 + 小时）。
  • 结果存储在表中，便于后续分析。
  • 对于大量站点/时间，考虑并行计算（使用 parfor 需要 Parallel Computing Toolbox）以加速。

注意事项和高级用法

  1. 时间范围：模型基于 20 年 ERA5 数据（2000-2020），但线性趋势允许外推到未来年份。对于历史数据，确保年份在合理范围内。
  2. 高度类型：选择正确类型。正高常见于地形数据，椭球高常见于 GNSS（参考文献 [1] Section 4）。
  3. 性能：单个调用快速，批量处理数千点可能需几分钟。网格文件 ~100 MB，内存友好。
  4. 验证：如文献 [2]，使用 radiosonde 数据验证 RH（RMSE ~5-10%）。对于 GNSS 处理，可输入到 GAMIT/GLOBK或PPP处理（参考文献 [2] Section 4）。
  5. 扩展：计算斜路径延迟需修改函数（参考 GitHub）。PWV 转换使用 Bevis 公式，精度依赖 Tm（参考文献 [1] Section 2.3）。
  6. 错误处理：如果 .bin 文件缺失，函数会报错。确保经纬度在 -180~180。
  7. 其他语言：Python 版本使用 NumPy 读取 .bin，类似调用 hgpt2(year, doy, hour, lat, lon, height, height_type)。
  8. 应用场景：GNSS 延迟校正（减少高度偏差）、PWV 估算（气候监测）、InSAR 大气校正（参考文献 [1] 引言）。

脚本获取地址

关注GNSS空间环境智能感知课题微信公众号，发送“HGPT2”，获取全部脚本下载地址。