文章目录
- 前言
- 一、数据准备
- 1. Sentinel-1数据
- 2. 精密轨道数据
- 3. DEM数据
- 二、处理流程
- 1. 提取SLC
- 2. 生成精配准所需hgt文件
- （1）外部DEM格式转换
- （2）生成初始地理编码查找表
- （3）查找表精化
- 3. 精配准
- 4. 去斜与镶嵌
- 5. 地理编码
- 6. 差分干涉与滤波
- 7. 相位解缠
- 8. 去趋势
- 9. mb计算（多基线计算相位时间序列）
- 10. 地表形变年形变速率计算
- 11. 时序形变计算
- 三、结果展示
前言
上一篇教程以日喀则市定日县地震为例，介绍了基于GAMMA软件开展两轨D-InSAR形变监测的基本流程。本篇将以河北燕郊地区为研究区，进一步介绍时序InSAR处理方法，重点讲解基于GAMMA软件的SBAS-InSAR数据处理流程，包括多景SAR影像配准、干涉对构建及时序形变反演等关键环节。本系列教程由地枢遥感整理，仅作 InSAR 技术学习与操作记录。
1️⃣S1原始数据；
2️⃣GAMMA处理过程数据；
3️⃣GAMMA处理成果数据。
本次使用数据链接(191GB大小 zip格式):https://pan.baidu.com/s/1ZuquRfiPrihbdo6UN0tDbg?pwd=d477
一、数据准备
1. Sentinel-1数据
本文收集了覆盖河北省燕郊镇的25景C波段开源Sentinel-1升轨SAR影像，平均每月一期（手动选择），开展SBAS-InSAR地表形变监测。Sentinel-1升轨Path69覆盖了燕郊镇。数据详细信息如表1所示。
| 轨道方向 | Path-Frame | 成像模式 | 极化方式 | 时间范围 | 影像数量 |
|---|---|---|---|---|---|
| 升轨 | 69-129 | IW | VV | 20240526-20260422 | 25 |
数据下载方法参考此前发布的专题推文：
InSAR处理全流程：Sentinel-1卫星数据获取指南
2. 精密轨道数据
精密轨道下载可使用由地枢遥感提供的轨道数据下载器，具体使用方法可参考此前发布的专题推文：
InSAR处理全流程：Sentinel-1精密轨道数据获取指南
3. DEM数据
DEM下载可使用由地枢遥感提供的DEM数据下载器，具体使用方法可参考此前发布的专题推文：
InSAR处理全流程：DEM数据获取指南
二、处理流程
1. 提取SLC
准备研究区的kml范围文件，调用read_S1_TOPS_SLC.py脚本进行SLC数据提取。本案例根据yj.kml几何边界文件自动过滤并剪切出研究区所在的Burst区域，生成子条带iw1的.slc、.slc.par和.tops_par文件。
read_S1_TOPS_SLC.py S1A_IW_SLC__1SDV_20240526T100608_20240526T100635_054041_069205_C250.zip--burst_selyj.kml--OPOD_dir./orbit--out_dir./SLC--polvv2. 生成精配准所需hgt文件
（1）外部DEM格式转换
使用srtm2dem命令将外部TIF转换为GAMMA可识别的.dem数据文件和dem.par参数文件。
srtm2dem DEM.tif yj.dem yj.dem.par1- -（2）生成初始地理编码查找表
运行gc_map命令计算初始查找表（lookup_table），该文件建立了地理坐标系和雷达坐标系之间的映射关系；同时，该命令会根据雷达成像的几何参数和DEM信息模拟出一个后向散射强度图（sim_sar）。
SLC_mosaic_S1_TOPS 20250521_slc120250521.rslc20250521.rslc.par82multi_look20250521.rslc20250521.rslc.par20250521.rmli20250521.rmli.par82gc_map20250521.rmli.par - yj.dem.par yj.dem seg.dem.par seg.dem20250521.lt1120250521.sim_sar uvinc psi pix ls_map82-（3）查找表精化
采用offset_pwrm计算局部互相关偏移值，随后运行offset_fitm计算偏移多项式。
offset_pwrm pix_sigma020250521.rmli20250521.diff_par20250521.offs20250521.ccp128128offsets264640.2offset_fitm20250521.offs20250521.ccp20250521.diff_par coffs coffsets0.256运行完拟合后，检查终端日志或日志文件，看精度是否符合要求。
拟合达标后，利用gc_map_fine对初始查找表进行优化，输出精细查找表lookup_table.fine。最后，运行geocode命令将DEM编码到SAR坐标系下，生成.hgt文件，用于精配准
3. 精配准
建立RSLC文件夹，调用S1_coreg_TOPS命令（图6），通过“强度匹配”与“谱分多样性”的迭代算法，实现时序影像千分之一像素的精配准。
S1_coreg_TOPS 20250521_slc12025052120240526_slc12024052620240526_rslc20250521.hgt82- -0.60.010.810配准结束后，检查生成的配准质量文件coreg_quality，确保方位向精度严格小于千分之一；同时查看生成的差分干涉图（*.diff.bmp），目视确认burst拼接处条纹连续、无错位。
4. 去斜与镶嵌
配准好的时序子条带数据仍带有相位斜坡，在整景合并前必须进行去斜处理。命令S1_deramp_TOPS_reference（专门为主影像去斜），并用SLC_mosaic_S1_TOPS将各个条带合并，生成最终的RSLC。对于辅影像，则调用S1_deramp_TOPS_slave严格对齐主影像进行去斜。
S1_deramp_TOPS_reference 20250521_slc1 SLC_mosaic_S1_TOPS 20250521_slc1.deramp20250521.slc20250521.slc.par82S1_deramp_TOPS_slave 20240526_slc12024052620250521_slc182-5. 地理编码
此处与3.2步骤一致，但此处是去斜后的，生成seg.hgt_sim。
6. 差分干涉与滤波
SBAS流程的核心在于构建短时空基线的自由组合网络，以最大程度抑制时空失相干。由于获取的燕郊镇Sentinel-1A影像为平均每月一期，因此，设置最大时空基线阈值分别为600m、96天，使用base_calc进行全组合自由构网，生成干涉对索引表itab，同时利用base_plot绘制时空基线网络连接图（如图5）。
调用mk_diff_2d命令，基于初始基线组合，引入雷达坐标系下的DEM扣除地形相位，生成各干涉对的差分干涉图（.diff）与相干性图（.cc）。
mk_diff_2d rslc_tab itab0DEM/seg.hgt_sim - mli/20250521.rmli mli diff_mb_1d825- - -原始干涉图存在大量斑点噪声，调用mk_adf_2d用Goldstein自适应滤波算法平滑相位，提高条纹相干性。自适应滤波值（滤波强度）设为0.5，滤波窗口大小设为64，滤波器步长设为16。
mk_adf_2d rslc_tab itab mli/20250521.rmli diff_mb_1d50.564167. 相位解缠
滤波前与滤波后均生成了相干性图，选择远离沉降漏斗区、相干性极高且地表判定绝对稳定的城市建筑区作为解缠参考点（本例中为2117 1218）。
调用mk_unw_2d，应用最小费用流（MCF）算法进行相位解缠，恢复出连续的形变相位场。这里设置解缠掩膜阈值为0.25，即相干性低于0.25的区域会直接被掩膜掉，不解缠这些区域。
mk_unw_2d rslc_tab itab mli/20250521.rmli diff_mb_1d0.250.051111211712181注：关于基线精化
可根据研究区处理情况灵活选择是否进行此步骤，相关命令为：mk_base_2d，之后进行二次差分干涉、滤波、解缠。
研究区范围较小且已使用Sentinel-1精密轨道，一般不容易出现显著的趋势相位；另外，由于研究区平坦，DEM模拟SAR图像无法正确反映平原区域结构特征，基线精化容易出错，故选择跳过。
8. 去趋势
查看部分干涉图仍有趋势向误差残留，使用mk_quad_2d命令去除趋势误差。adf.cc.ave.mask.bmp为mask（掩膜），避免把真正的形变当成轨道误差给滤掉。
mk_quad_2d rslc_tab itab mli/20250521.rmli diff_mb_1d diff_quad_1d101- -33adf.cc.ave.mask.bmp119. mb计算（多基线计算相位时间序列）
手动剔除因局部失相干或解缠跳变导致质量较差的干涉图，并保证干涉网络的连通性，通过编辑itab完成。
使用mb通过奇异值分解（SVD）算法来获取相位时间序列的最小二乘解，同时解算出高程残差hgt_out与时序相位标准差sigmal_ts。
mb diff_tab1 rmli_tab itab_1 - itab_ts diff1_ts/diff3_ts1diff1_ts/sigmal_ts1diff1_ts/hgt_out211712181616- -去趋势后，干涉图还残余有湍流大气成分。采用时空滤波法去除湍流大气成分，相关命令为tpf、fspf。
10. 地表形变年形变速率计算
使用ts_rate命令从求解出的已解缠相位时间序列中提取出长期的年线性形变速率（年下沉速率）。
ts_rate diff_ts/ts_tab rmli_tab itab_ts - SBAS_rate/sbas_rate SBAS_rate/sbas_const SBAS_rate/sbas_sigma_ts回归完成后，利用dispmap将雷达坐标系下的形变相位转换为标准雷达视线方向(LOS)的形变量速率文件(单位：m/year)；然后将其反向地理编码并导出为GeoTIFF成果。
11. 时序形变计算
同样地，为了获取燕郊地表沉降随时间推移的动态演变过程，使用dispmap将时序形变相位转为形变量（默认为LOS向，形变量单位为m），接着利用精细查找表lookup_table.fine运行geocode_back批量执行反向地理编码，将形变转为地理坐标系。
dispmap diff_ts/diff3_ts_001.diff.clean DEM/yj.rdc.sim_sar rslc/20260422.rslc.par diff_mb_1d/20260422_20260528.off diff_ts/diff3_ts_001.diff.clean.ldisp000geocode_back diff_ts/diff3_ts_001.diff.clean.ldisp2708DEM/yj.Fine_lookup_table diff_ts/diff3_ts_001.diff.clean.geo.ldisp21101-0011利用data2geotiff命令，导出GeoTIFF格式的时序累积形变成果。
data2geotiff DEM/yj.utm.par diff3_ts_001.diff.clean.ldisp2diff3_ts_001.diff.clean.ldisp.tif最后，根据影像日期将diff3_ts_*.diff.clean.ldisp.tif替换成相应的日期即可。
三、结果展示
将最终解算导出的年平均沉降速率GeoTIFF（sbas_los_rate_disp_utm.tif）与时序累积形变图直接导入QGIS或ArcGIS中，叠加高分辨率卫星光学底图，即可开展定量空间解译。
