哨兵1A数据DInSAR测量云南大理州漾濞县地表形变

    据中国地震台网中心测定，2021年5月21日21时48分在云南大理州漾濞县（北纬25.67度，东经99.87度）发生6.4级地震，震源深度8千米。

    本文以哨兵1A作为数据源，使用DInSAR的方法，对本次地震进行干涉测量处理。数据情况如下表所示：

卫星	Sentinel1A
获取模式	IW
数据级别	SLC
极化方式	VV
轨道方向	降轨
入射角	42.932
地面分辨率	5*20米
成像时间	2021年5月10日（震前）、2021年5月22日（震后）
参考DEM数据	ALOS World 3D DEM

本文使用的软件为SARscape5.6。

哨兵1A数据下载网址：https://scihub.copernicus.eu/dhus/#/home

1、数据导入

    处理数据之前，设置哨兵1数据的系统参数，打开/SARscape/Preferences/Preferences specific面板，选择Load Preferences->Sentinel TOPSAR(IW-EW)，在General Parameters面板中，设置Cartographic Grid Sze(m)为25米（本文使用25米作为制图分辨率）。

打开数据导入工具/SARscape/Import Data/SAR Spaceborne/Single Sensor/SENTINEL-1。在打开的面板中，

• 数据输入面板（Input Files）
• 输入两个时相的safe文件
• 可选文件面板（Optional Files）
• 输入所需处理地震范围的地理范围矢量文件
• 参数设置面板（Parameters）：主要参数（Principal Parameters）
• 生成强度数据（Make power QL）：True，数据导入时生成强度数据，以便查看；
• 极化方式（Polarization）：Copolarization Only，只输出同极化数据（DInSAR处理首选同极化数据）；
• 同轨道镶嵌（Make mosaic same track）：True，若输入同轨道相邻数据，自动进行镶嵌；
• 不进行子区域裁剪（Skip Sample Selection）：False，使用输入的地理范围进行子区域裁剪；（这一参数为6版本新增参数，使用该参数在数据导入时直接进行子区域裁剪）
• 对数据重命名（Rename the File Using Parameters）：True。软件会根据输入数据的文件名进行输出文件的自动命名。
• 数据输出面板（Output Files）：输出文件（Output file list）：自动读取ENVI默认的数据输出目录以及输入面板中的数据文件名。

注：1.如果要修改输出的路径，在右边单击文件夹图标选择输出文件夹目录。

2、如果要修改输出的路径，在输出文件名右键选择Edit菜单。

图 数据输入面板参数设置

输出的数据文件包括：整景图像的强度图数据（_pwr）、slc索引文件（.slc_list）带有地理坐标的外边框矢量文件（.shp）。

注：本次数据处理时，精密轨道文件（数据获取时间之后21天发布）还未发布，故本次数据导入没有使用精密轨道文件，轨道误差仅在轨道精炼中使用GCP拟合多项式去除。一般情况下，做InSAR处理，导入哨兵数据建议使用精密轨道文件。

2、基线估算

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

图 基线估算面板

基线估算结果：

Normal Baseline (m) = -48.988        Critical Baseline min - max(m) = [-6368.518] - [6368.518]

Range Shift (pixels)    = 0.128

Azimuth Shift (pixels)  = 0.392

Slant Range Distance (m)  = 916515.195

Absolute Time Baseline (Days) = 12

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

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

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

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

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

Master Incidence Angle = 42.932   Absolute Incidence Angle difference = 0.003

Pair potentially suited for Interferometry, check the precision plot

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

3、DInSAR处理

打开/SARscape/Interferometry/DInSAR Displacement Workflow工具，使用DInSAR流程化工具进行差分干涉处理，步骤包括：干数据输入、涉生成（包括复数据配准、干涉相位计算、去已知地形）、相干性计算、干涉图滤波、相位解缠、轨道精炼和重去平、相位转形变。

在进行DInSAR工作流之前，需要准备该区域的参考DEM，本文使用的是ALOS World 3D DEM。可使用/SARscape/General Tools/DEM Extraction/ALOS World 3D 30M工具自动下载。

图 DEM自动下载工具

一、     数据输入

DInSAR Displacement Workflow工具左侧第一步Input File面板，输入如下数据：

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

图 输入干涉像对的slc数据

图 输入参考DEM数据

图 设置制图分辨率

输入的数据文件设置好之后，点击Next。由于之前已经做了数据导入，故Import Generic SAR Date这步略过，直接点击Next到Interferogram Generation界面。

二、     干涉图生成

这一步对干涉像对进行配准，并计算干涉相位，生成干涉图，使用参考DEM对已知地形相位去除，参数说明如下：

• 主要参数（Principal Parameters）
  • 距离向视数（Range Looks）：7。
  • 方位向视数（Azimuth Looks）：2（自动计算得到的是1，由于SARscape自动计算视数的原则是取整而非四舍五入，在此手动调整为2）。
  • 制图分辨率（Grid Size for Suggested Looks）：25。
  • 配准时使用DEM（Coregistration With DEM）：Ture。
• 全局参数（Global）

生成TIFF数据（Generate Quick Look）：Ture，生成结果图的TIFF格式快视图，方便查看中间结果，其他步骤类似。

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

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

结果命名	内容
INTERF_out_master_pwr/ INTERF_out_slave_pwr	主影像/辅影像的多视强度图
INTERF_out_par	干涉像对配准时生成的偏移量数据
INTERF_out_par_orbit_off	配准时估算轨道偏移时用到的点矢量数据，包括在方位向和距离向 的点的位置坐标、测量到的偏移量、计算出的线性偏移量。
INTERF_out_par_winCC_off	配准时从强度数据的相差上估算交叉相关偏移量的点矢量数据。包含每个点上交叉相关的值（CC），范围是0-1。
INTERF_out_par_winCoh_off	配准时从相位数据的相差上估算相干性的点矢量数据，包含信噪比（SNR）和相干性，相干性值的范围是0-1
INTERF_out_srdem	地距转为斜距的DEM数据
INTERF_out_sint	地形合成相位
INTERF_out_int	干涉相位
INTERF_out_dint	去已知地形的差分干涉相位

图 去已知地形后差分干涉图_dint

三、     干涉图滤波和相干性计算

对上一步生成的差分干涉图进行滤波，对干涉图的噪声进行一定程度的抑制，计算干涉像对的相干系数。干涉图滤波有4种方法（Adaptive、Boxcar、Goldstein、Adaptive Non Local InSAR）可选，本次使用了Goldstein滤波方法，这种滤波方法的滤波器是可变的，提高了干涉条纹的清晰度、减少了由空间基线或时间基线引起的失相干的噪声。这种方法是最常用的方法。

• 主要参数（Principal Parameters）：
• 滤波方法（Filtering Method）：Goldstein
• 滤波参数（Filtering）：
• Goldstein Win Size：128
• Goldstein Min Alpha：1
• Goldstein Max Alpha：3

由于地震区域相干性比较低，默认滤波参数得到的干涉条纹还是存在较多噪声，本次处理增加了滤波窗口和Alpha最大最小值，效果是增强了滤波强度。

注：滤波参数的说明可以查看软件自带的help。

• 全局参数（Global）

生成TIFF数据（Generate Quick Look）：Ture，生成结果图的TIFF格式快视图，方便查看中间结果，其他步骤类似。

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

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

结果命名	内容
INTERF_out_fint	自适应滤波后的差分干涉图
INTERF_out_cc	相干系数图

图 滤波后的差分干涉图

图 相干系数图

    相干系数图是一副0-1的灰度图，值越大表示相干性越高，地表发生变化的区域相干性低，从相干系数图上可以明显的看到地震区域地表了发生了变化，该区域相干系数低。

四、     相位解缠

干涉图中的相位，从-π到π呈周期性变化，当真实相位值大于π时，相位会重新从-π开始，以2π为周期循环。相位解缠是将相位由差分相位恢复为真实值的过程。

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

图 相位解缠面板参数设置

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

结果命名	内容
INTERF_out_upha	相位解缠结果

图 相位解缠结果

五、     轨道精炼

使用一组代表稳定位置的GPC点，进行多项式拟合，对轨道误差进行去除，优化相位解缠的结果。

1、GCP选择

在GCP Slection面板，点击 按钮，自动打开流程化的控制点选择工具，软件自动在数据输入位置输入了相应的文件（Input file为_fint，reference file为_upha，DEM file为参考DEM文件）。点击Next，在控制点选择面板，鼠标变为选点的十字丝状态，单击鼠标左键就可以选点。在远离形变条纹、相干性高、解缠结果平滑的位置，选择若干控制点。

• GCPs：GCP列表，每个GCP点的坐标在该面板上显示
• Cartographic System：默认，与参考DEM坐标系一致
• Export：输出GCP文件，注意路径不要有中文字符。

图 选择控制点流程化工具

在控制点生成面板上，点击Finish，生成的控制点文件.xml，自动添加到面板中。点击Next按钮。

图 GCP文件生成

2、轨道精炼

    使用多项式模型，结合GCP点的相位，解算多项式系数，拟合多项式，对干涉图和解缠后的相位进行轨道精炼和重去平，得到优化的相位结果。

• 轨道精炼方法（Refinement Method）：Polynomial Refinement，软件默认使用多项式优化的方法来计算，该算法健壮性好。
• 多项式系数（Refinement Res Phase Poly Degree）：3，默认值3表示在距离和方位角方向上的相位残差与恒定的相位偏移将被校正。
• 使用DEM配准（Coregistration With DEM）：True

图 轨道精炼面板

点击Next，运行得到轨道精炼后的干涉图_reflat_fint和解缠相位图_reflat_upha，这一步处理之后生成的结果有：

结果命名	内容
INTERF_out_reflat_fint	优化后的滤波差分干涉图
INTERF_out_reflat_sint	优化后的合成相位图
INTERF_out_reflat_upha	优化后的解缠相位
INTERF_out_reflat_refinement.shp	轨道精炼用到的有效的控制点文件，斜距坐标系
INTERF_out_reflat_refinement_geo.shp	轨道精炼用到的有效控制点文件，地理坐标系
INTERF_out_reflat.txt	轨道精炼结果报告

同时输出了使用这组GCP点拟合的轨道精炼多项式及误差信息：

ESTIMATE A RESIDUAL RAMP

Points selected by the user = 16

Valid points found = 16

Extra constrains = 2

Polynomial Degree choose = 3

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

Polynomial Coefficients (radians) :

                   k0 = 12.2587151588

                   k1 = 0.0008717432

                   k2 = -0.0019063896

Root Mean Square error (m)= 68.2671075452

Mean difference after Remove Residual refinement (rad)  = 0.2129152330

Standard Deviation after Remove Residual refinement (rad)  = 1.4389082745

可以看出，本次使用了16个GCP点，全部位于有效像元，使用这16个GCP点，拟合的轨道精炼的多项式为：12.2587151588 + 0.0008717432*rg -0.0019063896*az，均方根误差RMSE为68.267米，去除残差后的平均差值为0.2129弧度，去除残差后的标准差为1.4389弧度。

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

图 轨道精炼和重去平前后的差分干涉图（左：原始结果，右：轨道精炼后）

结果命名	内容
INTERF_out_reflat_fint	优化后的滤波差分干涉图
INTERF_out_reflat_sint	优化后的合成相位图
INTERF_out_reflat_upha	优化后的解缠相位
INTERF_out_reflat_refinement.shp	轨道精炼用到的有效的控制点文件，斜距坐标系
INTERF_out_reflat_refinement_geo.shp	轨道精炼用到的有效控制点文件，地理坐标系
INTERF_out_reflat.txt	轨道精炼结果报告

六、     相位转形变与地理编码

将经过轨道精炼和重去平的解缠相位，转换为形变，并地理编码到地理坐标系。

• 主要参数（Principal Parameters）
  • 相干性阈值（Product Coherence Threshold）：3。相干性大于该值的相位转为形变值。
  • 垂直方向的形变（Vertical Displacement）：False。不计算垂直方向上的形变。
  • 斜坡形变（Slope Displacement）：False。不计算斜坡方向上的形变。
  • 用户自定义方向的形变（Displacement Custom Direction）：False。
  • 方位角（Azimuth Angle）：0。自定义方向的时候设置该角度。
  • 入射角（Inclination Angle）：0。自定义方向的时候设置该角度。
  • X方向上的水平分辨率（X Dimension（m）：25。X方向制图分辨率。
  • Y方向上的水平分辨率（Y Dimension（m）：25。Y方向制图分辨率。

图 相位转形变主要参数

点击Next按钮，进行相位转形变和地理编码处理。地理编码的坐标系与参考DEM一致。默认得到的是LOS方向的形变结果。这一步得到的结果有：

结果命名	内容
_dispwf_disp_precision	数据质量的估算结果图，代表相位转形变的精度
_dispwf_disp_ILOS	视线入射角，视线和水平面垂线的夹角
_dispwf_disp_dem	重采样到制图输出分辨率上的参考DEM数据，范围与形变结果一致
_dispwf_disp_cc_geo	地理编码的相干性系数图
_dispwf_disp_ALOS	视线方位角，正值是正北的顺时针方向，负值是正北的逆时针方向
_dispwf_disp	形变结果，单位为m，默认是LOS方向上的形变

点击Next。

七、     结果输出：

   在Output面板，选择是否删除中间结果，点击Finish，自动对形变结果进行分级设色，输出密度分割图。

注：1、软件此处有bug，不论勾选是否删除，都会删除中间结果，如需要保留，建议在此步骤之前手动备份中间结果文件夹。

2、在DInSAR流程化工具界面，点击Back和Next按钮，可切换至中间某一步查看参数或调整参数重新处理。

图 相位转形变的结果

形变结果是一个单波段的灰度图像，像元值就是形变量，单位为米，正值代表朝传感器方向运动（抬升），负值代表远离传感器方向移动（沉降）。

DInSAR结果显示，本次地震引起地表形变量级约为6cm，地震发生的位置，存在两个方向相反的形变区，西南侧沉降，东北侧抬升。将滤波差分干涉图进行地理编码，对应形变结果，在干涉图上也可以看出此特点。

注：1、干涉图的地理编码可使用/SARscape/Basic/Intensity Processing/Geocoding/Geocoding and Radiometric Calibration工具

2、在应急情况下，为了初步快速了解地震灾害的影响，直接对照干涉条纹判断形变区和估算形变量，也是常用的方法（一个相位变化周期对应二分之一的波长的形变）。

图 左：像元形变量的查询 右：地理坐标系的干涉条纹

4、结果渲染与表达

形变结果是一个单波段的灰度图像，像元值就是形变量，正值代表朝传感器方向运动（抬升），负值代表远离传感器方向移动（沉降）。

对形变结果进行彩色渲染有两种方式：颜色分级显示、色带显示。

• 颜色分级显示

使用ENVI的密度分割功能实现，在形变结果图上点击右键，点击New Raster Color Size，在密度分割面板，可根据制图需要设置各级DN值范围与颜色，

图 密度分割分级显示

• 色带显示

默认形变结果是灰度黑白显示，可在形变结果图上点击右键，选择Change Color Table，点击More，点击 ，可选择多种配色方案进行彩色渲染。 按钮可以对色带进行颜色调换。

图 多种彩色渲染方案可选

小技巧：使用彩色渲染形变结果，由于形变量较小，且在图像中属于弱信息（大部分像元是稳定的），可通过调整形变结果的直方图，凸显形变信息。点击 打开直方图，默认为2%线性拉伸，可手动调整直方图到最大最小位置，稳定位置为处在色度带中部的相近颜色，沉降和抬升在色带的末端，可凸显形变信息。

图 默认2%线性拉伸直方图彩色渲染效果

图 调整直方图全范围彩色渲染效果

• 叠加底图

可叠加雷达数据强度图或光学影像作为底图，形变图设置一定透明度显示。本次处理由于形变区域位于山区，可叠加山体阴影图作为底图，操作方法如下：

使用DEM数据生成山体阴影：打开ENVI工具箱/Terrain/Create Hill Shade Image，选择DInSAR工作流生成的DEM文件_dispwf_disp_dem（也可选择初始参考DEM），在山体阴影生成工具界面，点击Compute Elevation and Azimuth，输入该地区的经纬度（ENVI界面左下角自动显示鼠标所在位置的经纬度），输入时间（GMT时间），生成所设置时间和位置的山体阴影效果图，

图 使用DEM数据生成山体阴影图

形变结果图放置在山体阴影图上方，在ENVI工具栏中，使用 调整形变结果图的透明度，便可看到形变结果图叠加在山体阴影图上的效果。

图 形变结果渲染叠加山体阴影效果图

也可将地理编码后的差分干涉图叠加到山体阴影图上（干涉条纹的图例来自SARscape的Help，保存为图片，使用ENVI注记加入），效果如下：

图 差分干涉图叠加山体阴影效果图

ENVI工具栏Annotations下拉菜单，可以根据需要添加制图要素，如色度带、线段、符号、箭头标识等。

注：ENVI制图表达可参考相关博文。

图 地震形变结果制图表达

说明：

本次处理制图分辨率设置了25米，哨兵1数据制图分辨率默认为15米，分别尝试了15米、20米、25米，在形变导致失相干较为严重、植被覆盖区域相干性低等情况，低的制图分辨率，有利于得到较好的形变条纹。下面是不同制图分辨率下，形变条纹的对比结果：

图 不同制图分辨率得到的差分干涉图

本文所用的干涉像对原始数据：

链接：https://pan.baidu.com/s/1jjOyRXDDoI6FRs6kY6IPug 

提取码：envi 

同学们可自行下载练习，也可直接在欧空局官网下载。