SAR 影像与 InSAR 形变监测

引言

光学遥感依赖可见光与近红外，遇到云雨和夜间就无能为力；合成孔径雷达（SAR）用微波主动成像，全天时全天候工作，还能借干涉相位测量毫米级的地表形变。这两点让 SAR 成为地质灾害普查、城市地面沉降、矿区监测与基础设施健康诊断的主力数据源。

工程上的难点不在原理，而在误差链条太长。干涉相位里混着地形相位、平地相位、大气延迟、轨道误差、解缠误差与失相关噪声，任何一环处理不当，最终形变结果就会出现条带、跳变或整体偏移，而且往往肉眼难以察觉，必须靠残差图与相干性掩膜交叉验证才能发现。

另一个痛点是数据量与时序规模。Sentinel-1 单景 SLC 约 8GB，一个沉降监测区动辄几百景，加上配准、干涉、解缠、时序反演，单机根本跑不动，必须切块并行并接入任务编排。

此外，SAR 结果是视线向（LOS）的一维投影，不是垂直形变，也不是水平位移。很多人拿到速率图就直接当沉降量使用，忽略了入射角投影因子，误差可达百分之几十。理解几何投影是正确解读结果的前提。

本文按「成像几何 → 干涉相位 → 误差源 → 时序反演 → 工具链落地」的顺序展开。先讲清每个物理量从哪里来、怎么算，再给出可直接复用的参数取值与命令片段。SAR 影像的基础格式与波段约定可先看 遥感影像基础与波段组合 ，栅格读写与投影对齐问题见 遥感坐标参考系与投影 。

目录

1. SAR 成像几何与侧视观测
2. 散射机制与斑点噪声
3. 极化与干涉测量基础
4. 干涉图生成与去平地效应
5. 相干性来源与失相关分析
6. 相位解缠算法与质量评估
7. DInSAR 与形变时序反演
8. 大气延迟与轨道误差改正
9. 工具链与 Sentinel-1 数据流程
10. 权衡取舍
11. 常见坑清单
12. 小结

1. SAR 成像几何与侧视观测

SAR 是侧视雷达，天线斜向下照射地面，沿飞行方向称为方位向（azimuth），垂直于飞行方向称为距离向（range）。因为是斜距成像，同一条距离线上的地物会被压缩到同一斜距，于是产生三类几何畸变：透视收缩、叠掩与阴影。

透视收缩指迎坡在影像上被压缩，坡度越接近入射角余角，压缩越剧烈；当坡度大于入射角余角时进入叠掩，坡顶与坡底回波顺序颠倒，影像重叠无法解译；背坡则可能完全无回波形成阴影。这三类畸变在山区的形变监测中会直接造成假信号，必须在解译时用地形与入射角掩膜剔除。

分辨率与像元尺寸

方位向分辨率由合成孔径长度决定，卫星运动形成等效长天线，Sentinel-1 的方位向分辨率约 20m；距离向分辨率由脉冲带宽决定，IW 模式约 5m，因此单视 SLC 像元是 20m × 5m 的矩形。多视后通常重采样为 20m × 20m 的方形像元再进入干涉流程。

分辨率关系
方位向分辨率 = 天线长度 / 2
距离向分辨率 = c / (2 * 带宽)
Sentinel-1 IW: 方位向 ~20m, 距离向 ~5m, 多视 4x1 后约 20m x 20m

波段选择

波长决定穿透与相干性，这是选波段的第一依据。

轨道类型与重访周期

轨道类型决定了干涉对的可用性。太阳同步近极地轨道保证同一地方时过境，光照与电离层条件稳定，是大多数 SAR 卫星的选择。轨道高度决定重访周期：Sentinel-1 单星 12 天，双星 6 天；ALOS-2 约 14 天；TerraSAR-X 约 11 天。

重访周期越短，时间失相关越轻，时序分析的可用点密度越高。这也是为什么 Sentinel-1 双星恢复后，植被区形变监测的可行性显著提升。

观测几何与升降轨

侧视成像示意
      卫星轨道
        |
        | 斜距 R
        v
   [叠掩区]  <- 坡面向雷达倾斜
   [透视收缩]
   [阴影区]  <- 背坡无回波
入射角 theta 越大，几何畸变越严重，通常取 30 至 45 度最稳

升轨与降轨的观测方向不同，同一区域的形变在两轨中投影方向不同，联合升降轨可分解出东西与垂直分量，这也是为什么单轨结果只能给出 LOS 形变。若只做升降轨联合仍无法分解南北向，因为 SAR 对南北向形变几乎不敏感，这是侧视几何的固有盲区。

2. 散射机制与斑点噪声

SAR 记录的是后向散射强度，用归一化散射系数 σ0 表示，单位 dB。不同地物的散射机制差异明显：裸土与水面接近面散射，植被冠层是体散射，城市建筑与角反射器是二次散射，表现为极亮点。理解散射机制才能解释为什么同一场形变在不同地物上相干性天差地别。

斑点噪声的统计本质

单视 SAR 影像天然带斑点噪声。原因是每个分辨率单元内分布着大量随机相位的散射体，它们的回波相干叠加，强度服从指数分布，视觉上就是颗粒状噪点。斑点不是设备缺陷，而是相干成像的固有属性，用任何滤波都无法根除，只能抑制。

多视与自适应滤波

抑制手段是多视与滤波。多视在距离向与方位向平均若干像元，等效视数（ENL）提升 N 倍，斑点方差按 1/√N 下降，代价是空间分辨率下降。若要保持分辨率，则用自适应滤波。

import numpy as np  # 仅演示多视平均原理
block = np.ones((4, 4), dtype=float)  # 4x4 像元窗口
multilooked = block.mean()  # ENL 提升 16 倍，分辨率下降 4 倍

工程上常用 Refined Lee 滤波而非简单均值，它在同质区平滑、在边缘处保留，能在不模糊地物边界的前提下压低斑点。干涉流程中更推荐多视加相干性掩膜，而不是对强度图重滤波，因为滤波会改变像元的统计特性，影响后续相干性估计。

• Boxcar：最简单，边缘模糊严重，仅用于快速预览。
• Lee：基于局部统计的自适应，通用性好。
• Refined Lee：改进的边缘保持，生产环境常用。
• 非局部均值：质量最高但慢，适合小批量精细处理。

3. 极化与干涉测量基础

极化 SAR 用 HH、VV、HV、VH 四种收发组合描述地物，不同极化对地物结构敏感度不同。HH 对水平结构敏感，VV 对垂直结构敏感，交叉极化 HV 主要反映体散射，是植被生物量反演的重要输入。

极化分解

极化分解（如 Freeman-Durden 三分量、Yamaguchi 四分量）可分离表面散射、二次散射与体散射功率，常用于农作物分类与地表覆盖制图。分解结果以 RGB 合成显示，红色常代表二次散射，绿色代表体散射，蓝色代表表面散射，一眼就能区分城区、森林与裸地。

干涉相位构成

干涉测量用的是复数影像。两景配准后的 SLC，逐像元共轭相乘得到干涉相位：

phi_int = phi_master - phi_slave
        = phi_flat  平地相位，由椭球与轨道决定
        + phi_topo  地形相位，与高程成正比
        + phi_def   形变相位，我们真正要的
        + phi_atm   大气延迟相位
        + phi_orb   轨道误差相位
        + phi_noise 失相关与热噪声

四项已知项必须逐一扣除，剩下的才是形变。这就是干涉流程的全部逻辑：把可建模的相位剥离干净，让残余相位只反映地表位移。

相位缠绕与灵敏度

需要注意的是相位是缠绕的，值域只在 [-π, π]，真实相位差可能是若干倍的 2π。一个 C 波段干涉对，2π 对应约 2.8cm 的 LOS 位移，超出这个量级就会出现相位跳变，这也是快速形变区必须依赖时序方法的原因。L 波段 2π 对应约 11.8cm，容忍的形变梯度更大，这也是它在矿区等大形变场景更受青睐的原因之一。

干涉相干性的数学表达

干涉相干性定义在两个配准影像的局部窗口内，用复数互相关归一化得到：

gamma = | sum(s1 * conj(s2)) | / sqrt( sum(|s1|^2) * sum(|s2|^2) )
其中 s1, s2 为主辅影像的复数值，sum 在局部窗口内进行
gamma 取值 0 到 1，越大表示两次观测相位越一致

这个定义解释了为什么低后向散射区相干性差：分母的能量项太小，噪声占比高，归一化后的互相关被拉低。工程上遇到大面积水面或光滑裸地时，相干性图会呈现出与地物高度相关的斑块，这属于正常现象而非处理错误。

4. 干涉图生成与去平地效应

干涉流程第一步是配准。主辅影像要在方位向与距离向亚像元级对齐，通常要求配准精度优于 0.1 像元，否则相干性会急剧下降。Sentinel-1 精密轨道（POD）文件能把轨道误差降到厘米级，务必在配准前下载最新的 AUX_POEORB。

配准的三个阶段

• 粗配准：用轨道与影像元数据估算初始偏移，误差可达数十像元。
• 几何配准：用外部 DEM 与轨道做几何精配准，降到亚像元。
• 增强谱分集配准：用 ESD 方法消除 burst 间的方位向配准误差，是 Sentinel-1 TOPS 模式的关键步骤。

去平地与去地形

第二步是去平地效应。平地相位由参考椭球与轨道几何决定，可用轨道状态矢量精确建模后扣除，这一步不依赖外部 DEM，因此是最可靠的一环。

第三步是去地形相位，需要外部 DEM。常用 SRTM 30m 或 TanDEM-X 12m，DEM 的垂直精度直接决定残余地形相位的大小。若 DEM 与影像时相相差太大（如冰川消融区），DEM 本身的高程误差就会引入假形变。

gdal_translate -of GTiff -co COMPRESS=DEFLATE dem_30m.tif dem_utm.tif  # 转换并压缩 DEM
gdalwarp -t_srs EPSG:32650 -tr 30 30 dem_utm.tif dem_geo.tif  # 重投影到 UTM 50N
gdalinfo dem_geo.tif  # 检查范围与分辨率是否覆盖影像

去地形后得到差分干涉图，此时残余相位主要由形变、大气与噪声构成。若 DEM 存在系统性高程误差，会在干涉图上留下与地形相关的长波条纹，这是最常见的误判来源。

干涉流程的完整步骤

1. 下载主辅 SLC 与精密轨道文件
2. 粗配准  -> 几何精配准 -> ESD 增强谱分集配准
3. 生成干涉图，同时输出相干性图
4. 去平地相位（轨道建模）
5. 去地形相位（外部 DEM）
6. 多视（按目标分辨率选倍数）
7. 自适应滤波（可选，谨慎使用）
8. 相位解缠
9. 地理编码（重采样到地图投影）

每一步都应保留中间产物并检查，尤其第 3 步的相干性图与第 5 步的残余地形条纹，这两个是后续所有误差的源头。

5. 相干性来源与失相关分析

相干性 γ 取值 0 到 1，衡量两景影像局部相位的一致性，是判断干涉质量的核心指标。失相关有四个来源，理解它们才能设计合理的观测策略。

• 时间失相关：地表在两次观测间发生变化，植被区尤其严重，C 波段植被区通常 12 天就几乎完全失相干。
• 空间失相关：基线过长导致两次观测的视角差异过大，C 波段垂直基线一般要求小于 300m，超过 500m 基本不可用。
• 体失相关：微波穿透到冠层内部，不同深度的散射体随机移动。
• 热噪声失相关：低后向散射区（平静水面、光滑裸土）信噪比过低。

相干性估计窗口

相干性用局部窗口估计，窗口越大估计越稳但空间分辨率越低。常用 5×5 到 15×15 的窗口，多视后按多视倍数缩放。窗口太小会高估相干性，导致低质像元被误纳；窗口太大则模糊了相干性的空间变化，把真正的低相干区平滑掉。

工程上通常以 0.3 作为相干性掩膜阈值，低于此值的像元不参与解缠与时序反演，否则会污染整个解算。对时序分析，还要进一步要求像元在时间序列上保持稳定的高相干，而不是偶尔一次高相干。

6. 相位解缠算法与质量评估

解缠是把缠绕相位还原为连续相位的过程。它本质上是一个二维积分问题，理论上有无穷多解，必须引入假设：相邻像元真实相位差小于 π。这个假设在陡峭形变梯度区或低相干区会被打破，导致整片区域的解缠错误。

主流算法对比

• 枝切法（branch-cut）：识别残差点并连线阻断积分路径，速度快但对残差密集区效果差。
• 最小二乘法：全局平滑，抗噪但误差会扩散到整幅图。
• 最小费用流（MCF）：在精度与鲁棒性之间平衡，是 SNAPHU 的默认算法，也是工程首选。

snaphu.csh config.snaphu  # 使用 MCF 算法解缠
snaphu -s -f config.snaphu -o unwrap.out phase.in  # -s 表示统计输出
snaphu -M 200 -C conncomp.out unwrap.out  # 基于连通分量做解缠

解缠质量评估

解缠结果必须做质量评估。最直接的方法是残差图：把解缠相位重新缠绕回 [-π, π]，与原干涉图相减，残差应接近零，出现 2π 整倍数的条带就说明该处解缠错误。另一个指标是解缠前后的相干性加权残差，可用来生成置信度掩膜。

MintPy 提供的 bridging 方法能自动识别并改正孤立的解缠错误连通域，但前提是错误区域是孤立的；若错误连成大片，仍会传播。实际项目中，解缠往往是整个流程最容易出错的环节，建议保留中间产物并逐景目视检查，不要盲目批处理。

解缠误差的传播机制

解缠误差一旦发生，会以 2π 整数倍的偏移形式污染整条积分路径下游的所有像元。这意味着一个孤立的低相干斑块可能让它后面的整片区域整体抬升或下沉一个或多个 2π。这也是为什么相干性掩膜如此重要：宁可丢掉低相干区，也不能让错误传播。

一个实用的检查手段是把解缠结果按行或按列做差分，正常区域差分应平滑连续，出现阶跃就说明该处有解缠错误。对时序堆栈，还可以比较不同干涉对在同一像元解缠结果的一致性，不一致的像元标记为可疑。

7. DInSAR 与形变时序反演

两轨差分干涉（DInSAR）只用主辅两景，能捕捉一次显著形变事件，比如地震同震位移或火山喷发，精度可达厘米级，但受大气与失相关限制，无法测缓变累积形变。

时序方法通过成百上千景影像联合反演，把形变从噪声中分离出来。两条主流路线是 PS-InSAR 与 SBAS。

时序反演流程

MintPy 是当前开源时序反演的主流，输入为解缠后的干涉图堆栈，输出形变速率图（mm/yr）与累积形变时间序列。它的处理步骤依次是：加载干涉图堆栈、剔除低相干干涉对、估计并改正解缠错误、改正大气分层延迟、建立形变模型反演速率。

mintpy:
  mintpy.load.processor: isce
  mintpy.subset.lalo: "30.2:30.6,104.0:104.4"  # 监测区范围
  mintpy.network.coherenceBased: yes  # 按相干性剔除低质干涉对
  mintpy.unwrapError.method: bridging  # 解缠错误改正
  mintpy.reference.lalo: "30.40,104.20"  # 稳定参考点
  mintpy.troposphericDelay.method: height_correlation  # 分层延迟改正

参考点与模型选择

时序反演的关键是参考点选择。参考点必须位于已知稳定的基岩或深基础建筑上，一旦参考点本身在沉降，整幅速率图会整体偏移，且很难事后发现。

形变模型也要选对：线性模型适合稳定沉降区，若存在季节性波动（如冻胀、地下水开采季节变化），必须加入周期项，否则速率会被系统性高估或低估。对矿区等非线性沉降，还需引入指数或双曲模型。

8. 大气延迟与轨道误差改正

大气延迟是 C 波段时序分析的第二大误差源，可造成数厘米量级的假形变。它分两部分：湍流混合造成的短波随机延迟，以及地形相关的分层延迟（高程越高，对流层水汽越少，相位越滞后）。

改正顺序与方法

分层延迟可用 GACOS 数据或 ERA5 再分析资料建模后扣除，能消除大部分与地形相关的条纹。湍流延迟靠时空滤波抑制，因为它在时间上不相关、在空间上高频，而真实形变是缓慢累积的。

轨道误差表现为长波条纹，通常用多项式拟合残差相位后扣除。残余轨道误差在 Sentinel-1 精密轨道支持下一般小于 2cm，但若用了预报轨道，误差会明显放大，务必检查所用轨道类型。

• 分层延迟：用 GACOS 或 ERA5 逐景改正，效果稳定。
• 湍流延迟：用时间高通加空间低通滤波抑制。
• 轨道误差：用二次多项式拟合后扣除。
• 电离层延迟：L 波段才显著，Sentinel-1 基本可忽略。
• DEM 误差相位：用高程相关项建模，与分层延迟一同估计。

改正顺序也有讲究：先改轨道误差（长波、确定性最强），再改分层延迟，最后用滤波抑制湍流。顺序颠倒会把形变信号误当作误差扣除。

9. 工具链与 Sentinel-1 数据流程

开源工具链已相当成熟。ISCE2 负责配准与干涉，SNAPHU 负责解缠，MintPy 负责时序反演，这套组合被广泛验证。SNAP 图形界面友好，适合小批量与教学；GMTSAR 脚本简洁，适合熟悉 GMT 的团队；HyP3 与 ARIA 提供云端处理，可直接下载成品干涉图，适合快速验证而不适合深度定制。

Sentinel-1 数据从 ASF 或 Copernicus 下载，IW 模式 SLC 覆盖 250km 幅宽，2021 年 12 月后星座恢复双星，重访周期回到 6 天，这对提高时间采样率、缓解时间失相关非常关键。

python3 -m pip install asf_search  # 安装数据检索库
python3 -c "import asf_search as a; print(a.search(platform='Sentinel-1'))"  # 检索可用景
isce2_topsApp.py --reference ref.SAFE --secondary sec.SAFE --dem dem.tif  # 一步干涉

存储与并行

处理前务必检查：轨道文件是否为精密轨道、DEM 是否与影像时空匹配、子区是否跨越多个 IW burst 边界。跨 burst 的干涉图需要特殊拼接处理，否则会在 burst 接缝处出现相位不连续。

若监测区跨越大范围，建议切块并行。切片时相邻块保留 20 至 30 像元重叠，避免边缘解缠误差传播到相邻块，拼接时用重叠区做一致性检查。堆栈存储上，单精度复数干涉图可用 GDAL 的 COG 格式压缩存放，比原始格式省一半空间。

gdal_translate -of COG -co COMPRESS=DEFLATE -co PREDICTOR=2 intf.tif intf_cog.tif  # 压缩干涉图

权衡取舍

核心权衡是分辨率与相干性、精度与覆盖、算力与时效。没有普适最优解，必须针对监测区地物与形变速率量级来定。快速形变区要优先保证相位不缠绕，缓变区要优先保证长时序的稳定性。

精度与成本的量化对照

精度提升的代价是算力与人力成倍增长。很多项目其实只需要厘米级的趋势判断，用云端成品就足够，不必上全套时序反演。

常见坑清单

• 形变结果出现规则长波条纹：DEM 高程误差或轨道误差未扣除，换高精度 DEM 并检查轨道类型。
• 整幅速率图整体偏移：参考点选在了沉降区，重新选取基岩或深基础参考点。
• 解缠结果出现 2π 跳变条带：低相干区相位不连续，提高相干性掩膜阈值并启用 bridging 改正。
• 相干性整体偏低：时间或垂直基线过长，筛选更短基线组合重做。
• 结果只给 LOS 方向却当成垂直形变：侧视投影未做几何分解，需联合升降轨或注明方向。
• 跨 burst 接缝处相位不连续：未做 burst 拼接，需用 topsApp 的拼接模式或分段处理。
• 处理后发现用了预报轨道：预报轨道误差达分米级，务必下载精密轨道重跑。
• 山区形变被大气淹没：未做分层延迟改正，接入 GACOS 或 ERA5 数据。
• 季节性形变被当成线性趋势：模型未加周期项，改用带周期项的形变模型。
• 多视倍数设得过大：细节被抹平且相干性虚高，按目标分辨率反推多视倍数。

小结

SAR 与 InSAR 的价值在于全天候与毫米级灵敏度，但它的精度完全取决于误差链条的管理。把配准、去平地、去地形、解缠、大气改正、参考点选择这六个环节逐一做扎实，剩下的相位才是可信的形变信号。

对工程团队而言，建议先用 HyP3 或 ARIA 的成品干涉图快速验证区域可行性，再决定是否自建 ISCE 加 MintPy 的完整流水线。流水线一定要保留中间产物并做质量评估，形变监测的代价往往不在算力，而在误判。

还有一点值得强调：SAR 形变监测的结果必须与地面实测（水准、GNSS）做交叉验证。InSAR 给出的是相对参考点的 LOS 位移，GNSS 给出的是绝对三维位移，两者在垂直分量上应能互相印证。若差异超过 1cm/yr，说明处理链条中某个环节仍有系统偏差，需要回溯排查。

下一步可结合 变化检测方法 把形变与地表覆盖变化联合分析，或用 栅格数据读写与 GDAL 优化大规模 SLC 堆栈的存储与切片策略。