基于振幅的空间变迹滤波：一种用于大坝形变监测的PSI旁瓣抑制新方法

《Sensors》：Novel Amplitude-Based Approach for Reducing Sidelobes in Persistent Scatterer Interferometry Processing Using Spatially Variant Apodization

本文提出了一种基于振幅的空间变迹滤波（SVA）新方法，用于抑制合成孔径雷达（SAR）数据中的旁瓣。该方法在SNAP2StaMPS工作流中，仅对振幅进行滤波，同时保留原始干涉相位，从而在保证形变测量精度的前提下，有效减少了由角反射器（ECR-C）引起的虚假永久散射体（PS）探测，为大坝等基础设施的毫米级形变监测提供了更可靠的技术支撑。

　　

文章内容归纳总结

1. 引言

合成孔径雷达（SAR）是一种强大的遥感技术，它通过发射微波脉冲并记录从地球表面反向散射的回波，能够实现全天候、全天时的对地观测。干涉合成孔径雷达（InSAR）通过比较在不同时间或从略微不同位置获取的至少两幅SAR图像的相位信息，来测量地表形变和地形。在干涉测量方法中，永久散射体干涉测量（PSI）能够识别长期相干的散射体，即永久散射体（PS），这些散射体在时间上保持相位稳定性。这使得PSI在城市环境中特别有效，因为大坝等人造结构可以作为稳定的目标。通过分析同一区域多次SAR采集的相位差，PSI在形变测量中可以达到毫米级的精度。Sentinel-1（S-1）SAR任务以其6至12天的重访周期和视线（LOS）成像几何，特别适合使用PSI进行长期监测。

角反射器（CR）通过提供高度稳定的参考点来增强PSI应用。在研究地点，安装了一个C波段电子角反射器（ECR-C），用于与C波段SAR系统（如S-1）配合使用。这种有源反射器（AR）能够实现毫米级精度的精确形变监测，并通过响应发射脉冲和放大返回回波来主动与SAR信号交互。然而，尽管有这些优点，其强大的后向散射会引入旁瓣伪影，在SAR数据中表现为明亮的星状交叉图案。

旁瓣的产生是因为能量不仅辐射到SAR天线的主瓣中，也辐射到非预期的方向。它们由系统的冲激响应（IPR）描述，该响应表征了点目标在SAR图像中的表现。典型的IPR遵循sinc函数模式，具有一个尖锐的主瓣，周围环绕着逐渐减弱的旁瓣。sinc函数具有一个包含大部分SAR信号能量的主导中心峰（主瓣），以及从主瓣向外延伸的振荡（旁瓣）。旁瓣的产生是因为空间域中的矩形孔径在频域中变换为sinc函数。然而，旁瓣可能会扭曲或掩盖周围的信号信息，导致PSI导出的形变时间序列被误解。因此，旁瓣抑制是SAR图像分析中至关重要的预处理步骤。

2. 研究现状与创新方法

2.1. 旁瓣抑制与变迹技术

变迹是指对SAR信号应用加权函数以控制系统IPR形状的技术。各种技术已被开发用于减少SAR图像中的旁瓣伪影。线性变迹方法（如Hamming、Hanning和Blackman窗）应用固定的加权函数来减少旁瓣，但会通过主瓣展宽来折损图像分辨率。非线性变迹方法（如双变迹、复双变迹、三变迹和空间变迹）通过基于局部图像特征动态调整加权，而不是在整个图像上应用统一的函数，提供了更先进的旁瓣抑制，同时保留了分辨率。其中，空间变迹（SVA）是文献中研究最广泛的方法，因为它能够以逐像素为基础自适应滤波器权重，为每个像素分配一个单独的加权函数，以最优地减少旁瓣干扰。

然而，SVA在干涉测量背景下的应用相对较少被探索。现有研究表明，当应用于干涉图像对时，SVA通过增加采集间的相干性来提高图像分辨率和相位精度，但同时会降低相位信息。此外，SVA能够更针对性地选择PS候选点。一项进一步的发展涉及仅使用SVA来识别受旁瓣影响的像素，以创建旁瓣风险掩膜（SLRM）。SLRM标记受影响的区域，将其从进一步的滤波中排除，确保只有可靠的像素对最终的InSAR分析做出贡献。尽管如此，使用SLRM并没有利用SVA作为直接图像域滤波器的潜力，而将SVA直接应用于干涉数据则存在相位退化的风险，这对于PSI等相位敏感的应用至关重要。

2.2. 研究思路

为了解决这一研究空白，本文提出了一种基于振幅的SVA新方法，以保持相位完整性。关注振幅的动机在于PSI软件StaMPS（Stanford Method for Persistent Scatterers）使用振幅分散指数进行PS候选点的初始选择。通过在此阶段抑制与旁瓣相关的PS候选点，可以防止它们传播到后续的相位驱动分析中。SVA滤波器在PSI处理过程中进行了测试，测试对象是一个安装了ECR-C的大坝，该反射器会产生强烈的旁瓣，导致周围区域出现错误的PS探测。SVA被集成到已建立的S-1 PSI工作流中，通过SNAP2StaMPS（SentiNel Application Platform to Stanford Method for Persistent Scatterers）将SNAP中的预处理与StaMPS软件中的进一步PSI分析相结合。这种新颖的基于振幅的SVA方法仅应用于振幅，而原始相位信息保持不变。为了实现这一点，在滤波后重新组合了复SAR信号，确保旁瓣抑制不会损害相位相干性。

3. 材料与方法

3.1. 研究区域

索尔佩大坝（Sorpe Dam）是本次案例研究的研究区域。它是一个土石坝，于1926年至1935年间由鲁尔河谷水库协会在德国绍尔兰地区建造。该大坝的总库容约为7000万立方米，满库水位为海拔283.30米，坝顶高度为海拔285.60米。在满库容时，水库表面积约为3.3平方公里。为了防洪，大坝设有一个由带阶梯的溢流堰和一个泄洪闸组成的高水位泄洪系统，泄洪能力为每秒100立方米。

在鲁尔协会管理的所有大坝中，索尔佩大坝表现出最明显的旁瓣效应。这使其成为测试PSI分析中旁瓣抑制的理想场景，这归因于大坝相对于卫星过境的朝向以及ECR-C的工作原理的结合。旁瓣的影响以虚假PS点的形式可见，这些点从坝顶延伸至水中近350米。

3.2. 数据

本研究利用了20个S-1 Level-1产品，这些产品是在干涉宽幅（IW）模式下使用TOPSAR技术获取的。单视复（SLC）影像拍摄于2023年5月13日至2023年12月27日期间，在仅S-1A运行的下降轨道上，遵循12天的重访周期。选择下降轨道是因为在研究区域，下降轨道比上升轨道表现出更明显的旁瓣。这归因于特定的观测几何形状：S-1下降轨道从东向西移动，而坝顶从东南向西北方向延伸。这种几何排列与ECR-C的信号放大相结合，产生了强烈的后向散射，导致在SAR图像中观察到明显的旁瓣。

此外，本研究使用了来自航天飞机雷达地形测绘任务（SRTM）的1弧秒全球数字高程模型（DEM）。该数据集在PSI预处理期间用于配准，根据其轨道元数据和SRTM DEM对齐参考和次级的S-1 SLC分割产品。

3.3. 软件

本研究使用由欧空局（ESA）开发的SNAP，这是一个用于处理SAR和光学遥感数据的开源工具箱。SNAP（版本9）用于按照SNAP2StaMPS工作流对S-1数据进行预处理，以进行进一步的PSI处理。

SNAP2StaMPS工作流在SNAP中准备TOPSAR数据，以便与StaMPS PSI处理工作流无缝集成。该工作流围绕SNAP处理图和一组Python脚本构建，这些脚本充当包装器，允许用户通过配置文件配置关键处理参数并自动执行工作流。

StaMPS（版本StaMPS-4.1-beta）是一个基于MATLAB的软件包，用于PSI处理。StaMPS处理由SNAP生成的干涉图堆栈，并应用统计方法来识别用于位移监测的稳定PS。

3.4. SAR像素的笛卡尔和极坐标表示

图像形成后，最终的SAR产品表示为复SAR图像，包含两个分量：同相（I）分量代表信号的实部，而正交（Q）分量对应于虚部。笛卡尔（I, Q）和极坐标（振幅，相位）表示之间的关系由以下公式给出：

• •
  振幅 = √(I²+ Q²)
• •
  相位 = arctan(Q / I)

振幅和相位提供了SAR信号的直观理解，而I和Q分量简化了计算。这种关系对于理解本文提出的新颖的基于振幅的方法至关重要，该方法通过重新组合I和Q分量来保留原始相位，同时仅对振幅应用SVA滤波。

3.5. 空间变迹（SVA）

SVA通过识别主瓣峰值的像素并减少其周围的旁瓣来减少旁瓣。通常，SVA使用三点卷积，其中每个像素受其自身值及其两个相邻像素值的影响。SVA采用余弦加基座加权，其中余弦加基座参数w从0变化到0.5。加权函数的选择（例如Hamming、Hanning）直接控制主瓣宽度和旁瓣抑制之间的权衡。这些窗对应于w的固定值（例如，w=0表示矩形窗，w=0.5表示Hanning窗，w≈0.43表示Hamming窗）。相比之下，SVA不依赖于预定义的窗函数。相反，它评估w在0≤w≤0.5范围内的所有可能值，并为每个像素选择使输出能量最小化的值。

本研究应用的SVA方法将Wang等人（2012年）的现有方法扩展到二维（2D）处理，结合了方位向和距离向。因此，SVA分别应用于每个SAR图像的I和Q分量。计算出的加权函数w_u(m,n)用于像素(m,n)，使用其四个直接相邻像素计算。如果任何相邻像素包含NaN值，则将其从计算中排除。0.5的阈值直接来自规范1D SVA中使用的余弦加基座加权理论。所提出的2D实现将原始1D模板推广到最小的对称2D邻域（四个直接邻居），同时保留了w的理论约束。

使用计算出的加权函数，最终的SVA滤波按如下方式应用：

• •
  如果w_u(m,n) < 0，主瓣保持不变，意味着像素保留其原始值。
• •
  如果0 ≤ w_u(m,n) ≤ 1/2，通过将像素值设为零来减少旁瓣。
• •
  如果w_u(m,n) > 1/2，像素值被替换为其四个有效相邻值的平均值。

3.6. 新颖的基于振幅的方法

所提出的基于振幅的方法建立在SVA滤波之上。首先，如第3.5节所述，将SVA应用于I和Q分量，产生场景的旁瓣减少版本。在第二步中，通过将原始数据（SVA滤波前）的相位信息与SVA滤波数据的振幅相结合，重新计算新的I和Q分量。修改后的I和Q分量由以下公式给出：

• •
  I_modified= I_orig× (√(I_reduc²+ Q_reduc²) / √(I_orig²+ Q_orig²))
• •
  Q_modified= Q_orig× (√(I_reduc²+ Q_reduc²) / √(I_orig²+ Q_orig²))

其中，orig指原始未滤波数据，reduc代表旁瓣减少的对应数据。这些表达式源自第3.4节讨论的笛卡尔（I, Q）和极坐标（振幅，相位）表示之间的关系，确保在修改振幅以使用SVA滤波器减少旁瓣的同时，保留了原始相位。

4. 工作流

4.1. 使用SNAP2StaMPS进行标准预处理

传统SNAP2StaMPS工作流从整个时间序列中选择参考影像开始。然后，分别对参考产品和次级产品执行TOPSAR-Split和Apply-Orbit-File步骤。这些步骤使用SNAP2StaMPS内预定义的处理图自动执行。配准和干涉图计算使用第二个处理图执行，该图在单个过程中执行这两个步骤。Back-Geocoding算子使用其轨道元数据和SRTM 1弧秒全球DEM对齐两个S-1 SLC分割产品（参考和次级）。在Back-Geocoding之后，Enhanced-Spectral-Diversity步骤使用非相干互相关来估计整个子带的恒定距离偏移，以确定最终校正，确保完全配准。在此阶段，工作流分为两个分支，每个分支产生不同的中间输出。在第一个分支中，下一步是干涉图生成，计算复干涉图。随后是TOPSAR-Deburst步骤，将脉冲串合并为连续图像。第二个处理分支遵循相同的工作流，但在此阶段终止。在第一个处理分支中，附加步骤Topographic-Phase-Removal通过将SRTM 1弧秒全球DEM转换为SAR坐标并从干涉图中减去相应的相位来消除地形相位。该图完成后，生成两个中间输出：干涉图和配准图像对，从而结束工作流的这一阶段。SNAP2StaMPS工作流的最后一步，StaMPS导出，准备与StaMPS兼容格式的数据。

4.2. 使用StaMPS进行标准PSI处理

StaMPS工作流的方法由Hooper等人（2018年）描述。SNAP2StaMPS工作流生成的StaMPS导出产品作为此阶段的输入数据。在启动StaMPS的PSI处理之前，使用振幅分散指数执行PS候选点的初始选择。PSI处理在MATLAB环境中进行，遵循步骤1到7，并利用TRAIN软件，如Hooper等人（2018年）所述。在处理过程中估计时间相干性，并通过将观测数据与随机生成的相位数据进行比较，基于噪声特征概率性地选择像素。StaMPS的最终结果是一组具有坐标和形变时间序列的PS候选点。

4.3. 在SNAP2StaMPS工作流中实现SVA

SVA滤波在Back-Geocoding和Enhanced Spectral Diversity之后的配准之后应用。SVA在逐行基础上对复SAR数据进行操作。其有效性取决于精确对齐的输入。配准确保参考和次级采集在几何上一致。因此，SVA在配准后立即应用，允许滤波器在几何对齐且尽可能接近原始采集的数据上操作。由于SNAP本身不支持SVA，滤波通过snappy接口通过自定义Python脚本在外部实现。为了合并SVA滤波器，传统的SNAP2StaMPS工作流的第二个图被分成两部分，从而能够在配准和干涉图生成阶段之间在Python中外部应用SVA滤波器。

4.4. 在SNAP2StaMPS工作流中实现新颖的基于振幅的方法

新颖的基于振幅的方法的基本原理是，StaMPS最初基于振幅分散指数选择PS候选点，而相位特征仅在PSI处理期间被纳入。通过仅对振幅应用SVA滤波，减少了受旁瓣影响的像素并最小化了错误的PS候选点，同时保留了原始相位以确保基于相位信息的可靠PS选择。为了实现基于振幅的方法，传统的第二个SNAP2StaMPS处理图被重组为三步过程。第一步，执行标准处理直至旁瓣减少。然后通过snappy接口应用SVA滤波。在此阶段，存储每个场景的原始配准（未滤波）数据和SVA滤波数据，因为基于振幅的处理都需要这两者。第二步，SNAP中的Band Maths工具实现了基于振幅的方法：对于每个场景，Band Maths加载原始数据和SVA滤波数据，并通过重新计算I和Q分量，将原始相位与SVA滤波的振幅重新组合。选择Band Maths是因为它在SNAP中可用，最大限度地减少了将操作外包给外部软件的需要。最后，在Band Maths步骤之后，第二个处理图恢复，继续标准的SNAP2StaMPS工作流。

5. 结果与讨论

5.1. 结果

为了评估所提出方法的有效性，分析了三种处理变体：（1）原始处理数据，缩写为Original；（2）对整个场景应用了传统SVA滤波的数据，影响振幅和相位（SVA_entireScene）；（3）仅对振幅应用了SVA滤波的数据，保留原始相位（SVA_amplitudeBased）。为了进行比较，检查了PS点的空间分布、形变模式、IPR以及振幅和相位特征。

PS点的空间分布显示，SVA_entireScene和SVA_amplitudeBased都比Original方法产生了更少的PS点，空间分布只有微小差异。AOI1（红色）的结果表明，SVA_entireScene和SVA_amplitudeBased有效地移除了水上的PS点，而坝顶（AOI2，黄色）的大多数PS点被保留。比较SVA_entireScene和SVA_amplitudeBased，只观察到微小差异：SVA_amplitudeBased在AOI2中检测到两个额外的PS点，并在空间分布上表现出轻微的视觉变化。大多数使用SVA_entireScene和SVA_amplitudeBased检测到的PS也包含在Original选择中。所有来自SVA_amplitudeBased的PS点也被Original方法检测到，而SVA_entireScene没有识别出四个PS点。尽管如此，大多数PS点是两种方法共有的。

PSI结果通过绘制代表ECR-C信号的PS点的时间形变信号进一步分析。两种SVA变体SVA_entireScene（深蓝色）和SVA_amplitudeBased（浅蓝色）与Original（灰色）略有不同，在9月份偏差达到0.7毫米。Original和SVA_entireScene之间的均方根误差（RMSE）为0.388毫米，Original和SVA_amplitudeBased之间的RMSE为0.384毫米。SVA_amplitudeBased时间形变模式与SVA_entireScene的对齐度比Original信号更接近。为了进一步验证这些结果，分析了索尔佩大坝坝顶的另一个PS点。SVA_entireScene在第二个位置显示出更大的偏差，在9月份达到近3毫米，而SVA_amplitudeBased方法显示出最小的偏差，并保持在接近Original信号的位置。Original和SVA_entireScene之间的RMSE为0.947毫米，Original和SVA_amplitudeBased之间的RMSE为0.386毫米。

为了进一步验证ECR-C信号，分析了Original和SVA_entireScene版本的振幅信号的IPR。选择2023年5月15日的S-1影像，因为它是时间序列中的第一个。由于分析仅基于振幅，SVA_amplitudeBased方法的结果与SVA_entireScene相同。对于Original和SVA_entireScene，识别了最亮的像素（对应于ECR-C信号），并在方位向和距离向提取了150个像素的窗口。1D IPR轮廓使用三次插值过采样10倍。振幅存储并处理为32位浮点数组。轮廓归一化到其最大值并转换为对数刻度（dB）。由此产生的IPR轮廓，以ECR-C峰值（位置0）为中心，在距离向和方位向显示，突出了旁瓣抑制和主瓣特征的差异。在距离向，SVA_entireScene（深蓝色）在像素105附近显示出较低的旁瓣水平，而它主要遵循Original信号（灰色）的形状。与Original的平均偏差为0.91 dB。在方位向，差异变得更加明显。SVA_entireScene在几个位置偏离更多，其主瓣显得更清晰。这里的平均偏差为4.37 dB。

振幅图像、缠绕相位干涉图、逐像素相位差以及模拟IPR的综合概览显示，视觉上，两种SVA方法SVA_entireScene和SVA_amplitudeBased与Original振幅相比，都表现出明显的旁瓣强度降低，与模拟的旁瓣减少一致。SVA_entireScene和SVA_amplitudeBased的振幅图像看起来相同，因为两者都对振幅分量应用了相同的滤波。关于相位，SVA_amplitudeBased在很大程度上保留了Original干涉图的相位结构，仅显示出微小的偏差。相比之下，SVA_entireScene在相位模式中引入了更明显的改变。为了量化这些差异，计算了平均相位值（μ）和标准差（σ）。SVA_entireScene与Original平均值的偏差高达0.02弧度，而SVA