🛰️ Sentinel-1 时间序列 InSAR 数据处理完整流程 (基于GMTSAR)
🎯 本文将详细介绍如何使用 GMTSAR 软件对 Sentinel-1 数据进行时序 InSAR 处理，包括数据准备、预处理、配准、干涉图生成、相位解缠、最终时序分析等步骤。
🤔 文中没有提供的脚本，默认使用GMTSAR官方脚本，可直接调用。
🧭 处理流程概览
- 📁 数据准备
- 📚 图像配准
- 🌈 差分干涉
- 🌊 合并子条带
- 🗺️ 检查干涉对
- 🔓 相位解缠
- ☁️ GACOS大气校正
- 🛠️ SBAS解算
- 🧵 SBAS后处理
- 📈 形变时间序列提取
1. 📁 数据准备
首先按照子条带F1 F2 F3，对数据进行处理，文件目录设置如下： 
将相应的
*.tiff*.xml文件放置在raw/文件夹1
link_S1.csh Sentinel-1/data/path swathnumber
将轨道
*.EOF文件放置在raw/文件夹1
link_S1_orbits.csh Sentinel-1/orbit/path
此时会生成一个
data.in文件，其中具有数据文件名（无后缀）和轨道，用:分隔将DEM
dem.grd文件分别放置在raw/topo/文件夹
2. 📚 图像配准
- ⚙️ 首先进行粗配准，运行脚本
1
preproc_batch_tops.csh data.in dem.grd 1
完成后会生成一个
baseline_table.dat文件，其中包含数据的基线信息。可以将该文件移动至上层目录, 避免精密配准中，该文件重复生成。S1_20180201_ALL_F1 2018031.1874165505 1491 19.088688436880 56.434825714407 S1_20180207_ALL_F1 2018037.1879020806 1497 -55.518698696406 -34.273719507112 S1_20180213_ALL_F1 2018043.1874162052 1503 -41.389270891689 -22.026755662146 S1_20180219_ALL_F1 2018049.1878994363 1509 -50.211739839015 -25.142105633455 S1_20180225_ALL_F1 2018055.1874136333 1515 -20.370088939888 16.909205059498 S1_20180303_ALL_F1 2018061.1878994068 1521 -45.495652914617 -13.8650482657401
mv baseline_table.dat ../- 检查
baseline.ps文件，在基线图中间选择超级主影像，使所有影像在下一步的精密配准中与其配准。
- 检查
⚙️ 再进行精密配准 将
data.in中超级主影像所在行移动到文件第一行，然后进行精密配准。- 运行脚本，进行精密配准
1
preproc_batch_tops.csh data.in dem.grd 2
在其余子条带
F2F3中重复上述操作
- 运行脚本，进行精密配准
3. 🌈差分干涉
- 🧾 首先准备干涉列表，例如在
F1目录下运行：1
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
S1_20180508_ALL_F1:S1_20180607_ALL_F1
S1_20180508_ALL_F1:S1_20180520_ALL_F1同时也会生成时空基线图，检查是否有未连接的影像，手动添加干涉对，使时空基线图完整，没有单独的干涉对或影像。
手动调整后，可通过代码重新绘制新的时空基线图
1
plot_baseline.csh intf.in baseline_table.dat
🧾 修改配置文件
batch_tops.config"# 1 - start from make topo_ra"
"# 2 - start from make and filter interferograms, unwrap and geocode"
"proc_stage = 1"
🌟🌟🌟默认为1，从反向地理编码，生成topo_ra.grd开始。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，分辨率会降低，提高运算效率。- 🌈🌈🌈 差分干涉！！！
- 建议先将batch_tops.config中proc_stage改为1，对一个干涉对进行干涉，并检查。
1 2
head -1 intf.in > one.in intf_tops.csh one.in batch_tops.config
- 一旦完成，再次编辑batch_tops.config，将proc_stage改为2，跳过创建topo_ra文件。便可进行并行干涉处理。
1
intf_tops_parallel.csh intf.in batch_tops.config 10
所有干涉对都存放在
intf_all文件夹中，干涉图如下图所示：
剩余条带
F2F3的处理，仅需将F1目录中的intf.in和batch_tops.config复制过来，并修改为F2F3即可。
4. 🌊 合并子条带
- 首先制作merge.in文件，使用
create_merge_input.csh脚本1
create_merge_input.csh intf_list path mode
mode 0 是合并 3 个subswaths, mode 1 是合并 F1/F2, mode 2 is 合并 F2/F3
例如:create_merge_input.csh intflist .. 0 > merge.in ⚠️注意将含有超级主影像日期的行全部置顶，确保所有图像都使用相同的坐标处理，并且最终网格大小相同。
- 将配置文件
batch_tops.config复制到merge目录，dem.grd也链接过来。1 2
cp ../F2/batch_tops.config . ln -s ../topo/dem.grd .
- 合并子条带 运行脚本：
1
merge_batch.csh merge.in batch_tops.config
在合并完第一个干涉对后，会生成
trans.dat文件。- 💡 （可选）若干涉对数量巨大，可选择并行操作，需要在生成
trans.dat文件后，终止上个命令。运行并行脚本：1
merge_batch_parallel.sh inputfile config_file
- 💡 （可选）若干涉对数量巨大，可选择并行操作，需要在生成
5. 🗺️ 检查干涉对
由于数据本身、处理过程、操作失误等因素，部分干涉对出现错误，我们需要检查并剔除错误干涉对。
我们先将干涉图绘制出来，使用脚本：
1
./select_phasefilt.csh mode
mode 1将建立
phasefilt_all文件夹，汇总所有的干涉对以便查阅。错误干涉图在此删除即可。 mode 2将剔除的干涉图对应的文件夹移动到rm文件夹中。
6. 🔓 相位解缠
相位解缠是InSAR时间序列处理中最耗时的步骤之一。在平坦干燥地区，解缠速度较快；而在地形复杂，植被覆盖的地区，相位解缠准确度较低。
- 在
merge目录中，首先制作一个干涉图列表🧾1
ls -d 20* > intf.list
- 运行解缠脚本
1
unwrap_intf.csh intf.list
运行前需要先修改解缠参数
"snaphu.csh correlation_threshold maximum_discontinuity [<rng0>/<rngf>/<azi0>/<azif>]"
💡相干性阈值correlation_threshold一般选0.02 - 0.1，maximum_discontinuity地震相位跳跃，一般取0，[<rng0>/<rngf>/<azi0>/<azif>]为解缠范围，默认全部解缠。 🌟🌟🌟通常为提高运算效率，我们选择并行解缠：
1
./unwrap_parallel.csh intf.list 10
当前目录需要放置一个unwrap_intf.csh文件，注意与之前的不同，不需要循环语句，内容大致如下：
1 2 3
cd $1 snaphu[_interp].csh 0.02 0 cd ..
📢📢📢 如果在这一步，设置了解缠范围。需要将相干性文件
corr.grd进行对应的裁剪，使用脚本：1
corr_cut.csh intf.list
7. ☁️ GACOS大气校正 (选做)
GACOS大气校正能够较好的去除与地形相关的对流层大气误差, 通常我们用GACOS进行大气校正。
在
merge目录下，我们选择并行gacos校正：make_gacos_correction_parallel.csh
make_gacos_correction.csh1
make_gacos_correction_parallel.csh intflist Ncores
脚本会调用
make_gacos_correction.csh脚本进行大气校正，需要修改参数信息：
set gacos_path = /GACOS_path/ # GACOS file (*.ztd, *.rsc) path
set pixel_center = 41347
set line_center = 2422
💡这是参考点的雷达坐标，一定要修改，选择研究区内稳定参考点。地理坐标系转雷达坐标系（参考点）：
1
proj_ll2ra_ascii.csh trans.dat points_ll.txt points_ra.txt
至少需要三个点才能转换，三个点不能离得太近，不能在一条直线。
格式为xyname
👀 解缠和gacos校正后的干涉对仍需检查，需剔除（或重新处理）错误干涉对，可在
merge目录执行脚本：
1
gacos_select.csh intf.list
8. 🛠️ SBAS解算
- 准备输入文件，使用脚本自动生成，需要将
intf.in和baseline_table.dat复制到SBAS目录，运行：1
prep_sbas.csh intf.in baseline_table.dat ../merge unwrap.grd corr.grd
若采用了GACOS大气校正，
unwrap.grd应改为unwrap_gacos_corrected_detrended.grd脚本会生成intf.tab和scene.tab，以及一行简单的SBAS指令- 若在前几步删除了一些干涉对（yyyyddd_yyyyddd），此时对应的intf.in应重新生成，可使用脚本：
1
intflist2intfin.sh intf.list intf.in
我们在命令行输入sbas指令：
1
sbas intf.tab scene.tab N S xdim ydim [-atm ni] [-smooth sf] [-wavelength wl] [-incidence inc] [-range -rng] [-rms] [-dem]
输入参数:
intf.tab -- 干涉对列表
scene.tab -- SAR数据列表（天数从2014年1月1日起算）
N -- 干涉对数量
S -- SAR数据的数量
xdim and ydim -- 干涉图的尺寸
-smooth sf -- 平滑因子，默认为0，一般选3
-atm ni -- 大气校正（CPS方法）迭代次数，默认为0(跳过校正)，一般选3
-wavelength wl -- 雷达波长
-incidence theta -- 入射角，在xml中可以查看（incidence）
-range rng -- 雷达到干涉图中心的距离，在PRM文件中可以查看（near range）
-rms -- 形变速率的均方根误差 rms.grd
-dem -- 输出的dem误差 dem_err.grd
-mmap -- 使用硬盘代替运行内存，降速明显
-robust -- only work with -atm turnned on, estimate velocity with records that has atm correction
输出参数:
disp_##.grd -- 累积形变量 (mm)
vel.grd -- 形变速率 (mm/yr)
9. 🧵 SBAS后处理
求出的
disp_##.grd和vel.grd需要进一步处理才能变为我们最终的结果。
首先修改文件名，将
disp_##.grd修改为日期命名，运行脚本：1
post_sbas.csh
然后进行掩膜和地理编码，这两步同时进行。此时的数据还是雷达坐标系，需转为地理坐标系。
先进入merge目录计算干涉对平均相干性，操作如下：
1 2 3
ls */corr.grd > corr.list stack_corr.csh corr.list meancorr.grd mask_meancorr.csh meancorr.grd thresholds
thresholds一般取0.12
再将相干性阈值文件
mask_file和地理编码文件trans.dat链接到SBAS目录，运行脚本：1
proj_disp_ra2ll.csh disp_ra.list mask_file
生成
yyyymmdd_mask_ll.grd与vel_mask_ll.grd文件
接着进行参考点校正，注意修改脚本里的参考点坐标。
1
make_reference.csh
生成
yyyymmdd_mask_ll_referenced.grd与vel_mask_ll_referenced.grd文件绘制累积形变时间按序列图：
1 2
ls 20*referenced.grd > disp_referenced.list plot_disp_ll_zxj.csh disp_referenced.list -300 100
绘制形变速率图：
1
plot_grd.csh vel_mask_ll_referenced.grd -60 20该脚本适用于绘制任何WGS84坐标的栅格数据
grd文件转kml, 用于在Googleearth中查看：
1
grd2kml.csh vel_file cpt_flie
10. 📈 形变时间序列提取
时间序列分析是InSAR的关键环节，我们在这里介绍如何提取一些点位的时间序列。
首先需要创建一个
points.list, 里面每一行包含一个点的信息, 格式为xyname
同时需要存在一个grd_list, 即我们要提取的累计形变栅格文件，一般为disp_referenced.list，执行脚本即可生成特征点的时间序列。1
extract_ts_std.csh
脚本中，
gmt grdclip mask.grd -Sa0.2/NaN -Sb0.2/1 -Gmask.grd，其中0.2即200米圆形范围内取平均，可根据需求修改。若需要快速提取特征点的时间序列，且不需要范围内取平均，可采用以下脚本，能够大大提高提取速度。
1
extract_ts_track.csh
使用以下脚本可以快速绘制时间序列图像：
1
plot_timeseries.csh name_ts_disp.txt