GMTSAR做InSAR时序分析：从环境配置到SBAS反演实战指南

简介：这份PDF教程聚焦GMTSAR开源InSAR处理软件的实际应用，面向熟悉GMT、希望掌握Sentinel-1干涉测量与SBAS时序分析的测绘、遥感方向学生和科研人员。教程以华北平原为实验区，从SAR数据与精密轨道下载、DEM拼接起步，覆盖GACOS大气改正、陆态网GNSS位移对比、干涉对选取、参数配置、相位解缠与融合生成时间序列的完整流程，并配有常用命令和脚本片段，便于读者快速建立处理链路。压缩包仅含1个PDF文件，大小约622KB，内容高度浓缩；目前已有3898人学习下载。通过学习该教程，可以理解GMTSAR的关键参数含义与处理逻辑，掌握结合GACOS削弱大气延迟、联合陆态网数据生成地表位移时间序列的方法，为地质灾害监测和地壳形变研究提供直接参考。文中还给出重要参数的含义与常见问题处理思路，能减少自学摸索时间，对InSAR初学者尤其实用。

1. 为什么这年头还有人坚持用 GMTSAR 做 InSAR 时序

GNSS 都躺进数据中心了，InSAR 处理软件却一直是付费商业软件和开源工具两条腿走路。GMTSAR 这个名字听起来像 GMT 的一个插件，实际上它是一套完整的 InSAR 处理框架，底层调用 GMT 做绘图和网格插值，核心计算用 C 和 Fortran 实现。这两年 Sentinel-1 数据免费开放，SBAS 时序分析需求暴增，GMTSAR 因为完全开源、命令行可控、批处理能力强，成了很多课题组和中小公司做地表形变监测的首选工具。它不需要 license，跑在 Linux 服务器上，配合 Python 脚本可以一口气处理上百景 SAR 数据。这篇笔记会把 GMTSAR 做 SBAS 时序的完整链路拆开，从数据准备到最终形变速率图，每一步参数怎么调、坑在哪里都写清楚，适合想摆脱商业软件束缚、自己掌控处理流程的从业者。

2. 从零搭建 GMTSAR 环境：装对了才能少折腾

2.1 安装前的依赖考虑

GMTSAR 对系统环境不算挑，但依赖项有讲究。Ubuntu 22.04 和 20.04 我都实测过，CentOS 7 也能跑但编译容易出问题。核心依赖包括 GMT（版本 5.4 或 6.x）、gcc/gfortran、libtiff、libhdf5、netcdf，还有用于读取 Sentinel-1 数据的 S1TBX 库（现在叫 snap-engine 的 Python 接口，但 GMTSAR 用的是自己的 S1TOOLS）。

最省事的方式是直接用发行版自带的包管理器装 GMT，然后单独编译 GMTSAR 本体：

sudo apt-get install gmt gmt-dcw gmt-gshhg libgmt-dev \ libnetcdf-dev libhdf5-dev libtiff5-dev libfftw3-dev \ gfortran make cmake git git clone https://github.com/gmtsar/gmtsar.git cd gmtsar autoremp -i ./configure --prefix=/usr/local/GMTSAR make -j4 sudo make install

编译完成之后必须设置环境变量，这一步漏掉会导致命令找不到：

export GMTSAR=/usr/local/GMTSAR export PATH=$PATH:$GMTSAR/bin export GMTHOME=/usr/lib (GMTHOME 根本不用设，但很多人误配)

参数说明：--prefix指定安装路径，建议装在独立目录方便卸载；make -j4是四核并行编译，机器内存小于 8GB 时不要开-j4，容易编译到一半被 OOM 杀掉。环境变量这里有个常见的误操作是export GMTHOME=/usr/lib，GMTSAR 源码里确实会找 G M T H O M E，但新版 GMT6 的库路径是自动检测的，手动指定反而会干扰运行。

2.2 数据目录结构与预处理脚本

装好软件之后要做的事是规划数据目录。GMTSAR 处理 Sentinel-1 数据的标准工作流要求把原始数据、轨道文件、DEM 按固定结构摆放。我一般用下面的目录骨架：

mkdir -p ~/insar/{raw,prep,proc,dem,orbits}

raw放从 ASF 下载的 Sentinel-1 原始数据（zip 包），orbits放精密轨道文件（POEORB），dem放研究区域的 SRTM 或 Copernicus DEM，prep存放预处理中间结果，proc放干涉处理结果。要注意的是不同版本 GMTSAR 对目录名有硬编码要求，比如topo/目录名在proc下是默认找 DEM 的地方，不要随便改。

cd ~/insar/prep alias s1_prep="$GMTSAR/bin/prep_s1"

预处理脚本会自动检测 zip 文件的脉冲总数和方位向零多普勒时刻，并生成raw_*后缀的临时目录。这一步常见报错是S1TOOLS路径未配置，需要在~/.bashrc里加上export S1TOOLS=/usr/local/GMTSAR/s1tools，否则预处理直接失败。

3. 单景干涉图生成流程：从 SLC 到 wrapped 相位

3.1 主影像选择与基线计算

做 SBAS 时序之前，先得把每一对干涉对各自处理出来。GMTSAR 的做法是先选定一景主影像，其余影像一一配对做干涉。主影像的选择不是随便挑一景，要看基线分布。用sat_baseline命令计算每景影像相对参考时间的垂直基线：

cd ~/insar/raw $GMTSAR/bin/sat_baseline -s1 -sw swath1 -r 20240101 -p 20240101

这条命令生成baseline_table.dat，第二列是垂直基线数值（单位米），第三列是时间基线（单位天）。我选择主影像的原则是让所有影像相对主影像的时间基线居中，同时垂直基线绝对值尽量小。从baseline_table.dat读出每景的基线后，画一个简单的分布图，挑居中且斜率不过大的那景。

常见做法是直接用plot_baseline脚本出图，但脚本依赖 GMT 的psbasemap，GMT6 版本下需要把-JX参数改成-JX新版写法。这一步不涉及干涉计算，但基线表错了后面全错，必须在进入下一步之前检查baseline_table.dat的内容是否与轨道状态向量一致。

3.2 配准与干涉图生成命令解析

配准用的是align_s1脚本，基于幅度互相关做粗配准，再基于相位做精配准。批量处理时可以写一个循环，逐对跑：

for i in $(seq -w 1 20); do $GMTSAR/bin/align_s1 -r 20240101 -m 20240201 -sw swath1 \ -rl 4 -ra 16 -refine -min_nr 32 done

-rl 4 -ra 16分别表示距离向和方位向的过采样倍数，4 和 16 是经验值，低相干区域可以提高到-rl 8 -ra 32，代价是计算时间几乎翻倍。-refine参数启用相位精配准，-min_nr 32是精配准时最少使用的窗口数。这一步产出的 SLC 偏移量文件会被interferogram命令读走。

干涉图生成用make_interferogram：

$GMTSAR/bin/make_interferogram -S 20240101 -M 20240201 \ -sw swath1 -rl 4 -ra 16 -filter 2 -pbase 0

-filter 2是滤波强度参数，取值范围 0 到 8，数值越大滤波越强。处理城市区域时我习惯-filter 2，山区或者植被覆盖大的区域用-filter 4到6，但滤波太强会抹掉真实形变信号，这个参数需要根据研究区反复试。-pbase 0表示不限制垂直基线，时序分析场景下保留所有干涉对，后续在 SBAS 反演时再做筛选。

3.3 滤波与相位解缠的参数经验

干涉图生成之后是一系列后处理：去平、滤波、解缠。GMTSAR 推荐的做法是用snaphu做相位解缠，snaphu不是 GMTSAR 自带的，需要单独安装。装好之后用snaphu_interf命令：

cd ~/insar/proc/20240101_20240201 $GMTSAR/bin/snaphu_interf -d -f filt_20240101_20240201.int \ -c correlation_20240101_20240201.grd -s 0.01

-s 0.01是 snaphu 的成本函数阈值，这个值决定了解缠相位的空间连续性约束强度。阈值取 0.01 是比较保守的数，相干性好的区域可以放到 0.05，相干性差的区域建议降到 0.005。解缠后的相位是 SBAS 输入的关键数据，格式为unw_phase_*.grd。

注意：解缠失败是最常见的翻车点，表现为解缠相位图出现明显断层或呈椒盐状。原因通常是-s值设置过高导致相位跳跃被过度惩罚，或滤波参数过大把真实形变条纹滤掉了。我处理山地数据时用-filter 4加-s 0.005的组合，成功率明显提高。

4. SBAS 时序反演：把干涉对变成形变时间序列

4.1 干涉对网络设计

SBAS 的核心思想是把多景影像组成的小基线集干涉对联合反演，得到形变速率和累积形变量。反演之前要先生成干涉对列表，GMTSAR 用sbas脚本自动完成配对和反演两步。干涉对列表文件格式是每一行两个日期，中间用空格或制表符分隔，注意日期格式必须是YYYYMMDD：

20240101 20240113 20240101 20240125 20240113 20240125

我一般用sbas_prep脚本基于空间基线和时间基线阈值自动生成配对表。阈值设置参考：空间基线小于 150 米，时间基线小于 60 天。如果是哨兵一号数据，时间基线还可以放宽到 120 天，因为 Sentinel-1 的轨道控制精度高，回归周期短，干涉对质量整体较好。配对表生成后要人工检查一遍——把每一对的空间基线和时间基线画成网络图，发现有孤立节点说明该影像与所有邻景基线都超限，这个时间点的形变信息在反演时会被约束得很差。

4.2 反演命令与输出文件解读

生成干涉对列表后，运行 SBAS 反演命令：

cd ~/insar/proc $GMTSAR/bin/sbas -f list.txt -d unwrap_list.txt \ -r 20240101 -z 0.01

-f list.txt指定干涉对列表，-d unwrap_list.txt指定解缠相位文件路径的列表，-z 0.01是正则化参数。正则化参数的反演逻辑是：当干涉对网络不连通或个别影像只在很少的干涉对中出现时，最小二乘解是病态的，需要加一个阻尼项让时间序列解不至于发散。-z 0.01表示阻尼权重为基线差的 1%，这个值是我从默认的 0.1 降下来的，因为 0.1 在多数场景下会让形变速率被严重平滑，真实的地表位移幅度会被低估。

反演结果的核心输出是：

velocity.grd # 平均形变速率 temporal_coherence.grd # 时间相干性 demeaned_ts.grd # 去除平均后的形变时序

时间相干性是对每个像素的反演残差做统计得到的，取值 0 到 1。质量评估时我以temporal_coherence.grd是否大于 0.3 作为可信区域的筛选线，低于 0.3 的像素形变信号基本可视为噪声。平均形变速率图可以直接叠加到光学影像上出结果图，但要注意 GMTSAR 输出的速率单位是米/年，不是毫米/年，出图前要乘以 1000，这是无数人栽过的坑。

4.3 反演精度的验证方法

反演完了不能直接拿去写报告，对照独立观测数据验证精度是必须的。GMTSAR 本身不提供验证模块，需要把形变时间序列导出成文本格式，和 GPS 或水准数据比较：

$GMTSAR/bin/grd2xyz velocity.grd -s | awk '{print $1, $2, $3*1000}' > velocity_mmyr.txt

grd2xyz -s不加-s会输出所有像素，数据量巨大；加了-s只输出有效网格。导出的第三列是形变速率（米/年），把它乘以 1000 转换成毫米/年。和 GPS 数据对比时，要先把 GPS 的速率插值到 InSAR 像素位置，直接用最近邻会导致比较结果偏差。我一般是用 kriging 插值把 GPS 速率场网格化，再在对应像素上计算残差的 RMS。实测下来，城市区域 RMS 在 2 到 3 毫米/年属于正常水平，超过 5 毫米/年就要回头检查配准和解缠环节是否出了问题。

5. 参数调优与避坑指南：SBAS 处理中那些血泪教训

5.1 参数该怎么调

GMTSAR 的默认参数不是为所有场景设计的。make_interferogram里的滤波参数对城区和山区要求完全不同，而 SBAS 反演的正则化参数和基线阈值之间的相互制约也经常被忽视。我整理了一张参数调整优先表：

5.2 避坑记录

坑 1：配准失败但程序不报错 现象：干涉图看起来有信号，但相位条纹完全杂乱无章，和真实地形不吻合。 原因：align_s1在低相干区域（植被覆盖区）配准窗口没能收敛，输出偏移量是局部极小值，后续干涉处理把错位当成形变信号。 解决：处理山区时，把-refine改成-refine -iter 4，同时把-rl和-ra从 4/16 降到 8/32，这组参数能显著提升配准成功率。

坑 2：时间相干性极高但形变速率全是零 现象：velocity.grd所有像素值接近 0，时间相干性却大于 0.8。 原因：反演输入的解缠相位文件路径列表写错了。unwrap_list.txt里写的是空的或全 0 的相位文件路径，SBAS 反演把全 0 序列完美拟合了。 解决：在运行sbas之前，手动用grdinfo检查任意一个解缠相位文件的数值范围，非 0 的相位波动是基本前提。

坑 3：SAR 图像拼接后出现条带差异 现象：多帧数据拼接后，相邻帧重叠区域形变速率相差 5 毫米/年以上。 原因：每一帧数据各自做了主影像选择，帧与帧之间的主影像基准不同，导致拼接处出现系统性偏差。 解决：多帧数据做 SBAS 时，要把所有帧的影像合并到一起选主影像，而不是分别处理后再拼接。具体做法是预处理时把所有帧的数据放进同一个目录，用统一的时间参考点做主影像。

坑 4：解缠相位图出现同心圆纹理 现象：解缠相位呈现规则的环形条纹，和地形无关，且各干涉对都出现。 原因：频率估计没做好，干涉图的空间参考相位没有去干净，导致每一对干涉图都有一个固定的斜坡相位。 解决：回看make_interferogram的日志，检查-pbase参数是否是 0。如果是 0，改用实际基线写入，例如-pbase 80，一般能消除该现象。

坑 5：SBAS 形变时间序列在最后一个时间点突然大幅跳变 现象：前 80% 时间点的形变趋势平稳，最后一景反演出的位移突然增加或减少 10 毫米以上。 原因：最后一景影像只和极少干涉对连接，反演网络在该节点约束不足，导致解被拉到噪声方向。 解决：在sbas反演时增加一个约束，或直接删除该时间节点重新反演。也可以通过调整配对阈值把该影像和更多相邻影像配对，再反演。

6. 一个加速批量处理的技巧：利用并行计算

SBAS 处理的瓶颈通常在干涉图生成环节，几十景数据产生的上百个干涉对要逐个执行align_s1和make_interferogram。干涉对之间没有依赖关系，天然适合并行。这里有个很实用的经验：不用再装额外的并行框架，GNU Parallel 就能解决问题。

cat list.txt | parallel -j 8 --colsep ' ' \ "cd ~/insar/proc/{1}_{2} && \ $GMTSAR/bin/make_interferogram -S {1} -M {2} \ -sw swath1 -filter 2 -pbase 0"

-j 8开启 8 个并行任务，适合 16 核服务器，核数少的机器建议-j 4。执行并行时要注意，make_interferogram会在当前目录生成临时文件，所以用cd切到独立目录再跑。--colsep ' '把 list.txt 的列传给{1}和{2}占位符。

并行处理最大的限制是磁盘 IO。如果干涉对数量上百，建议把不同干涉对的工作目录放在不同磁盘分区上，避免读写冲突。我在一台 8 核机器上处理 50 景数据，串行要跑 9 个小时，并行开 8 任务缩短到 1.5 小时，几乎线性加速。

解缠和反演阶段也可以并行，但解缠需要小心：snaphu 的重启动模式适合在内存受限时使用。snaphu_interf命令加上-t 2可以启用双线程处理，实测对提升大场景解缠速度有帮助。反演阶段 SBAS 内部默认使用单线程，但写好的输出文件多，不建议强行并行，顺序处理反而稳定。

一个小习惯是我每次处理完一批数据后，都会把干涉对列表、解缠相位列表和处理日志归档到单独目录，方便复盘。GMTSAR 的默认行为是覆盖写临时文件，如果过程中参数调错了想重跑，日志被覆盖就很难回溯当时用了什么参数。从那以后我每次处理数据都强制走一遍归档命令，把list.txt、unwrap_list.txt和baseline_table.dat复制一份带日期的版本。这个习惯帮我在返回去重处理旧数据时省了大量时间，希望帮到你。

延伸话题：GMTSAR 之外的 InSAR 处理工具对比、快速上手 GMTSAR 的开源数据集推荐（比如某区域模拟数据）、GMTSAR 与 ISCE 的转换接口哪个更实用。

本文还有配套的精品资源，点击获取