Sentinel-1A地面沉降监测：PS-InSAR与SBAS流程选型与避坑

简介：这份PPT案例剖析面向遥感、测绘与地质灾害监测方向的学习者，围绕Sentinel-1A降轨VV极化影像在天津地区的地面沉降监测展开，重点解决PSInSAR与SBAS两条技术路线的流程理解与参数配置问题。内容以32景Level 1 SLC影像（2015年5月至2016年7月）和90米SRTM DEM为数据基础，串联数据导入、参考DEM下载、连接图生成、干涉工作流、第一次与第二次反演、大气相位去除、轨道精炼、相位解缠及地理编码等环节，并给出TOPSAR系统参数与GPU加速设置思路。压缩包仅含1个pptx文件，大小约36.93MB，便于集中查阅和课堂演示。已有505人学习，适合希望把SAR沉降监测理论落到软件操作与案例复盘的读者，可从中获取完整技术流程、关键参数说明、结果展示框架及常见处理要点。

1. 从一份 32 景 Sentinel-1A 的案例 PPT 说起：PS 和 SBAS 该怎么选

有人拿一份沉降监测的案例 PPT 来问：同样 32 景 Sentinel-1A 降轨 SLC，为什么前后做了 PS 和 SBAS 两套流程，结果还都对得上？这问题问到了点子上。Sentinel-1A 的 C 波段 VV 数据，26×20 平方公里的建成区，2015 年 5 月到 2016 年 7 月，一年多的跨度，做地面沉降监测时在 PS-InSAR 和 SBAS 两条路之间选一条，其实取决于你要回答什么问题：是要一批稳定点上的年速率，还是要一片面上的累计形变场。这份案例把两套流程的参数、节点、取舍都摆了出来，适合刚上手 SAR 干涉、手里已经有几景数据但卡在连接图和反演那一步的人。下面按数据落地、流程拆解、翻车点、成果验证的顺序,把这份案例里能直接抄的部分拆开。

2. 数据落盘与预处理：32 景降轨 SLC 进工作流前的四道关

2.1 先定 TOPSAR 系统参数和 OpenCL 平台，别急着导数据

Sentinel-1A 用的是 TOPS 成像模式，跟传统条带模式最大的差别在于：一幅 SLC 里包含三个子条带，每个子条带又由若干 burst 拼成，burst 之间在方位向有重叠。这意味着配准不能按整幅影像做一次多项式拟合了事，必须按 burst 做方位向配准，精度要求压到千分之一像元量级。这个参数如果一开始没设对，后面干涉图上会出现周期性的相位跳变，看起来像条纹，实际是配准残差。

我一般会在动手之前先把两件事配好。第一件是干涉处理系统参数选SENTINEL_TOPSAR，这是告诉软件后面所有配准、干涉、解缠都按 TOPS 的几何来处理，别用默认的条带模式参数。第二件是 OpenCL 平台，路径在/SARscape/Preferences下的系统参数里。32 景 SLC 做全组合干涉，计算量不小，单靠 CPU 跑连接图那一步能耗掉大半天；把 GPU 挂上去之后，同样的活能压到可接受的范围。这里有个前提：显卡驱动和 OpenCL 运行时版本要匹配，选错平台在连接图阶段不会立刻报错，而是在干涉工作流跑到一半时抛出一个内存或设备异常，排查起来很费劲。

提示：OpenCL 平台列表里如果出现多个设备，优先选独显那块，别选到集显或者 CPU 模拟设备上，否则加速效果接近于零，还多占内存。

2.2 导入、轨道、DEM：一条不能省的准备链

数据格式是 Level 1 SLC，从官方开放数据门户按成像日期检索，注册后免费下载，一个 zip 里带影像和标注文件。导入走/SARscape/Import Data/SAR Spaceborne/SENTINEL 1，这一步会把 zip 解成软件内部的 SLC 格式，同时把轨道信息、成像几何、burst 边界一起读进去。

紧接着是精密轨道。Sentinel-1A 的预报轨道精度在米级，做干涉够用但不够好；精密轨道文件一般在成像后二十天左右发布，精度到厘米级。这个案例的 32 景跨越一年多，如果精密轨道不全，会出现一个问题：部分像对的基线估计偏差，导致去平地相位不干净，残余的平地条纹叠加在形变相位上。常见做法是把所有精密轨道文件一次性放进轨道目录，导入时软件按成像时间自动匹配；匹配不上的景单独列出来，缺的用预报轨道顶上，并在后面的连接图里降低这几景的权重或者干脆剔除。

DEM 用的是 90 米采样间隔的 SRTM DEM V4。90 米对沉降监测来说勉强够用：它决定的是地形相位的去除精度，在平原城市里地形起伏小，残余地形相位对速率的影响通常在毫米每年量级；如果研究区里有高层建筑密集区或者明显的坡地，90 米 DEM 的高程误差会被放大到几十米的水平分辨率上，形成局部假形变。这种地方我一般会换用 30 米 DEM 复核一遍，看差异点是否集中在同一片区。

2.3 用脚本核对文件清单与轨道配对

导完数据别急着往下走，先花两分钟把清单核一遍。血泪经验是：32 景里混进一景成像时间差了几天的、或者某景的轨道文件没对上，前面所有干涉图都会带着这个问题跑下去。

import os, glob, re from datetime import datetime slc_dir = "/data/s1a/slc" # SLC 解压后的目录 orb_dir = "/data/s1a/orbit" # 精密轨道 .EOF 目录 # Sentinel-1 文件名里带成像起始时间，形如 S1A_IW_SLC__1SDV_20150512T... pat = re.compile(r"(\d{8}T\d{6})") slc = {} for f in glob.glob(os.path.join(slc_dir, "*.zip")): m = pat.search(os.path.basename(f)) if m: slc[datetime.strptime(m.group(1), "%Y%m%dT%H%M%S")] = f orb = {} for f in glob.glob(os.path.join(orb_dir, "*.EOF")): # 轨道文件头里有 VALIDITY_START，这里简化成取文件名中的日期 m = re.search(r"(\d{8}T\d{6})", os.path.basename(f)) if m: orb[datetime.strptime(m.group(1), "%Y%m%dT%H%M%S")] = f print("SLC 景数:", len(slc), " 轨道文件数:", len(orb)) # 每景找时间差最小的轨道，超过 24 小时视为缺失 for t, f in sorted(slc.items()): if not orb: print(t, "无轨道") continue best = min(orb, key=lambda o: abs((o - t).total_seconds())) gap_h = abs((best - t).total_seconds()) / 3600 if gap_h > 24: print(t, "轨道缺口 %.1f 小时" % gap_h)

这段脚本干的事很朴素：从文件名里抠出成像时间，跟轨道目录里的时间做最近邻匹配，超过 24 小时就报缺口。参数上要改的地方只有两个目录路径和那个 24 小时阈值；如果你的数据跨了不同年份，建议把阈值收到 12 小时，Sentinel-1 的精密轨道是按天发布的，12 小时足以覆盖单景的整个采集弧段。脚本跑完输出干净了，再进下一步，比后面在干涉图上一张张找异常省事得多。

2.4 裁剪与坐标系：26×20 平方公里的边界怎么定

整景 IW SLC 覆盖范围远大于研究区，直接全幅处理会让连接图的组合数翻几倍，内存和时间都吃不住。这个案例的研究区是 26×20 平方公里，做法是先按边界裁一个子区出来，所有后续处理都在这块上跑。

裁剪要注意两点。一是裁剪边界不要切在 burst 重叠区上，否则子条带拼接会在裁边位置产生相位不连续；稳妥的做法是留出几百米余量，处理完再按研究区边界裁成果。二是坐标系，SLC 本身是斜距-方位向的雷达坐标，所有地理相关的操作（DEM 引入、地理编码、成果出图）都涉及投影转换。我一般统一用 UTM，研究区跨带的话按中心经度选带号，别中途混用地理坐标和投影坐标，混用之后地理编码出来的速率图会有整体偏移，这种偏移在图上看着像沉降，其实是坐标错了。

3. PS-InSAR 链路拆解：从振幅离差指数到第二次反演的六步

3.1 主影像与连接图：PS 的骨架

PS 流程的第一步是选主影像。逻辑是所有从影像都跟主影像做干涉，主影像选得偏，部分像对的空间基线或时间基线就会拉长，相干性掉下来，后面 PS 点选得稀。判据一般是看基线组合：空间基线在临界基线的百分之几以内、时间基线尽量居中。32 景的案例里，覆盖时间从 2015 年 5 月到 2016 年 7 月，主影像落在中间位置，两头的时间基线都在半年出头，这是个均衡的选择。

连接图生成这一步要填的参数包括最大时间基线和最大空间基线。填得宽松，像对多、冗余大、计算重；填得紧，像对少、某些时段的数据被孤立。经验值是按时间跨度分档：一年半的跨度，最大时间基线放到 300 天左右比较合适，能保证每景至少有三到四个像对；空间基线按 Sentinel-1 的轨道管特性，放到 150 米以内通常不会损失相干性。

注意：连接图生成后一定要看一眼像对网络图，有没有哪一景只连了一两个像对。这种孤立景在第一次反演里会被当作低质量数据剔除，白导进来了。

3.2 候选点选择：振幅离差指数 0.25 不是拍脑袋

PS 候选点的筛选靠两个指标：振幅离差指数和相干系数。振幅离差指数衡量的是像素在时间序列上的幅度稳定性，值越小越稳；相干系数衡量的是干涉相位的质量。

阈值怎么定？建成区里人工建筑、路面、裸露岩石的振幅离差指数普遍在 0.3 以下，植被、水面、裸土则在 0.4 以上。案例里用的 0.25 是个偏严的取值，好处是选出来的点可信度高，坏处是在郊区和新开发片区会丢点，密度可能降到每平方公里几十个。如果研究区一半是建成区一半是农田，我一般会把阈值放宽到 0.35，靠后面的相干系数和解缠相干性再筛一道，两级筛选比一级卡死更实用。

相干系数阈值同理，常用 0.7 起步。这两个阈值是联合作用的：振幅离差指数筛的是"点本身稳不稳"，相干系数筛的是"这个点在干涉图里相位干不干净"，两者都过才进候选集。

3.3 两次反演的分工与参数表

PS 的第一次反演解的是位移速率和高程残差，用的是线性模型加相位解缠；第二次反演解的是大气相位，在去除大气之后再估一次速率和高程残差。这两步的分工很多人搞混，以为第二次只是"精化一下"，实际上第二次才是拿到可用速率的关键——大气相位在单幅干涉图上的量级能到几厘米，直接换算成年速率就是几十毫米每年的假信号。

第二次反演的大气滤波是时空联合的：时间维做高通，把缓慢的形变留住、把快变的大气滤掉；空间维做低通，因为大气在空间上是平滑的。这两个窗口的取值直接决定成果的"软硬"：空间窗口给小了，局部的真实沉降漏斗会被当成大气抹平；给大了，大尺度的大气残差会漏进速率图。案例的研究区是 26×20 平方公里，空间窗口放在一到两公里之间比较合适，正好比单个建筑群大、比整个研究区小。

3.4 地理编码与成果导出

地理编码是把雷达坐标下的速率转到地理坐标，需要DEM、需要控制点。控制点（GCP）的作用是校正轨道残余误差，一个 GCP 只能做平移校正，多个 GCP 才能拟合出旋转和尺度项。案例里明确列了"一个 GCP"和"多个 GCP"两种地理编码，实际做事时只要研究区跨度超过十公里，就用多个 GCP，均匀分布在四角和中心，选在已知稳定的区域。

导出的成果一般包括平均位移速率、累计形变、相干系数、高程残差几张栅格，外加一份点文本。拿到点文本后我习惯先做一遍统计，看看速率分布有没有异常长尾：

import numpy as np # SARscape 导出的点文本，列序随版本略有差异，先看头两行确认 raw = np.loadtxt("ps_velocity.txt", skiprows=1, usecols=(0, 1, 2, 3, 4)) lon, lat, vel, hres, coh = raw.T print("点数:", len(vel)) print("速率均值 %.2f mm/a, 标准差 %.2f" % (vel.mean(), vel.std())) print("P1/P99: %.2f / %.2f" % (np.percentile(vel, 1), np.percentile(vel, 99))) # 高程残差超过 30 米的点大概率是解缠错误或建筑相位中心偏移 bad = np.abs(hres) > 30 print("高程残差异常点数:", bad.sum())

这里usecols的列序一定要按实际文件改，不同导出设置下列的顺序会变，照抄容易读串列，把高程残差当成了速率。统计出来的标准差如果在 5 毫米每年以内，说明整体解算稳定；P1 和 P99 的差值如果超过 200 毫米每年，八成是有局部解缠错误，得回到第二次反演之前去查那几景的干涉图。

4. SBAS 链路拆解：小基线连接图、3D 解缠与连接图编辑的取舍

4.1 连接图的三个阈值怎么定

SBAS 跟 PS 最大的区别在于：PS 是单主影像、全组合；SBAS 是多主影像、只保留小基线组合。这个差别决定了 SBAS 在植被覆盖区、非城市区的表现更好，因为小基线像对的去相干程度低。

连接图生成的三个阈值是空间基线、时间基线、多普勒中心频率差。空间基线的常用判据是临界基线的百分之二左右，Sentinel-1 的临界基线在几公里的量级，取百分之二大概是一百多米；时间基线按地表类型定，建成区可以放到 200 天以上，农田和植被区最好收到 60 到 90 天。多普勒差这个参数容易被忽略，它控制的是方位向频谱重叠度，差得太多干涉相位会被频谱失配污染，一般控制在几十赫兹以内。

三个阈值联合起来决定的是连接图的连通性。SBAS 最怕的不是像对少，而是连接图断成几个不连通的子集——一旦断开，解算时只能分块处理，块与块之间没有公共参考，拼起来会有整体偏移。生成完连接图，先看连通分量，理想情况是所有景归到一个连通集里。

4.2 干涉工作流里的滤波与解缠

干涉工作流这一步要跑：配准、干涉图生成、去平地、去地形、自适应滤波、相干估计、相位解缠。多视参数是要想清楚的，案例的 SLC 分辨率是 5 米乘 20 米，如果按距离向 4、方位向 1 做多视，输出格网大致是 20 米乘 20 米，正好接近正方形的像素，便于后续出图和统计。

自适应滤波的目的是压噪，但滤波强度过头会把小尺度的形变信号一起抹掉。判断滤波是否过头有个土办法：看滤波前后的干涉图，如果条纹的连续性明显改善但边缘变得模糊，说明滤波力度刚好；如果条纹被抹成一片均匀的灰，那就是过头了。

相位解缠这一步，SBAS 允许用三维解缠，也就是在时间维和空间维同时约束。三维解缠的结果里会带_del_upha后缀的像对，这些像对的解缠结果是被当作参考基准使用的。案例里明确写了"建议保留"，原因是删掉它们之后，后续的反演会失去时间维的约束，残余相位会被错误地分配到位移速率上。这个坑我见过不止一次：有人嫌中间产物占地方，批量清理临时文件时把_del_upha一起删了，结果第二次反演出来的速率整体偏移几毫米每年，怎么也找不到原因。

4.3 连接图编辑：哪些像对该保留、哪些必须剔除

连接图编辑是 SBAS 特有的、也是最考验经验的一步。它不是可选项，案例里标了"可选"但实际操作中几乎必做。判据是三类结果的质量：相干系数图、滤波后的干涉图、相位解缠图。

这里有个取舍原则值得说清楚：如果三张图里只有一张质量好，保留；三张都不行，剔除。很多人看到干涉图条纹乱就一刀切全删，结果某些时段的数据被清空，时间序列出现断档，速率拟合的斜率直接失真。反过来，什么都不删、全丢进去反演，大气和失相干的相位会被硬塞进模型，成果图上出现成片的假沉降区，看着像真的一样。

4.4 轨道精炼与重去平：GCP 与一个 GCP 的差别

轨道精炼和重去平这一步解决的是残余的轨道误差和相位常数。做法是选地面控制点，用多项式拟合残余相位并扣除。

一个 GCP 和多个 GCP 的差别在这里体现得最直接：一个 GCP 只能定一个常数偏移，等于把所有像对的相位整体平移一下；多 GCP 才能拟合出随空间变化的趋势项。案例里研究区是 26×20 平方公里，这个尺度上轨道残余主要表现为沿方位向的线性或二次趋势，至少需要五到六个 GCP 才能约束住。

GCP 选在哪里？必须是已知稳定的区域，而且要在整幅图上分布开。选在沉降漏斗中心的"稳定点"，是整个流程里最隐蔽的错误——它会把整片区域的速率基准抬起来，图上所有点的沉降量都被抵消掉一部分。判断 GCP 是否可靠有个交叉验证法：先用一批候选点做精炼，看残差分布，把残差最大的几个点剔掉再跑一遍，两遍结果一致才用。

精炼完之后是第一次反演和第二次反演，逻辑跟 PS 类似，但 SBAS 的第一次反演还要做精度估计，输出速率的标准差。第二次反演去大气、重估速率和高程残差，然后地理编码。整条链路走下来，SBAS 的成果是一片连续的面，PS 的成果是一组离散点，这也是开头那个问题的答案：要面就用 SBAS，要点就用 PS，两者不冲突，可以对同一批数据各跑一遍互相校验。

5. 避坑清单：五个把 mm/a 结论带偏的翻车现场

5.1 现象：干涉图上出现周期性条纹，像连在一起的短横线

原因基本是 TOPS 配准没做方位向精配，或者配准精度不够。Sentinel-1 的 burst 之间方位向有重叠，如果只用整幅影像做一次多项式配准，重叠区的相位会错位，形成规则的条纹。解决方式是在干涉工作流里启用增强谱分集配准，把方位向配准精度压到千分之一像元以内。判断有没有做对，看配准报告里的方位向残差，超过 0.001 像元就得重跑。

5.2 现象：速率图上出现一条笔直的边界，两边速率差几毫米

这是裁剪边界切在了 burst 重叠区上，子条带拼接的相位不连续被当成了形变。解决方式是把裁剪范围外扩几百米，处理完再按研究区边界裁成果；已经切坏了的话，回到裁剪那一步重来，别在后面用插值修补，修补出来的边界在时间序列上会露出破绽。

5.3 现象：整片区域速率统一偏了 5 到 10 毫米每年

九成是参考点选在了沉降区。PS 的参考点是速率零点，选错等于给整个成果加了个常数偏移。SBAS 里对应的是 GCP 选在漏斗中心。解决方式是换参考点：找一块长期稳定的基岩区、老城区核心，或者干脆用外部水准数据标定一个点，把整幅速率平移到那个基准上。这种偏移很危险，因为图上形变纹理是对的，只是数值全错，不跟外部数据比根本看不出来。

5.4 现象：时间序列上出现规律的正弦波动，周期约一年

典型的大气相位残留，尤其是水汽的季节性变化。第一次反演里它会被当成形变，拟合出虚假的季节性沉降。解决方式是在第二次反演里加强大气滤波，时间窗口按季节周期设置；如果残留还在，考虑在时间序列分析阶段加一个季节项一起解算，把大气和真实形变分开。判断是真季节形变还是大气，看空间分布：大气是平滑的大片，真形变通常跟着建筑群或地下水漏斗走。

5.5 现象：_del_upha像对找不到，第二次反演报错或成果异常

这是清理中间文件时误删了三维解缠的参考像对。SBAS 的三维解缠把这些像对作为时间维的锚点，删了之后解算缺少约束，残余相位被错误分配到速率上。解决方式只能重跑解缠；预防的办法是在清理脚本里把带_del_upha的文件名加进白名单，别用通配符一把梭。

6. 精度验证与成果交付：形变速率图怎么变成能签字的结论

6.1 内部精度：先看速率标准差和时间序列残差

反演结果自带一组精度指标，别跳过。PS 的输出里有每个点的速率标准差，SBAS 的第一次反演里有精度估计。这两个量的物理含义是拟合残差，反映的是模型解释观测的程度。建成区的速率标准差一般在 1 到 3 毫米每年，如果某片区域普遍高到 8 毫米以上，说明那块地方的相位质量有问题，要么失相干，要么解缠错了。

时间序列残差也值得看。把每个点的观测相位减去拟合的线性模型，剩下的残差如果呈现系统性趋势，说明形变不是线性的，用线性速率描述会失真；如果残差是随机噪声，那线性模型成立。案例的覆盖时间只有一年多，线性假设基本可用，但如果换成三五年跨度的数据，非线性形变（比如地下水回灌后的反弹）就必须用分段模型或者时间序列分解来解。

6.2 外部验证：水准点、角反射器、跨技术比对

内部精度只能说明解算自洽，不能说明结果正确。外部验证的路子有三条。

第一条是水准点比对。研究区里如果有历史水准测量点，把雷达成果在点位附近做邻域平均，跟水准年速率比。注意坐标系和时间基准要对齐，水准的年速率是按什么时段算的，雷达成果就得用同一时段的数据重新拟合，直接拿全时段的平均速率去比，误差能到几毫米。

第二条是角反射器。人工布设的角反射器在雷达图上信噪比极高，相位中心稳定，可以作为绝对基准点，也可以用来验证几何定位精度。案例里没有提角反射器，但如果做长期监测，布几个是值得的投入。

第三条是跨技术比对。PS 和 SBAS 各跑一遍，两条链路的结果如果在一个毫米每年以内吻合，可信度就上来了。这也是这份案例同时给两套流程的价值所在：不是为了炫技，是为了互检。

6.3 出图与交付的几处硬要求

成果图有三条底线。一是色标范围要固定，别让每张图自动缩放，自动缩放会让读者误判不同区域的沉降量级；同一批成果统一用一个色标，比如正负一百毫米每年。二是要标参考点位置，图上不标参考点，别人没法复现你的基准。三是要带分辨率说明，Sentinel-1 的 SLC 是 5 米乘 20 米，多视之后格网变了，成果图上的空间分辨率必须写清楚，不然拿去跟 30 米格网的成果叠加会出偏差。

交付的数据包里我一般会放这几样：速率栅格（带投影信息）、累计形变的时间序列、相干系数图、精度指标图、点位文本、处理参数记录。参数记录这一项最容易被省略，但它是唯一能让别人复现你结果的东西——主影像日期、连接图阈值、滤波窗口、GCP 坐标，一条条写清楚。

做这行久了有个习惯：每次拿到一份新的沉降成果，不管图多漂亮，先问参考点在哪、GCP 怎么选的、_del_upha还在不在。这三件事对上，数值才敢往外报；对不上，图再好看也只是张图。从那以后我每次跑完反演都强制走一遍参考点和精度指标的复核，宁可多花半天，也不想到交付那天再返工。希望帮到你。

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