PS-InSAR技术实战：基于StaMPS的永久散射体形变监测全流程解析

1. 从“点”开始：理解PS-InSAR的核心价值

在InSAR（合成孔径雷达干涉测量）的世界里，我们通常处理的是整片区域的形变信息，比如一幅SAR影像覆盖的几十上百平方公里。但很多时候，我们真正关心的，是这片区域内那些“不动”的点——比如高楼上的角反射器、裸露的岩石、或者长期稳定的建筑物。这些点，我们称之为永久散射体（Permanent Scatterers， PS）。PS-InSAR技术，就是专门从海量像元中，把这些“钉子户”找出来，并精确计算它们长时间序列形变的方法。如果说前两篇我们搭建了环境、跑通了流程，是在铺路和造车，那么这篇PS处理，就是教你如何在这条路上，精准地找到并追踪那些最值得信赖的“路标”。

为什么PS如此重要？想象一下，你用InSAR监测一个城市的地面沉降。整幅影像里，有农田、树林、湖泊，这些地方因为植被生长、水位变化，雷达信号回波强度随时间剧烈变化，信噪比很低，计算出的形变结果可能跳来跳去，不可信。但城市里那些钢筋混凝土的建筑屋顶、桥梁的特定结构，它们对雷达信号的反射非常稳定，几乎不受时间和天气影响。这些点就是PS点。通过分析这些稳定点的相位变化，我们能够剥离掉大气延迟、轨道误差等公共误差，得到毫米级甚至亚毫米级精度的形变时间序列。这对于监测大坝、桥梁、高铁沿线、矿区沉降、滑坡体蠕动等，意义重大。

ISCE2负责生成精确的干涉图堆栈，而StaMPS则是PS-InSAR处理的王牌工具。本篇，我们就聚焦于StaMPS 4.1的PS处理核心流程。网络上关于StaMPS的教程不少，但大多停留在命令罗列。我将结合我处理上百景哨兵数据积累的经验，重点拆解每个步骤背后的物理意义和参数设置的“所以然”，并分享那些手册里不会写、但实际操作中一定会遇到的“坑”。

2. 启程前的最后检查：数据与环境确认

在正式启动StaMPS的PS流程之前，我们必须确保从ISCE2产出的“原料”是合格且准备就绪的。很多人在这一步栽跟头，不是因为StaMPS命令复杂，而是输入数据本身就有问题。

2.1 干涉图堆栈的完整性验证

首先，回到你的ISCE2处理目录。假设你的项目目录是SenDT/，里面应该有一个merged/文件夹，存放着所有配准后的SLC（单视复数影像）和生成的干涉图对。

你需要检查两个关键文件：

1. baselines文件：这个文件记录了每一对干涉图的时空基线（垂直基线和时间基线）。用cat baselines命令查看。确保行数等于你生成的干涉图数量，并且没有出现NaN（非数字）或异常大的值。一个常见的坑是，如果配准步骤有某景影像失败，但流程仍继续，可能导致基线计算错误，进而影响后续相位解缠。
2. date_list.txt或类似的主影像日期列表文件：这个文件列出了所有SLC影像的日期。StaMPS需要根据这个列表来组织时间序列。确保日期格式正确（通常是YYYYMMDD），并且顺序与baselines文件中的主影像日期对应。

注意：ISCE2的stackSentinel.py脚本通常会自动生成这些文件。但如果你是自己手动组干涉对，务必确保baselines文件的格式是StaMPS可读的（通常是：主影像日期、从影像日期、垂直基线、时间基线）。

2.2 StaMPS工作目录的初始化

StaMPS处理需要在独立的工作目录中进行，避免污染ISCE2的原始数据。我通常的做法是：

cd /path/to/your/area mkdir StamPS_PS cd StamPS_PS # 将ISCE2的关键输出链接过来 ln -s ../SenDT/merged/merged . ln -s ../SenDT/merged/baselines . ln -s ../SenDT/merged/date_list.txt .

这里，merged文件夹的链接是关键。StaMPS会从这个文件夹里读取所有配准后的SLC（*/*.slc.full）和干涉图（*/*.int）。使用软链接而不是复制，可以节省大量磁盘空间。

2.3 环境变量与Matlab路径配置

StaMPS 4.1运行依赖于Matlab。你需要确保两件事：

1. Matlab可执行文件路径：在终端中，which matlab应该能返回路径。如果没有，需要在你的shell配置文件（如~/.bashrc）中添加export PATH=/path/to/matlab/bin:$PATH。
2. StaMPS的Matlab工具箱路径：这是最容易出错的地方。你需要在Matlab的启动脚本（startup.m）或者直接在StaMPS的配置文件里，添加StaMPS和其依赖工具（如snaphu）的路径。

一个更稳妥的方法是，在运行StaMPS的Matlab脚本前，在终端里临时设置Matlab路径：

export MATLABPATH=/path/to/StaMPS:/path/to/StaMPS/matlab:/path/to/snaphu

然后通过matlab -nodesktop -nosplash -r “stamps(1,1)”这样的命令启动，并在Matlab命令行里再次用addpath确认路径已添加。我个人的习惯是在StamPS工作目录下创建一个小的setup_env.m脚本，里面写好所有addpath命令，每次启动Matlab后先运行它。

3. 核心第一步：相位校正与噪声估计

一切就绪，我们开始运行StaMPS。第一步通常是stamps(1,1)。这个步骤看似简单，实则包含了多个关键操作。

3.1 相位校正的物理意义

从ISCE2生成的干涉图，其相位包含了几何相位（地形相位）、形变相位、大气相位、轨道误差相位和噪声。stamps(1,1)首先会利用外部DEM（来自ISCE2处理）去除地形相位。这一步之后，干涉图中剩余的相位主要就是形变、大气和噪声了。

但这里有一个至关重要的细节：多普勒质心频率差异校正。哨兵数据是TOPS（Terrain Observation with Progressive Scans）模式，不同时刻获取的影像，即使经过配准，其多普勒质心也可能有微小差异。这个差异会引入一个与距离向坐标成线性关系的相位项。如果不校正，它会污染后续的大气相位估计，尤其是在像幅边缘。StaMPS的这一步会自动估计并移除这个线性相位趋势。你需要关注日志输出，看校正量是否在合理范围内（通常很小）。如果发现校正量异常大，可能预示着配准质量有问题。

3.2 噪声估计与像素初选

校正后，StaMPS会开始估计每个像素点的相位噪声水平。它使用一个基于空间相关性的模型。简单理解，如果一个像素点周围的像素相位都很杂乱（不相干），那么这个点本身的相位噪声就大；反之，如果周围像素相位平滑一致，噪声就小。

基于这个噪声估计，StaMPS会进行第一轮像素筛选。它会计算每个像素的“相位稳定性”指标，并设定一个阈值（如默认的0.3）。高于这个阈值的点，被认为是潜在的“候选PS点”。这一步会淘汰掉绝大部分像元（可能95%以上），只留下那些相位相对稳定的点进入后续处理。

实操心得：这个阈值（weed_standard_dev）不要轻易改动。调低它会纳入更多点，但也会引入更多噪声，增加后续计算负担和误判风险。除非你处理的是特别贫瘠的岩石山区，信号整体都很差，否则保持默认是稳妥的选择。你可以通过后续步骤查看候选点的密度图来评估初选效果。

4. 相位解缠：从缠绕相位到绝对形变

这是PS-InSAR中最核心、也最考验算法功力的步骤，对应stamps(2,2)。经过第一步，我们得到了许多候选PS点，但它们的相位值是被“缠绕”在[-π, π]区间内的。我们需要把这些缠绕的相位“解开”，恢复其真实的、连续的相位值。

4.1 三维相位解缠的挑战

传统的二维相位解缠（比如处理单幅干涉图）已经很难。PS-InSAR是三维相位解缠：在空间（x,y）和时间（t）三个维度上同时进行。难点在于：

1. 空间不连续：PS点是离散分布的，不像连续区域那样有明确的相邻关系。
2. 时间基线网络复杂：干涉图对之间形成的是一个复杂的网络，有的时间间隔长，有的短，解缠需要在时间维度上保持一致性。
3. 高噪声：尽管经过了筛选，候选PS点的相位仍包含噪声，可能在某些干涉对上出现“残差”。

StaMPS采用了一种非常聪明的方法：它不直接对每个PS点进行三维解缠，而是先利用所有干涉图的信息，估计出一个“最可能”的相位时间序列模型（包括线性形变速率和非线性形变），然后基于这个模型去指导每个点的二维空间解缠。

4.2 关键参数解析与设置

运行stamps(2,2)前，通常需要修改parms结构体中的一些参数。在Matlab命令行中操作：

% 加载参数 load(‘parms.mat’) % 查看当前参数 parms % 修改关键参数 parms.llook = 20; % 多视比（距离向）。应与ISCE2生成干涉图时的一致！ parms.n_win = 32; % 空间滤波窗口大小。用于估计空间相关性，默认32通常够用。 parms.grid_size = 100; % 解缠用的网格大小（米）。城市区域可设小点（如50），山区可设大点（200）。 parms.unwrap_method = ‘3D’; % 解缠方法。‘3D’是推荐的核心算法。 % 保存修改 save(‘parms.mat’， ‘parms’)

重点解释parms.llook：这个参数必须与你在ISCE2的stackSentinel.py中设置的多视比完全一致！如果ISCE2你用了20:4（距离向:方位向），那么parms.llook应该等于20。如果不一致，StaMPS在读取干涉图时会误判像素位置，导致所有后续处理都是错的。这是我踩过的最大的坑之一，现象是解缠后的相位图一片混乱，PS点位置完全不对。

关于parms.grid_size：StaMPS会将研究区域划分成一个个网格，在每个网格内分别进行相位解缠，然后再拼接起来。网格尺寸越小，计算越精细，但耗时越长，且在小网格内可能因PS点太少而解缠失败。对于城市区域，PS点密集，可以设置较小的网格（如50米）以获得更细节的解缠结果；对于山区或乡村，PS点稀疏，需要设置较大的网格（如100-200米）以保证每个网格内有足够的点进行可靠解缠。

4.3 解缠过程监控与常见问题

运行stamps(2,2)后，控制台会输出大量信息。你需要关注几点：

• “Percentage of pixels unwrapped”：成功解缠的像素百分比。理想情况下应该在90%以上。如果过低（比如低于70%），说明很多点的相位噪声太大，或者参数设置（特别是grid_size）不合适。
• 解缠迭代次数：StaMPS会迭代优化解缠结果。通常迭代几次后就会收敛。如果迭代次数非常多（>10次），可能意味着数据质量有问题，或者存在强烈的非线性形变信号。
• 程序会生成很多中间图，如phase_std.ps（相位标准差图）。用ps_plot(‘v-doi’…)等命令查看这些图，可以帮助你直观判断解缠质量。好的解缠结果，PS点的相位标准差应该比较低，且空间分布均匀。

常见问题与排查：

• 解缠结果出现条带状或块状异常：这通常是parms.llook设置错误导致的。立即检查并修正此参数，然后从stamps(1,1)重新开始。
• 大量PS点解缠失败：首先检查候选PS点密度图（ps_plot(‘d’））。如果密度本身就很低，可能是第一步的噪声阈值weed_standard_dev设得太高，或者研究区域本身缺乏稳定散射体（如茂密森林、水域）。如果密度正常但解缠失败，尝试增大parms.grid_size，或者检查干涉图堆栈中是否存在质量极差的干涉对（可通过查看ph_disp.ps等图辅助判断）。

5. 大气相位屏估计与剔除

成功解缠后，我们得到了每个PS点“绝对”的相位时间序列。但这个相位里还混着我们需要的大气延迟相位和形变相位。stamps(3,3)和stamps(4,4)就是用来分离它们的。

5.1 大气相位的时空特性

大气延迟（主要是对流层水汽）引起的相位误差，在空间上是低频变化的（平滑的“屏”），在时间上是高频变化的（与天气快速相关）。而地表形变，在空间上可以是高频的（单个建筑物沉降），在时间上通常是低频的（缓慢持续沉降）。

StaMPS利用这种特性差异来分离两者。stamps(3,3)首先会用一个高通时间滤波（比如滤除周期长于1年的信号）和一个低通空间滤波（比如滤除尺度小于1公里的变化），从解缠后的相位中初步估计出大气相位屏（APS）。

5.2 迭代优化与非线性形变估计

stamps(4,4)是一个迭代过程。它用估计出的APS去校正原始相位，然后重新估计形变（包括线性速率和非线性部分），再用残差相位更新APS估计，如此反复，直到收敛。

这里的关键是滤波器的设置（在parms中）：

• parms.filter_time：时间滤波器的截止周期。例如设为365，意味着认为周期大于365天（一年）的信号是形变，小于的是大气噪声。这个值需要根据你的研究区域形变特征来定。对于缓慢沉降，这个值可以设大一点；对于季节性形变明显的区域，要小心设置。
• parms.filter_space：空间滤波器的窗口大小（单位：米）。例如设为1000，意味着认为空间尺度大于1公里的变化是大气相关的。这个值通常与大气扰动的典型尺度有关，默认值（如1000-2000米）在多数情况下是合理的。

实操心得：分离大气和形变是PS-InSAR的精华，也是难点。没有绝对正确的参数。我的建议是：先用默认参数跑一遍全程。然后，重点分析结果中那些明显的、大范围的、与地形高度相关的相位图案。如果发现这样的图案，它很可能是残余的大气误差（因为水汽分布常与地形相关）。这时，你可以尝试略微减小parms.filter_space，让空间滤波更“激进”地移除大尺度信号，再重新运行stamps(4,4)。观察形变速率图是否变得更合理（比如，山区本应无显著形变的地方，速率值是否接近零了）。这是一个需要反复调试和验证的过程。

6. 最终产品生成与可视化解读

经过上述步骤，我们终于得到了“干净”的形变相位时间序列。stamps(5,5)会将这些相位转换为实际的地表形变量（单位：毫米），并生成最终的结果文件。

6.1 结果文件解读

处理完成后，工作目录下会生成几个关键文件：

• ps_plot_v-doi.eps：形变速率图（平均每年形变毫米数）。这是最常用的成果图。暖色（红、黄）通常表示远离卫星的形变（如沉降），冷色（蓝）表示靠近卫星的形变（如抬升）。
• ps_plot_v-doi.mat：包含所有PS点经纬度、形变速率、高程误差等数据的Matlab文件。
• ts_params.mat：包含时间序列形变数据。
• 一系列以日期命名的.mat文件：每个文件包含该日期所有PS点相对于参考日期的累积形变量。

6.2 使用MATLAB进行深度分析与制图

StaMPS自带了很多绘图函数，但为了发表或报告，我们通常需要更精美的定制化图表。这里分享一段我常用的MATLAB代码片段，用于提取单个PS点的时间序列并绘图：

load(‘ts_params.mat’) % 加载时间序列参数 load(‘ps_plot_v-doi.mat’) % 加载PS点信息 % 假设你想查看某个特定位置的点（例如经纬度 lon0, lat0） [~， idx] = min(abs(lon_mat-lon0) + abs(lat_mat-lat0)); % 找到最近点的索引 % 提取该点的形变时间序列（单位：毫米） d_cum = ts_params.d_cum(idx， :); % 累积形变 d = ts_params.d(idx， :); % 单个日期对的形变？这里需要注意，ts_params结构可能版本不同 % 更通用的方法是使用ph_disp（相位）进行转换 ph = ts_params.ph(idx， :); % 相位值 wavelength = 0.0555; % 哨兵1号C波段波长，单位米 d_cum_mm = -ph * wavelength / (4*pi) * 1000; % 转换为毫米，负号取决于相位符号约定 % 获取日期 date_list = ts_params.day; % 日期序列，可能是相对于参考日期的天数 % 需要将天数转换为实际日期 master_date = ‘20180101’; % 你的主影像日期，需自行替换 master_datenum = datenum(master_date， ‘yyyymmdd’); actual_dates = master_datenum + date_list; % 绘图 figure(‘Position’， [100， 100， 800， 400]) plot(actual_dates， d_cum_mm， ‘b-o’， ‘LineWidth’， 1.5， ‘MarkerFaceColor’， ‘b’) datetick(‘x’， ‘yyyy-mm’， ‘keepticks’) xlabel(‘Date’) ylabel(‘Cumulative Deformation (mm)’) title([‘PS Point at (‘， num2str(lon_mat(idx)， ‘%.4f’)， ‘， ‘， num2str(lat_mat(idx)， ‘%.4f’)， ‘)’]) grid on

这段代码能帮你深入分析特定点的形变过程，比如判断形变是匀速、加速还是存在突变。

6.3 结果验证与误差分析

得到形变图后，切勿直接下结论。必须进行交叉验证：

1. 与已知事实对照：研究区域内是否有已知的沉降区、滑坡点？你的结果是否与之吻合？
2. 检查空间模式：形变速率图是否显示出与地质构造、地下水开采区、重大工程活动相关的空间格局？如果形变图案杂乱无章，或与地形高度重合，可能暗示大气相位剔除不净。
3. 分析时间序列：随机选取一些PS点，绘制其时间序列。曲线应该是相对平滑的，符合物理过程。如果出现剧烈的、无规律的跳动，可能是该点本身不稳定，或者在处理环节（如相位解缠）出了问题。
4. 定量评估：计算整个区域PS点形变速率的统计值（均值、标准差）。在理论上稳定的区域（如基岩出露区），形变速率应接近于0，其标准差可以视为本次监测的精度水平。如果能控制在每年1-2毫米以内，说明处理质量很高。

最后，PS-InSAR的结果是相对形变，即每个点相对于一个“参考点”的形变。这个参考点通常是处理过程中自动或手动选择的一个假设稳定的点。你需要在成果中明确说明参考点的位置及其稳定性假设。如果可能，用现场水准测量或GPS数据对几个关键PS点进行绝对验证，是提升成果可信度的最佳方式。

整个PS处理流程，从数据检查、参数调试到结果验证，是一个需要耐心和经验的循环。它不像流水线点击按钮就能出完美结果。每一个参数背后都有其地球物理或数学意义，每一次调整都需要结合对研究区域的先验认知和对中间结果的细致判读。这份工作一半是科学，一半是艺术。当你第一次看到清晰的、符合预期空间格局的形变图从杂乱的数据中浮现出来时，那种成就感，正是我们从事技术工作的乐趣所在。