SARscape5.7中DS-InSAR：E-PS操作说明

PS-InSAR方法测量PS点的形变，SBAS-InSAR方法测量分布式散射体DS（Distributed Scatters）的形变。随着技术的发展，在这一领域取得了许多研究进展，SARscape5.7版本开始，提供了能够同时提取PS和DS测量值的新方法，即增强型永久散射体（E-PS）和增强型短基线（E-SBAS）。

E-PS和E-SBAS这两种方法都可以同时测量PS点和DS像元的位移，可以得到覆盖更多的干涉叠加处理结果，但遵循不同的处理理念。E-PS使用自适应滤波器以保留区分PS技术的细节水平。E-SBAS处理可以在PS和DS像元时间序列上得到非线性运动。这两种技术在绝对精度、管理非连续或非线性时间序列的能力以及覆盖范围方面各有特点。

本文介绍SARscape5.7中E-PS方法处理流程，所使用的例子数据集如下：

  • Sentinel1A SLC数据，IW模式，VV极化方式，覆盖区为临海地区。
  • SRTM30米DEM。
  • 数据在SARscape软件中进行了导入处理，导入时使用了工程区地理坐标系的矢量文件进行裁剪处理。

说明：本教程假设您对SARscape中的PS-InSAR有一定的处理经验。

在数据处理之前，打开/SARscape/Preferences/Preferences specific面板，选择Load Preferences->Sentinel TOPSAR(IW-EW)，设置哨兵数据处理系统参数。选择/SARscape/Interferometric Stacking/PS & E-PS/PS Overview，浏览PS和E-PS数据处理说明。

1、E-PS软件处理流程

在SARscape软件中，PS处理包含五个步骤，E-PS是在PS处理的第二步干涉工作流之后增加一步滤波操作，软件处理流程图如下：

图1：E-PS处理流程

  • 连接图生成

对输入的SLC数据集，软件根据设定的基线阈值自动选择一个主影像，其他影像跟该景组成干涉像对。

  • 干涉工作流

    这一步是根据像对的连接关系，对每一对像对进行干涉工作流处理。所有图像都要配准和重采样到参考文件，这涉及到4个因素（至少2个）的过采样，避免在大基线的情况下干涉条纹过密而混淆。不同于标准的InSAR处理，因为PS方法是找点目标，所以不进行光谱移动和普通的多普勒带宽滤波。会生成每个从影像和参考主影像的干涉图。

  • E-PS：自适应滤波

进行滤波处理：通过自适应滤波器生成滤波后的干涉图。

  • PS第一次反演

    第一次生成的形变结果，包括速率和高程改正值，这个结果没有去除大气影响的相位。

  • PS第二次反演

    使用大气校正步骤，生成最终的形变结果。

  • 地理编码

所有PS相关结果（形变速率、高度残差、形变序列、KML、矢量文件等）投影到地图系统中。

2、E-PS数据处理

在SARscape5.7软件中，E-PS处理工具如下：

图2：E-PS处理工具

一、连接图生成

（1）在Toolbox中，双击/SARscape/Interferometric Stacking/PS & E-PS/1 - Connection Graph工具。打开在PS Connection Graph工具面板。

（2）输入数据（Input Files），单击，选择多时相的_slc_list数据。

（3）可选文件（Optional File），此处不选择，系统自动选择主影像。

（4）数据输出（Output Files）面板，单击按钮，在输出路径下，输入文件根名称。

注：系统会自动创建：根名称为_ps_processing的文件目录。

（5）单击Exec按钮开始执行。

这一步生成的结果有：

结果内容

说明

auxiliary.sml

工程文件，文件记录相关的处理信息，包括处理步骤和生成的文件记录等

_pwr以及相应的.hdr和.sml文件

主影像的强度图

CG_report

像对的基线等情况

Reference_selection.txt

报告文件，其中显示输入选择参数，如所选主参数、基线和理论精度值以及阈值等

CG_baseline.txt

CG_position.txt

存储图形可视化的文件（CG_baseline.txt和CG_position.txt）。

注：可使用/SARscape/Interferometric Stacking/Stacking Tools/Plot Viewer工具，输入auxiliary.sml文件，显示连接图。

二、PS干涉工作流

PS干涉处理是对每一对干涉像对进行配准、和干涉处理， PS干涉测量过程是在斜距几何下进行的，所以无法对其进行直观查看。

（1）在Toolbox中，双击/SARscape/Interferometric Stacking/PS & E-PS/2 - Interferometric Process工具，打开PS Interferometric Process面板。

（2）数据输入（Input Files）：点击，选择sml文件；

（3）坐标系设置（DEM/Cartographic System）：输入参考DEM文件；

（4）参数设置（Parameters）： Principal Parameters参数选项，按照默认参数。

  • Generate Dint Multilooked for Quick View：默认设置为False。设置为True以生成多视的差分干涉图的快视图。
  • Az Looks for Quick View：1
  • Rg Looks for Quick View：4
  • Rebuild All：默认为False。设置为True可从头开始整个过程。
  • Atmosphere External Sensors：用于进行大气校正的传感器数据
    • Not SELECTED：不进行大气校正
    • GACOS：使用GACOS数据对干涉结果进行大气校正
    • MERIS：使用MERIS数据对干涉结果进行大气校正，用于ENVISAT ASAR的数据
  • 使用DEM进行配准：默认为True。配准时使用DEM

（5）单击Exec按钮进行处理。

图3：2 - Interferometric Process工具

三、E-PS：自适应滤波

此工具的目的是通过在干涉相位上执行自适应滤波，来增强分析非城市地区标准永久散射体（PS）。

这种方法参考了Ferretti（2011）和Fornaro（2015）的工作成果。可以独立执行PS和分布式散射体（DS）的联合处理。

为实现这个目的，需要进行两个主要步骤。首先，对统计上相似的像素集合进行识别。这步骤基于配准和定标的堆栈SAR数据，对其振幅信息进行Kolmogorov-Smirnov (KS) and Anderson-Darling (AD)检验。

然后，对经过统计检验识别的所有DS点，利用相似像素集合来估计协方差矩阵。采用高斯性假设下的最大似然估计(MLE)或利用协方差矩阵的最大主成分对对应于DS的SLC相位进行最优加权。相干性高于特定阈值的DS与PS联合处理，用于时序位移的最终估计。

（1）在Toolbox中，/SARscape/Interferometric Stacking/PS & E-PS/ E-PS /Adaptative Filtering工具，打开Adaptative Filtering面板。

（2）数据输入（Input Files）：点击，选择sml文件；

（3）参数设置（Parameters）： Principal Parameters参数选项，主要有以下参数。

  • 从相似像素点识别DS点的方法（Shp Map Method）：提供两种方法KS和AD。
  • Win Az Size (px)、Win Rg Size (px)：自适应滤波在方位向和距离向上的窗口大小。
  • 自适应滤波方法（Adaptive Filtering Method）：提供两种方法MLE和PCA。
  • DS threshold：设置区分PS和DS像素的阈值。
  • 是否使用滤波批处理模式（Activate Filtering batch Mode）：True。
  • 批处理数量（Batch size (acq number) ）：默认10。
  • 是否删除原始的去平后干涉图（Delete orig dint files）：False。
  • 重新构建（Rebuild All）：False。

注：上述参数的说明参考帮助文件，暂时缺乏设置经验。

（4）按照默认参数，单击Exec按钮进行处理。

图4：Adaptative Filtering参数面板

这一步生成的结果有：

结果内容

说明

_rsp_slc

滤波和配准后的SLC数据

cc_ds

DS像素多时相相干系数

cc_ds_avg

多时相相干系数平均值，设置使用滤波批处理模式时才会生成。

KS/AD_shp_neighbors

存放在Map文件夹内。统计相似像素的数量，每个值的范围是1到方位向和距离向上的窗口大小。

work_interferogram_stacking文件夹中还包括滤波后的干涉图_fil_DS_dint。

四、PS第一次反演

第一次反演会自动选择一个或者多个PS参考点，即永久散射体，之后重点分析这些可靠的单一目标（一个像素）的历史相位，从而估算形变速率和残余高程。

SARscape默认25 km2选择一个PS参考点，当处理区域小于25 km2时，只使用一个PS参考点处理整个区域；当处理面积较大时，每25 km2划分为一个子区域，相邻子区域之间有30%重叠，整个区域将分成若干个子区域，每个子区域单独处理，最后将每个子区域结果合并。

第一次模型反演得到位移速率和残余高程。这一步还没有去除大气中的任何相位分量。

（1）在Toolbox中，双击/SARscape/Interferometric Stacking/PS & E-PS/3 - Inversion: First Step工具，打开First Step面板。

（2）数据输入（Input Files）：点击，选择sml文件。

（3）参数设置（Parameters）：Principal Parameters参数选项：

  • 形变采样频率（Displacement Sampling (mm/year)）：设置形变速率采样频率，默认为1
  • 最小形变速率（Min Displacement Velocity (mm/year)）：预设最小形变速率，默认-25；
  • 最大形变速率（Max Displacement Velocity (mm/year)）：预设最大形变速率， 默认25；

注：以上两个参数，用户根据已经掌握的资料对此区域沉降的预估，设置一个合理的形变范围区间，如果未知可先按照默认设置，对第一次反演的位移速率进行分析，进一步设置最小和最大形变速率

  • 残余高度采样频率（Residual Height Sampling (m)）：设置残余高度采样频率，默认为2；
  • 最小残余高度（Min Residual Height (m)）：预设最小地形残余高度，默认-70米；
  • 最大残余高度（Max Residual Height (m)）：预设最大地形残余高度，默认-70米；

注：以上两个参数取决于研究区是否有高层建筑，最高达到多少米，用来去除高程残差，设置个高层建筑高度的正负范围。

  • 选择一个参考点的最大面积（Area For Single Reference Point (sqkm)）：25
  • 子区重叠区比例（Overlap for SubArea (%)）：30
  • 候选点个数（Number of Candidates）：300，默认一个子区域选择300个候选点，之后选出一个作为参考GCP点
  • 重新构建（Rebuild All）：False

（4）单击Exec按钮执行处理。

图5：First Step参数面板

五、PS第二次反演

这个步骤是对第一次线性模型反演的结果估算并去除大气相位成分，得到最终形变速率。

（1）在Toolbox中，双击/SARscape/Interferometric Stacking/PS & E-PS/4 - Inversion: Second Step工具，打开Second Step面板。

（2）数据输入（Input Files），点击，输入sml文件

（3）参数设置（Parameters）：Principal Parameters参数选项。

  • Atmosphere High Pass Size：输入以天为单位的窗口大小。应用于时间分布相关滤波。默认365。大气变化的空间分布比例。它通过一个正方形窗口实现：大的窗口适合校正大范围的变化，而小窗口适合校正孤立的局部变化。窗口越小，过滤效果越强。
  • Atmosphere Low Pass Size：输入以米为单位的窗口大小。应用于空间分布相关滤波。默认1200。大气变化的时间分布比例。它通过一个瞬时窗口实现：大的窗口适合校正频率慢的大气变化，小的窗口适合校正频率快的大气变化。窗口越大，过滤效果越强。
  • 重新构建（Rebuild All）：False。

（4）单击Exec按钮执行处理。

图6：Inversion: Second Step中参数设置

六、地理编码

PS产品地理编码处理后，默认输出为PS矢量点文件，也可以选择输出栅格格式的形变速率结果。

（1）在Toolbox中，点击/SARscape/Interferometric Stacking/PS & E-PS/5 - Geocoding工具，打开PS Geocoding界面。

（2）数据输入（Input Files）：点击，输入sml文件。

（3）可选文件（Optional File）：Refinement GCP File。

优化GCP点文件。为了获得更可靠的形变结果，可以输入一个或者多个GCP（可以来自GPS或者其他途径），这些GCP是用来优化形变趋势评估。如果选择一个GCP，修正系数是一个平均速率常数偏移量；如果是多个GCP，修正系数是一个从所有GCP通过最佳拟合计算得到的。此处不选择

（4）坐标系设置（DEM/Cartographic System）：点击，输入参考dem文件。

（5）参数设置（Parameters）：Principal Parameters参数选项。

  • 产品相干性阈值（Product Coherence Threshold）：设置相干系数阈值，相干系数小于这个阈值的像素不会保留在PS结果中。这里默认设置为75
  • 生成KML文件（Generate KML）：False。激活该项将会生成KML文件。
  • KML文件最大刻度（Upper Limit KML Scaling）：生成KML文件参数激活后，该参数可设置，KML文件中预期最大刻度，默认为10。
  • KML文件最小刻度（Lower Limit KML Scaling）：生成KML文件参数激活后，该参数可设置，KML文件中预期最小刻度，默认为-10。
  • 生成矢量文件（Make Geocoded Shape）：True，生成PS点矢量文件
  • 生成栅格文件（Make Geocoded Raster）：False，不生成PS栅格结果
  • 生成斜距坐标的矢量文件（Make Slant Shape）：False，不生成斜距坐标的矢量文件
  • 重新处理（Rebuild All）：False。
  • 优化数据集（Refinement Stacking）：可使用输入的refinement GCP点对以下数据进行优化：
    • 速率形变优化（Stack Velocity Disp Refinement）
    • 残余高程优化（Stack Residual Height Refinement）
    • 所有产品优化（Stack All Products Refinement）
  • 形变优化的多项式阶数（Refinement Displacement Poly Degree）：多项式的阶数，用于估计位移斜坡，在重新去地形处理时，将从输入位移日期逐日删除。如果这个值高于输入的地面控制点的数量，它将自动减少。默认值3表示在距离和方位角方向上的位移斜坡加上一个恒定的位移偏移将被修正。如果只需要位移偏移校正，则将多项式度设置为1。
  • 垂直形变（Vertical Displacement）：当激活并设置为True时，形变和速率将投影到垂直方向上。
  • 坡度形变（Slope Displacement）：当激活并设置为True时，形变和速率将投影到最大坡度方向上。
  • 自定义方向形变（Displacement Custom Direction）：当激活并设置为True时，形变和速率将投影到自定方向上，需要设置方位角（Azimuth Angle：以度为单位，从北顺时针方向）和倾斜角（Inclination Angle：以度为单位，从水平面开始计算）。
  • X分辨率（X Dimension）(m)：输出像元X方向大小，以米为单位，默认15米。
  • Y分辨率（Y Dimension）(m)：输出像元Y方向大小，以米为单位，默认15米。

注：投影坐标系为经纬度坐标时，当设置像元值大于0.2，软件内部会粗糙的自动转换为度的单位，当设置像元值小于0.2，软件自动识别以度为单位。

  • Other Parameters参数项，常做如下设置：
  • Max Points in Shape：每个矢量文件中最大值点数量，默认为100000。如果PS点非常多，会生成多个分块矢量文件。当处理区域较大时，可设置为1000000。

注：由于shapefile格式的限制，它不能超过2GB，因此最多允许约7000万个点要素，所以该值也不宜设置过大。当点太多时，PS点存放在多个shapefile文件。

  • Geocoding参数项，常进行如下设置：
  • Dummy Removal：True，去除无效值

（6）单击Exec按钮执行处理。

图7：地理编码数据输入及参数选项

PS处理结果中，得到的矢量文件中，包含每个PS点的信息，存储在geocoding/E-ps_EPS_75_0.shp文件中,这里文件名75标识相干性阈值为0.75。

七、结果浏览

PS默认处理结果为Shp点矢量文件，可以使用GIS软件进行浏览和制图，本例使用SARscape中提供的工具进行简单的浏览分析。

（1）打开地理编码后的SAR平均强度图\geocoding\mean_geo，作为底图。

（2）打开E-PS点矢量文件shp，叠加在SAR平均强度图上显示。

（3）在Toolbox中，点击/SARscape/General Tools/Time Series Analyzer/Vector工具。打开TS Vector Analyzer面板，设置一个显示区间，本例中设置-30和30，单击Color Apply，对PS点进行彩色显示。

（4）同时叠加PS结果。

注：如果有多个PS矢量图层，在所要渲染和分析的图层上点击右键设置“Set as active layer”，使该图层处于激活状态，在TS Vector Analyzer面板上对该图层进行操作。

如下图所示，E-PS相比PS在相同地表覆盖区域得到的结果点密度有了显著提高，其中海上人工建筑也得到了相应的结果，相比PS结果，这些海上人工建筑上大部分区域没有得到结果。

图8：E-PS（左图）和PS结果（右图）

posted @ 2024-06-05 10:11  ENVI-IDL技术殿堂  阅读(1484)  评论(0)  编辑  收藏  举报