标签：1A 相位 DInSAR 漾濞县 DEM INTERF 形变 干涉 out

    据中国地震台网中心测定，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 

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

标签：1A,相位,DInSAR,漾濞县,DEM,INTERF,形变,干涉,out
来源： https://www.cnblogs.com/enviidl/p/16266590.html

本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享；
2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关；
3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关；
4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除；
5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。