基于 GMTSAR 的 Sentinel-1 时序 InSAR 处理全流程详解

本文旨在详细阐述使用 GMTSAR 软件对 Sentinel-1 数据进行时序干涉测量（InSAR）处理的标准作业流程。涵盖内容包括数据准备、预处理、图像配准、干涉图生成、相位解缠、SBAS 解算以及最终的时序分析等关键步骤。

说明：文中未显式提供的脚本默认指代 GMTSAR 官方工具包中的标准脚本，用户可直接调用。

处理流程概览

1. 数据准备：整理原始影像、轨道文件及 DEM 数据。
2. 图像配准：执行粗配准与精密配准。
3. 差分干涉：生成干涉图并进行质量初筛。
4. 合并子条带：将分条带处理的数据进行拼接。
5. 干涉对检查：可视化检查并剔除低质量干涉对。
6. 相位解缠：对相位进行解缠处理。
7. 大气校正 (GACOS)：利用 GACOS 数据去除对流层延迟误差（可选）。
8. SBAS 解算：进行短基线集时序反演。
9. SBAS 后处理：地理编码、参考点校正及绘图。
10. 时序提取：提取特定点位的形变时间序列。

1. 数据准备 (Data Preparation)

处理的第一步是按照 Sentinel-1 的子条带（Subswaths: F1, F2, F3）建立标准化的文件目录结构，如下图所示：

1.1 链接数据文件

将对应的 *.tiff 和 *.xml 文件链接至 raw/ 目录。在 raw 目录下执行以下命令：

link_S1.csh

link_S1.csh Sentinel-1/data/path swathnumber

1.2 链接轨道文件

将精密轨道 *.EOF 文件链接至 raw/ 目录。在 raw 目录下执行以下命令：

link_S1_orbits.csh

link_S1_orbits.csh Sentinel-1/orbit/path

注意：执行完毕后会生成 data.in 文件，内容为数据文件名（无后缀）与轨道文件的对应列表，以冒号分隔。

1.3 准备 DEM

将数字高程模型文件 dem.grd 分别链接至 raw/ 和 topo/ 目录。

2. 图像配准 (Image Registration)

2.1 粗配准 (Coarse Registration)

首先运行批处理脚本进行粗配准：

preproc_batch_tops.csh data.in dem.grd 1

说明：运行完成后会生成 baseline_table.dat 文件，包含所有影像的时空基线信息。建议将该文件移动至上级目录，避免在精密配准步骤中被覆盖。

基线表内容示例： S1_20180201_ALL_F1 2018031.1874165505 1491 19.088688436880 56.434825714407 ...

mv baseline_table.dat ../

检查生成的 baseline.ps 图表，根据时空分布选择一幅位于中心位置的影像作为超级主影像 (Super Master)，所有从属影像将在下一步中与其进行精密配准。

2.2 精密配准 (Fine Registration)

编辑 data.in 文件，将选定的超级主影像所在行移动至文件第一行。

运行以下命令进行精密配准：

preproc_batch_tops.csh data.in dem.grd 2

提示：请在其余子条带（F2, F3）中重复上述配准操作。

3. 差分干涉 (Differential Interferometry)

3.1 生成干涉对列表

以 F1 目录为例，根据时空基线阈值筛选干涉对：

select_pairs.csh baseline_table.dat 50 100

参数说明：

• 50：时间基线阈值（天），干涉对时间间隔需小于 50 天。
• 100：垂直基线阈值（米），干涉对垂直基线需小于 100 米。

生成的 intf.in 文件格式如下： S1_20180508_ALL_F1:S1_20180514_ALL_F1 S1_20180508_ALL_F1:S1_20180526_ALL_F1 …

质量检查：检查生成的时空基线图，确保网络连通性，避免出现孤立的影像或干涉对。如有必要，需手动在 intf.in 中添加干涉对。

手动调整后，可使用以下脚本重新绘制基线图以确认： plot_baseline.csh

plot_baseline.csh intf.in baseline_table.dat

3.2 配置参数文件

修改配置文件 batch_tops.config：

• proc_stage = 1
  • 1：从反向地理编码开始，生成 topo_ra.grd。
  • 2：跳过地形生成，直接开始干涉图制作、滤波、解缠和地理编码（适用于并行处理）。
• master_image = S1_20180426_ALL_F1
  • 重要：必须修改为实际选定的超级主影像名称，且 F1/F2/F3 需保持一致。
• range_dec = 8 / azimuth_dec = 2
  • 多视。对于小区域研究，可调整为 4:1 以提高分辨率。
• dec_factor = 1
  • 小范围建议设为 1；大范围处理可设为 2 以降低分辨率并提高效率。

3.3 执行差分干涉处理

建议采用分步策略：

1. 测试运行：将 proc_stage 设为 1，提取 intf.in 中的第一行进行单对测试，确保参数设置无误。

   head -1 intf.in > one.in
   intf_tops.csh one.in batch_tops.config
   

2. 并行处理：测试通过后，将 batch_tops.config 中的 proc_stage 修改为 2（跳过 topo_ra 重复生成），执行并行处理。

   intf_tops_parallel.csh intf.in batch_tops.config 10
   

   生成的干涉图将存放于 intf_all 文件夹中。

   

操作提示：F2 和 F3 条带的处理仅需复制 F1 中的 intf.in 和 batch_tops.config 并修改对应条带名称即可。

4. 合并子条带 (Merge Subswaths)

4.1 生成合并列表

使用 create_merge_input.csh 脚本生成 merge.in 文件：

create_merge_input.csh intf_list path mode

参数 mode 说明：0 (合并3个子条带), 1 (合并F1/F2), 2 (合并F2/F3)。 示例：create_merge_input.csh intflist .. 0 > merge.in

注意：请务必检查 merge.in，将包含超级主影像日期的行置于文件顶部，以确保所有图像基于统一的坐标系和网格处理。

4.2 配置合并环境

将配置文件与 DEM 链接至 merge 目录：

cp ../F2/batch_tops.config .
ln -s ../topo/dem.grd .

4.3 执行合并

运行合并脚本：

merge_batch.csh merge.in batch_tops.config

首个干涉对合并完成后会生成 trans.dat 查找表文件。

并行优化：若数据量较大，可在 trans.dat 生成后终止上述命令，改用并行脚本：

merge_batch_parallel.sh inputfile config_file

5. 干涉对检查 (Quality Control)

为剔除因数据质量或处理异常导致的错误干涉对，需进行目视检查。

运行绘图脚本： select_phasefilt.csh

./select_phasefilt.csh mode

• mode 1：在 phasefilt_all 文件夹中汇总所有干涉图供预览。发现错误图件后，可直接在该文件夹中删除对应文件。
• mode 2：根据 phasefilt_all 中剩余的文件，自动将已剔除的干涉对原始文件夹移动至 rm 目录存档。

6. 相位解缠 (Phase Unwrapping)

相位解缠是时序处理中的计算瓶颈，也是决定结果精度的关键步骤。

6.1 准备解缠列表

在 merge 目录下生成列表：

ls -d 20* > intf.list

6.2 执行解缠

单线程解缠命令：

unwrap_intf.csh intf.list

参数配置：需预先编辑脚本中的参数： snaphu.csh correlation_threshold maximum_discontinuity [<rng0>/<rngf>/<azi0>/<azif>]

• correlation_threshold：相干性阈值，通常设为 0.02 - 0.1。
• maximum_discontinuity：相位跳变阈值，一般设为 0。
• [...]：解缠范围，留空则默认全图解缠。

6.3 并行解缠 (推荐)

为提高效率，建议使用并行解缠脚本：

unwrap_parallel.csh

./unwrap_parallel.csh intf.list 10

该操作依赖于当前目录下的 unwrap_intf.csh 单元脚本（内容应剔除循环语句，仅保留单次执行逻辑）：

unwrap_intf.csh

cd $1
snaphu[_interp].csh 0.02 0
cd ..

注意：若设置了特定解缠范围，需同步裁剪相干性文件 corr.grd： corr_cut.csh

corr_cut.csh intf.list

7. GACOS 大气校正 (Atmospheric Correction - Optional)

GACOS (Generic Atmospheric Correction Online Service) 校正可有效去除与地形相关的对流层延迟误差。

7.1 使用 GMTSAR 内置脚本 (v2025+)

GMTSAR 2025年7月2日更新版本已集成 GACOS 校正功能。推荐使用并行版本：

make_gacos_correction_parallel.csh intflist GACOS_path ref_range ref_azimuth dem.grd Ncores

7.2 获取参考点雷达坐标

需将选定的参考点经纬度转换为雷达坐标：

proj_ll2ra_ascii.csh trans.dat points_ll.txt points_ra.txt

• points_ll.txt 格式：lon lat name。
• 注意：至少需要三个非共线且间距适中的点才能完成坐标转换。

7.3 质量复查

校正后需再次检查干涉图质量： gacos_select.csh

gacos_select.csh intf.list

8. SBAS 解算 (SBAS Processing)

8.1 准备 SBAS 输入数据

将 intf.in 和 baseline_table.dat 复制至 SBAS 工作目录，运行准备脚本：

prep_sbas.csh intf.in baseline_table.dat ../merge unwrap.grd corr.grd

提示：若进行了 GACOS 校正，输入文件应替换为 unwrap_gacos_corrected_detrended.grd。该步骤将生成 intf.tab 和 scene.tab。

若前期剔除了部分干涉对，需更新 intf.in： intflist2intfin.sh

intflist2intfin.sh intf.list intf.in

8.2 执行 SBAS 反演

运行核心解算命令：

sbas intf.tab scene.tab N S xdim ydim [-atm ni] [-smooth sf] [-wavelength wl] [-incidence inc] [-range -rng] [-rms] [-dem]

关键参数详解：

输出结果：

• disp_##.grd：累积形变量 (mm)
• vel.grd：形变速率 (mm/yr)

9. SBAS 后处理 (Post-processing)

9.1 格式转换与地理编码

将 disp_##.grd 重命名为日期格式： post_sbas.csh

post_sbas.csh

计算平均相干性并生成掩膜： mask_meancorr.csh

ls */corr.grd > corr.list
stack_corr.csh corr.list meancorr.grd
mask_meancorr.csh meancorr.grd 0.12

将雷达坐标系结果转为地理坐标系（WGS84）： proj_disp_ra2ll.csh

proj_disp_ra2ll.csh disp_ra.list mask_file

9.2 参考点校正

选定稳定区域作为参考点，校正所有形变图： make_reference.csh

make_reference.csh

注意：运行前需编辑脚本，填入准确的参考点经纬度。

9.3 绘图与输出

绘制累积形变时间序列图： plot_disp_ll_zxj.csh

ls 20*referenced.grd > disp_referenced.list
plot_disp_ll_zxj.csh disp_referenced.list -300 100

绘制年形变速率图： plot_grd.csh

plot_grd.csh vel_mask_ll_referenced.grd -60 20

生成 KML 文件以在 Google Earth 中查看：

grd2kml.csh vel_file cpt_flie

10. 时序数据提取 (Time Series Extraction)

10.1 基于坐标列表提取

创建 points.list 文件，格式为 x y name。

方法 A：标准提取（区域平均） 计算指定半径（如 200m）范围内的平均形变值，平抑噪声： extract_ts_std.csh

extract_ts_std.csh

核心命令参考：gmt grdclip mask.grd -Sa0.2/NaN -Sb0.2/1 -Gmask.grd

方法 B：快速轨迹提取（单点） 不进行空间平均，直接提取像素值，速度较快： extract_ts_track.csh

extract_ts_track.csh

10.2 绘制时序曲线

使用以下脚本快速生成单点时序图： plot_timeseries.csh

plot_timeseries.csh name_ts_disp.txt