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单视复数影像和生成的干涉图对。你需要检查两个关键文件baselines文件这个文件记录了每一对干涉图的时空基线垂直基线和时间基线。用cat baselines命令查看。确保行数等于你生成的干涉图数量并且没有出现NaN非数字或异常大的值。一个常见的坑是如果配准步骤有某景影像失败但流程仍继续可能导致基线计算错误进而影响后续相位解缠。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。你需要确保两件事Matlab可执行文件路径在终端中which matlab应该能返回路径。如果没有需要在你的shell配置文件如~/.bashrc中添加export PATH/path/to/matlab/bin:$PATH。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处理去除地形相位。这一步之后干涉图中剩余的相位主要就是形变、大气和噪声了。但这里有一个至关重要的细节多普勒质心频率差异校正。哨兵数据是TOPSTerrain 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三个维度上同时进行。难点在于空间不连续PS点是离散分布的不像连续区域那样有明确的相邻关系。时间基线网络复杂干涉图对之间形成的是一个复杂的网络有的时间间隔长有的短解缠需要在时间维度上保持一致性。高噪声尽管经过了筛选候选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_sizeStaMPS会将研究区域划分成一个个网格在每个网格内分别进行相位解缠然后再拼接起来。网格尺寸越小计算越精细但耗时越长且在小网格内可能因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 结果验证与误差分析得到形变图后切勿直接下结论。必须进行交叉验证与已知事实对照研究区域内是否有已知的沉降区、滑坡点你的结果是否与之吻合检查空间模式形变速率图是否显示出与地质构造、地下水开采区、重大工程活动相关的空间格局如果形变图案杂乱无章或与地形高度重合可能暗示大气相位剔除不净。分析时间序列随机选取一些PS点绘制其时间序列。曲线应该是相对平滑的符合物理过程。如果出现剧烈的、无规律的跳动可能是该点本身不稳定或者在处理环节如相位解缠出了问题。定量评估计算整个区域PS点形变速率的统计值均值、标准差。在理论上稳定的区域如基岩出露区形变速率应接近于0其标准差可以视为本次监测的精度水平。如果能控制在每年1-2毫米以内说明处理质量很高。最后PS-InSAR的结果是相对形变即每个点相对于一个“参考点”的形变。这个参考点通常是处理过程中自动或手动选择的一个假设稳定的点。你需要在成果中明确说明参考点的位置及其稳定性假设。如果可能用现场水准测量或GPS数据对几个关键PS点进行绝对验证是提升成果可信度的最佳方式。整个PS处理流程从数据检查、参数调试到结果验证是一个需要耐心和经验的循环。它不像流水线点击按钮就能出完美结果。每一个参数背后都有其地球物理或数学意义每一次调整都需要结合对研究区域的先验认知和对中间结果的细致判读。这份工作一半是科学一半是艺术。当你第一次看到清晰的、符合预期空间格局的形变图从杂乱的数据中浮现出来时那种成就感正是我们从事技术工作的乐趣所在。
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单视复数影像和生成的干涉图对。你需要检查两个关键文件baselines文件这个文件记录了每一对干涉图的时空基线垂直基线和时间基线。用cat baselines命令查看。确保行数等于你生成的干涉图数量并且没有出现NaN非数字或异常大的值。一个常见的坑是如果配准步骤有某景影像失败但流程仍继续可能导致基线计算错误进而影响后续相位解缠。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。你需要确保两件事Matlab可执行文件路径在终端中which matlab应该能返回路径。如果没有需要在你的shell配置文件如~/.bashrc中添加export PATH/path/to/matlab/bin:$PATH。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处理去除地形相位。这一步之后干涉图中剩余的相位主要就是形变、大气和噪声了。但这里有一个至关重要的细节多普勒质心频率差异校正。哨兵数据是TOPSTerrain 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三个维度上同时进行。难点在于空间不连续PS点是离散分布的不像连续区域那样有明确的相邻关系。时间基线网络复杂干涉图对之间形成的是一个复杂的网络有的时间间隔长有的短解缠需要在时间维度上保持一致性。高噪声尽管经过了筛选候选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_sizeStaMPS会将研究区域划分成一个个网格在每个网格内分别进行相位解缠然后再拼接起来。网格尺寸越小计算越精细但耗时越长且在小网格内可能因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 结果验证与误差分析得到形变图后切勿直接下结论。必须进行交叉验证与已知事实对照研究区域内是否有已知的沉降区、滑坡点你的结果是否与之吻合检查空间模式形变速率图是否显示出与地质构造、地下水开采区、重大工程活动相关的空间格局如果形变图案杂乱无章或与地形高度重合可能暗示大气相位剔除不净。分析时间序列随机选取一些PS点绘制其时间序列。曲线应该是相对平滑的符合物理过程。如果出现剧烈的、无规律的跳动可能是该点本身不稳定或者在处理环节如相位解缠出了问题。定量评估计算整个区域PS点形变速率的统计值均值、标准差。在理论上稳定的区域如基岩出露区形变速率应接近于0其标准差可以视为本次监测的精度水平。如果能控制在每年1-2毫米以内说明处理质量很高。最后PS-InSAR的结果是相对形变即每个点相对于一个“参考点”的形变。这个参考点通常是处理过程中自动或手动选择的一个假设稳定的点。你需要在成果中明确说明参考点的位置及其稳定性假设。如果可能用现场水准测量或GPS数据对几个关键PS点进行绝对验证是提升成果可信度的最佳方式。整个PS处理流程从数据检查、参数调试到结果验证是一个需要耐心和经验的循环。它不像流水线点击按钮就能出完美结果。每一个参数背后都有其地球物理或数学意义每一次调整都需要结合对研究区域的先验认知和对中间结果的细致判读。这份工作一半是科学一半是艺术。当你第一次看到清晰的、符合预期空间格局的形变图从杂乱的数据中浮现出来时那种成就感正是我们从事技术工作的乐趣所在。
)
