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

1.项目简介

本项目是对Ghobadi-Far等人（2022）提出的沿轨分析方法的工程化实现与扩展，其核心创新在于利用奥地利格拉茨技术大学（TU Graz）发布的ITSG-Grace2018高频日解重力场模型，替代传统的月平均解，从而在亚月（Sub-monthly）时间尺度上捕捉瞬时地球物理质量迁移信号 。

GRACE（Gravity Recovery and Climate Experiment）及其后续任务GRACE-FO（Follow-On）通过测量两颗低轨卫星之间的微米级（K波段微波测距，KBR）乃至纳米级（激光干涉测距，LRI）距离变化，为全球水储量、冰盖消融及固体地球形变提供了前所未有的观测窗口。然而，传统的重力场反演通常以“月”为时间单位，生成球谐系数（Level-2产品）或质量浓度（Mascon，Level-3产品）。这种时间平均化处理虽然有效降低了观测噪声，但也不可避免地造成了时间混叠（Aliasing）效应，导致洪水爆发、风暴潮引发的海洋质量重新分布等高频瞬变信号被平滑甚至丢失 。

本项目通过引入ITSG-Grace2018日解数据，结合严格的轨道动力学正演模拟，构建了一套从球谐系数空间到卫星视线方向加速度空间的转换体系。系统包含正演建模、信号频域滤波、轨迹筛选及交叉点分析等模块，旨在验证高频重力场模型的精度，并为水文学、海洋学及冰冻圈科学提供更高时间分辨率的观测数据支持。

GitHub仓库：https://github.com/SingyuTang/SHC2LGD

2. 科学原理与理论框架

2.1 卫星重力测量的基本观测量

GRACE-FO任务的核心在于测量两颗卫星（GRACE-C与GRACE-D）之间的距离变化率（Range Rate, ）和距离变率的变化率（Range Acceleration, ）。在理想情况下，这一物理量直接反映了卫星所受引力势差沿视线方向的投影。然而，传统的处理流程倾向于将这些沿轨观测值通过复杂的最小二乘平差反演为全球均匀格网或球谐系数。这一过程是一个不仅病态（Ill-posed）而且高度依赖正则化的数学逆问题。本项目采用“正演法”（Forward Modeling），即直接利用已有的重力场模型计算理论上的卫星受力，并将其与实际观测进行对比。这种方法避免了逆问题的数值不稳定性，使得对特定区域、特定时刻的重力异常分析更加直观和精确。

2.2 视线重力差（LGD）的数学定义

视线重力差（LGD）是连接全球重力模型与局部卫星观测的关键物理量。它定义为两颗卫星所在位置的引力矢量差在视线方向上的投影。设 和 分别为时刻 两颗卫星在地固系（ECEF）中的位置矢量， 为该位置的地球引力加速度矢量。两卫星间的视线单位矢量 定义为：

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

此时，LGD标量 可表示为：$\delta g_{12}^{LOS} = (\mathbf{g}(\mathbf{r}_2) - \mathbf{g}(\mathbf{r}_1)) \cdot \mathbf{e}_{12}$

该物理量的单位通常为 。LGD之所以重要，是因为它与卫星搭载的LRI激光干涉仪直接测量的距离变率 高度相关。Ghobadi-Far等人（2018）的研究表明，在频率高于1 mHz的频段内，LGD与经过非保守力修正后的距离变率残差几乎完全线性相关。因此，计算合成LGD是验证重力场模型是否准确捕捉了高频信号的最佳手段。

2.3 基于球谐函数的引力梯度计算

本项目中的核心算法 sh2lgd.py 实现了从球谐系数（Spherical Harmonic Coefficients, SHC）到引力矢量的转换。地球引力势 在球坐标系 下展开为：$V(r, \theta, \lambda) = \frac{GM}{R} \sum_{l=0}^{L} \left(\frac{R}{r}\right)^{l+1} \sum_{m=0}^{l} [C_{lm} cos(m\lambda) + S_{lm} sin(m\lambda)]\bar{P}_{lm}(\cos\theta)$
其中： 为地心引力常数。 为地球参考半径。 为完全归一化的勒让德函数。 为球谐系数，包含了地球质量分布的信息。

为了得到加速度 ，必须计算引力势在直角坐标系下的偏导数。这涉及到勒让德函数及其导数的复杂递归计算。代码中采用了Cunningham或Pines方法的变体（具体体现为 _recursive_EF_vectorized 函数），通过递归关系直接计算梯度分量，极大地提高了在大规模时间序列上的计算效率。

注意区分【Ghobadi-Far2022】中的公式（1）描述的为地球引力势变化在球坐标系下的展开式

$$\delta V(r,\varphi,\lambda) = \frac{GM}{R} \sum_{l=0}^{L} \left(\frac{R}{r}\right)^{l+1} \sum_{m=0}^{l} [\delta \bar{C}_{lm} cos(m\lambda) + \delta \bar{S}_{lm} sin(m\lambda)] \overline{P}_{lm}(\sin\varphi)$$

其中：：完全归一化的异常球谐系数（即随时间变化的重力场系数，也就是L2产品或ITSG产品中的球协系数）。：完全归一化的勒让德函数（Associated Legendre Function）。

2.4 从月解到日解的范式转移

标准的GRACE Level-2产品是月平均解。这意味着，如果在一个月内亚马逊流域发生了一次持续3天的洪水脉冲，或者阿根廷海盆出现了一个周期为20天的海洋涡旋，这些信号在月平均场中会被平均化，导致幅值衰减甚至信号混淆（Aliasing）。

本项目特异性地采用了 ITSG-Grace2018 日解数据。ITSG（格拉茨技术大学大地测量研究所）利用卡尔曼滤波（Kalman Filter）技术处理GRACE/GRACE-FO数据。

• 正则化策略：由于单日的卫星覆盖极其稀疏（无法覆盖全球），直接解算每日重力场是不可能的。ITSG通过引入时间相关性约束（通常是三阶自回归过程 AR-3），利用前后的观测信息来约束当天的解 。

• 信号保留：这种平滑处理在抑制噪声的同时，能够保留约500公里尺度以上的亚月变化信号。

• 截断阶数：由于高阶项受噪声影响严重，ITSG日解通常截断在40阶（Degree 40），这在本项目代码配置 MAX_DEGREE = 40 中得到了体现 。

3. 数据架构与输入规范

本系统的运行依赖于特定格式的输入数据，理解这些数据的物理意义与结构是正确使用软件的前提。

3.1 动态重力场模型：ITSG-Grace2018

• 数据来源：ITSG (Institute of Geodesy at Graz University of Technology) 。
• 文件格式：.gfc (Gravity Field Coefficient)。这是一种标准的ASCII格式，包含头文件（元数据）和数据块。
• 卡尔曼滤波日解（Kalman Smoothed Daily Solutions）：本项目的核心输入，最高40阶。
• 命名规范：代码 sh2lgd.py 期望的文件名模式需严格匹配，例如 ITSG-Grace2018_n40_YYYY-MM-DD.gfc。
• 物理含义：这些系数代表了相对于某个静态场的全量重力场或异常重力场（取决于具体发布版本，通常包含水文、冰冻圈及GIA信号，但大气海洋信号AOD1B已被扣除或需额外处理）。
• 下载链接：https://ftp.tugraz.at/outgoing/ITSG/GRACE/ITSG-Grace_operational/daily_kalman/daily_n40/

3.2 卫星精密轨道星历：GNV1B

为了计算沿轨LGD，必须知道卫星在每一时刻的精确位置。
数据来源：GRACE-FO Level-1B Navigation Data (GNV1B)。
坐标系：通常为地固系（ECEF, Earth-Centered Earth-Fixed），与重力场模型的旋转坐标系一致。
文件格式：通常为ASCII或二进制，包含GPS时间标签及位置速度矢量。
采样率：标准发布频率为1Hz，但重力场反演通常降采样至5秒（0.2Hz）。代码中的 OrbitLoader 类负责处理这种降采样与插值。

3.3 背景重力场模型

LGD通常表现为相对于静态地球重力场的“异常值”。

常用背景场模型：GGM05C（GRACE与GOCE结合模型）、GOCO06s（卫星组合模型）、重力场组合模型（如前96阶使用L2产品，97阶后使用GOCO06s）或Level-2月度重力场模型（如单独使用JPL某个月份的重力场模型）。

作用：在 sh2lgd.py 中，日解系数 会减去背景场系数 ，得到 。$\Delta C_{lm}(t) = C_{lm}^{daily}(t) - C_{lm}^{background}$这确保了计算出的LGD仅反映随时间变化的质量迁移（如水循环），而非地球本身巨大的静止质量引力。

3.4 辅助数据与依赖库

Python环境：

• numpy: 核心矩阵运算库，用于处理球谐递归算法中的张量运算。
• scipy: 用于读取和保存MATLAB格式文件（.mat）。matplotlib: 用于绘图。

MATLAB环境

• Signal Processing Toolbox: 提供 cwt, icwt 等小波分析函数。
• Wavelet Toolbox: 必须包含“Morse”小波支持。

4. 核心功能模块说明

本项目由5个核心脚本组成，涵盖了从数据计算、信号处理到可视化的全流程。以下是各文件的功能简介。

4.1 正演计算引擎：sh2lgd.py

这是整个数据处理流的起点。

• 主要功能：将抽象的球谐系数（SHC）转换为具体的沿轨重力加速度差（LGD）。

• 工作流程：

  1. 读取指定日期的ITSG日解文件和静态场文件，计算系数差值。

  2. 加载当天的卫星轨道数据（GNV1B），确定卫星位置。

  3. 正演计算：根据卫星位置和重力场系数，算出两颗卫星在视线方向上的引力差。

  4. 区域筛选：支持仅计算特定经纬度区域（如亚马逊流域）的数据以节省时间。

  5. 输出：生成包含时间戳和原始LGD数值的 .mat 文件。

4.2 信号处理与滤波：S03cwt_fliter_lgd_ra.m & process_cwt_fliter_signal_data.m

原始的LGD数据包含各种频率成分，这两个MATLAB脚本负责提取我们感兴趣的信号。

• S03cwt_fliter_lgd_ra.m (批处理驱动)：

  • 这是用户操作的主入口。

  • 负责定义处理的时间范围（如2020年6月-8月），自动遍历每一天的数据文件。

  • 调用核心处理函数，实现批量自动化作业。

• process_cwt_fliter_signal_data.m (算法实现)：

  • 功能：执行连续小波变换（CWT）。

  • 目的：利用小波的时频特性，提取特定频段（通常为 1-12 mHz）的信号。这一步能有效去除低频的轨道误差和高频的测量噪声，保留真实的地球物理信号（如水文质量变化）。

  • 输出：生成滤波后的 cwt_time-lgd-*.mat 文件。

4.3 单日结果可视化：plot_single_lgd.py

用于微观层面的数据验证。

• 功能：读取某一天的处理结果，绘制该日卫星轨道经过特定区域时的重力变化曲线。

• 用途：用户可以用它来对比“原始信号”与“滤波后信号”，检查滤波器是否工作正常，或者查看某次特定洪水事件在当天的信号强度。

4.4 时变序列分析：plot_multi_lgd_cross_over_area.py

该脚本用于宏观层面的时变分析。

• 瀑布图（Waterfall Plot）：通过在Y轴（纬度）和X轴（LGD值+偏移量）上堆叠多天的轨道数据，形成随时间演变的视觉效果。可以直观地看到重力异常随时间的演变过程。
• 物理意义：如果亚马逊河的洪水波向下游传播，或者土壤湿度随季节增加，用户将观察到LGD曲线的峰值在不同日期的轨道上发生位移或幅值增长。这是验证时变重力场模型捕捉动态过程能力的强有力工具。

5. 用户指南与操作流程

本指南假定用户具备基本的Linux或Windows命令行操作能力，并已安装Python (3.8+) 和 MATLAB (R2020a+)。

5.1 环境部署

• Python依赖安装：

pip install numpy scipy matplotlib

• MATLAB工具箱检查： 在MATLAB命令窗口输入 ver，确认已安装 “Signal Processing Toolbox” 和 “Wavelet Toolbox”。

5.2 数据准备

在运行代码前，请确保以下数据已下载并按目录结构存放：

1. ITSG Daily Solutions (.gfc):
   • 来源: ITSG-Grace_operational_Kalman_daily_n40
   • 存放路径示例: BASE_DIR/itsg_dataset/YYYY-MM/
2. Static Background Model:
   • 用于扣除平均场（如GOCO06s或GRACE-FO Level2 GSM产品）。
   • 存放路径示例: 'BASE_DIR\grace_products\GSM-2_2020122-2020152_GRFO_UTCSR_BB01_0603'
3. GRACE-FO Orbit Data (GNV1B):
   • 包含卫星的位置。代码中通过 OrbitLoader 类加载 GROOPS 工作区数据，你需要根据实际情况修改 load_orbit 函数以适配你的轨道文件格式（如 GNV1B等）。
   • 如果按照博客详细介绍利用GRACE1B多日数据计算LGD工作流程二_基于LRI1B多日数据已经创建了工作根目录workdir并进行了相关处理，只需要将本项目文件全部拷贝到该目录下即可运行，因为本项目中的 load_orbit 函数默认读取workdir/gracefo_dataset路径下的GNV1B轨道文件。

注：BASE_DIR为GROOPS工作目录，也就是上面所提到的workdir，如 BASE_DIR = r'G:\GROOPS\PNAS2020Workspace'

5.3 运行步骤

步骤一：LGD正演计算 (Python)

编辑 sh2lgd.py：

# 配置示例
BASE_DIR = r'G:\GROOPS\PNAS2020Workspace'		# 工作空间根目录
STATIC_FILE_REL = r'grace_products\GSM-2_2020122-2020152_GRFO_UTCSR_BB01_0603'		# 背景场文件 (相对于 BASE_DIR 或绝对路径)
START_DATE = '2020-06-01'	# 处理时间范围起始时间
END_DATE = '2020-08-30'		# 处理时间范围终止时间
MAX_DEGREE = 40  # 根据ITSG模型的阶数设定
REGION_CFG = {
        'lon': (88, 92),  # 感兴趣区域经度范围
        'lat': (22, 26),  # 感兴趣区域纬度范围
        'lat_limit': (-80.0, 80.0),  # 轨道纬度限制
        'dir': 'asc'  # 轨道方向 'asc' (升轨) 或 'desc' (降轨)
    }

预期输出：results/ 目录下生成一系列 time-lgd-YYYY-MM-DD.mat 文件。

步骤二：信号滤波 (MATLAB)

在 MATLAB 中打开并运行 S03cwt_fliter_lgd_ra.m。

配置：设置 start_datenum, end_datenum 和频率范围（f_high，例如 30 mHz）。
功能：脚本会自动读取上一步生成的 .mat 文件，应用 CWT 滤波保留高频/亚月信号。
输出：results/cwt_time-lgd-YYYY-MM-DD.mat
注意：这一步必须运行，但是结果不会使用滤波后的数据，因为作者偷懒没有修改结果可视化部分的代码（从其他项目复制过来的代码逻辑），所以包含本部分运行结果的的逻辑部分。

步骤三：结果可视化 (Python)

1. 多日序列区域分析 (plot_multi_lgd_cross_over_area.py)

用于展示特定区域（如孟加拉湾、亚马逊流域等）随时间的重力变化。

# 修改日期列表和区域范围
date_list = ['2020-06-04', '2020-06-10', ...]
lon_range = (88, 92)
lat_range = (22, 26)
lat_limit = (-80.0, 80.0)  # 轨道延申纬度范围
orbit_direction = 'asc'  # 轨道方向

运行：

python plot_multi_lgd_cross_over_area.py

2. 单日对比分析 (plot_single_lgd.py)

查看特定日期的原始信号与滤波后信号的差异。

date_str = '2020-07-07'
lon_range = (88, 92)  # 经度范围（度）
lat_range = (22, 26)  # 纬度范围（度）
lat_limit = (-80.0, 80.0)  # 轨道延申纬度范围
orbit_direction = 'asc'  # 轨道方向

运行：

python plot_single_lgd.py

6. 结果解释与科学意义

6.1 数据的物理量级

LGD的典型单位是 。

• 静态场异常：由山脉、海沟引起的静态重力异常在LGD中可达数千 。
• 时变信号：由水文（如亚马逊洪水）引起的信号通常在 量级 1。
• 噪声水平：LRI的仪器噪声极低（），但ITSG日解由于截断和混叠，其模型误差可能在 左右。因此，能够从模型中提取出连续、具有物理规律的 信号是模型有效的有力证明。

6.2 滤波器的作用

CWT滤波器的通带设置（如1-12 mHz）对应特定的空间波长。

• GRACE轨道速度 。
• 频率 与空间波长 的关系为 。
• （大陆尺度）。
• （大型流域或海盆尺度）。通过仅保留这一频段，用户实际上是在“聚焦”观察特定尺度的质量异常，排除了深部地幔对流（超长波）和测量噪声（短波）的干扰。

6.3 潜在应用领域

1. 极端水文事件监测：追踪洪水波峰在河流系统中的传播速度。
2. 海洋动力学：验证高频海洋模型（如AOD1B）中未被捕捉的涡旋或波动信号。
3. 冰川质量平衡：监测格陵兰岛或南极冰盖在融雪季节的快速质量损失。

7. 总结

项目通过融合ITSG-Grace2018的高频解算能力与严密的轨道正演算法，使得研究人员能够突破月平均场的限制，深入探索地球系统在亚月尺度上的动态变化。通过对LGD这一直接观测量的分析，本工具为理解气候变化驱动下的快速质量迁移提供了强有力的量化手段。

8. 参考文章

1.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.

2.Mayer-Gürr, Torsten; Behzadpur, Saniya; Ellmer, Matthias; Kvas, Andreas; Klinger, Beate; Strasser, Sebastian; Zehentner, Norbert (2018): ITSG-Grace2018 - Monthly, Daily and Static Gravity Field Solutions from GRACE. GFZ Data Services. doi.org: 10.5880/ICGEM.2018.003.

3.Ghobadi-Far, Khosro, Shin-Chan Han, Steven Weller, Bryant D. Loomis, Scott B. Luthcke, Torsten Mayer-Gürr, and Saniya Behzadpour. 2018. “A Transfer Function Between Line-of-Sight Gravity Difference and GRACE Intersatellite Ranging Data and an Application to Hydrological Surface Mass Variation.” Journal of Geophysical Research: Solid Earth 123(10):9186–9201. doi:10.1029/2018JB016088.