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

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-SBAS方法处理流程，所使用的例子数据集如下：

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

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

在数据处理之前，打开/SARscape/Preferences/Preferences specific面板，选择Load Preferences->Sentinel TOPSAR(IW-EW)，设置哨兵数据处理系统参数。

注：查看SARscape/Interferometric Stacking/SBAS & E-SBAS/E-SBAS/E-SBAS Overview，浏览E-SBAS数据处理说明。

E-SBAS处理流程

在进行E-SBAS之前，需要先完成SBAS处理流程，再将SBAS处理结果文件中的辅助文件作为输入文件，进行E-SBAS处理流程。软件处理流程如下：

图1：E-SBAS处理技术流程

E-SBAS数据处理

    在SARscape5.7软件中，E-SBAS所需的处理工具如下：

图2：E-SBAS处理工具

第一步：SBAS处理

首先对数据集进行SBAS处理，此处不再赘述SBAS处理流程，可参考教程：《SARscape雷达图像处理入门与提高2023版》

第二步：E-SBAS处理

1、连接图

（1）在Toolbox中，双击/SARscape/Interferometric Stacking/SBAS & E-SBAS/1 - Connection Graph，打开E-SBAS Connection Graph面板；

  • Input Files：输入SBAS地理编码处理之后，*_SBAS_processing路径下的sml文件
  • Output Files：输入根目录名称，软件运行之后会自动添加_ESBAS_processing后缀

图3：E-SBAS Connection Graph参数面板

这一步运行完之后，自动生成E-SBAS工程路径，同时生成连接图，如下图所示。

图4：E-SBAS连接图

注：在E-SBAS处理完之前，slc数据集和第一步做完的_SBAS_processing文件夹不能挪动位置，E-SBAS处理过程会用到这些数据。

2、干涉处理

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

（2）数据输入（Input Files）：单击，选择*_ESBAS_processing路径下的auxiliary.sml文件。

（3）可选文件（Optional Files）：有以下两个可选的输入文件。

  • GeoMetry GCP File：几何控制点文件；
  • Optional Water Vapour File List：如果Parameters面板参数设置里使用GACOS数据对干涉结果进行大气校正，则这里需要输入每个时相的水汽文件，该文件是从GACOS网站免费下载的。

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

参数设置（Parameters）：Principal Parameters（主要参数）。

  • 是否生成多视快视图（Generate Dint Multilooked for Quick View）：设置是否生成干涉图快视图。
  • 快视图距离向视数（Az Looks for Quick View）：1。
  • 快视图方位向视数（Rg Looks for Quick View）：4。
  • 重新计算（Rebuild All）：False。
  • 大气传感器（Atmosphere External Sensors）：用于进行大气校正的传感器数据
    • Not SELECTED：不进行大气校正
    • GACOS：使用GACOS数据对干涉结果进行大气校正
    • MERIS：使用MERIS数据对干涉结果进行大气校正，用于ENVISAT ASAR的数据
  • 使用DEM进行配准（Coregistration With DEM）：True。

（5）单击Exec按钮，运行E-SBAS干涉工作流。

图5：E-SBAS Interferometric Process面板

该步骤运行完之后，会生成链接图对应像对的干涉结果，存放在_ESBAS_processing\work\work_interferogram_stacking路径下。如果选择输出快视图，结果存放位置：_ESBAS_processing\interferogram_stacking\tiff_multilooked_dint。同时生成相应的_meta索引文件，方便打开浏览多个结果文件。

3、反演

反演步骤利用SBAS第二次反演产品初始化E-SBAS反演，估算大气相位。通过模型反演，去除来自SBAS的低通分量和大气相位部分信号，最终拟合高通位移和速率，得到逐日位移。

（1）在Toolbox中，双击/SARscape/Interferometric Stacking/SBAS & E-SBAS/E-SBAS/3 – Inversion，打开E-SBAS Inversion面板。

（2）数据输入（Input Files）：单击，选择*_ESBAS_processing路径下的auxiliary.sml文件。

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

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

（4）单击Exec按钮，运行E-SBAS反演。

图6：E-SBAS Inversion面板

反演步骤生成了以下产品。

产品

说明

Height

大气校正的DEM校正值（米）

Height_lp

来自SBAS处理和重采样为E-SBAS干涉图相同分辨率的DEM校正值（米）

Velocity

大气校正之后得到的平均形变速率（毫米/年）

precision_vel

估算形变速率时相应的平均精度（毫米/年）

precision_height

估算残余高度（residual height）时相应的平均精度（米）

cc

每个数据对的相干系数，它显示有多少形变趋势符合所选模型

slant_atm_meta

在斜距几何上每期大气相关成分

slant_dint_reflat_meta

大气校正后，在斜距几何上每期差分干涉图

slant_disp

大气校正后，在斜距几何上的每期形变量

slant_disp_full_meta

PS和DS的大气校正后逐日形变结果

 

生成两个GCP文件：Ref_GCP.shp和Ref_GCP_geo.shp。

4、地理编码

（1）在Toolbox中，双击/SARscape/Interferometric Stacking/SBAS & E-SBAS/E-SBAS/4 – Geocoding，打开E-SBAS Geocoding面板。

（2）数据输入（Input Files）：单击，选择*_ESBAS_processing路径下的auxiliary.sml文件。

（3）可选文件（Optional Files）：

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

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

（5）参数设置（Parameters）：

  • Principal Parameters参数选项：
  • 产品相干性阈值（Product Coherence Threshold）：设置相干系数阈值，相干系数小于这个阈值的像素不会保留在结果中。这里默认设置为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，生成矢量结果文件。
  • 生成栅格文件（Make Geocoded Raster）：False，不生成栅格结果文件。
  • 生成斜距坐标的矢量文件（Make Slant Shape）：False，不生成斜距坐标的矢量文件
  • 重新处理（Rebuild All）：False。
  • 大地水准面类型（Geoid Type）：提供两种类型EGM96和EGM2008。
  • 优化半径（Refinement Radius (m)）：GCP在斜距上解缠相位上有效像素半径。
  • 形变优化的多项式阶数（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。如果结果点非常多，会生成多个分块矢量文件。当处理区域较大时，可设置大一些，如1000000。

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

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

（6）单击Exec按钮执行E-SBAS地理编码处理。

图7：E-SBAS Geocoding面板

"geocoded" 文件夹里包括了地理编码生成的结果，主要包括：

结果内容

说明

*_PS_DS__XX_Y.shp

带有相关文件（.shx、.dbf和Google Earth.kml）的点矢量结果。如果点的数量大于地理编码参数中定义的“Shape Max Nr of points”，则会创建多个矢量文件。命名方式如：根名称_ PS_XX_Y.shp其中XX是“Product Coherence Threshold”（0.XX）的小数部分，Y是多个shapefile情况下的增量数字。

mean_geo

SAR平均强度图像和附带文件(.sml, .hdr)

 

shapefile属性表中的字段说明如下：

-Velocity: 每个PS点的平均形变速率mm/y。

-Coherence: 多时相的相干性，是干涉相位在时序上最佳拟合的指标。

-MuSigma: 是一个质量指标，是均值/标准差得到的值，均值是在时间上的均值。

-Scatterer：标识这个点的类型：PS、DS。

-HPrecisio:它表示了模型估计的高程测量的平均精度，单位是米。它的计算考虑了数据的空间基线和多时相相干性。

-VPrecisio: 它表示了形变速率测量的平均精度，单位是毫米/年。计算时考虑了数据的时间基线和多时相相干性。

-Range:点距离向的像元坐标。

-Azimuth:点方位向的像元坐标。

-SubArea ID：子区ID号。

-Lon/ Lat: 点的地理坐标。

-xpos/zpos:基于DEM参考系统的地理坐标。

-zpos: 基于DEM椭球体的高程。

-Z: 基于大地水准面的高程。

-ALOS: 方位向上的视线入射角。

-ILOS: 垂直方向上的视线入射角。

-Hcorrectio: 高程相对于输入的参考DEM的校正量，单位为米。

-D_date: 该时相相对于第一个时相的形变量，单位为mm。

浏览结果

浏览输出的Shp点矢量文件，可以使用GIS软件进行浏览和制图，本例使用SARscape中提供的工具进行简单的浏览分析。

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

(2)打开E-SBAS点矢量文件*_PS_DS_75_0.shp，叠加在SAR平均强度图上显示。

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

(4)同时叠加SBAS点矢量结果。

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

如下图所示，E-SBAS相比SBAS在相同地表覆盖区域得到的结果点密度有了显著提高，其中海上人工建筑也得到了更多的结果。

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

posted @ 2024-06-12 08:55  ENVI-IDL技术殿堂  阅读(1762)  评论(2)  编辑  收藏  举报