Sentinel-1A SLC到沉降速率图:PS-InSAR与SBAS避坑

发布时间:2026/11/3 14:04:32
Sentinel-1A SLC到沉降速率图:PS-InSAR与SBAS避坑

简介:这份PPT以Sentinel-1A卫星影像的PS/SBAS地面沉降监测为主题,面向遥感、测绘、地质灾害监测方向的学习者与技术人员,通过天津案例串联数据预处理、PSInSAR与SBAS处理流程、参数设置、SRTM DEM辅助数据使用及沉降结果展示,帮助读者理解InSAR地面沉降监测从数据到产品的完整链路。包内为1个pptx文件,压缩包约36.93MB,内容以技术流程、参数配置和案例结果为主,适合课堂汇报、项目复盘或自学查阅。当前已有505人学习,可作为SAR/InSAR入门与案例实践的参考材料,也能为地面沉降监测方案设计、流程梳理和结果解读提供较具体的示例。

1. 从 32 景 SLC 到一张沉降速率图:PS 与 SBAS 到底在解什么反问题

Sentinel-1A 的 32 景 Level 1 SLC 影像,跨越 2015 年 5 月到 2016 年 7 月、约 26×20 平方公里的地面沉降监测区,最后被压缩成一张色标 ±100 mm/a 的速率图——整个过程中,没有任何一个地面点被真正测量过。这就是 SAR 差分干涉最反直觉的地方:它测的不是沉降本身,而是两次成像之间电磁波往返相位的差,再把相位差翻译成视线向的毫米级位移。PS 和 SBAS 是两套把这笔相位账算清楚的策略:一个紧盯少数永久散射体,时间基线拉得很长;另一个把短基线像对连成网,覆盖面更宽。搞懂它们的差别,本质上是在搞懂误差假设不同、结论就会不同。

2. 干涉相位里混着六种信号:为什么必须先扣地形再谈形变

2.1 从往返相位到缠相位:差分干涉的记账方式

卫星发射的是一段连续载波,回波被记录下来时,携带的是从天线到地面再返回的整段路程对应的相位。两次成像之间,如果地面沿视线向移动了 $d$,这段路程就变化了 $2d$,相位变化量是

$$\varphi_{defo} = -\frac{4\pi}{\lambda} d$$

Sentinel-1 工作在 C 波段,波长约 5.55 厘米,代入上式可知:一个完整的 $2\pi$ 相位周期,对应视线向位移 $\lambda/2 \approx 27.7$ 毫米。反过来说,想拿到 1 毫米的位移精度,相位测量精度得到 0.227 弧度,约 13 度。单看一景图像当然做不到,这就是为什么必须用多时相。

麻烦在于,两景影像做差得到的那张干涉图,相位里混进的东西远不止形变。一个像素的干涉相位可以拆成六项:

分量物理来源是否可扣
平地相位 $\varphi_{flat}$参考椭球面上斜距随位置变化用轨道几何精确扣除
地形相位 $\varphi_{topo}$地面高程引起的斜距差用外部 DEM 扣除,残差留在相位里
形变相位 $\varphi_{defo}$想要的信号保留
大气相位 $\varphi_{atm}$对流层折射率随高度、湿度变化部分可估,难完全扣
轨道误差精密轨道残余、基线误差用多项式拟合去除
噪声 $\varphi_{noise}$热噪声、时间去相干靠滤波和多视抑制

地形相位的量级可以用一个数感受一下:垂直基线 120 米、斜距 850 公里、入射角 38 度时,高程模糊度是

$$h_{amb} = \frac{\lambda R \sin\theta}{2 B_\perp} \approx 121\ \text{米}$$

也就是说,高程每变化 121 米,相位就绕一整圈。90 米采样间隔的 SRTM DEM 在城区高层建筑区,绝对高程误差到十几米是常事,换算成相位就是几十度——比我们想测的形变信号大一到两个数量级。所以,"先去地形相位再谈形变"不是流程上的礼貌,而是信噪比上的硬约束。

import numpy as np lam = 0.0555 # C 波段波长,m theta = np.deg2rad(38.0) # 入射角 R = 850e3 # 斜距,m B_perp = 120.0 # 垂直基线,m h_err = 15.0 # DEM 高程误差,m v = -0.035 # 视线向形变速率,m/a(负号表示远离卫星) dt = 0.4 # 时间基线,a # 地形相位:高程误差被基线放大成相位误差 phi_topo = 4 * np.pi * B_perp * h_err / (lam * R * np.sin(theta)) # 形变相位:双程,所以是 4*pi/lambda 而不是 2*pi/lambda phi_defo = -4 * np.pi * v * dt / lam phi_atm = np.deg2rad(25.0) # 一次大气延迟的量级示意 phi_int = phi_topo + phi_defo + phi_atm # 任何干涉图拿到的都只是主值,落进 [-pi, pi] phi_wrapped = np.angle(np.exp(1j * phi_int)) print(f"地形相位 {np.rad2deg(phi_topo):.1f}°, " f"形变相位 {np.rad2deg(phi_defo):.1f}°, " f"缠绕后 {np.rad2deg(phi_wrapped):.1f}°")

这段代码算的是一个像素的相位合成。三个参数值得记住:B_perp决定地形相位的敏感度,基线越长,同样的高程误差被放大得越厉害;h_err是外部 DEM 的绝对精度,它是 PS 和 SBAS 都要在反演里重新估一个"高程残差"的直接原因;phi_wrapped用np.angle取出主值,是后面所有解缠算法的输入。那一步np.exp(1j*phi)再取角,就是模拟相位在复平面上的绕圈——超过 $2\pi$ 的部分被折叠回去了,解缠要做的事就是把它猜回来。

2.2 32 景、26×20 平方公里、5 米 × 20 米:这组参数定了什么

数据处理前先要固定的表,往往比算法本身更能决定成败:

参数取值影响
传感器 / 波段Sentinel-1A / C 波段波长 5.55 cm,形变敏感度约 27.7 mm/2π
极化方式VV城区建筑、裸地回波强,利于 PS 点提取
轨道方向降轨视线向朝东偏,对东西向形变不敏感
数据格式Level 1 SLC保留相位,是干涉的唯一可用级别
分辨率5 m × 20 m单像素约 100 m²,城区点密度高
影像数量32 景决定组网规模与时间采样率
覆盖时间2015-05 ~ 2016-07约 14 个月,跨一个完整水文年
辅助 DEM90 m 采样 SRTM V4提供初值,残差进反演

32 景这个数量很关键。如果按 SBAS 全组网,理论上有 $C_{32}^{2} = 496$ 个干涉对;如果按 PS 的单主影像策略,只有 31 对。同样的数据,两种策略的干涉图数量差一个数量级,计算量、内存占用、误差传播路径全都不一样。这也是为什么处理平台里会把干涉处理系统参数单独设成SENTINEL_TOPSAR模式,并挂上 OpenCL 的 GPU 平台做加速——逐像素生成干涉图、算相干性、算振幅离差,是典型的可并行的密集浮点运算,32 景在 CPU 上跑一遍足够让人喝掉半壶茶。

还有一个容易被忽略的物理约束:Sentinel-1 的 TOPSAR 成像模式采用方位向扫描,一景由若干 burst 拼接而成。相邻 burst 之间本来就有相位跳变,必须在配准阶段用增强谱分集把方位向配准精度压到千分之一像素量级,否则拼接缝上会出现周期性条纹。这个坑在后面单独说。

3. PS-InSAR:为什么"挑永远不变的点"反而更准

3.1 振幅离差指数:把"稳定"变成一个能算的数

PS 的核心假设只有一句:存在一类地面目标,它们在所有影像里的后向散射特性几乎不变。这类目标通常是建筑物的墙角、桥梁的金属构件、裸露的岩石——尺寸远大于波长、以单个散射体为主,不产生盖尔伯特(Gel'bert)随机散射。

"不变"怎么量化?用振幅的时间统计量。对每一个像素,沿着时间轴取 $N$ 景的振幅 $A_i$,算均值 $\mu$ 和标准差 $\sigma$,两者之比就是振幅离差指数:

$$D_A = \frac{\sigma_A}{\mu_A}$$

这个量的物理含义很直接:如果像素的回波由一个稳定散射体主导,振幅只受热噪声扰动,$D_A$ 就很小;如果像素里是一堆随机分布的小散射体,回波相位与振幅都服从 Rayleigh 统计,$D_A$ 就会明显偏大。工程上常用的阈值在 0.25 到 0.40 之间,阈值越低,选出的点越可靠但越稀疏。

import numpy as np def amplitude_dispersion(amp, thresh=0.25): """ amp: 形状 [n_image, n_line, n_pixel] 的振幅栈 thresh: 常用 0.25 ~ 0.40,越低越保守 """ mu = amp.mean(axis=0) # ddof=1 用无偏估计;影像数 32 时两者差别可忽略,但写清楚更好 sig = amp.std(axis=0, ddof=1) da = np.divide(sig, mu, out=np.full_like(mu, np.inf), where=mu > 0) return da < thresh, da

where=mu > 0那一句不是多余的防御,水体、阴影区域的均值为零或极小,直接除会产生 inf 和 nan,这些值一旦漏进后续的相干系数图或者反演矩阵,整个最小二乘会直接崩掉。这类小节在文档里不会写,但排查的时候会让人找一整天。

3.2 单主影像连接图:用时间基线换相干性

PS 的连接图策略与 SBAS 完全不同:选一景主影像,其余全部与它配准成对。选主影像的准则通常是让所有像对的平均相干性最大,可以从"时间基线居中的影像"里挑,也可以根据已有的相干性估计迭代挑。

这里有个反直觉的设计取舍。单一主影像意味着部分像对的时间基线可以拉到一两年,长相干时间去相干本应非常严重,但 PS 点恰恰是那些不受时间去相干影响的目标——它们的相位在时间上是"硬"的。所以 PS 敢用长基线,因为它的选点机制已经把脆弱的像素排除掉了。同时,单主影像把垂直基线整体压小,地形相位残差被约束住。用一句话概括:PS 用空间基线换时间采样,SBAS 用时间基线换空间覆盖。

干涉工作流的步骤是固定的:SLC 配准 → 生成差分干涉图 → 去地形相位 → 生成相干系数图 → 生成振幅离差指数。这几步的产物同时交给后面的候选点选择和反演。

3.3 两次反演的分工:第一次解参数,第二次解大气

PS 的相位模型通常写成线性叠加:

$$\varphi_i = -\frac{4\pi}{\lambda} v \cdot t_i + \frac{4\pi B_{\perp,i}}{\lambda R \sin\theta}\Delta h + \varphi_{atm,i} + \varphi_{noise,i}$$

以 $v$(形变速率)、$\Delta h$(高程残差)和常数项为未知数,每个 PS 点用 31 个观测做最小二乘,这是第一次反演的核心。

import numpy as np def ps_first_inversion(phi, t, b_perp, lam, R, theta): """ phi: 形状 [n_pair] 的解缠相位 t: 时间基线,a b_perp: 垂直基线,m 返回: 速率(m/a)、高程残差(m) """ coef = 4 * np.pi * b_perp / (lam * R * np.sin(theta)) A = np.column_stack([t, coef, np.ones_like(t)]) # [时间, 高程, 常数] x, *_ = np.linalg.lstsq(A, phi, rcond=None) v = -x[0] * lam / (4 * np.pi) # 还原成速率 dh = x[1] # 高程残差 resid = phi - A @ x # 留给第二次反演 return v, dh, resid

注意A的三列结构:第一列是时间,第二列是基线系数,第三列是常数。常数项吸收的是参考相位偏差和部分轨道残差,不给它留位置,$v$ 和 $\Delta h$ 就会被系统性偏置。resid是第一次反演剩下的相位,它里面主要就是大气和噪声,这是第二次反演的原料。

第二次反演要干的事是:对残差相位在空间上低通、时间上高通滤波,估出每个时刻的大气延迟屏,扣掉,再用扣完的相位重估一次 $v$ 和 $\Delta h$。为什么要两遍?因为大气相位在时间上也可能表现为缓变趋势,跟线性形变在时域上不容易分开;但大气在空间上是分层的、随高程变化,形变在空间上往往是局部的漏斗状。利用这种"时域相似、空域不同"的差异,才能把它们剥离开。这是典型的玄学环节,参数调不好,扣完大气形变也跟着被削掉一块。

3.4 地理编码:一个 GCP 和多个 GCP 差在哪

反演出来的东西还在斜距-方位坐标系里,得转到地图坐标上。这一步要轨道、要 DEM、要地面控制点。控制点数量直接决定系统误差的形态:

一个 GCP 只能解一个平移量,等于把整幅图整体挪一挪,轨道残余、基线误差造成的斜坡完全没纠正;多个 GCP 可以做仿射甚至二次多项式拟合,把残余的趋势面吸收掉。城区做 PS 时,如果只给一个控制点,很容易看到整幅速率图带一个方向一致的小斜坡,看着像区域性的构造沉降,实际是几何校正没做干净。多个 GCP 也有讲究:要选在稳定、远离形变梯度大的地方,否则控制点本身就在动,等于用一个错误的基准去校准。

4. SBAS:短基线组网之后,误差从哪来

4.1 连接图与连通分量:SBAS 会自然产生"孤岛"

SBAS 的思路是不挑点,而是挑像对。给时间基线和空间基线各设一个阈值,只保留两个维度都短的像对,这样每张干涉图的相干性都还不错,解缠更容易成功。

import numpy as np from scipy.sparse.csgraph import connected_components n = 32 t_days = np.arange(n) * 36.0 # 约 36 天一景 b_perp = np.random.normal(0, 40, n) # 垂直基线,示意 t_max, b_max = 180.0, 150.0 # 双阈值,按实际几何调整 adj = np.zeros((n, n), dtype=int) for i in range(n): for j in range(i + 1, n): if abs(t_days[i] - t_days[j]) <= t_max and abs(b_perp[i] - b_perp[j]) <= b_max: adj[i, j] = adj[j, i] = 1 k, label = connected_components(adj, directed=False) print(f"连通分量个数 {k}") # 大于 1 就说明组网出现了孤岛

这段代码点出了 SBAS 最本质的结构问题。阈值一收紧,某些影像可能因为基线条件不满足,跟任何其他影像都连不上,整张网络被切成几个互不相连的子集——这就是"孤岛"。孤岛内部可以相对定标,孤岛之间没有观测约束,只能靠先验或者共同参考点强行接起来。connected_components返回的k大于 1,就得回到连接图去调阈值或补像对。

组网之后的方程是超定的,但秩亏——因为所有干涉相位都是相对量,整体有一个未知常数。解法通常是奇异值分解求最小范数解,把形变的时间序列拆成若干段线性叠加。这也解释了 SBAS 为什么对形变模型的假设是"在相邻短基线之间近似线性":它不需要形变长期线性,只需要局部线性。

4.2 自适应滤波与相位解缠:SBAS 最脆弱的一环

SBAS 的干涉工作流比 PS 多了两步:自适应滤波和相位解缠。滤波是为了压噪声、提相干性,代价是空间分辨率损失,滤波器窗口设大了,小尺度的形变梯度就被抹平了。

解缠是整条链上最容易失控的地方。缠相位是 $[-\pi, \pi]$ 的主值,解缠要给每个像素加上整数倍的 $2\pi$,本质上是一个全局优化问题。常用的最小费用流类算法,在相干性好的区域几乎不会出错,但在低相干区域——水体、植被、施工工地——一个错误的 $2\pi$ 跳变会像裂纹一样向周围蔓延,最后在速率图上表现为一条突兀的界面,两边差出几十毫米每年。

这也是为什么 SBAS 的第一次反演除了估速率,还要做精度估计:用残差相位的方差推断每个像素速率的标准差,把不可信的区域标出来。精度产品不是装饰,它是判断"这张图能不能用"的唯一依据。

4.3 轨道精炼与重去平:解缠完了为什么还要回头

SBAS 流程里有个容易让人困惑的步骤:第一次反演之后,做轨道精炼和重去平,然后再解缠一次。都已经解缠过了,为什么还要再来一遍?

原因是第一次解缠时的输入里带着系统性的相位趋势——残余的轨道误差、DEM 误差、还有解缠本身引入的整周偏差。这些趋势在解缠阶段是"看不见"的,因为解缠只保证相邻像素的相位连续,不保证整体没有斜坡。第一次反演把这些趋势估出来(通常是空间上的低阶多项式),从相位里扣掉,得到"重去平"的干涉图,再解缠一遍,质量会明显好于第一遍。

连接图编辑这一步也在这里。三条经验规则值得记住:三维解缠的结果里带_del_upha后缀的像对建议保留;相干系数图、滤波后干涉图、相位解缠图三者中哪一个质量好就保留哪一个,不必强求三件套齐全;含有明显大气相位的像对建议剔除。第三条尤其重要,一张被强对流层延迟污染的干涉图,进了反演就是往整个时间序列里掺沙子,宁可少一对。

4.4 第二次反演与精度产品

第二次反演做的事和 PS 类似:估大气相位、扣除、重估速率和 DEM 残差。区别在于 SBAS 的大气估计是在时间序列上做高通、空间上做低通,利用的是大气在空间上相关、时间上不相关的特性。

出来的产品有两类:平均形变速率图和累计形变图。前者是斜率,后者是积分。累计形变对解缠误差极其敏感——某一景解缠错了一个整周,从那之后所有时刻的累计值都平移 27.7 毫米,图上表现为一个台阶。所以看累计形变曲线时,如果看到某个时刻突然跳变然后保持平移,八成不是地面真的塌了一下,是解缠错了。

5. 避坑与排查:TOPSAR 处理链上最容易翻车的五处

5.1 burst 边界出现周期性横条纹

现象:差分干涉图在方位向上每隔固定距离出现一条相位不连续的横线,密集时像栅栏。

原因:TOPSAR 的方位向扫描使相邻 burst 之间存在多普勒中心差异,配准精度不够时,拼接处的相位被引入一个常量偏移。方位向配准要求达到千分之一像素量级,常规的强度相关配准根本达不到。

解决:改用增强谱分集做方位向精配准,配准后检查 burst 拼接处的相位是否连续;如果仍有残余,在拼接缝附近做局部相位解缠验证。这条几乎没有取巧空间,配准不过关,后面所有步骤都是白做。

5.2 速率图上出现与地形高度相关的条纹

现象:把速率图和 DEM 叠在一起看,速率的高低起伏跟地形起伏明显相关,山区或高楼区尤其明显。

原因:地形相位没扣干净,残差进了形变项。90 米 SRTM 在高楼密集区的高程误差可以到十几米,按前面算的高程模糊度,相位残差能到几十度,被误当成形变。

解决:在 PS 反演里显式估计高程残差并输出,检查残差量级是否合理;在 SBAS 里检查重去平后的残余相位是否还有地形相关结构。如果残差普遍偏大,考虑换更高精度的 DEM,或者提高空间基线阈值把长基线像对剔掉。

5.3 局部孤立"沉降漏斗"其实是大气

现象:某一小块区域在所有时相上速率都异常,形状跟地形、跟已知工程活动都对不上,换个时间段重跑,位置还会漂移。

原因:湍流型大气延迟在空间上呈斑块状,尺度几百米到几公里,量级能到几厘米,在单张干涉图上完全可能伪装成一个沉降中心。

解决:看时间序列而不是只看平均速率——真实沉降是单调的,大气是随机的;用第二次反演的空时滤波压一遍;跨季节分组处理,如果"漏斗"只在夏季出现,基本可以判定是大气。

5.4 解缠错误在低相干区蔓延成大面积极值

现象:速率图上出现大片 ±100 mm/a 的极端值,色标被撑爆,正常区域全糊成一片灰。

原因:低相干区解缠失败,整周误差沿着路径传播,反演时这些错误被当成真实形变。相干性阈值设得太低、滤波窗口太小,都会放大这个问题。

解决:提高相干性阈值,宁可让低相干区留空;用三维解缠把时间维的连续性也用上;对解缠结果做人工检查,把明显错误的像对在连接图里剔掉。图上大片极值几乎都不是地质现象,别急着解释。

5.5 全图整体偏移,像区域沉降其实是基准问题

现象:整幅图的速率都朝一个方向偏,看着像全区均匀沉降,但幅度不大,几毫米每年。

原因:参考点或控制点选在了形变区,或者只给了一个 GCP 导致轨道残余没被吸收。基准错了,所有相对量都被整体平移。

解决:控制点选在稳定基岩或已知不动点上,尽量给多个 GCP 做多项式校正;处理后检查整幅图的速率均值,如果显著不为零,先怀疑基准而不是地质。

6. 取舍:PS 与 SBAS 是在选择不同的误差假设

把两条流程并排放,会发现它们并不是精度高低的关系,而是对"什么样的像素可信"给出了不同答案。

维度PS-InSARSBAS
选点方式振幅离差指数筛永久散射体全像素,靠相干性加权
连接图单主影像,31 对多主影像短基线组网,可达数百对
时间基线可以很长受阈值约束,普遍较短
空间覆盖城区密集,郊野稀疏覆盖均匀,低相干区留空
主要误差选点漏检、大气与形变混叠解缠错误传播
适合场景城市建筑、基础设施、单体形变区域面状沉降、大范围普查

真正的取舍在误差结构上。PS 把风险前置到选点:只要点选对了,后面每一点的时序都是干净的,代价是空间覆盖残缺,郊野和农田几乎拿不到点。SBAS 把风险后置到解缠:覆盖面广,但整周误差一旦进去就会沿着网络扩散,累计形变图上留下永久的台阶。所以在实际工程里,我一般会两条链并行跑同一批数据,用 PS 的高置信点去校验 SBAS 的面状结果,两套结果在重叠区对不上的地方,往往就是问题所在——要么 PS 点被建筑自身的非线性形变污染,要么 SBAS 在那一带解缠翻了车。

还有一层常被忽略的取舍是视线向的几何约束。降轨、VV 极化、约 38 度入射角,这个组合对垂直向形变敏感,对东西向几乎不敏感,南北向完全测不到。速率图上那个负号代表远离卫星,至于地面到底是下沉还是水平移动,单靠一个轨道的几何是分不开的,得上升轨降轨联合解算。这一点在写报告时最容易被含糊过去,也是我看别人结果时第一个会问的问题。

血泪经验是:这套流程里最贵的不是算力,是判断。参数能调,算法能换,但"这个漏斗是真的还是大气""这条条纹是配准还是解缠",只能靠一遍遍对比、看时间序列、看空间形态慢慢积累。从那以后我每次出图之前,都会强制走一遍三件事:检查连通分量个数、检查残差相位有没有地形相关、检查累计形变曲线有没有台阶。这三步花不了半小时,但能省掉后面返工好几天。希望帮到你。

本文还有配套的精品资源,点击获取