尧图网页设计 
文章目录前言一、数据准备1. Sentinel-1数据2. 精密轨道数据3. DEM数据二、处理流程1. 提取SLC2. 生成精配准所需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?pwdd477一、数据准备1. Sentinel-1数据本文收集了覆盖河北省燕郊镇的25景C波段开源Sentinel-1升轨SAR影像平均每月一期手动选择开展SBAS-InSAR地表形变监测。Sentinel-1升轨Path69覆盖了燕郊镇。数据详细信息如表1所示。表1 SAR影像信息轨道方向Path-Frame成像模式极化方式时间范围影像数量升轨69-129IWVV20240526-2026042225数据下载方法参考此前发布的专题推文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--polvv图1 read_S1_TOPS_SLC.py参数2. 生成精配准所需hgt文件1外部DEM格式转换使用srtm2dem命令将外部TIF转换为GAMMA可识别的.dem数据文件和dem.par参数文件。srtm2dem DEM.tif yj.dem yj.dem.par1- -图2 srtm2dem 参数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 精化查找表后模拟的SAR影像左和真实的SAR影像右3. 精配准建立RSLC文件夹调用S1_coreg_TOPS命令图6通过“强度匹配”与“谱分多样性”的迭代算法实现时序影像千分之一像素的精配准。图4 配准参数和命令举例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。图5 时空基线图调用mk_diff_2d命令基于初始基线组合引入雷达坐标系下的DEM扣除地形相位生成各干涉对的差分干涉图.diff与相干性图.cc。mk_diff_2d rslc_tab itab0DEM/seg.hgt_sim - mli/20250521.rmli mli diff_mb_1d825- - -图6 差分干涉图部分原始干涉图存在大量斑点噪声调用mk_adf_2d用Goldstein自适应滤波算法平滑相位提高条纹相干性。自适应滤波值滤波强度设为0.5滤波窗口大小设为64滤波器步长设为16。mk_adf_2d rslc_tab itab mli/20250521.rmli diff_mb_1d50.56416图7 滤波前左后右局部细节对比7. 相位解缠滤波前与滤波后均生成了相干性图选择远离沉降漏斗区、相干性极高且地表判定绝对稳定的城市建筑区作为解缠参考点本例中为2117 1218。调用mk_unw_2d应用最小费用流MCF算法进行相位解缠恢复出连续的形变相位场。这里设置解缠掩膜阈值为0.25即相干性低于0.25的区域会直接被掩膜掉不解缠这些区域。mk_unw_2d rslc_tab itab mli/20250521.rmli diff_mb_1d0.250.051111211712181图8 相位解缠结果注关于基线精化可根据研究区处理情况灵活选择是否进行此步骤相关命令为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- -图9 mb参数去趋势后干涉图还残余有湍流大气成分。采用时空滤波法去除湍流大气成分相关命令为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图10 ts_rate参数回归完成后利用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替换成相应的日期即可。三、结果展示将最终解算导出的年平均沉降速率GeoTIFFsbas_los_rate_disp_utm.tif与时序累积形变图直接导入QGIS或ArcGIS中叠加高分辨率卫星光学底图即可开展定量空间解译。图11 地表年平均沉降速率图图12 时间序列形变图 Sentinel-1的GAMMA SBAS-InSAR教程：以河北省燕郊镇地面沉降为例
文章目录前言一、数据准备1. Sentinel-1数据2. 精密轨道数据3. DEM数据二、处理流程1. 提取SLC2. 生成精配准所需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?pwdd477一、数据准备1. Sentinel-1数据本文收集了覆盖河北省燕郊镇的25景C波段开源Sentinel-1升轨SAR影像平均每月一期手动选择开展SBAS-InSAR地表形变监测。Sentinel-1升轨Path69覆盖了燕郊镇。数据详细信息如表1所示。表1 SAR影像信息轨道方向Path-Frame成像模式极化方式时间范围影像数量升轨69-129IWVV20240526-2026042225数据下载方法参考此前发布的专题推文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--polvv图1 read_S1_TOPS_SLC.py参数2. 生成精配准所需hgt文件1外部DEM格式转换使用srtm2dem命令将外部TIF转换为GAMMA可识别的.dem数据文件和dem.par参数文件。srtm2dem DEM.tif yj.dem yj.dem.par1- -图2 srtm2dem 参数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 精化查找表后模拟的SAR影像左和真实的SAR影像右3. 精配准建立RSLC文件夹调用S1_coreg_TOPS命令图6通过“强度匹配”与“谱分多样性”的迭代算法实现时序影像千分之一像素的精配准。图4 配准参数和命令举例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。图5 时空基线图调用mk_diff_2d命令基于初始基线组合引入雷达坐标系下的DEM扣除地形相位生成各干涉对的差分干涉图.diff与相干性图.cc。mk_diff_2d rslc_tab itab0DEM/seg.hgt_sim - mli/20250521.rmli mli diff_mb_1d825- - -图6 差分干涉图部分原始干涉图存在大量斑点噪声调用mk_adf_2d用Goldstein自适应滤波算法平滑相位提高条纹相干性。自适应滤波值滤波强度设为0.5滤波窗口大小设为64滤波器步长设为16。mk_adf_2d rslc_tab itab mli/20250521.rmli diff_mb_1d50.56416图7 滤波前左后右局部细节对比7. 相位解缠滤波前与滤波后均生成了相干性图选择远离沉降漏斗区、相干性极高且地表判定绝对稳定的城市建筑区作为解缠参考点本例中为2117 1218。调用mk_unw_2d应用最小费用流MCF算法进行相位解缠恢复出连续的形变相位场。这里设置解缠掩膜阈值为0.25即相干性低于0.25的区域会直接被掩膜掉不解缠这些区域。mk_unw_2d rslc_tab itab mli/20250521.rmli diff_mb_1d0.250.051111211712181图8 相位解缠结果注关于基线精化可根据研究区处理情况灵活选择是否进行此步骤相关命令为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- -图9 mb参数去趋势后干涉图还残余有湍流大气成分。采用时空滤波法去除湍流大气成分相关命令为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图10 ts_rate参数回归完成后利用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替换成相应的日期即可。三、结果展示将最终解算导出的年平均沉降速率GeoTIFFsbas_los_rate_disp_utm.tif与时序累积形变图直接导入QGIS或ArcGIS中叠加高分辨率卫星光学底图即可开展定量空间解译。图11 地表年平均沉降速率图图12 时间序列形变图