用哨兵1A干涉数据对DInSAR测量台湾地震引起的地表形变信息

2016年2月6日早晨3时57分，台湾高雄发生6.7级地震，震源深度15千米。震中位于高雄市美浓区，震源深度16.7公里，云林县草岭的最大震度达6级，高雄、屏东、台南、嘉义等地的最大震度均达5级。本节以哨兵1A slc作为数据源，使用DInSAR的方法，对2016年2月6日台湾南部的地震事件进行干涉测量。

本文所使用软件版本：ENVI5.3.1+SARscape5.2.1

数据情况：

  • Sentinel1A IW slc像对
  • 极化方式：VV极化
  • 入射角：39.14
  • 成像时间：2016年1月9日（震前），2016年2月24日（震后）
  • 成像范围：台湾南部地区
  • 地面分辨率：20米
  • 辅助数据：90米分辨率的SRTM DEM数据
  • 数据来源：https://scihub.esa.int/

1、数据准备

第一步：设置系统参数

/SARscape/Preferences，设置Load Preferences为SENTINEL_TOPSAR，这套参数是专门针对哨兵1A数据的TOPSAR (IW)模式做InSAR处理时的系统参数。

图 设置SARscape系统参数

如果计算机支持GPU，在General Parameters面板，设置Opencl Platform Name的类型为相应的设备，如NVIDIA CUDA显卡。

设置Opencl的类型

为了便于后面操作的数据读取和输出，设置ENVI的系统参数中的输入输出路径。ENVI->File->Preferences，在Directories中设置默认的输入路径和输出路径。

设置ENVI系统参数

第二步：数据导入

（1）数据文件输入

网站下载的数据解压，打开数据导入工具/SARscape/Import Data/SAR Spaceborne/SENTINEL 1，在Input File List分别输入两景哨兵1A数据的元数据文件manifest.safe。

Optional Input Orbit File List，轨道文件，该文件是可选文件，在这里没有轨道文件，就不输入。

图 数据输入面板

（2）参数设置

切换到Parameters面板，主要参数就是对输出数据命名设置，推荐选择Rename the File Using Parameters：True，可以对输出的数据自动按照数据类型进行命名。如果要输出镶嵌后的SLC数据文件，可以在Other Parameters设置Generate SLC mosaic为True，否则默认是输出各个条带的slc数据。数据类型软件自动识别，不用做设置。

图 数据输入面板参数设置

（3）输出设置

之前设置过默认的输出路径，这里直接按照默认即可，如果要改输出路径，在数据输出路径上点击右键，选择Change Output Directries。

图 数据输出面板的输出设置

设置好参数后，点击Exec执行，完成后弹出对话框，点击End。

输出的数据文件包括：整景图像的强度图数据（_pwr）、slc索引文件（.slc_list）带有地理坐标的外边框矢量文件（.shp）、GoogleEarth文件（.kml）。打开矢量文件，ENVI主界面view->Reference Map Link，可以自动链接Arcgis Online的在线底图，查看数据的地理范围。

图 数据的地理范围

第三步：参考DEM数据的下载

打开工具/SARscape/General Tools/Digital Elevation Model Extraction/SRTM-3 Version 4，在Input File输入两个时相的slc_list文件，DEM/Cartographic System设置GEO GLOBLE WGS84，其他参数按照默认。

注：1）这一步也可以在InSAR界面下进行。

2）如果网速许可，最好下载SRTM Version4版本的DEM数据。SRTM Version2版本的数据有空洞，会对最终的结果有所影响。

图 根据数据范围下载DEM

自动下载的SRTM DEM数据

2、基线估算

打开工具/SARscape/Interferometry/Interferometric Tools/Baseline Estimation，输入主从影像的_slc_list文件，点击Exec执行，输出基线估算的结果。

基线估算面板

基线估算结果：

Normal Baseline (m) = -64.213    Critical Baseline min - max(m) = [-6474.081] - [6474.081]

Range Shift (pixels) = 13.012

Azimuth Shift (pixels) = 4.669

Slant Range Distance (m) = 879019.561

Absolute Time Baseline (Days) = 36

Doppler Centroid diff. (Hz) = 17.918    Critical min-max (Hz) = [-486.486] - [486.486]

2 PI Ambiguity height (InSAR) (m) = 242.615

2 PI Ambiguity displacement (DInSAR) (m) = 0.028

1 Pixel Shift Ambiguity height (Stereo Radargrammetry) (m) = 20379.668

1 Pixel Shift Ambiguity displacement (Amplitude Tracking) (m) = 2.330

Master Incidence Angle = 39.723    Absolute Incidence Angle difference = 0.004

Pair potentially suited for Interferometry, check the precision plot

基线估算的结果显示，这两景数据的空间基线为64.213米，远小于临界基线6474.081米，时间基线36天，做DInSAR的一个相位变化周期代表的地形变化为0.028米。

3、DInSAR工作流

打开/SARscape/Interferometry/DInSAR Displacement Workflow工具。

（1）文件输入

  • Input File面板，输入20160109的VV极化方式的_slc_list作为主影像，20160214的VV极化方式的_slc_list作为从影像；
  • DEM/Cartographic System面板，输入参考DEM文件；
  • Parameters面板，Grid Size：20

注：数据有两种极化方式：VH和VV，选择同极化方式的做差分干涉处理。

图 数据输入面板

输入的数据文件设置好之后，点击Next。弹出自动计算视数的面板，算出来的视数为5：1，点击确定。

图 自动计算的视数

（2）干涉图生成

  • 主要参数
    • 距离向视数（Range Looks）：5。自动添加。
    • 方位向视数（Azimuth Looks）：1。自动添加。
    • 制图分辨率（Grid Size for Suggested Looks）：20。根据之前的设置自动添加。
    • 配准时使用DEM（Coregistration With DEM）：Ture。
  • 全局参数（Global）

生成TIFF数据（Make TIFF）：Ture，生成TIFF格式的中间结果，如果需要使用中间结果，如写文章时候作为插图，可以设置为True，其他步骤类似。

图 干涉图生成参数设置面板

处理完成之后，ENVI视窗中自动加载了去平后的干涉图，以及主从影像的强度图，这一步生成的数据结果存放在ENVI默认的临时文件路径下，默认路径为：C:\Users\<计算机用户名>\AppData\Local\Temp，自动生成了一个名为SARsTmpDirXXXXXXXXX的文件夹，这一步生成的结果有：

  • INTERF_out_dint：去平后的干涉图
  • INTERF_out_int：干涉图
  • INTERF_out_master_pwr：主影像强度图
  • INTERF_out_slave_pwr：从影像强度图
  • INTERF_out_par：文本文件保存的配准时的偏移量数据，
  • INTERF_out_sint：合成的相位
  • INTERF_out_srdem：斜距几何下的参考DEM
  • INTERF_out_par_orbit_off：估算轨道偏移时用到的点矢量数据，包括在方位向和距离向 的点的位置坐标、测量到的偏移量、计算出的线性偏移量。
  • INTERF_out_par_winCC_off：从强度数据的相差上估算交叉相关偏移量的点矢量数据。包含每个点上交叉相关的值（CC），范围是0-1。
  • INTERF_out_par_winCoh_off：从相位数据的相差上估算相干性的点矢量数据，包含信噪比（SNR）和相干性，相干性值的范围是0-1

图 去平后的干涉图_dint

（3）滤波和相干性计算（Adaptive Filter and Coherenace Generation）

  • 主要参数（Principal Parameters）：
    • 滤波方法（Filtering Method）：Goldstein

提供了三种滤波方法：

  • Adaptive

这种方法适用于高分辨率的数据（如TerraSAR-X或COSMO-SkyMed）

  • Boxcar

使用局部干涉条纹的频率来优化滤波器，该方法尽可能的保留了微小的干涉条纹。

  • Goldstein

    这种滤波方法的滤波器是可变的，提高了干涉条纹的清晰度、减少了由空间基线或时间基线引起的失相干的噪声。这种方法是最常用的方法。

图 滤波和相干性生成参数设置面板

点击Next按钮，进行干涉图滤波和相干性生成处理，处理完成之后，自动加载滤波后的干涉图_fint和相干性系数图_cc。这一步处理之后生成的结果有：

  • INTERF_out _fint：滤波后的干涉图
  • INTERF_out _cc：相干性系数图

图 滤波后的干涉图和相干系数图

（4）相位解缠（Phase Unwrapping）：相位的变化是以2π为周期的，所以只要相位变化超过了2π，相位就会重新开始和循环。相位解缠是对去平和滤波后的相位进行解缠处理，使之与线性变化的地形信息对应，解决2π模糊的问题。

  • 主要参数（Principal Parameters）
    • 解缠方法（Unwrapping Method Type）：Minimum Cost Flow
    • 解缠分解等级（Unwrapping Decomposition Level）：1
    • 解缠最小相干性阈值（Unwrapping Coherence Threshold）：0.2。相干系数大于该阈值的像元进行解缠。

图 相位解缠面板参数设置

说明：

1）解缠方法（Unwrapping Method Type）提供了三种：

  • 区域增长法(Region Growing)：选这种方法，则不要设置过高的相干性阈值（0.15-0.2是比较好的）以便留下足够的自由增长空间，相位突变部分在解缠后的图像上以解缠孤岛存在，这种方法降低了由相位突变引起的误差，
  • 最小费用流(Minimum Cost Flow)：默认的解缠方法，当有大面积的低相干或是其他限制增长的因素而使解缠困难时，最小费用流算法可以取得比区域增长法更好的结果。这种方法采用正方形的格网，考虑了图像上所有的像元，对相干性小于阈值的像元做了掩膜处理。

    Delaunay MCF：和最小费用流法的不同在于，这种方法不是考虑了图像上所有的像元，而是仅考虑了相干性大于阈值的部分，而且不是用正方形的格网而是用了德罗尼三角形格网。只有对相干性高的部分进行解缠，不受低相干像元的影响。对于有大量相干性低的地物存在的时候，如影像上存在大量水体、浓密植被等，推荐使用该方法，最小化相位突变的影响。

    2）分解等级（Unwrapping Decomposition Level）：该分解是为了用迭代的方法对数据做多视和疏采样：干涉图以一个较低的分辨率被解缠然后被重采样成原来的分辨率。使用了分解可减少解缠错误，提高处理效率。用户可以指定迭代的次数，每个迭代相当于3次采样过疏。这里可以填的数有：-1、0、1、2、3。其中，-1和0代表不执行分解，用原始的像素采样，当形变很大或是地形很陡峭的情况下（多相位/高密度分布），分解可引起交迭效应，建议设置为-1或0；1代表最小的分解等级。3：最大的分解等级。建议这个次数不要超过3。通常情况设置原始的像素采样（如-1）或者最小的分解等级（1）。

    点击Next按钮，进行干涉图滤波和相干性生成处理，处理完成之后，自动加载滤波后的相位解缠结果图_upha。这一步处理之后生成的结果有：

    • INTERF_out _upha：相位解缠结果

    图 相位解缠结果

    （5）控制点选择（GCP Selection）：输入用于轨道精炼的控制点文件，可以用已有的文件也可在此选择控制点。

    在Refinement GCP File（Mandatory）项中，点击按钮，自动打开流程化的控制点选择工具，软件自动输入了相应的文件。

    图 选择控制点流程化工具

    在控制点生成面板上，点击Next，打开控制点选择工具，鼠标变为选点状态，单击鼠标左键就可以选点。在图像周边远离形变区域的地方、相位好的地方选择若干控制点。

    图 控制点选择

    单击Cartographic System按钮，查看控制点的参考坐标系统，该坐标系是从参考DEM上自动读取的。

    图 设置控制点的坐标系

    单击Export按钮，查看控制点的存放路径和文件名。生成的控制点文件为INTERF_out_upha_gcp.xml。

    注：如果默认输出的路径有中文字符的文件夹或者文件名，需要改到只有英文字符的文件路径中。

    图 控制点文件输出

    在控制点生成面板上点击Finish，生成了控制点文件，并自动添加到InSAR流程化处理面板的Refinement GCP File（Mandatory）项。

    图 生成控制点文件

    在面板上点击Next按钮，进入下一步。

    （6）轨道精炼和重去平（Refinement and Reflattening）进行轨道精炼和相位偏移的计算，消除可能的斜坡相位，对卫星轨道和相位偏移进行纠正。这一步对解缠后的相位是否能正确转化为形变值很关键。

    • 主要参数（Principal Paremeters）
      • 轨道精炼方法（Refinement Method）：Polynomial Refinement，软件先使用多项式优化的方法来计算，该算法的健壮性更好，即使是在基线小的情况下也可以使用。
      • 轨道精炼的多项式次数（Refinement Res Phase Poly Degree）：3。在重去平的过程中用到的估算相位斜坡的多项式次数，若输入的控制点个数较少，次数会自动降低。默认的3表示在距离向和方位向加上一个恒定的相位偏移的相位斜坡会被修正，如果仅需要相位偏移校正，这个次数可以设置为1。
      • 配准时是否考虑到地形因素（Coregistration With DEM）：Ture。默认参数，不需要用户设置。

    图 轨道精炼和重去平参数设置面板

    点击Next按钮，进行轨道精炼和重去平处理，处理完成之后，将优化的结果显示在Refinement Results的面板，内容如下：

    ESTIMATE A RESIDUAL RAMP

     

    Points selected by the user = 64

    Valid points found = 64

    Extra constrains = 2

    Polynomial Degree choose = 3

    Polynomial Type : = k0 + k1*rg + k2*az

    Polynomial Coefficients (radians) :

            k0 = 8.9544091183

            k1 = -0.0045068910

            k2 = -0.0018570690

    Root Mean Square error (m)= 64.1891189212

    Mean difference after Remove Residual refinement (rad) = -0.0431859477

    Standard Deviation after Remove Residual refinement (rad) = 1.9514668017

    这一步得到的结果有：

    • INTERF_out_reflat_fint：重去平后的干涉图
    • INTERF_out_reflat_sint：重去平后的合成相位图
    • INTERF_out_reflat_upha：重去平后的解缠结果
    • INTERF_out_reflat_srdem：重去平后的斜距DEM
    • INTERF_out_reflat.txt：轨道精炼所用的轨道修正参数
    • INTERF_out_reflat_refinement.shp：轨道精炼用到的有效的控制点文件，斜距坐标系。
    • INTERF_out_reflat_refinement_geo.shp：轨道精炼用到的有效控制点文件，地理坐标系。

    （7）相位转形变以及地理编码（Phase to Displacement Conversion and Geocoding）：将经过绝对校准和解缠的相位，结合合成相位，转换为形变数据以及地理编码到制图坐标系统，默认得到的是LOS方向的形变。

    • 主要参数（Principal Parameters）
      • 相干性阈值（Product Coherence Threshold）：0.2。相干性大于该值的相位转为形变值。
      • 垂直方向的形变（Vertical Displacement）：False。不计算垂直方向上的形变。
      • 斜坡形变（Slope Displacement）：False。不计算斜坡方向上的形变，在监测如滑坡形变的时候可激活该选项。
      • 用户自定义方向的形变（Displacement Custom Direction）：False。
      • 方位角（Azimuth Angle）：0。自定义方向的时候设置该角度。
      • 入射角（Inclination Angle）：0。自定义方向的时候设置该角度。
      • X方向上的水平分辨率（X Dimension（m））：20。制图分辨率，X方向，按之前的设置默认为20。
      • Y方向上的水平分辨率（Y Dimension（m））：20。制图分辨率，y方向，按之前的设置默认为20。

    图 相位转形变主要参数

    • 地理编码参数（Geocoding）
      • 去除图像外的无用值（Dummy Removal）：True。对图像外的区域做掩膜处理

    图 地理编码参数

    点击Next按钮，进行相位转形变和地理编码处理。地理编码的坐标系是以参考DEM的坐标系为准。这一步得到的结果有：

    • _slc_out_disp_dem：重采样到制图输出分辨率上的参考DEM数据，经过地理编码，范围和输出的SAR产品一致。
    • _slc_out_disp_cc_geo：地理编码的相干性系数图
    • _slc_out_disp_precision：数据质量的估算结果图，代表形变的精度。
    • _slc_out_disp_ALOS：视线方位角，正值是正北的顺时针方向，负值是正北的逆时针方向。
    • _slc_out_disp_ILOS：视线入射角，视线和水平面垂线的夹角。
    • _slc_out_disp：传感器观测方向的形变，即LOS方向上的形变。

    图 相位转形变的结果

    （8）结果输出（Output）：

    结果默认输出在ENVI的默认输出路径下，文件名命名为sentinel1_69_20160109_100020274_IW_SIW1_A_VV_slc_list_out_disp。若想保留中间结果便于查看，将Delete Temporary Files不勾选。

    点击Finish，输出结果。结束DInSAR Diaplcement处理的工作流。生产的形变数据自动进行密度分割配色展示。

    图 形变结果自动配色显示

    注：Back和Next按钮，可切换至中间某一步查看参数或调整参数重新处理。

    4、结果分析

    结果显示，本次地震发生的位置，导致了一个区域明显的抬升。如果对整景数据做DInSAR处理，右边的大部分没有形变的区域，都参与到处理中，容易引入一些误差（能看出右侧的稳定区域，结果也有一些抬升的情况）。在这种情况下，可以在_fint的数据上，对地震形变区域手动绘制一个矢量区域，用SARscape的裁剪工具对原始的SLC数据进行裁剪，然后对子区域做DInSAR处理，能得到更精确的结果。

    ENVI中手动绘制矢量范围：打开_reflat_fint数据，ENVI主菜单File->New-> Vector Layer，选择_fint数据，绘制矢量区域，如下图所示：

    图 绘制子区域的矢量范围

    在矢量图层上点击右键，选择Save As，保存为subarea.shp文件。

    点击/SARscape/General Tools/Sample Selections/Sample Selection SAR Geometry Data工具，打开Sample Selection SAR Geometry Data面板。

    • 数据输入项（Input Files）：点击Browse按钮，选择需要裁剪的两景数据的_VV_slc_list文件：

    图 子区裁剪的数据输入面板

    • 可选文件项（Optional File）：
      • 矢量文件（Vector File）：可以输入地理坐标下的子区范围的矢量文件裁剪，用的是矢量文件的外接多边形。点击，选择上一步绘制生成的subarea.shp矢量文件。
      • DEM文件（DEM File）：当用坐标范围选择子区时，需要输入DEM数据，这里范围是斜距范围，所以不用设置DEM。
      • 参考文件（Input Reference File）：如果对哨兵1A数据_slc_list文件进行裁剪的话，需要输入多视后的_list_pwr强度数据作为参考文件。点击，选择一景数据的_slc_list_pwr数据作为参考文件。

    图 输入裁剪矢量范围及参考的强度数据

    • 参数项（Parameters）：
      • 配准（Make Coregistration）：False。不对输入的待裁剪的数据进行配准处理
      • 地理范围（Geographical Region）：False。选择True，代表输入的范围是地理坐标下的范围。选择False，代表输入的范围是斜距坐标下的范围
      • West/First Column：最西边/第一列的坐标，可输入地理坐标（经纬度格式如29.30）。也可输入行列号。
      • North/First Row：最北边/第一行的坐标。可输入地理坐标（经纬度格式如29.30）。也可输入行列号。
      • East/Last Column：最东边/最后一列的坐标。可输入地理坐标（经纬度格式如29.30）。也可输入行列号。
      • South/Last Row：最南边/最后一行的坐标。可输入地理坐标（经纬度格式如29.30）。也可输入行列号。
      • 使用最大和最小坐标（Use Min and Max Coordinates）：False。激活该项目的话只用到输入的矢量文件的角点坐标进行裁剪。

    图 子区裁剪的参数设置面板

    • 输出文件项（Output Files）：默认的文件名中添加了_cut。

    图 裁剪结果输出面板

    图 裁剪结果强度图

    然后对裁剪后的slc数据，进行DInSAR工作流处理（参数和之前相同），得到地震区域的形变结果。在子区域DInSAR处理中，在轨道精炼一步，选择的控制点如下图所示：

    子区域DInSAR处理中轨道精炼的GCP

    得到的轨道优化的结果报表：

     

    ESTIMATE A RESIDUAL RAMP

     

    Points selected by the user = 27

    Valid points found = 27

    Extra constrains = 2

    Polynomial Degree choose = 3

    Polynomial Type : = k0 + k1*rg + k2*az

    Polynomial Coefficients (radians) :

            k0 = 8.6332973041

            k1 = -0.0025818939

            k2 = -0.0017582238

    Root Mean Square error (m)= 16.5580109079

    Mean difference after Remove Residual refinement (rad) = -0.1456119011

    Standard Deviation after Remove Residual refinement (rad) = 0.7900878373

    图 裁剪后的SLC做DInSAR处理得到的结果

posted @ 2022-05-19 17:14  ENVI-IDL技术殿堂  阅读(3665)  评论(6)  编辑  收藏  举报