简介这份PPT资料面向地质灾害监测、InSAR遥感处理方向的学习者与工程技术人员围绕Sentinel-1A卫星影像系统讲解PS与SBAS两类时序InSAR地面沉降监测方法。内容以天津沉降区为案例覆盖数据预处理、PSInSAR与SBAS技术流程、参数设置、辅助DEM数据使用、地理编码及结果展示等环节并给出传感器、波段、极化方式、轨道方向、影像数量与覆盖时间等关键参数便于对照SARscape操作流程理解完整技术链路。资源包共1个pptx文件约36.93MB以图文幻灯片形式呈现适合作为课堂汇报、项目复盘或自学参考。目前已有500人学习可帮助读者快速建立从SLC数据导入到形变速率产品输出的整体认识掌握PS与SBAS处理中的参数配置与流程差异为实际沉降监测项目提供可借鉴的案例框架与排错思路。1. 从一份 PPT 标题说起Sentinel-1A 影像做 PSSBAS 地面沉降监测到底在干什么一份名为「案例剖析基于Sentinel-1A卫星影像PSSBAS地面沉降监测」的 PPT 标题信息量其实很密。Sentinel-1A 是欧空局哥白尼计划里的 C 波段合成孔径雷达卫星重访周期 12 天免费开放是当前做 InSAR 形变监测最常用的数据源之一。PSSBAS 是 Persistent Scatterer Small Baseline Subset 的缩写中文一般叫「永久散射体小基线集」本质是把 PS-InSAR 和 SBAS 两条技术路线做融合既保留 PS 方法对点状稳定目标的敏感性又用短基线干涉对组合提升分布式散射体区域的点密度。地面沉降监测是它的典型落地场景——城市地下水超采、矿区采空、软土固结、地铁施工影响这些缓慢累积的毫米级形变用传统水准测量覆盖成本极高而 InSAR 一次能覆盖上千平方公里。这份 PPT 面向的是想用免费 SAR 数据做区域形变普查的从业者测绘、地质、城市地下空间管理、矿山安全监测方向的技术人员。它要回答的核心问题是——拿到 Sentinel-1A 的 SLC 数据之后怎么一步步走到可用的沉降速率图中间哪些参数决定成败。下面我按自己实际跑通这套流程的顺序把选型理由、命令、参数和踩过的坑讲清楚。2. 数据准备与预处理从 Sentinel-1A SLC 到可用的干涉对2.1 为什么选 Sentinel-1A 而不是其他 SAR 数据做地面沉降数据源选择直接决定后面所有环节的天花板。Sentinel-1A 的 C 波段波长约 5.6 cm在形变监测里是个折中点波长比 X 波段长抗失相干能力强适合植被稀疏的城市和农田比 L 波段短对毫米级形变更敏感但穿透植被能力弱。免费、全球覆盖、12 天重访、IW 模式幅宽 250 km这几个条件叠加起来让它成为区域沉降普查的首选。我一般会优先选 IW 模式的 SLC 产品升轨或降轨取决于研究区的地形和形变方向。升轨对东西向形变敏感降轨对垂直向更敏感做沉降监测通常优先降轨。如果研究区有大量建筑物PS 点密度会很高如果是农田或裸地就得靠 SBAS 的分布式散射体来补点。PSSBAS 的价值就在这里——不用在 PS 和 SBAS 之间二选一。2.2 用 SNAP 做轨道校正和配准的最小命令预处理第一步是轨道校正和配准。Sentinel-1A 的精密轨道文件POD有 20 天左右的延迟做历史数据时直接用精密轨道做近实时监测时先用快速轨道后期再替换。下面是用 SNAP 命令行工具做轨道校正和配准的典型流程# 1. 应用精密轨道文件EOF 格式 gpt Apply-Orbit-File -Ssource./S1A_IW_SLC__1SDV_20230101.SAFE \ -PorbitTypeSentinel Precise \(Auto Download\) \ -PpolyDegree3 \ -t ./S1A_20230101_orb.dim # 2. 拆分子条带IW 模式有 3 个子条带 gpt Split -Ssource./S1A_20230101_orb.dim \ -PsubswathIW2 \ -PselectedPolarisationsVV \ -t ./S1A_20230101_IW2.dim # 3. 配准到主影像以 20230101 为主 gpt Back-Geocoding -Ssource./S1A_20230101_IW2.dim \ -Ssource./S1A_20230113_IW2.dim \ -PdemNameSRTM 1Sec HGT \ -PdemResamplingMethodBILINEAR_INTERPOLATION \ -PresamplingTypeBILINEAR_INTERPOLATION \ -t ./S1A_20230101_20230113_IW2.dimApply-Orbit-File里的polyDegree3是轨道多项式拟合阶数精密轨道用 3 阶足够快速轨道可以降到 2 阶。Split只保留 IW2 是因为沉降监测通常只用一个子条带就够三个子条带全处理会让数据量翻三倍。Back-Geocoding的demName用 SRTM 1Sec如果研究区在高纬度或 SRTM 有空洞换成 ASTER GDEM。配准精度直接决定干涉相位质量配准误差超过 0.1 个像素干涉图就会出现明显的条纹噪声。2.3 干涉图生成与相干性估计的参数设置配准之后生成干涉图这一步的核心参数是多视数和滤波。Sentinel-1A IW 模式的距离向和方位向分辨率不同通常做 4:1 多视距离向 4方位向 1让像元接近正方形。滤波用 Goldstein 或 BoxcarGoldstein 对低相干区域更友好但会引入相位偏差做高精度形变时我倾向用 Boxcar 加自适应滤波。# 生成干涉图多视比 4:1 gpt Interferogram -Ssource./S1A_20230101_20230113_IW2.dim \ -PsubswathIW2 \ -PselectedPolarisationsVV \ -PcohWinAz3 -PcohWinRg10 \ -t ./ifg_20230101_20230113.dim # Goldstein 相位滤波 gpt GoldsteinPhaseFiltering -Ssource./ifg_20230101_20230113.dim \ -Palpha1.0 \ -t ./ifg_filtered.dimcohWinAz3和cohWinRg10是相干性估计窗口方位向 3、距离向 10 是经验值窗口太大会平滑掉细节太小相干性估计噪声大。alpha1.0是 Goldstein 滤波的过平滑因子取值 0 到 1越大滤波越强但相位细节损失越多。做沉降监测时如果研究区主要是城市alpha 可以设 0.6 左右如果是农田设 1.0 保证相干性。2.4 去平地效应和相位解缠的边界条件去平地效应用外部 DEMSRTM 1Sec 在大部分地区够用但在地形起伏大的区域DEM 误差会残留成相位斜坡。相位解缠是 InSAR 里最玄学的一步PSSBAS 用的是三维解缠时间空间比二维解缠稳定但对相干性要求更高。# 去平地效应和地形相位 gpt TopographicPhaseRemoval -Ssource./ifg_filtered.dim \ -PdemNameSRTM 1Sec HGT \ -PdemResamplingMethodBILINEAR_INTERPOLATION \ -t ./ifg_topo_removed.dim # 导出为 Gamma 或 SNAPHU 格式做解缠 gpt SnaphuExport -Ssource./ifg_topo_removed.dim \ -PstatCostModeDEFO \ -PinitMethodMCF \ -t ./snaphu_inputstatCostModeDEFO是形变模式的统计代价适合地面沉降这种连续形变场。initMethodMCF是最小费用流初始化比区域增长法稳定。解缠后的相位如果出现大面积跳变先回头查相干性图相干性低于 0.3 的区域解缠结果基本不可信后面做时序分析时要掩掉。3. PSSBAS 时序处理PS 点选取与 SBAS 组合策略3.1 PS 点候选提取的振幅离散指数阈值PSSBAS 的第一步是选 PS 候选点。经典 PS-InSAR 用振幅离散指数ADI来筛选ADI 小于 0.25 的像元才进入候选。ADI 的计算方式是振幅标准差除以振幅均值反映像元在时间序列上的稳定性。建筑物、裸露岩石、桥梁这些硬目标的 ADI 通常很低植被和农田的 ADI 很高。import numpy as np from osgeo import gdal # 读取振幅时间序列堆栈 amp_stack [] for i in range(1, 25): ds gdal.Open(f./amp_{i:02d}.tif) amp_stack.append(ds.ReadAsArray()) amp_stack np.array(amp_stack) # 计算振幅离散指数 amp_mean np.mean(amp_stack, axis0) amp_std np.std(amp_stack, axis0) adi np.where(amp_mean 0, amp_std / amp_mean, np.inf) # 阈值筛选ADI 0.25 为 PS 候选 ps_candidates adi 0.25 print(fPS 候选点数: {np.sum(ps_candidates)})amp_mean 0的判断是为了避免除零水体像元的振幅均值接近零ADI 会变成无穷大。阈值 0.25 是经验值城市区域可以放宽到 0.3 增加点密度但会引入更多噪声点。如果研究区有大量低矮建筑ADI 分布会偏大这时候可以结合相干性做二次筛选相干性均值大于 0.7 的像元优先保留。3.2 小基线集组合的时空基线阈值怎么定SBAS 的核心是干涉对组合不是所有影像两两组合而是按时间和空间基线阈值筛选。时间基线一般设 60 到 120 天空间基线设 150 到 300 米。时间基线太短干涉对数量不够解算方程秩亏太长失相干严重。空间基线太长地形相位残差大去平地效应不干净。from itertools import combinations import pandas as pd # 读取影像元数据 meta pd.read_csv(./slc_metadata.csv) # 列date, bperp垂直基线 max_temporal 90 # 天 max_perp 200 # 米 pairs [] for (i, row_i), (j, row_j) in combinations(meta.iterrows(), 2): dt abs((row_j[date] - row_i[date]).days) db abs(row_j[bperp] - row_i[bperp]) if dt max_temporal and db max_perp: pairs.append((row_i[date], row_j[date], dt, db)) print(f满足条件的干涉对: {len(pairs)})max_temporal90和max_perp200是我在多个项目里试出来的折中值。如果研究区植被覆盖高时间基线要压到 60 天以内如果地形平坦、城市密集可以放宽到 120 天。空间基线 200 米对应 Sentinel-1A 的临界基线约 5000 米200 米以内地形相位残差可控。干涉对数量至少要 30 对以上否则时序解算不稳定。3.3 时序形变解算从相位到毫米级速率时序解算分两步先估计大气相位和轨道误差再解算形变速率。PSSBAS 用迭代最小二乘先假设形变是线性解出速率后再看残差如果残差有季节性周期加入周期项重新解算。import numpy as np from scipy.linalg import lstsq # 相位矩阵行是干涉对列是 PS 点 # 设计矩阵每行对应一个干涉对主影像为 -1从影像为 1 # 形变模型v * dt c线性速率 常数项 n_ifg len(pairs) n_ps np.sum(ps_candidates) A np.zeros((n_ifg, 2)) b np.zeros(n_ifg) for k, (d1, d2, dt, db) in enumerate(pairs): A[k, 0] dt A[k, 1] 1.0 # 逐点解算 velocity np.zeros(n_ps) for p in range(n_ps): phase_vec phase_stack[:, p] # 该点的干涉相位序列 x, res, rank, sv lstsq(A, phase_vec) velocity[p] x[0] * 365.25 / (4 * np.pi) * 0.056 # 转成年速率单位米/年 print(f平均沉降速率: {np.mean(velocity)*1000:.2f} mm/年)4 * np.pi是相位到视向位移的转换系数Sentinel-1A 的波长 0.056 米相位变化 2π 对应半个波长的位移。x[0] * 365.25是把日速率转成年速率。解算出来的速率是视向LOS的不是垂直向要做垂直向沉降分析需要升轨和降轨联合解算或者假设水平形变可忽略用入射角余弦做投影。3.4 大气相位校正GACOS 还是线性估计大气相位是 InSAR 时序分析里最大的误差源之一尤其是对流层延迟在山区能造成几厘米的假形变。常用做法有两种用外部气象数据GACOS做校正或者在时序解算里把大气相位当成低频信号估计掉。GACOS 需要下载对应区域的延迟图精度高但依赖外部数据线性估计不需要额外数据但对长波长大气信号抑制有限。# 下载 GACOS 数据需注册账号这里只示意命令结构 wget http://www.gacos.net/download?date20230101lat30.5lon114.3 \ -O ./gacos_20230101.ztd # 用 GACOS 校正干涉相位 gpt GACOSCorrection -Ssource./ifg_topo_removed.dim \ -PgacosFile./gacos_20230101.ztd \ -t ./ifg_gacos.dim如果研究区地形平坦、大气稳定线性估计够用如果地形起伏超过 500 米强烈建议用 GACOS。我踩过的坑是在山区做沉降监测没用 GACOS结果解算出来的沉降中心跟实际采空区对不上后来加了 GACOS 校正形变场才合理。4. 避坑与排查PSSBAS 落地时最容易翻车的五个地方4.1 相干性整体偏低PS 点稀疏到无法解算现象干涉图上看不到明显条纹相干性均值低于 0.3PS 候选点不到总像元的 1%。原因最常见的是时间基线太长或者研究区植被覆盖高C 波段在植被区失相干严重。另一个容易被忽略的原因是配准精度不够配准误差超过 0.2 像素相干性会断崖式下降。解决先把时间基线压到 60 天以内重新组合干涉对如果还是不行检查配准步骤用Back-Geocoding时把demResamplingMethod改成BICUBIC_INTERPOLATION比双线性更锐利。植被区可以考虑用 Sentinel-1B 的数据补充但 1B 已经退役只能靠 1A 的 12 天重访这时候 SBAS 的分布式散射体比 PS 更靠谱。4.2 解缠相位出现大面积跳变形变场不连续现象解缠后的相位图上有明显的 2π 跳变带形变场被切成几块速率图上有不自然的断层。原因解缠算法在低相干区域失败或者去平地效应不干净残留的地形相位被解缠算法误判成形变。SNAPHU 的DEFO模式对连续形变场友好但如果相干性低于 0.2任何解缠算法都救不回来。解决先掩掉相干性低于 0.25 的区域再解缠。如果跳变带出现在地形起伏大的地方换更高精度的 DEM比如 TanDEM-X 90 米。解缠后做一次相位闭合检查闭合相位超过 1 弧度的干涉对直接剔除。4.3 沉降速率出现季节性波动误判为地下水开采现象解算出来的形变时间序列有明显的年周期夏季抬升、冬季沉降幅度 5 到 10 毫米。原因这是典型的冻胀或土壤水分变化引起的季节性形变不是地下水开采。如果研究区在季节性冻土区或者浅层土壤含水量变化大这个信号会很强。解决在形变模型里加入周期项把年周期和半年周期一起解算分离出长期趋势。如果周期项幅度大于趋势项说明这个点的形变主要由季节性因素控制做沉降分析时要单独标注。# 加入周期项的形变模型 # v * dt a * sin(2*pi*dt/365.25) b * cos(2*pi*dt/365.25) c A np.zeros((n_ifg, 4)) for k, (d1, d2, dt, db) in enumerate(pairs): A[k, 0] dt A[k, 1] np.sin(2 * np.pi * dt / 365.25) A[k, 2] np.cos(2 * np.pi * dt / 365.25) A[k, 3] 1.04.4 升轨降轨结果对不上垂直向沉降算错现象升轨解算的 LOS 速率和降轨解算的 LOS 速率差异很大投影到垂直向之后沉降量差了一倍。原因升轨和降轨的入射角、方位角不同对水平形变的敏感方向也不同。如果研究区有显著的水平形变比如滑坡、地铁施工引起的侧向位移只用单轨数据投影到垂直向会引入很大误差。解决用升轨和降轨联合解算解出东西向和垂直向两个分量。如果只有单轨数据在报告里明确说明这是 LOS 向速率不是垂直向沉降。我一般会在成果图里同时标注 LOS 向和垂直向避免用户误读。4.5 参考点选错整个形变场整体偏移现象解算出来的形变场整体抬升或整体沉降没有明显的沉降中心或者沉降中心跟已知采空区位置偏移。原因时序解算需要一个参考点假设参考点形变为零。如果参考点本身在沉降整个形变场会整体偏移。参考点选在基岩出露区或者已知稳定的 GPS 站附近最可靠。解决先用 GPS 或水准数据验证参考点的稳定性如果没有外部数据选在形变场边缘、相干性高、ADI 低的点。解算后检查参考点附近的形变速率如果参考点周围 1 公里内有明显形变梯度说明参考点选得不对换点重算。5. 从速率图到可交付成果验证方法与一个实用技巧5.1 用交叉验证判断沉降速率是否可信PSSBAS 解算出来的速率图不能直接交付得先做验证。最可靠的方法是用独立数据交叉验证GPS 站的时间序列、水准测量点、或者另一套 InSAR 处理流程的结果。我一般会留出 20% 的 PS 点不参与解算作为验证集看预测速率和实际相位的残差。# 交叉验证留出 20% 干涉对做验证 np.random.seed(42) n_val int(0.2 * n_ifg) val_idx np.random.choice(n_ifg, n_val, replaceFalse) train_idx np.setdiff1d(np.arange(n_ifg), val_idx) # 用训练集解算速率 A_train A[train_idx] # 用验证集计算残差 residual phase_stack[val_idx] - A[val_idx] x rmse np.sqrt(np.mean(residual**2)) print(f验证集 RMSE: {rmse:.3f} 弧度)RMSE 小于 0.5 弧度说明形变模型拟合得不错大于 1 弧度说明模型漏掉了某些形变信号比如非线性沉降或大气残差。这时候要回头检查干涉对质量或者加入更多形变模型项。5.2 沉降速率图的出图参数与色标设置速率图是最终交付物色标设置直接影响可读性。沉降速率通常在 -50 到 10 毫米/年之间色标用发散型零值附近用白色或浅灰负值沉降用暖色正值抬升用冷色。如果研究区沉降速率超过 100 毫米/年色标要分段否则小沉降区域会被压成一片。import matplotlib.pyplot as plt from matplotlib.colors import TwoSlopeNorm fig, ax plt.subplots(figsize(10, 8)) norm TwoSlopeNorm(vmin-50, vcenter0, vmax10) im ax.imshow(velocity * 1000, cmapRdBu, normnorm) plt.colorbar(im, labelLOS 向速率 (mm/年)) ax.set_title(Sentinel-1A PSSBAS 地面沉降速率图) plt.savefig(./subsidence_rate.png, dpi300, bbox_inchestight)TwoSlopeNorm让零值居中负值和正值对称分布比线性色标更直观。vmin-50和vmax10是根据研究区实际速率范围调的如果沉降中心超过 -50把 vmin 调到 -100但要在图例里注明。5.3 一个实用技巧用时间序列聚类识别沉降模式速率图只能看平均趋势看不出沉降的时间演化。我习惯对 PS 点的时间序列做聚类把沉降模式相同的点归成一类这样能区分「线性沉降」「加速沉降」「季节性波动」三种典型模式。用 KMeans 对标准化后的时间序列聚类k 取 3 到 5聚类结果叠加到速率图上能直接看出哪些区域的沉降在加速。from sklearn.cluster import KMeans from sklearn.preprocessing import StandardScaler # 提取每个 PS 点的时间序列 ts np.zeros((n_ps, n_ifg)) for p in range(n_ps): ts[p] phase_stack[:, p] # 标准化后聚类 scaler StandardScaler() ts_norm scaler.fit_transform(ts) kmeans KMeans(n_clusters4, random_state42, n_init10) labels kmeans.fit_predict(ts_norm) # 统计每类的平均速率 for k in range(4): mask labels k print(f类别 {k}: 点数 {np.sum(mask)}, 平均速率 {np.mean(velocity[mask])*1000:.2f} mm/年)这个技巧帮我省了很多解释成本用户看到速率图只知道哪里沉看到聚类图能知道哪里在加速沉。加速沉降的区域才是真正需要现场核查和干预的。做这套流程最大的教训是不要跳过相干性检查直接解算。我早期跑一个矿区项目急着出图没看相干性就解缠结果解出来的沉降中心跟采空区偏了 2 公里后来回头查相干性图发现采空区上方植被茂密相干性只有 0.15解缠结果全是噪声。从那以后我养成了一个习惯每生成一景干涉图先看相干性均值低于 0.3 的直接标记为低质量不进入后续解算。这个习惯比任何参数调优都管用。希望帮到你。本文还有配套的精品资源点击获取