这几年做地表形变监测应该绕不开InSAR和ISCE。ISCE这个软件包圈内习惯叫它科学计算环境专门用来处理干涉合成孔径雷达数据尤其是哨兵一号的TOPS模式处理起来很顺手。不过大多数人熟练的是用ISCE出视线向形变图也就是LOS向的形变真正要绘制“三维形变”——把东向、北向、垂直向三个分量都标在图上很多新手会卡在流程设计这一步。本文就从我实际跑过的数据出发聊聊怎么把ISCE当成前端的干涉处理工具再用脚本把三维形变场解算并绘制出来内容覆盖原理、流程、参数和踩坑适合正在学习InSAR或者准备做区域形变监测的朋友参考。我自己第一次接触三维形变时也以为ISCE能直接出一个矢量图结果发现根本不是这么回事。三维形变解算的难点不在于“画”而在于如何把多个轨道、多种几何的观测拼在一起构造出一个可解的方程组。ISCE擅长的是给你精确的干涉图、解缠相位和雷达几何信息后面的矩阵解算和可视化需要自己补上。这篇文章会把这条链路完整拆开每一步都说明为什么这么干以及有哪些地方容易翻车。1. 为什么单一视线向形变满足不了工程需求1.1 InSAR到底测的是什么形变InSAR通过两次雷达观测的相位差提取出地表沿雷达视线方向的位移。这个视线向也就是LOS方向是一个从地面指向卫星的方向。雷达从侧视方向照射地面所以它对“朝着卫星靠近”或者“远离卫星”的位移特别敏感但对沿卫星飞行方向的分量基本不敏感。你可以把它想象成手电筒照着一根立着的棍子棍子在光束方向上的长短变化最明显而左右水平移动几乎看不出影子变化。InSAR测量到的LOS位移就是地表东、北、垂直三维位移在视线方向上的投影。问题在于这个投影是一个标量不是矢量我们拿到的LOS形变图中往往把垂直沉降和水平移动混在一起。比如地面发生滑坡时坡体既往下沉又向河谷方向滑动单看LOS图很容易低估真实移动量甚至符号判断错误。更重要的是很多灾害监测场景需要区分东西、南北和垂直分量。矿山开采沉陷、冻土区地面变形、城市地面沉降、滑坡长期蠕变这些过程的力学机制完全不同。垂直分量适合判断水量、压实、构造应力水平分量则与剪切、滑移、侧向流动更相关。如果只提供一张LOS图地质工作者、安全监测人员很难直接使用。这也是三维形变近些年越来越受重视的原因。1.2 三维形变解算的数学基础要解出三个未知数理论上至少需要三个独立的观测方程。但一条轨道的InSAR只提供一个LOS向观测值所以单轨道数据根本无法实现三维解算。常规思路是使用不同轨道几何的数据比如升轨和降轨的LOS形变再结合GPS、光学像素偏移或方位向偏移跟踪等额外信息构建超定方程组。假设某像素的三维位移为 d_e、d_n、d_u对应的LOS形变为 d_los那么存在这样的线性关系d_los a_e·d_e a_n·d_n a_u·d_u其中 a_e、a_n、a_u 是雷达视线单位矢量在东、北、垂直方向上的投影系数它们由卫星航向角和入射角共同决定。ISCE在干涉处理时会输出描述每个像素雷达几何的文件包括入射角、航向角这些数据可以用来计算上面的系数。很多教程喜欢直接给出固定公式但不同处理软件、不同坐标约定下符号定义有差异建议从ISCE的geometry文件中读取角度而不是抄公式了事。升轨和降轨的卫星飞行方向一南一北视线方向差别明显因此两个LOS方程可以提供一定的约束。但这两个方程在南北方向的敏感度都不高法矩阵在这个方向上是病态的。也就是说即便解出了结果南北向分量的误差也会被放大到离谱的水平。解决问题的常见办法有两种一是假设南北向位移近似为零但这只在某些构造背景下成立很容易引入系统性偏差二是引入额外观测比如用GPS网形变数据作为弱约束或者在ISCE里做方位向偏移跟踪提取轨道方向的位移场然后一起参与解算。后面这种方式更稳健也是我比较推荐的方案。2. ISCE在整个流程里的角色与工作流设计2.1 ISCE的看家本领ISCE的核心能力集中在干涉处理链路从原始SLC数据开始完成配准、干涉、去平地、滤波、相位解缠和地理编码。它对哨兵一号TOPS模式的支持尤其好Burts拼接和相位残留校正都做得比较细处理结果在科研和工程项目中认可度很高。ISCE还有一套完整的几何产品输出例如每个像素对应的入射角、航向角、高程和雷达成像几何信息这些非影像数据恰恰是后续三维解算最需要的东西。很多人把ISCE想象成“一键式”软件其实它是模块化的。你可以通过配置文件控制每一步是否执行也可以把中间结果输出成标准GeoTIFF方便和GDAL、GMT、Python等工具衔接。我个人的体感是ISCE不是最傻瓜的软件但它的透明度和可定制性是最好的。尤其是在做三维形变时你需要清楚自己的数据经过了怎样的算法处理而ISCE每一步都有日志和输出文件排查问题相对容易。2.2 为什么三维解算不能靠ISCE一键完成到目前为止ISCE官方并没有提供一个“输入多条轨道数据、输出三维形变场”的终极按钮。原因也很务实三维解算需要的输入千差万别有人手里只有升降轨数据有人额外有GPS有人还想把像素偏移结果融进去不同案例的权重策略不一样解算结果的正负号约定也不统一。把这些全部塞进一个通用模块维护成本高还容易误用。所以实际工作流通常是这样的用ISCE批量生成每个轨道几何的LOS形变图、入射角和航向角然后统一到公共地理网格再用Python或其他脚本进行多方程联合解算最后绘图。ISCE在整个流程中的角色更像“精密的测量仪器”三维解算则是基于测量数据的“后期建模”。这种拆分方式有个好处你可以在任何一步替换或增加数据源而不需要推倒重来。举个例子我常用的流程是脚本化地跑ISCE每个干涉对输出到一个独立目录目录里统一放置los_ll.grd、unwrap_ll.grd和geometry文件。之后用一个名为“三维解算”的Python脚本读取这些文件自动构建设计矩阵生成东、北、垂直形变图。修改某个轨道的输入时只需要重新处理那一个目录解算脚本不用变。你如果硬要把所有步骤粘在一起改起来反而束手束脚。3. 核心操作流程从原始数据到三维形变场3.1 数据准备选什么数据怎么组织三维形变解算最常用的数据源是哨兵一号重访周期短、数据免费、TOPS模式覆盖范围大ISCE对它的支持也最完善。数据准备阶段就要考虑几何多样性尽量同时准备升轨和降轨的数据轨道方向差得越明显对东西和垂直分量的约束越好。如果有多个不同视角的降轨或升轨也可以一并使用但要保证这些干涉对在空间和时间上尽量接近同一时段否则形变本身随时间变化解算出来就成了“混合时间窗口”。下载原始SLC后还需要准备精密轨道数据和对应的DEM。ISCE在干涉处理时对轨道的精度很敏感不要偷懒使用原始播发轨道尽量下载精密轨道文件。DEM方面一般用覆盖研究区的SRTM或ASTER数据事先切成与研究区匹配的区域。把所有文件放在有规律命名的目录中比如asc/desc分开放每个干涉对单独一个子目录。这个习惯会在后面批量处理时省下很多时间。个人建议维护一个数据清单表第一列是干涉对名称第二列是轨道号第三列是主影像日期第四列是辅影像日期第五列是时间基线第六列是空间基线。每次处理之前先扫一眼像升轨、降轨、不同入射角的数据是否都覆盖到了避免处理到最后发现某个方向缺数据白跑一遍流程。3.2 用ISCE运行干涉处理ISCE处理TOPS数据的入口是topsApp.py。你需要在配置文件中写入主影像、辅影像、DEM、极化方式、是否滤波、是否解缠等设置。不同版本配置略有差异但大体思路一致。我粗略总结一个典型的配置逻辑reference/secondary两个组件填写原始SLC数据路径DEM字段填写裁好的高程文件unwrap选项设为Truefiltering设置为适度坐标系和输出路径指定清楚运行命令时可以分步执行便于检查中间结果。比如先跑配准和干涉查看强度图和干涉图是否正常再跑滤波和解缠。一次跑完整条流水线确实省事但一旦中间某个步骤出错整条链路都得重来。我建议新手拆成几步执行topsApp.py --dostep start topsApp.py --dostep unwrap topsApp.py --dostep geocode每一步都可以在对应目录下检查输出文件。比如干涉阶段重点看coherence和interferogram相干性如果大面积低于0.3就要考虑是不是基线太长、时间去相干或者数据本身有问题不要急着解缠。ISCE处理完成后主要输出包括解缠相位图、相干性图、视线向单位矢量信息、地理编码后的相位图层等。需要记住这里的解缠相位单位是弧度要换算成以米为单位的LOS位移公式大致是d_los -λ · φ_unw / (4π)其中λ是雷达波长φ_unw是解缠相位。符号取决于解缠算法的约定不同ISCE版本可能不一样。最稳妥的方法是用一个已知地面位移的区域进行验证比如某个有GPS连续站的区域对比LOS投影方向上的位移确认正负号正确后再进入批量流程。3.3 提取LOS形变并统一网格在ISCE生成结果后需要把雷达坐标系下的数据统一转到公共的地理网格上。ISCE本身有地理编码功能可以输出GEOCoded的干涉图和解缠相位。但注意不同轨道、不同框架的地理编码结果可能网格大小、范围不完全一致需要用一个统一的参考网格把所有数据重采样到一起。最简单的方法是选定一个网格基准比如研究区中心点坐标和0.0001度分辨率使用GDAL或Python的栅格重采样工具将所有LOS形变图、入射角、航向角都插值到这个网格上。重采样时注意选择合适的插值方法形变数据用双线性插值即可相干性或者标志位数据用最近邻插值避免引入伪值。同时ISCE输出的los文件通常包含视线向量在雷达坐标或地理坐标下的分量以及入射角和航向角。这些角度最终会用于构建设计矩阵。我建议在这一步把所有栅格统一命名例如los_e.tif、los_n.tif、los_u.tif、inc_asc.tif、az_asc.tif这样。整理好命名后面写脚本就不用每个文件单独判断含义了。这里有个容易忽略的地方入射角和航向角可能以弧度为单位也可能以度为单位务必统一。我踩过一次坑把弧度的入射角当成角度直接代入公式结果解算出来的垂直形变被放大一度让人怀疑是算法写错了最后才发现是单位问题。建议在读取文件时显式检查数据范围度数应该在0到90之间弧度则更小。3.4 构建三维解算方程与加权最小二乘完成多轨道数据的统一网格化后就可以进入解算环节。假设每个像素上有N个观测值每个观测值对应一个LOS方向。注意这里每个观测值可能是升轨LOS、降轨LOS、GPS垂直分量或者方位向偏移跟踪结果。对每个像素构建观测矩阵A第一列该观测对应的东向投影系数第二列该观测对应的北向投影系数第三列该观测对应的垂直向投影系数对应观测向量d例如升轨LOS位移、降轨LOS位移、GPS垂直位移。然后求解最小二乘问题。考虑到不同观测的精度差异可以引入权重矩阵W通常取观测误差方差的倒数作为权重。Python伪代码如下import numpy as np # 假设已有数组 # los每个观测的形变值 # a_e, a_n, a_u投影系数矩阵形状为 (N,) # weights每个观测的权重 A np.stack([a_e, a_n, a_u], axis1) W np.diag(weights) # 加权最小二乘 x np.linalg.inv(A.T W A) A.T W los如果矩阵出现接近奇异的情况还可以加入正则化项例如岭估计lam 1e-3 M A.T W A lam * np.eye(3) x np.linalg.inv(M) A.T W los这一步需要特别注意GPS垂直位移作为伪观测时对应的A行是[0, 0, 1]观测值为垂直位移。如果只用升降轨LOS没有外部约束法矩阵往往是奇异的解算出的南北分量会严重震荡。因此我建议要么加入GPS作为弱约束要么至少用“南北向位移近似为零”的约束。弱约束的意思是把GPS观测当成方差有限的方程放进正规方程而不是硬性赋值。实际操作中可以把权重设置为GPS误差倒数平方比如5毫米的GPS误差权重就是200的平方注意这里权重应该是1/σ²需要仔细。此外不同轨道的LOS观测噪声其实分布不均匀最直观的是看相干性。相干性越低相位噪声越高。我建议以相干性作为权重依据相干性高于阈值时给予高权重低于阈值时直接掩膜。这样可以显著提升解算效果。3.5 结果绘制与输出解算完成后就得到东向、北向、垂直向三个形变分量栅格。绘制时基础选项包括使用统一的裁剪范围、统一色标、标注正负号含义。我个人习惯用Python的matplotlib或者GMT出图色标用Blue-Red对比明显。因为垂直和水平分量的量级可能相差较大不要强行共用同一个色标建议分开设置以免大值把细节压没。如果关心形变矢量场还可以在垂直形变背景图上叠加水平方向箭头箭头长度代表东西南北的合成水平位移方向就是位移方位角。这一步看似简单但需要注意地图投影和坐标轴尺度比例尤其在高纬度地区经纬度每度的实际距离差异很大直接用经纬度坐标画箭头会得到歪曲的矢量。更稳妥的做法是先使用 pyproj 把经纬度转换成等距离投影坐标再计算箭头分量。输出文件建议保存为GeoTIFF同时附带一个元数据文件记录解算所用到的观测轨道、时间范围、解算方式、采样网格和权重设置。三维形变结果如果用于报告或论文元数据可以帮你快速回溯流程也便于后续复算。别相信自己的记忆隔三个月再看数据你一定会忘记当时怎么处理的了。4. 参数设置与精度控制4.1 滤波和解缠参数怎么选ISCE干涉处理中滤波参数和解缠算法对最终形变场影响很大。滤波太强形变信号会被平滑掉滤波太弱相位噪声把真实信号淹没。就哨兵一号陆地场景而言我常用Goldstein滤波滤波窗口在32到64之间滤波系数alpha设为0.4到0.6。如果你研究的是范围大但形变梯度缓的沉降场可以适当加大滤波窗口如果研究滑坡、断层等形变梯度大的区域滤波就要保守一点。相位解缠算法方面ISCE内置多种算法常用的是最小费用流或者区域生长。对三维解算来说重要的是解缠结果不能有大面积跳变。检查方法很简单把解缠相位图在空间上做剖面看是否连续或者生成解缠质量图观察低可信度区域是否连片。如果解缠结果出现长条状异常多半是森林、陡峭山区等低相干区域引起的这时可以设置相干性掩膜只解缠相干性较高的像素但注意不要造成数据空洞。另外对多轨数据解缠时要保持同一套掩膜标准。如果升轨数据在某个区域能解缠降轨数据却被掩膜掉最后解算时那个区域就只剩下单方向观测三维解算可信度会大幅下降。这种情况下要么采集新的数据要么在解算脚本中将观测数低于2的像素识别出来并置空不要硬插值。4.2 入射角和航向角设计矩阵的关键设计矩阵的投影系数由两个角度决定入射角和航向角。入射角是雷达波束与垂直方向的夹角通常在20度到45度之间航向角是卫星飞行方向与北方的夹角升轨和降轨约相差180度。需要注意不同ISCE版本输出的角度定义可能有所不同有时los.rdr文件里直接给出三个方向的单位向量分量有时则给出入射角、航向角需要自己换算。我建议从geometry中直接读取单位向量分量而不是自己用角度算。因为软件已经考虑了地形、地球曲率等修正比手动计算更精确。很多第三方教程喜欢展示一个统一的公式但公式的前提是雷达视线几何标准化遇到山区地形起伏时入射角会受局部坡面影响单位向量每个像素都不同。这恰恰是ISCE的优势所在——它可以逐像素输出真正的视线单位矢量直接在三维解算中使用。验证设计矩阵是否正确可以用一个简单的办法对单个轨道数据把单位向量和LOS形变做点乘还原看能不能得到原始的LOS形变值。还原结果如果相差很大说明单位向量或LOS形变的符号、单位有问题需要回到上一步修正。这一步我每次都会跑三分钟能节省后面三天排查。4.3 加权方案与病态问题三维解算的病态问题主要来自南北向弱观测。即便有升轨和降轨数据南北向的投影系数仍然很小直接最小二乘会让南北分量变得离谱。加入GPS约束后方程条件数会改善但也有GPS站点稀疏、插值误差大的情况。这时候需要合理设置权重GPS约束权重不能太高也不能太低。太高会导致结果完全被GPS插值控制失去InSAR本身的空间分辨率太低则起不到约束作用。一个工程上的实用做法是把GPS垂直位移当作“软约束”的同时对南北方分量施加空间平滑正则化。比如在解算方程后面添加一个平滑惩罚项让南北向形变在空间上不过度突变。这样既能利用InSAR的大范围空间信息又不会让异常单点噪声污染结果。具体平滑项可以用Laplacian算子不过要考虑分辨率尺度算起来不复杂。误差传播分析也很重要。解算得到的东、北、垂直分量各自带有不确定度最好输出不确定度栅格。不确定度可以通过协方差矩阵对角线计算。很多三维形变图只画形变值不画误差结果看似精美实际上分量的精度完全未知。我习惯在出图时把不确定度大于阈值的区域设为透明比画满色彩更有说服力。5. 常见问题与排查实录5.1 升轨和降轨解算的垂直形变符号相反这个问题我遇到过不止一次。现象是垂直分量图看起来像是“反转”了沉降区变成抬升区。多数情况下不是解算矩阵写错了而是某个轨道的数据符号约定错了。ISCE输出的LOS形变可能定义为正值表示远离卫星也可能反过来。不同版本、不同处理链路的符号并不统一。排查时不要直接改符号先用一个已知垂直形变明显的区域检查如果该区域确实在下沉那视线向位移应该是负值若约定朝向卫星为正。对照ISCE的输出看看符号能否对应。若不能对应就在读取该轨道数据时统一调整符号。切忌每轨道各调各的否则解算结果会更乱。5.2 解缠跳变导致三维形变呈条带状三维形变图出现条带状、沿方位向分布的异常往往是某个干涉对的解缠相位出现了条带跳变。这种情况在山区和城市边缘的高层建筑区特别明显。跳变的值通常接近整周数所以看似连续但投影到三维解算中会变成大面积假信号。排查方法是先不急着做三维解算把每个轨道的LOS形变单独成图沿轨道方向做差分。差分值如果出现规律性跳变说明解缠有问题。可以回到ISCE重新设置掩膜或者在结果后处理中检测相位跳边并剔除以该像素为核心的观测。三维解算中遇到条带状明显的数据宁可用观测数少但干净的像素不要保留含有解缠误差的伪观测。5.3 地理编码后不同轨道数据对不齐在把升轨和降轨数据重采样到公共网格后经常发现边缘或山区数据有几十个像素的漂移。这可能是不同轨道地理编码使用的DEM精度不一致也可能是重采样网格投影参数不一致。虽然看起来只是零点几个像元的差但在高分辨率形变场中这种错位会直接引进大到厘米级的误差。解决的办法有两个方向一是所有数据统一到同一版本DEM和同一投影参数不要混用不同地理编码结果二是在重采样前先做一次互相关配准利用强度图计算轨道间偏移量手动校正后再重采样。ISCE对于每一条轨道的地理编码定位通常已经比较准确但如果你发现两轨数据的河流、道路边缘错开就需要校正后再进入解算。5.4 快速排查速查表现象可能原因建议解决方式垂直分量符号反转LOS或单位向量符号不统一用已知形变点验证符号南北分量出现大尺度条纹缺少额外约束引入GPS或方位向偏移约束形变场边缘出现断裂不同轨道解缠掩膜不一致统一相干性阈值与掩膜处理东西分量与水平GPS不符入射角、航向角单位错误检查角度单位和投影系数整幅形变场噪声过大未加权或权重不合理按相干性或误差倒数加权三维解算结果空洞太多升降轨重叠范围小降低遮挡要求或补充轨道数据每个人跑三维形变解算时都会遇到自己数据特有的一堆怪问题。遇到问题不要一张张图手工改而是把中间结果全部堆在同一个目录下用批处理脚本统一生成检查图。把检查图按照“观测值-投影系数-中间解算结果”三个层次组织一旦某个结果不对很快就能定位到是输入、几何还是解算的问题。我后来把这条经验固化成了固定模板稳定跑了很多实验大幅减少了重复劳动。最后再分享一个心得刚开始研究三维形变时我总想一下子把南北、东西、垂直三个分量全做出来结果解出来的南北分量一团糟垂直分量也不对劲。后来我改成了“先求垂直再求水平”的两步策略先垂直确定下来再求解东西和南北。这个思路虽然在严格数学意义上只是变换了解算顺序但实际操作中很有帮助你能更早发现哪个方向的分量不可信。另外绘制的三维形变图要记得标注清楚投影坐标系和时间基线不然成果交出去之后别人根本不知道这个形变对应的是哪两幅影像之间的变化。地形变监测这件事细节决定结果到底能不能用希望这篇内容能帮你少走一些弯路。
