做InSAR的人应该都有印象2015年4月25日尼泊尔廓尔喀地震发生之后全球很多研究组在几天内就用哨兵数据拿到了同震形变图。Sentinel-1A当时刚进入业务化运行不久ISCE也正好是少数能把TOPS模式处理得比较干净的软件。这篇博文我把自己的完整处理流程写一遍包括数据准备、ISCE命令、参数调整和踩坑记录目标是用这套流程复现2015年尼泊尔Mw 7.8地震的LOS向同震形变场。我尽量不把文章写成操作手册而是把每一步背后的逻辑也讲清楚。毕竟InSAR数据处理这活参数随便改一下结果的差别就很大只看命令很难真正上手。如果你已经有点基础想用ISCE干一票形变遥感这篇文章可以直接照着走。1. 一个“大地震的形变场”到底怎么测1.1 2015年尼泊尔地震为何适合被InSAR观测2015年4月25日尼泊尔发生Mw 7.8地震震中位于廓尔喀地区是喜马拉雅山主逆冲断裂上的一次典型破裂事件。地震造成的地表形变范围大、信号强又发生在山地和高原交界地带这种场景非常适合用雷达干涉测量来观测。为什么不说GPS也不说光学影像因为InSAR可以直接测量地表在雷达视线方向LOS的厘米级位移大面积连续成像不需要地面布点。尼泊尔地震的震中区域山高谷深GPS站点稀疏光学影像虽然能看到地表破坏但没有办法直接给出空间连续的位移场。而Sentinel-1作为C波段SAR卫星在多云多雨地区仍有穿透能力能够获取到可靠的面状形变信号。另外尼泊尔地震发生在2015年4月Sentinel-1A卫星在2014年发射后已经完成了初期定标全球覆盖的干涉宽幅模式IW数据积累稳定。时间上刚好赶上哨兵数据在亚洲大型地震中第一次成规模用于地震形变分析。对于做InSAR数据处理的人来说这是一个练手和验证流程的绝佳案例。1.2 选哨兵数据和ISCE软件的理由哨兵数据指的是欧洲空间局的Sentinel-1雷达卫星。它默认使用TOPS模式成像幅宽250公里重访周期12天单景覆盖范围大特别适合区域尺度的形变监测。和早期的ENVISAT ASAR或ALOS PALSAR数据相比哨兵数据免费开放分辨率适中而且一直稳定运行这些年已经成为InSAR应用的主流数据源。TOPS模式有个特点卫星在成像过程中主动进行方位向扫描导致每一景数据被分成了若干个burst。处理的时候必须对相邻burst之间的重叠区域做精确配准否则干涉图会出现明显的相位跳变。这比传统条带模式复杂得多用一般商业雷达软件容易出问题。ISCE是美国JPL开发的开源InSAR处理框架全称是Interferometric synthetic aperture radar Scientific Computing Environment。它最大的优势就是对TOPS模式也就是Sentinel-1数据支持得很完整而且把InSAR处理的各个环节封装成了可以单独调用的脚本。如果你以后还要处理ALOS-2、TerraSAR-X、甚至未来的NISAR数据ISCE这套架构同样适用。我一直建议刚入门的人别只围着商业软件转。ISCE虽然装起来有点门槛但一旦跑通了你对每一步处理到底干了什么会理解得更深排查问题也更有底。这篇文章的核心流程全部用ISCE完成。1.3 整体数据处理管线速览在进入细节之前先把整个得到同震形变场的流程串一遍。拿到两景Sentinel-1 SLC数据后第一步是准备轨道数据和DEM第二步用ISCE的topsApp建立干涉工程第三步生成干涉图做多视和滤波第四步做相位解缠把包裹相位还原成真实位移第五步地理编码输出GeoTIFF格式的LOS向形变图。整个流程看起来简单但每一步都藏着不少坑。比如主影像和辅影像的时相搭配如果时间基线太长或者垂直基线太大相干性会很差又比如相位解缠的时候如果区域过大、噪声过多解缠结果会出现大片跳变。这些我都会在后面的章节展开。还有一个必须提前说的事情ISCE处理本身需要比较多的计算资源。尼泊尔地震这种范围单景SLC大概十几GB处理完中间产品可能到几十GB。如果你的电脑内存不足16GB或者磁盘空间小于100GB建议先别开始不然中间很容易因为磁盘写满或内存溢出中断。2. 数据准备SLC、轨道、DEM一个都不能少2.1 Sentinel-1 SLC数据的选择与下载ISCE处理InSAR需要的数据是Level-1的SLC产品也就是单视复数数据不是大家平时下载的GRD强度影像。SLC保留了相位信息是干涉处理的基础。选择数据时要特别注意三点同一轨道、同一模式、同一极化。以2015年尼泊尔地震为例我用的是Sentinel-1A的IW模式VV极化SLC数据。轨道号选择覆盖震中区域的同一条轨道确保主影像和辅影像的空间几何一致。下载前最好先查一下欧洲空间局的数据目录确认你想用的两景数据都覆盖了目标区域且震前和震后都有数据。时间上建议选地震前后各一景比如主影像选2015年4月13日辅影像选2015年4月25日或者4月29日。有人会问地震当天能不能直接用来作辅影像可以但如果震中区域发生严重地表去相干干涉图的信噪比会明显下降。稳妥的做法是震前和震后各留一景形变场只需要差分掉非形变相位时间间隔短一些能减少噪声。下载的时候建议直接搜SLC产品并下载原始的zip包或者SAFE格式。ISCE能直接读取SAFE目录下的manifest.safe文件。注意核对产品文件名里的绝对轨道号比如Path 142和切片号比如Frame 3072两条数据需要完全对应。不要小看这一步我见过不少人把不同轨道的两景数据放在一起跑结果生成了一堆噪声。2.2 精轨数据与DEM的准备轨道文件是InSAR处理中的关键辅助数据。哨兵卫星有预报轨道和精轨两套数据处理InSAR一定要用精轨。精轨数据文件名通常以AUX_POEORB开头代表了厘米级的轨道精度能显著改善干涉图中由轨道误差引起的条纹。下载轨道时要注意卫星编号、时间范围和绝对轨道号要和SLC数据的获取时间对齐。还有一个容易忽略的点轨道文件下载后ISCE并不会自动找到它需要在配置文件里指定轨道目录。我建议单独建一个orbits文件夹把下载的精轨文件按日期重命名方便ISCE扫描。比如这样组织project/ ├── SLC/ │ ├── S1A_IW_SLC__1SDV_20150413TXXXXXX/ │ └── S1A_IW_SLC__1SDV_20150425TXXXXXX/ ├── orbits/ │ └── S1A_OPER_AUX_POEORB_OPOD_20150425TXXXXXX_V20150413TXXXXXX_20150413TXXXXXX.EOF └── dem/ └── srtm1.dem再说DEM。ISCE做地形相位去除和地理编码都依赖DEM。推荐使用SRTM 1弧秒数据大约30米分辨率在尼泊尔山区基本够用。ISCE自带的dem.py脚本可以自动从网络获取SRTM数据但我更倾向于自己先把覆盖目标区域的tif文件裁切好再转成ISCE需要的dem格式。原因很简单自动下载如果断网或者中途失败耗一整晚也不是没可能。用dem.py的常见方式是这样的dem.py --bbox 80.0 26.0 88.0 29.0 -r -s 1 -c -f其中--bbox是经纬度范围-s 1表示1弧秒-c表示裁切-f表示填充。也可以用gdal来手工处理。关键是最后生成的dem文件需要是一个带地理参考的单波段高程栅格ISCE会把它转成自己内部使用的格式。如果DEM范围没完全覆盖SAR影像后面的处理会报错所以宁可扩大一点范围比如在目标区域四周多留0.2度。2.3 主辅影像时相组合和基线检查在正式运行ISCE之前最好先估算一下主辅影像之间的时间和空间基线。时间基线太长会导致时间去相干垂直基线太大会导致地形相位过于密集尤其是山区容易超过奈奎斯特采样极限。对于尼泊尔这种地形起伏大的区域垂直基线最好控制在100米以内时间基线就是哨兵的12天或24天。ISCE提供了一种快速检查基线的方法在配置好topsApp.xml之后可以直接运行startup步骤它会输出主辅影像的轨道状态和估算的空间基线。如果垂直基线上百甚至几百米建议换一景数据或者调整主影像的选择。做同震形变场并不是说主影像一定要严格在震前有时候为了压低时间和垂直基线也会选震后更晚的影像做辅影像再用其他方法隔离形变。但这里直接用震前震后各一景最简单。有一个实操细节在处理前检查数据覆盖的时候把震中附近区域用地图工具打开确认两景SLC的成像时间都是白天还是夜间。哨兵数据通常默认降轨白天或者升轨夜间这不影响干涉但会影响后续形变图的时间解释。尼泊尔地区覆盖角度很多选择时只要保证同一轨道的同一帧即可。3. ISCE处理实战把SLC变成干涉图3.1 topsApp工程的目录组织与配置ISCE处理TOPS数据的主程序是topsApp.py。它通过一个XML格式的配置文件控制整个流程。一般工作流程是先建一个工程目录把SLC数据、轨道、DEM放好然后创建topsApp.xml并填写参数。topsApp.xml需要花点心思。核心部分大致包括主影像参数、辅影像参数、干涉参数、滤波参数和解缠参数。以我常用的配置为例topsApp component nameisce property namesensor name valueSentinel-1/ property nameorbit directory value./orbits/ property namedem filename value./dem/srtm1.dem/ property nameregion of interest value84.5 27.5 87.5 29.0/ property nameswaths value1/ property namepolarization valueVV/ property nameazimuth looks value2/ property namerange looks value8/ /component /topsApp这里的swaths默认可以填1也可以填123表示处理全部的IW子测绘带。尼泊尔地震的形变场范围不算特别大我建议处理整个swath 1到3都打开后面通过roi来裁切。region of interest是经纬度范围ISCE会在这个范围内生成干涉图能够省掉很多不必要的数据量。另外还有一个非常关键的参数master和slave的SLC路径。在topsApp.xml中要用完整的SAFE文件路径或者SLC目录路径。我习惯把两个SLC目录放在SLC文件夹下配置文件里直接用相对路径这样工程可迁移性更好。3.2 从startup到interferogram的分步执行topsApp.py支持一次跑完全部步骤也可以逐步运行。我强烈建议分步执行因为中间任何一步出错你都能准确定位到是哪里出了问题。示范命令如下topsApp.py --startup --endstartup topsApp.py --startup --endbursts topsApp.py --startup --endmerge topsApp.py --startup --endfilter topsApp.py --startup --endunwrap topsApp.py --startup --endgeocode第一步startup主要做数据解包、轨道状态计算和基线估计。如果这里就报错多半是SLC路径写错、SAFE格式不对或轨道文件缺失。这是整个流程里最容易出现问题的地方也是最好排查的地方。第二步bursts阶段会把TOPS模式的burst数据逐段处理包括方位向配准和共轭相乘生成单burst干涉图。这一步计算量很大而且对配准精度要求极高。如果这个阶段输出很多“burst misregistration”之类的警告就要回去检查轨道文件或配准参数。第三步merge把所有burst的干涉图拼接成完整的干涉图。拼接前需要检查相邻burst重叠区域有没有相位跳变。如果跳变严重说明burst级配准不够好需要回头调粗配准窗口大小或提升过采样率。第四步filter是滤波。ISCE默认用Goldstein滤波主要目的是压制干涉图中的残余噪声同时保留形变条纹。滤波强度通常设置在0.4到0.8之间。对于山区形变场我一般取0.6。滤波太强会把细节磨掉太弱则噪声大对后面解缠不友好。第五步unwrap和第六步geocode后面专门讲。受限于篇幅我在这里只加了最常用的几行命令实际ISCE还会生成很多中间产品包括coherence、wrapped phase、filtered phase等。每次跑完filter之后记得先打开幅度图和滤波后的干涉图看看确认条纹是否连续、相干性是否足够再决定要不要继续解缠。这一步肉眼检查能省掉不少返工时间。3.3 多视、滤波和干涉质量核验多视是InSAR处理里提高信噪比的常规操作。Sentinel-1的原始分辨率和像素间隔在距离向和方位向差别很大如果直接做干涉会得到一副横向很长、纵向很窄的图像。多视的目的有两个一是让方位向和距离向的地面分辨率相近二是降低相位噪声。对于Sentinel-1 IW模式一个常见的多视组合是range looks8、azimuth looks2这样大约能把像元变成几十米量级。如果你想要更高的形变场空间分辨率也可以用41的组合但噪声会明显增加。尼泊尔地震这种地表形变信号强的案例我用的是82处理速度快形变场平滑度也很好。多视参数在topsApp.xml里设置。但要注意多视参数会影响干涉图分辨率进而影响相位解缠的可靠性。一般经验是目标区域形变梯度大需要保持较多条纹细节就少做多视如果目标是区域尺度的整体形变多视可以多一点。滤波后必须做一次干涉质量核验。打开coherence图看目标区域相干性是否整体大于0.3。尼泊尔山区植被覆盖比较多如果相干性太低说明两景之间的时间去相干严重这时候需要尝试缩短时间基线或者换一个轨道方向的数据。刚才提到的多视和滤波参数也需要根据相干性来做微调。有一个小技巧在merge之后、filter之前先看一下未滤波的干涉图和幅度叠加图确认形变条纹的走向和区域大致符合地质构造背景。尼泊尔地震是逆冲型地震视线向形变场会根据卫星观测几何出现明显的抬升区和下沉区。如果你看到的条纹杂乱无规律大概率不是地震信号而是处理问题。4. 从干涉图到同震形变场解缠与地理编码4.1 相位解缠的原理和关键参数InSAR的干涉相位是主辅影像相位差经过模2π包裹的结果如果形变量超过一个雷达波长的一半相位就会在2π周期内来回跳变。Sentinel-1的C波段波长约5.6厘米同震形变场往往有几分米到几米的位移所以必须把缠绕的相位解缠成连续的绝对相位这一步就是相位解缠。ISCE默认支持两种解缠方法snaphu和最小费用流。snaphu是比较经典的选择适合大多数场景。在topsApp.xml里设置解缠方法component nameunwrap property nameunwrapping method valuesnaphu/ property nameuse coherence valuetrue/ property namesnaphu cost mode valuesmooth/ /component这里有两点经验。第一使用coherence加权能让解缠算法避开低相干区域我建议开启。第二如果目标区域有很明显的低相干区比如水体或强植被区可以先用相干性阈值生成掩膜让解缠算法忽略这些区域否则它们会像断路一样把空间连续相位割裂开。解缠是个计算密集过程如果区域范围很大内存占用会比较高。我通常会把region of interest缩小到震中周边2到3度范围内比如经纬度范围84°E到88°E、27°N到30°N这样既能覆盖主要形变带又不会让snaphu陷入几小时的等待。解缠完成后会有unwrapped phase的GeoTIFF输出。肉眼检查时正常解缠结果应该是一个平滑的相位面没有明显的“台阶”。如果你看到一大片区域相位值异常跳变很可能就是解缠失败。4.2 地理编码与LOS向形变场输出解缠得到的是雷达坐标系下的相位图。为了把它和实际地面位置对应起来必须做地理编码。ISCE会根据DEM信息把每个雷达像元映射到经纬度坐标同时生成GeoTIFF。这一步可以理解为把雷达坐标系下的结果投影到标准地图坐标系中。地理编码时ISCE会用到之前准备的DEM。拿到结果后通常还会把它转成标准的位移量单位是米。相位到位移的转换公式很简单就是相位除以相位到形变的换算系数import numpy as np import gdal phase gdal.Open(filt_topophase.unw.geo.tif) phase_data phase.GetRasterBand(1).ReadAsArray() disp_los phase_data * 0.028 # Sentinel-1 C波段相位到形变的系数约2.8cm即2π对应2.8cm这里要注意0.028这个系数是针对单路径干涉的近似值。更严格的做法是根据雷达入射角逐个像元计算但在误差允许范围内很多处理流程直接用这个近似系数。如果后续要做形变场解释最好还是把LOS向位移投影到垂直和水平方向那就需要更精细的几何参数。地理编码后的结果通常会保存成unw_filt_phase_geo.tif或类似的名字。为了后续可视化我一般会再转成标准的uint16或float32的GeoTIFF并设置缺失值为NaN。这样导入GMT或Matplotlib时不会因为NoData值导致色带错乱。4.3 形变结果的整饰与解读拿到LOS向形变场首先要做定性检查。在地图上叠加震中位置和主要断层迹线看形变场符号是否合理。对于喜马拉雅主逆冲断裂地表抬升区和下沉区的空间分布应该与断层倾角以及卫星侧视方向一致。我常用GMT来做后处理gmt makecpt -Cpolar -T-1.0 1.0 -H los.cpt gmt grdimage los_displacement.grd -Clos.cpt -R84/88/27/30 -JM15c -Ba1 -P nepal_los.ps如果没有GMT用Python画也是一样的。关键是色表选择应该让零位移区域落在中性色上抬升和下沉用明显对立的颜色。形变场图的纵轴最好标注为“LOS displacementm”并注明卫星轨道方向因为升轨和降轨看到的形变是不完全一样的。从物理角度讲同震形变场可以反映破裂几何。尼泊尔地震的InSAR形变场在震中以南应该有明显抬升在北部有相对下沉或者不显著这与逆冲型地震的模型预测是吻合的。如果你的结果趋势和这个大致一致说明处理流程基本靠谱。如果你手头有地表位移积分等额外资料可以把InSAR结果和GPS做一个沿视线向的对比。这一步能快速验证数厘米到分米量级的形变场是否准确。2015年尼泊尔地震周围虽然GPS站点少但仍有几个永久站可以辅助验证。5. 实战排错那些ISCE处理中常踩的坑5.1 下载失败和路径问题做InSAR第一步就容易被下载和数据整理绊住。哨兵SLC数据动辄几个GB从欧空局的数据平台下载时如果网络不稳定经常会断。建议用支持断点续传的下载工具或者直接在ASF数据搜索页下载它可以自动匹配轨道。ISCE对文件路径特别敏感。路径里最好不要有中文、空格和特殊符号。我曾经因为工程目录名字里带了个“-”导致GNU make解析出错找了好久才发现。更常见的是把轨道文件目录路径写错ISCE解包时找不到对应的EOF文件直接报错退出。还有一个坑下载的SLC数据如果是zip包ISCE一般能自动解包但最好先检查zip包完整性和校验和。如果用ASF的bulk download有时候会下载成小体积损坏文件一跑就报CRC错误。先把这个关卡住能节省不少时间。5.2 配准与干涉条纹异常TOPS模式的干涉处理最关键的是burst之间的配准。如果配准误差超过千分之一个像素merge后的干涉图就会在burst边界出现明显的条带。ISCE在测试配准时会输出偏移估计值如果你看到某个burst的偏移和相邻burst差异很大要考虑是不是轨道辅助文件太粗或者DEM误差太大。干涉条纹异常还有一种情况你在滤波后的干涉图上看到很多密集平行条纹但这些条纹和地震形变方向无关。那可能是残余地形相位。如果DEM有洞或者高程误差大需要补充DEM数据或者在做干涉时用双差法剔除地形相位。ISCE中的topophase步骤就是负责去除地形相位的如果这一步有误后续解缠就会输出错误形变。对于尼泊尔这种高地形起伏区我建议在跑topsApp之前先把DEM和雷达模拟图像叠一下检查DEM覆盖和坐标是否与SAR影像匹配。如果DEM数据范围不对生成的合成相位图也会错干涉相位就会在地形陡峭区域出现大量残差条纹。5.3 解缠跳变与后处理修复解缠跳变非常常见尤其是在低相干区和形变梯度大的地方。具体表现是unwrapped phase图上出现一片一片的相位“孤岛”相邻区域相位的整数周期对不上。这种情况即使重跑snaphu也不一定有效更适合用掩膜和加权策略来解决。一个比较稳妥的做法是生成相干性掩膜把相干性低于0.2的区域直接排除在解缠范围之外gdal_calc.py -A correlation.geo.tif --outfilecoh_mask.tif \ --calcA 0.2 --NoData0然后在配置文件中启用掩膜。这样解缠算法会跳过峡谷、水体、强植被等低相干区域减少大片跳变的概率。代价是掩膜区得不到形变值但很多时候这些区域本身不可靠留白反而更诚实。如果解缠结果中只在局部有个别跳变也可以通过后期中值滤波或者图像修复法修正。实际项目中我更喜欢重新调整多视和滤波参数再解缠一次。把时间花在源头参数上比后期一遍遍修图强得多。6. 一点实战心得与建议做ISCE处理最大的感受就是不要被命令行吓住。topsApp.py整个流程虽然长但每一步都有明确的数据产品和肉眼可查的中间结果这让我在调试时能很清晰地判断问题是出在配准、地形去除还是解缠。相比黑盒式的商业软件ISCE给了我更多掌控感。最后分享一个我自己常用的习惯每跑完一步我会把关键中间产品复制到一个固定的checkpoints目录里并写下这个环节的参数和备注。比如“2015-04-13 / 2015-04-25IW VV8:2多视滤波0.6snaphu”。这样过几个月回来看还能立刻回忆起当时的处理状态。尤其当你要处理多轨道、多时段数据时这个习惯能帮你减少大量重复劳动。如果你也是第一次跑2015年尼泊尔地震的InSAR流程我建议先用一个小一点的裁切区域试通全流程再扩大到完整形变场。先把流程跑通比一次追求完美结果重要得多。等到整条链路熟练了再回头细调参数你会发现自己对InSAR的理解已经上了一个台阶。