其他分享
首页 > 其他分享> > 哨兵1A数据DInSAR测量云南大理州漾濞县地表形变

哨兵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。在打开的面板中,

注: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面板,输入如下数据:

图 输入干涉像对的slc数据

图 输入参考DEM数据

图 设置制图分辨率

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

二、     干涉图生成

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

生成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滤波方法,这种滤波方法的滤波器是可变的,提高了干涉条纹的清晰度、减少了由空间基线或时间基线引起的失相干的噪声。这种方法是最常用的方法。

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

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

生成TIFF数据(Generate Quick Look):Ture,生成结果图的TIFF格式快视图,方便查看中间结果,其他步骤类似。

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

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

结果命名

内容

INTERF_out_fint

自适应滤波后的差分干涉图

INTERF_out_cc

相干系数图

图 滤波后的差分干涉图

图 相干系数图

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

四、     相位解缠

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

图 相位解缠面板参数设置

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

结果命名

内容

INTERF_out_upha

相位解缠结果

图 相位解缠结果

五、     轨道精炼

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

1、GCP选择

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

图 选择控制点流程化工具

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

图 GCP文件生成

2、轨道精炼

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

图 轨道精炼面板

点击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

轨道精炼结果报告

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

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

图 相位转形变主要参数

点击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