基于GLDAS_NOAH_025_3H土壤湿度数据计算LGD

1.系统概述与科学背景

传统的重力反演方法通常生成月平均的全球重力场模型（Level-2 球谐系数或 Level-3 Mascon产品）。然而，这种时空平均化处理往往会模糊甚至消除高频的地球物理信号，例如短期的洪水事件、海啸传播或快速变化的洋流系统。为了充分利用GRACE-FO搭载的激光测距干涉仪（Laser Ranging Interferometer, LRI）所提供的前所未有的纳米级测量精度，本程序采用了一种“正演建模（Forward Modeling）”的方法。它利用高时空分辨率的水文模型（如GLDAS）作为输入，通过积分格林函数法，直接计算卫星轨道高度上的瞬时重力摄动。

本程序旨在实现 Ghobadi-Far et al. (2022) 在 Journal of Geophysical Research: Solid Earth 发表的论文 “Along-Orbit Analysis of GRACE Follow-On Inter-Satellite Laser Ranging Measurements for Sub-Monthly Surface Mass Variations” 中的部分核心算法。
主要功能是利用 GLDAS (Global Land Data Assimilation System) 水文模型数据（土壤湿度），计算沿 GRACE-FO 卫星轨道的视线方向重力差 (Line-of-Sight Gravity Difference, LGD) 模拟值。这可以用于验证 GRACE-FO LRI (Laser Ranging Interferometer) 观测到的高频（亚月级）质量变化信号。

核心功能

数据读取: 处理 GLDAS NOAH 0.25度分辨率的土壤湿度数据。
LGD 计算: 基于牛顿积分法（Point-mass integration）和球谐函数理论，将地表质量变化 (SMC) 转换为卫星高度处的 LGD 信号。
轨道分析: 筛选特定区域（如亚马逊流域、孟加拉湾等）的卫星过境轨迹。
可视化: 绘制不同日期的 LGD 沿轨变化图。

Github仓库： https://github.com/SingyuTang/GLDAS2LGD

2. 理论背景与原理

2.1 什么是 LGD?

传统的 GRACE/GRACE-FO 数据处理通常生成月平均重力场模型（Level-2 产品），这会平滑掉亚月级（Sub-monthly）的快速变化信号（如洪水、风暴潮）。
LGD (Line-of-Sight Gravity Difference) 定义为两颗卫星位置处的重力矢量差在视线方向（Line-of-Sight）上的投影。它是直接反映瞬时质量变化的观测量，能够捕捉高频信号。

设和分别为GRACE-FO C星（后随卫星）和D星（前导卫星）在地球中心地球固定坐标系（ECEF）中的位置矢量。星间视线方向的单位矢量定义为：

$\mathbf{e}_{12}^{LOS} = \frac{\mathbf{r}_2 - \mathbf{r}_1}{|\mathbf{r}_2 - \mathbf{r}_1|}$

设为由地表质量异常引起的重力扰动矢量。则两颗卫星处的重力扰动分别为和。LGD 计算公式为： $LGD = (\delta \mathbf{g}_2 - \delta \mathbf{g}_1) \cdot \mathbf{e}_{12}^{LOS}$

2.2 表面质量变化 (SMC) 与负荷效应

程序通过处理GLDAS水文模型数据来模拟地表质量变化。在物理上，土壤湿度、积雪、地表水等被视为附着在地球表面的薄层质量负载（Mass Load，计算使用去背景场后的质量异常）。这种负载会产生两种引力效应：

直接引力效应 (Direct Effect): 水体本身的质量异常对卫星产生的万有引力。

计算某一时刻的SMC相对于背景场（如月平均值）的偏差：

$\text{Direct Mass Potential}=\text{Current SMC−Mean SMC}$

弹性负荷效应 (Elastic Loading Effect): 地球并非刚体。地表水体的重量会压迫地壳，导致固体地球发生形变。这种形变不仅改变了地表的几何形状（引起观测点位置变化，虽然对卫星重力影响较小但物理上存在），更重要的是它重新分布了地球内部的质量（地幔物质的移动），从而产生了附加的重力场变化。

负荷Love数 (Load Love Numbers):为了描述这种弹性响应，引入了负荷Love数。对于球谐展开的第阶，总的引力位势变化由下式给出： $\delta V_{lm} \propto (1 + k_l) \times \text{Direct Mass Potential}$ 其中：代表直接的牛顿引力效应。代表固体地球变形产生的附加引力效应。

2.3 积分格林函数法 (Integral Green’s Function Approach)

虽然球谐函数法在处理全球尺度问题时效率很高，但在处理高频、局部化的水文信号（如亚马逊洪水）并将其映射到具体的卫星轨道点时，空间域的数值积分方法（格林函数法）更为直观且易于处理局部窗口。本程序采用的是空间域积分方法。对于轨道上的每一个点，程序计算其周围一定半径（cutoff_deg）内的所有地表质量单元对该点的引力贡献。引力矢量的计算核心 (_calculate_gravity_contribution):代码中并未直接调用球谐系数库，而是构建了一个物理内核，针对每一个质量网格点和卫星位置，计算引力矢量。在局部东北天（North-East-Up, NEU）坐标系下，由质量引起的引力分量可以通过包含Love数的格林函数卷积得到。数学上，这等价于对球谐级数的求和公式进行变换，得到关于球面距离的闭合解析式或截断级数形式的核函数：

$\delta \mathbf{g}_{NEU} = \sum_{j \in \text{Region}} \mathcal{G}(\psi_{ij}, k_l) \cdot dm_j$

其中是卫星与地面质量点之间的球面角距离。由于引力随距离的平方衰减，且 LGD 对高频信号敏感，为了提高计算效率，代码引入了截断距离 cutoff_deg（默认25度）。这意味着只计算卫星星下点周围约 2700 公里范围内的质量贡献，忽略地球背面的微小影响。这在保证精度的同时极大地加速了计算。

2.4 坐标系转换 (Coordinate Transformation)

地球物理模型（如GLDAS）的数据是在地理坐标系（经度、纬度）下定义的，其产生的引力矢量自然表达在局部水平坐标系（NEU: North, East, Up）中。然而，卫星轨道和视线矢量是在 ECEF（地心地固直角坐标系）中定义的。为了进行点积运算，必须将 NEU 坐标系下的重力矢量转换到 ECEF 坐标系。gldas2lgd.py 中的 neu_to_ecef 函数实现了这一旋转变换： $\mathbf{\delta}g_{ECEF} = \mathbf{R} \cdot \mathbf{\delta}g_{NEU}$ 旋转矩阵由卫星的地心纬度和经度决定： $\mathbf{R} = \begin{bmatrix} -\sin\phi \cos\lambda & -\sin\lambda & \cos\phi \cos\lambda \\ -\sin\phi \sin\lambda & \cos\lambda & \cos\phi \sin\lambda \\ \cos\phi & 0 & \sin\phi \end{bmatrix}$ 这一步是连接地面水文模型与空间观测几何的关键环节。

3. 环境依赖

本程序基于 Python 3 开发，依赖以下库：

Numpy: 数值计算
Scipy: 科学计算 (用于勒让德多项式 lpn)
NetCDF4: 读取 GLDAS .nc4 文件
Matplotlib: 绘图
Cartopy: 地图投影与绘制
自定义依赖 (需确保主要脚本同级目录下存在以下模块):
— S02compute_grace_lgd: 用于加载卫星轨道 (OrbitLoader)。
— S05plot_lgd_ra_cwt_filter: 用于筛选特定区域的轨道 (filter_complete_tracks_passing_region)。
— read_love_numbers: 用于读取负荷勒夫数。

依赖S02compute_grace_lgd和S05plot_lgd_ra_cwt_filter可到仓库进行下载。

安装基础依赖：

pip install numpy scipy netCDF4 matplotlib cartopy

4. 数据准备

本系统的运行高度依赖于外部数据输入。用户需要准备两大类数据：地表质量变化数据（GLDAS）和卫星轨道数据（GRACE-FO GNV1B）。

4.1 GLDAS 水文数据

全球陆地数据同化系统（GLDAS）提供了高时空分辨率的陆地表面状态变量。本程序默认配置使用的是 GLDAS Noah 模型 Level 4 产品。

数据产品名称: GLDAS Noah Land Surface Model L4 3 hourly 0.25 x 0.25 degree V2.1
短名称: GLDAS_NOAH025_3H
版本: 2.1
空间分辨率:
时间分辨率: 3小时 (3-hourly)
文件格式: NetCDF4 (.nc4)
数据提供方: NASA GES DISC

下载地址：NASA Earthdata (GES DISC)
文件格式： GLDAS_CLSM025_DA1_D.A20200501.022.nc4
存放结构： 建议按年或月存放，I:\LGD\GLDAS_CLSM025_D_2.2_2020

关键变量说明：GLDAS模型将土壤水分分为四层。为了获得总的陆地水储量变化（TWS），gldas2lgd.py 会读取并累加以下四个变量：

SoilMoi0_10cm_inst: 0-10 cm 土层土壤水分 ()
SoilMoi10_40cm_inst: 10-40 cm 土层土壤水分 ()
SoilMoi40_100cm_inst: 40-100 cm 土层土壤水分 ()
SoilMoi100_200cm_inst: 100-200 cm 土层土壤水分 ()

4.2 GRACE-FO 轨道数据

程序需要 GRACE-FO 的精密轨道数据（GNV1B）来确定卫星在特定时间的位置。代码中通过 OrbitLoader 类加载 GROOPS 工作区数据，你需要根据实际情况修改 load_orbit 函数以适配你的轨道文件格式（如 GNV1B等）。

如果按照博客详细介绍利用GRACE1B多日数据计算LGD工作流程二_基于LRI1B多日数据已经创建了工作根目录workdir并进行了相关处理，只需要将本项目文件全部拷贝到该目录下即可运行，因为本项目中的 load_orbit 函数默认读取workdir/gracefo_dataset路径下的GNV1B轨道文件。

4.3 勒夫数 (Love Numbers)

用于计算负荷效应。需准备 PREM 模型或其他模型的勒夫数文件，并通过 read_love_numbers 模块读取。本程序已提供勒夫数文件和读取函数。

6. 使用指南

第一步：配置路径与参数

打开 lgd_processor.py，修改 main 函数中的 CONFIG 字典与待处理日期列表 date_list：

CONFIG = {
        # 路径设置
        'gldas_dir': r"I:\LGD\GLDAS_NOAH025_3H_2.1_2020",
        'groops_workspace': r'G:\GROOPS\PNAS2020Workspace',
        'output_dir': r"./results_lgd",

        # 背景场设置 (用于计算距平)
        'bgd_start': '2020-05-01',
        'bgd_end': '2020-05-31',

        # 区域设置
        'region_lon': (88, 92),
        'region_lat': (22, 26),
        'track_lat_limit': (-80.0, 80.0),

        # 轨道与计算参数
        'orbit_direction': 'asc',  # 'asc' or 'desc'
        'orbit_interval': 5,  # 秒
        'cutoff_deg': 20.0,  # 积分截断距离
        'n_max': 200,  # 阶数

        # 绘图参数
        'results_dir': r"./results_lgd",  # 上一步代码保存结果的文件夹
        'track_lat_limit': (-80.0, 80.0)  # 绘图时的纬度截断范围
    }

date_list = [
        '2020-06-04', '2020-06-10', '2020-06-15', '2020-06-21', '2020-06-26',
        '2020-07-02', '2020-07-07', '2020-07-13', '2020-07-18', '2020-07-24', '2020-07-29',
        '2020-08-04', '2020-08-09', '2020-08-15', '2020-08-20'
    ]

第二步：运行计算程序

运行 lgd_processor.py。程序将执行以下操作：

读取 GLDAS 数据计算月平均背景场。
针对 date_list 中的每一天，提取轨道并计算瞬时 LGD。
打印进度并在 output_dir 生成 .npz 文件。

1 2	`python lgd_processor.py`

第三步：运行绘图程序

计算完成后，打开 lgd_plot.py，确保 CONFIG 中的 results_dir 指向刚才的输出目录，然后运行：

1 2	`python lgd_plot.py`

程序将生成 PNG 图片，展示目标区域内 LGD 随时间的变化。

8. 参考文章

Ghobadi-Far, Khosro, Shin-Chan Han, Christopher M. McCullough, David N. Wiese, Richard D. Ray, Jeanne Sauber, Linus Shihora, and Henryk Dobslaw. 2022. “Along-Orbit Analysis of GRACE Follow-On Inter-Satellite Laser Ranging Measurements for Sub-Monthly Surface Mass Variations.” Journal of Geophysical Research: Solid Earth 127(2):e2021JB022983. doi:10.1029/2021JB022983.

博客文章：详细介绍利用GRACE1B多日数据计算LGD工作流程二_基于LRI1B多日数据

学习

#GRACE #GRACE-FO #Python #LGD #GLDAS #NOAH #Soil Moisture

基于GLDAS_NOAH_025_3H土壤湿度数据计算LGD

https://singyutang.github.io/2025/12/10/基于GLDAS-NOAH-025-3H土壤湿度数据和GLDAS-CLSM-D数据计算LGD/

作者

SingyuTang

发布于

2025年12月10日

许可协议

基于ITSG日解数据的GRACE-FO沿轨视线重力差（LGD）反演与分析上一篇

详细介绍利用GRACE1B多日数据计算LGD工作流程二_基于KBR1B多日数据下一篇