工具:ENVI5.6 ; SarScape5.6 ;
数据详情:
![](https://i-blog.csdnimg.cn/blog_migrate/8e22c9f88f4951fdda1bb72394855555.png)
哨兵一号数据:
![](https://i-blog.csdnimg.cn/blog_migrate/94a149fcf90f28cad15d5d041973ee80.png)
参数设置
![](https://i-blog.csdnimg.cn/blog_migrate/c522ddc4deab0633214b7b4ff6f2963c.png)
导入哨兵固定参数:
![](https://i-blog.csdnimg.cn/blog_migrate/4689628ead240eada078a7c3615d81a0.png)
数据导入
(1)safe文件:
![](https://i-blog.csdnimg.cn/blog_migrate/32517d5e7eb83be2d7a0298484b34352.png)
(2)矢量文件:
![](https://i-blog.csdnimg.cn/blog_migrate/5f85b62276213c2f4ce8303ed4c894a8.png)
在QGIS中打开矢量文件,主图层为光学卫星影像esri satellite,研究区域如下:
![](https://i-blog.csdnimg.cn/blog_migrate/b1bf8b2bea3b1ffef4152a8f61042271.png)
(3)参数
由于slc数据的极化方式是VV,选择copolarzation ,若是VH或HV极化,则选crosspolarzation
(4)设置输出路径、运行
(5)结果
得到_pwr强度图tif文件:
![](https://i-blog.csdnimg.cn/blog_migrate/01bec5b944b183222c96ef8d6769b31d.png)
![](https://i-blog.csdnimg.cn/blog_migrate/fde6a43c6579219bb1d66071da5837e2.png)
得到的SAR影像和光学卫星影像左右颠倒,这是因为轨道方向为 降轨。(升轨为上下颠倒)
(注:本次未采用精轨文件,精轨文件可在欧空局官网:ESA 进行下载)
基线估计
![](https://i-blog.csdnimg.cn/blog_migrate/0311955ac4acbb8feecf0616bcfd1198.png)
基线估算的结果显示,这两景数据的空间基线为50.069米,远小于临界基线6368.418米,时间基线12天,做DInSAR一个相位变化周期代表的地形变化为0.028米。
DInSAR形变工作流
(1)输入数据
![](https://i-blog.csdnimg.cn/blog_migrate/573993e43fb363293f1ac17acff9651d.png)
制图分辨率选择了25m(哨兵1数据制图分辨率默认为15米),在形变导致失相干较为严重、植被覆盖区域相干性低等情况,低的制图分辨率,有利于得到较好的形变条纹。
(2)导入基本SAR数据
第二部分已导入过数据,这一步直接默认,next
(3)生成干涉图
![](https://i-blog.csdnimg.cn/blog_migrate/1e60876d3d5868ff007ad903a1bc462c.png)
运行结束后,自动加载了去已知地形的差分干涉图 _dint :
![](https://i-blog.csdnimg.cn/blog_migrate/966bb8f303702f716ae2c82057168154.png)
然后,在默认输出路径的文件夹下:
![](https://i-blog.csdnimg.cn/blog_migrate/69c361d30e39c7d6b9eb1d4aebe77f3f.png)
找出 去已知地形的差分干涉图 _dint_ql.tif (预览图):
![](https://i-blog.csdnimg.cn/blog_migrate/98106ec205082c08828dd593de7350c8.png)
(注:关于Quick Look图,又称快视图、预览图,在文件中有_ql后缀,个人将理解其为干涉图和_pwr多视强度图的叠加影像图,其配色和干涉图有所不同)
(4)自适应滤波和相干性生成
对上一步生成的差分干涉图进行滤波,对干涉图的噪声进行一定程度的抑制,计算干涉像对的相干系数。干涉图滤波有4种方法(Adaptive、Boxcar、Goldstein、Adaptive Non Local InSAR)可选,本次使用了Goldstein滤波方法,这种滤波方法的滤波器是可变的,提高了干涉条纹的清晰度、减少了由空间基线或时间基线引起的失相干的噪声,这种方法是最常用的方法。
由于地震区域相干性比较低,默认滤波参数得到的干涉条纹还是存在较多噪声,本次处理增加了滤波窗口和Alpha最大最小值,效果是增强了滤波强度。
![](https://i-blog.csdnimg.cn/blog_migrate/eb4cd62a8b39b34ab29ed6bcae19cd40.png)
得到结果:
![](https://i-blog.csdnimg.cn/blog_migrate/94ee1211ce9ed810cc18b9578999795d.png)
![](https://i-blog.csdnimg.cn/blog_migrate/6ed071dc3b4b957c0da5301d1c51b2de.png)
自适应滤波后的差分干涉图_ql.tif(预览图)
![](https://i-blog.csdnimg.cn/blog_migrate/044a3b08b84ca28a595e93300568c160.png)
相干系数图_ql.tif(预览图)
相干系数图是一副0-1的灰度图(图中 小图为相干系数分布直方图,其横坐标是相干系数,纵坐标是像素点个数),值越大颜色越白,表示相干性越高,地表发生变化的区域相干性低。
将相干系数分布直方图进行拉伸,只显示相干性低于0.3的像素点,可明显看到如下地震区域相干系数较低,表明其地表了发生了变化。
![](https://i-blog.csdnimg.cn/blog_migrate/853b65368803e86c4104d04ee900b551.png)
(5)相位解缠
干涉图中的相位,从-π到π呈周期性变化,当真实相位值大于π时,相位会重新从-π开始,以2π为周期循环。相位解缠是将相位由差分相位恢复为真实值的过程。
默认参数,点击next,得到相位解缠结果:
![](https://i-blog.csdnimg.cn/blog_migrate/d197749d78e48ef227b4acef30f1dd31.png)
相位解缠结果(预览图:INTERF_out_upha_ql.tif)
(注:每次要生成图片时,在Global设置中将Generate Quick Look设为True,以此来生成TIF)
(6)轨道精炼和重去平
使用一组代表稳定位置的GPC点,进行多项式拟合,对轨道误差进行去除,优化相位解缠的结果。
![](https://i-blog.csdnimg.cn/blog_migrate/583ee9da1bb9185bad6cebe1806c74ab.png)
在远离形变条纹、相干性高、解缠结果平滑的位置,选择若干控制点:
![](https://i-blog.csdnimg.cn/blog_migrate/05d655b4d83378776816167ab27cf33c.png)
cc
单击finish,下一步
![](https://i-blog.csdnimg.cn/blog_migrate/5a6c8d7718bdaa8111e94463570be92c.png)
使用多项式模型,结合GCP点的相位,解算多项式系数,拟合多项式,对干涉图和解缠后的相位进行轨道精炼和重去平:
![](https://i-blog.csdnimg.cn/blog_migrate/05831f824d04d53d2d31276e9e30d7c8.png)
运行得到轨道精炼后的干涉图_reflat_fint和解缠相位图_reflat_upha:
![](https://i-blog.csdnimg.cn/blog_migrate/85d53b122e12d0a2ae4b401536407246.png)
_reflat_fint
![](https://i-blog.csdnimg.cn/blog_migrate/cd120752a90f2ccc0aecd572288fc4ea.png)
_reflat_upha
轨道精炼报告如下:
ESTIMATE A RESIDUAL RAMP
Points selected by the user = 11
Valid points found = 11
Extra constrains = 2
Polynomial Degree choose = 3
Polynomial Type : = k0 + k1*rg + k2*az
Polynomial Coefficients (radians) :
k0 = -1.2727723347
k1 = -0.0018806371
k2 = 0.0004458157
Root Mean Square error (m)= 148.8535827091
Mean difference after Remove Residual refinement (rad) = -1.5275151390
Standard Deviation after Remove Residual refinement (rad) = 3.2334229996
可以看出,本次使用了11个GCP点,全部位于有效像元,使用这11个GCP点,拟合的轨道精炼的多项式为:-1.2727723347 - 0.0018806371 * rg + 0.0004458157 * az ,均方根误差RMSE为148.85米,去除残差后的平均差值为-1.5275弧度,去除残差后的标准差为3.233弧度。
(7)相位转形变与地理编码
将经过轨道精炼和重去平的解缠相位,转换为形变,并地理编码到地理坐标系。
![](https://i-blog.csdnimg.cn/blog_migrate/f28b6bcdc62a46a3b8f3c0dae6dd306a.png)
得到结果如下:
![](https://i-blog.csdnimg.cn/blog_migrate/784dccf4b3b8148151e7b9f902bc7eb7.png)
(8)输出结果
得到形变结果如下,是一个单波段的灰度图像,像元值就是形变量,单位为米,正值代表朝传感器方向运动(抬升),负值代表远离传感器方向移动(沉降)
如下所示,东北侧像元值为0.0615,表明本次地震引起地表形变量级约为6cm,地震发生的位置,存在两个方向相反的形变区,西南侧沉降,东北侧抬升
![](https://i-blog.csdnimg.cn/blog_migrate/f53eab6a3692e3e53b21a6e6a620dac4.png)
东北侧
结果渲染与表达
(1)颜色分级显示
使用ENVI的密度分割功能实现,在形变结果图上点击右键,点击New Raster Color Size,在密度分割面板,可根据制图需要设置各级DN值范围与颜色
![](https://i-blog.csdnimg.cn/blog_migrate/26679a6424f5cccc3e1b907b2b1f00fd.png)
(2)叠加底图
利用山体阴影生成工具,可叠加雷达数据强度图或光学影像作为底图,形变图设置一定透明度显示。
使用DEM数据生成山体阴影:打开ENVI工具箱/Terrain/Create Hill Shade Image,选择下图中的文件,其为Dinsar工作流生成的DEM模型:_cut_dispwf_disp_dem (该DEM经过工作流处理,与原DEM相比,已被裁剪)。
点击Compute Elevation and Azimuth,输入该地区的经纬度(ENVI界面左下角自动显示鼠标所在位置的经纬度),输入时间(GMT时间)。
通过调节透明度开关,调节形变图或山体阴影的透明度,得到最终结果如下!
![](https://i-blog.csdnimg.cn/blog_migrate/f55cc5cf6254a0783a2e4c59a672c4a8.png)
形变图叠加
也可与颜色分级图叠加:
![](https://i-blog.csdnimg.cn/blog_migrate/c19b38bc5480e5b22e6de171f15dd9b4.png)
颜色分级形变图叠加
ps:中间步骤太多过不了审,私信可答疑、提供数据,欢迎交流哦😋