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

1. 简介

文章详细介绍利用GRACE1B单日数据计算LGD工作流程二(基于LRI1B单日数据)介绍了如何基于GROOPS软件计算单日动力学轨道并结合单日的LRI数据计算单日RA和LGD，计算结果与PNAS论文GRACE Follow-On revealed Bangladesh was flooded early in the 2020 monsoon season due to premature soil saturation文末提供的LGD数据集基本吻合，但是与论文中Fig. 1结果不同，本文更新的代码对此进行了优化，能够与该篇论文的Fig. 1结果一致，并提供动力学轨道和LRI数据批量处理代码和工作流程。

1.1. GOCO06s静态重力场模型与月度Level2 产品结合

之前进行轨道积分使用的重力场模型为GOCO06s（0-200阶），对此进行更新，使用2019年GRACE-FO Level2月解产品计算的平均场作替换GOCO06s的前96阶，第97-200阶球协数据仍使用GOCO06s重力场模型。使用该方法计算的的重力场模型能够准确计算卫星动力学轨道，得出RA与LGD，详细代码与数据见附件代码中的程序A02~A04。

1.2. 利用GRACE LRI1B多日数据批量计算LGD

本次更新的代码可基于GROOPS软件、Matlab代码及Python代码进行动力学轨道计算与LGD（或RA）的批量计算，所有代码与数据见附件，工作流程在详细介绍利用GRACE1B单日数据计算LGD工作流程一(基于KBR1B单日数据)一文中已有详细介绍，这里结合新增的功能再次进行介绍。

2. 工作流程

工作目录（workdir）：G:\GROOPS\PNAS2020Workspace
将下载附件中的代码拷贝到工作目录中即可运行代码。

2.1 创建目录

在工作目录下创建四个文件夹，分别命名为input、output、image、gracefo_dataset如下，包含GROOPS处理脚本，目录结构如下。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
workdir/
├── input/
├── output/
├── image/
├── gracefo_dataset/
├── A00check_orbit_coverage.py
├── A01get_lri_kbr_records_num.py
├── A02read_grace_fo_shc_file.py
├── A03read_goco06s_file.py
├── A04combine_earth_gravity_field.py
├── process_cwt_filter_signal_data.m
├── S00GroopsIntegrateOrbitLRI-202007.xml
├── S01SatelliteTrackingFromLRI-202007.xml
├── S02compute_grace_lgd.py
├── S02compute_grace_lgd_gapfilled.py
├── S03cwt_filter_lgd_ra.m
├── S04lgd_ra_ipsd.py
├── S05plot_lgd_ra_cwt_filter.py
├── S05plot_lgd_ra_cwt_filter_gapfilled.py
├── S06compare_pnas_data.py
├── S07plot_lgd_ra_spatial_map.py
├── S08convert_twsgrid_format_from_fengwei_toolbox.m
├── S09plot_GRACE_L2_monthly_tws_from_matlab.py
├── S10plot_lgd_ra_multi_time_series_cross_over_area.py
└── main_hybrid_pipeline.py

2.2 数据准备

数据下载包含AOD1B RL06数据和GRACE-FO Level 1B RL04 .noLRI数据和.LRI数据。

AOD1B RL06：ftp://isdcftp.gfz.de/grace-fo/Level-1B/GFZ/AOD/RL06/2020/AOD1B_2020-06-02_X_06.asc.gz （通过filezilla软件下载）

GRACE-FO数据下载链接：https://podaac.jpl.nasa.gov/dataset/GRACEFO_L1B_ASCII_GRAV_JPL_RL04 或 https://search.earthdata.nasa.gov/search/granules?p=C2036882118-POCLOUD&pg[0][v]=f&tl=1743148786.317!3!!

GRACE-FO数据数据下载完成后，解压后将.noLRI 和 gracefo_1B_2020-07-07_RL04.ascii*.LRI所有文件夹拷贝到workdir/gracefo_dataset/下。AOD1B数据下载完成后，解压后将文件拷贝到workdir/gracefo_dataset/下。

*注意：AOD1B数据要在实际处理日期上多添加一天的文件。

2.3 程序运行

  1. 运行S00GroopsIntegrateOrbitLRI-202007.xml ，进行动力学轨道积分，参数修改如下
1
2
3
4
5
6
7
8
9
10
11
# global 模块参数修改
workspace → "G:\GROOPS\PNAS2020Workspace"
resampleInterval → 2
instrument → lri
background → ggm05b
timeStart → 58970 (2020-05-01 00:00:00)
timeEnd → 59091+23/24+59/1440+59/86400 (2020-08-30 23:59:59)
timeIntervalsDaily(uniformSampling)
├── timeStart → 58970 (2020-05-01 00:00:00)
├── timeEnd → 59091+23/24+59/1440+59/86400 (2020-08-30 23:59:59)
└── sampling → 2/86400

参数修改后在GROOPS运行S00GroopsIntegrateOrbitLRI-202007.xml脚本。

1.运行S01SatelliteTrackingFromLRI-202007.xml，进行LRI1B数据提取，参数修改如下

1
2
3
4
5
# global 模块参数修改
workspace → "G:\GROOPS\PNAS2020Workspace"
instrument → lri
timeStart → 58970 (2020-05-01 00:00:00)
timeEnd → 59091+23/24+59/1440+59/86400 (2020-08-30 23:59:59)

参数修改后在GROOPS运行S01SatelliteTrackingFromLRI-202007.xml脚本，运行后轨道结果保存在workdir/output路径下。

2.打开main_hybrid_pipeline.py，计算lgd和ra，先在main函数中初始化lgd处理器，lgd处理器初始化及计算lgd和ra的执行代码如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

if __name__ == "__main__":
    # 创建处理器实例
    processor = LgdProcessor(
        start_date='2020-08-01',    # 计算的起始日期
        end_date='2020-08-30',      # 计算的终止日期
        groops_workspace='G:\GROOPS\PNAS2020Workspace',     # 工作目录，即上文中的workdir
        instrument='LRI',           # 传感器类型，'LRI' 或 'KBR'
        lon_range=(88, 92),         # 研究区经度范围，东经为正，西经为负
        lat_range=(22, 26),         # 研究区纬度范围，北纬为正，南纬为负
        lat_limit=(-80.0, 80.0),    # 绘制穿过研究区轨道时轨道的延申纬度
        direction='asc'             # 限制穿过研究区范围的轨道类型，升轨（'asc'）、降轨（'desc'）和全部（'both'）
    )

    #  计算lgd和ra
    processor.execute(2)

运行后lgd和ra结果保存在workdir/results路径下。

3.Matlab中运行S03cwt_fliter_lgd_ra.m，计算lgd和ra小波分析重构后的结果，若只计算一天，则执行代码中的单日计算，代码如下

1
2
3
4
5
6
7
%% 单日计算
clear;clc;
date_str = '2020-06-15';
data_type = 'lgd';  % 'lgd' or 'ra'
freq = 0.5;     % 采样频率
process_cwt_fliter_signal_data(date_str, data_type, fs=freq, show_plots=true);

若批量计算，则执行批量计算部分代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
%% 批量计算
clear;clc;
freq = 0.5;     % 采样频率

start_datenum = datenum(2020, 8, 1);
end_datenum = datenum(2020, 8, 30);

date_nums = start_datenum:end_datenum;

date_list = cellstr(datestr(date_nums, 'yyyy-mm-dd'));

for i=1: length(date_list)
    process_cwt_fliter_signal_data(date_list{i}, "ra", fs=freq, show_plots=false);
    process_cwt_fliter_signal_data(date_list{i}, "lgd", fs=freq, show_plots=false);
end

运行后lgd和ra结果保存在workdir/results路径下。

4.打开main_hybrid_pipeline.py，计算lgd和ra的ASD曲线，在main函数中执行以下代码（这里不再展示初始化lgd处理器的代码）

1
2
#  计算lgd和ra的ASD图
processor.execute(4)

运行后lgd和ra结果保存在workdir/results路径下。

5.打开main_hybrid_pipeline.py，绘制lgd和ra的重构前后结果对比图，在main函数中执行以下代码（这里不再展示初始化lgd处理器的代码）

1
2
#  绘制lgd和ra的重构前后结果对比图
processor.execute(5)

运行后绘图结果保存在workdir/results路径下。

6.如果研究区域与上文介绍的PNAS论文一致，则可以进行对比，打开main_hybrid_pipeline.py，在main函数中执行以下代码即可（这里不再展示初始化lgd处理器的代码）。执行第该步骤必须更新PNAS数据目录（pnas_data_dir）。

1
2
3
# 执行第6步前必须更新PNAS数据目录（pnas_data_dir）
processor.update_config(pnas_data_dir='I:\LGD\GRACE-Follow-On-line-of-sight-gravity-processing-main\GRACE_data')
processor.execute(6)

运行后绘图结果保存在workdir/results路径下。

7.打开main_hybrid_pipeline.py，绘制lgd和ra的重构后的空间结果图，在main函数中执行以下代码（这里不再展示初始化lgd处理器的代码）

1
2
#  绘制lgd和ra的重构前后结果对比图
processor.execute(7)

运行后绘图结果保存在workdir/results路径下。

8.Matlab中运行S08convert_twsgrid_format_from_fengwei_toolbox.m，将fengwei工具包生成的月度tws结果转化成python能提取的格式，执行代码如下。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

clear;clc;

file_path = 'G:\GRACE_processing2\GRACE_results\gird_025_GSM_GFZ_RL06_DUAN_flt300_2002_2024_leakagefree.mat';
output_dir = 'grid_tws';

if ~exist(output_dir, 'dir')
    mkdir(output_dir);
end

load(file_path);
str_year = str2double(reshape(str_year, 233, 1));
str_month = str2double(reshape(str_month, 233, 1));

[path_str, name, ext] = fileparts(file_path);
output_file_name = [name, ext];

output_file_path = fullfile(output_dir, output_file_name);
save(output_file_path, 'grid_data', 'str_year', 'str_month', 'time', '-v7.3');

fprintf('已保存文件: %s\n', output_file_path);

运行后格式转换后的结果保存在workdir/grid_tws/路径下。

9.打开main_hybrid_pipeline.py，绘制grid tws空间分布图，在main函数中执行以下代码（这里不再展示初始化lgd处理器的代码）

1
2
3
4
5
6
7
8
9

# 执行第9步前必须更新TWS grid文件路径（tws_grid_filepath）、绘制的日期列表（date_spec)，同时可指定
# 绘图显示范围(tws_extent，不指定则使用self.lon_range和self.lat_range，指定None为全球)
processor.update_config(tws_grid_filepath=r'.\grid_tws\gird_025_GSM_GFZ_RL06_DUAN_flt300_2002_2024_leakagefree.mat')    # 上一步骤生成的转化格式后的grid tws文件
processor.update_config(date_spec=["2020-05", "2020-06", "2021-07", "2021-08"])     # 需要绘制的月份
processor.update_config(tws_extent=[80, 100, 10, 30])   # 绘图显示的范围
# processor.update_config(tws_extent=None)                # 绘图显示的范围，tws_extent=None为绘制全球
processor.execute(9)

注意：上面代码中指定的tws_grid_filepath需要经过第8步进行转换，相当于第8步转换后tws格网文件用于该步骤。
运行后绘图结果保存在workdir/results/tws_grid_plots路径下。

10.打开main_hybrid_pipeline.py，绘制穿过指定区域的多日lgd曲线图，在main函数中执行以下代码（这里不再展示初始化lgd处理器的代码）

1
2
3
4
5
6
7
8
9
10

# 执行第10步时必须更新日期列表(date_list)
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'
]
processor.update_config(date_list=date_list)
processor.execute(10)

运行后绘图结果保存在workdir/results路径下。

附件

上述中涉及到的所有程序脚本都打包到压缩包中，解压后文件目录如下。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
Scripts_grace_lri1b_lgd_20251110/
├── A00check_orbit_coverage.py
├── A01get_lri_kbr_records_num.py
├── A02read_grace_fo_shc_file.py
├── A03read_goco06s_file.py
├── A04combine_earth_gravity_field.py
├── process_cwt_filter_signal_data.m
├── S00GroopsIntegrateOrbitLRI-202007.xml
├── S01SatelliteTrackingFromLRI-202007.xml
├── S02compute_grace_lgd.py
├── S02compute_grace_lgd_gapfilled.py
├── S03cwt_filter_lgd_ra.m
├── S04lgd_ra_ipsd.py
├── S05plot_lgd_ra_cwt_filter.py
├── S05plot_lgd_ra_cwt_filter_gapfilled.py
├── S06compare_pnas_data.py
├── S07plot_lgd_ra_spatial_map.py
├── S08convert_twsgrid_format_from_fengwei_toolbox.m
├── S09plot_GRACE_L2_monthly_tws_from_matlab.py
├── S10plot_lgd_ra_multi_time_series_cross_over_area.py
└── main_hybrid_pipeline.py

下载链接：GRACE-FO LRI1B多日数据反演LGD和RA脚本

GitHub仓库：链接

学习
#GRACE #GRACE-FO #Python #Matlab #LGD #ASD #GROOPS #LRI1B
详细介绍利用GRACE1B多日数据计算LGD工作流程二_基于LRI1B多日数据
https://singyutang.github.io/2025/11/10/详细介绍利用GRACE1B多日数据计算LGD工作流程二-基于LRI1B多日数据/
作者
SingyuTang
发布于
2025年11月10日
许可协议