介绍
2017年九寨沟地震是一场发生于2017年8月8日的地震,其震中位于中华人民共和国四川省阿坝州九寨沟县漳扎镇(北纬33.20度,东经103.82度),震级为Ms7.0,震源深度为12千米。维基百科详细介绍:(https://zh.wikipedia.org/wiki/2017%E5%B9%B4%E4%B9%9D%E5%AF%A8%E6%B2%9F%E5%9C%B0%E9%9C%87),地震给当地人民造成了巨大的生命和财产损失。愿逝者安息!
哨兵1号(Sentinel-1)卫星是欧洲航天局哥白尼计划(GMES)中的地球观测卫星,由两颗卫星组成,载有C波段合成孔径雷达,可提供连续图像(白天、夜晚和各种天气)。最难能可贵的是Sentinel-1开放数据,结合InSAR技术,我们能够及时看到地震的巨大破坏力。震区的差分干涉图如下,最后导出部分也有结果展示,在滤波后的差分干涉图中能够看到大约8个干涉条纹,因此粗略估计的视线向地面形变量为22.4cm。
Sentinel-1数据准备
首先需要收集数据,数据下载网址(https://scihub.copernicus.eu/dhus/)。登录之后,先在网络上查明地震发生的位置。
在Sentinel-1数据网站上,将图层切换到卫星地图,找到对应的位置。
在Sentinel-1数据网站上框选出感兴趣的区域,然后依次选择Sentintel-1,IW,SLC模式。
点击搜索那个放大镜,进行数据搜索。
搜索出来的数据比较多,选择InSAR处理要用到的数据。
SAR数据选择对干涉来说还是很重要的,覆盖范围、升降轨数据、日期、基线各种因素都会影响最后的结果。利用InSAR来研究地震的形变场主要是使用DInSAR技术,需要使用两景SAR数据和数字高程模型。所以选择震前震后两景数据,日期分别是20170811和20170730(注意SAR影像的时间是UTC时间,我们这里是东八区,所以实际卫星采集的数据的时刻还要在这个基础上再加八个小时),选择好数据之后点击下载。
数据下载完成之后打开Preview检查一下,确保数据覆盖地震发生区域。
Sentinel-1数据处理
Skysense软件启动界面
1. 新建工程
从软件的【文件】菜单打开【新建工程】,弹出新建工程对话框,设置好工程名和路径,选择使用的卫星数据。
现有主流的星载数据基本都支持,选择“Sentinel IW” (Interferometric Wides wath)数据。
可以看到软件既支持单个subswath的处理也支持三个subswath一起处理。相比于最常用的条带数据,TOPS数据的Burst和Subswath的重合区域的拼接是一个难题,这里使用三个subswath联合处理。
2. 数据导入
新建好工程之后,可以在【文件】菜单栏选择【导入SLC数据】。
也可以在“工程树栏“里右键单击然后【导入SLC数据】,软件可以自动下载轨道数据和DEM数据,节省了操作人员很多时间。
导入的数据会提示没有轨道数据,点击【下载轨道数据】就可以自动下载。
由于精密轨道数据要在数据生产出来之后20天才能拿到,所以这里默认下载了粗轨(预测轨道)数据。
还可以自动下载DEM数据,可以选择ESD(Enhanced Spectral Diversity)提高配准精度。Sentinel TOPS数据配准要求达到千分之三像素,主要用基于地形的配准方法和ESD方法才能达到。导入的使用到了DEM,进行了基于地形的配准配准操作, ESD算法辅助配准可以进一步提高配准的精度。
点击导入窗口下方的【开始】就可以进行数据导入的操作。
3. 干涉处理
在Skysense中进行干涉处理要选取感兴趣的区域(Area Of Interest, AOI)。
在工程树栏的空白处右键单击,选择【新建AOI】,然后给新建的AOI命名。
Skysense中侧重于时序分析,只能选择PS和SBAS方法。DInSAR并不是软件的主要处理方法,但是PS操作使用两景图像就是DInSAR处理,所以这里选择【CuPS方法】。
3.1 配准
在工程树上的【AOI节点】上右键单击,选择【配准】,弹出配准窗口。
在配准窗口的【选择主副影像】选项卡中选择分析的数据,需要将全部数据从候选数据区选中移到分析数据中,虽然有点奇怪,但是这么设计应该是有原因的。
单击右下角【查看基线】,Sentinel-1数据的轨道管半径很小(3Sigma 150m),这两景图像的时间间隔为12天,垂直基线不到40米,时空失相干都应该不大,期待有一个好的处理结果。
在配准窗口的【选取区域】选项卡进行AOI的裁剪,裁剪出想要的区域,可以选择全景,也可以单独截出自己感兴趣的区域。这里选择全景,图像尺寸68465*13120。
3.2 生成干涉图
在【配准影像】上右键单击,选择【生成干涉图】操作。
用【全选】按钮将所有数据选中,然后选择多视数。Sentinel-1数据分辨率为距离向5m方位向20m,所以这里使用了距离向为4方位向为1的多视数。
生成完干涉图之后,可以使用鼠标滚轮进行放大缩小操作,从干涉图中可以看到这两景图像的相干性还是非常好的,有非常明亮清晰的干涉条纹。
3.3 去平地
在【干涉图】上右键选择【去平地相位】操作,弹出去平地窗口。
同样地,点击【全选】,然后点【开始】。
去平地之后数据的地形条纹很明显。
3.4 模拟地形相位
由于去平地之后的数据存在明显的地形相位,所以要通过使用DEM进行地形相位的模拟。
在【去平地干涉图】节点上右键单击,选择【模拟地形相位】。参数设置如下。
DEM可以自己下载或者自动下载,因为刚才导入数据的时候已经下载了DEM,所以这里只要选中之前下载好的DEM就行。SRTM分辨率为90m,4*1多视之后的SAR图像分辨率为20m, (90/20)向上取整就是5,方位距离方向分辨率相同,所以都是5。
由于图像很大,这一步的时间比较久,由于Skysensen软件使用了并行计算,前面的操作步骤处理速度都很快,到这步等待的时间略长。
3.5 生成差分干涉图
得到模拟的地形相位之后生成差分干涉图(这是我后来补的图,所以出现了后续结果),在【去平地干涉图】节点上右键单击选择【生成差分干涉图】。
从差分干涉图中可以看到地震产生的干涉条纹。
3.6 滤波
粗略来看,差分干涉图里的条纹数可以反映地震在震中产生的地面形变,但是为了提高信噪比压制噪声,就要对差分干涉图进行一下滤波,在【差分干涉图】-【未滤波结果】节点上选择【滤波】-【GoldStein滤波】。
滤波完的结果条纹更加清楚和明亮。
沿着图中直线的方向,从最中间的紫色圆环向外,从紫色到红色为一个条纹,一共能够观察到8个条纹,粗略估计视线向地表形变为22.4cm。
3.7 地理编码和导出
滤波完之后应该对干涉图数据进行相位解缠,然后转化为形变量,由于我对于滤波的参数还不是十分熟悉,处理结果中噪声也比较大,而且DInSAR产生的形变量其实里边有很多大气、地形残余之类的误差没有去除,所以这次的DInSAR处理只看条纹。
这里就把滤波后的结果直接【地理编码】,地理编码选择DEM,选择利用强度和轨道进行配准,地理编码结果输出的分辨率设置为和SAR图像一样的20m。
地理编码的结果如图。但这样还是只能看看,要想知道地震的干涉条纹在哪里,还有最后一步,导出到Google Earth。
导出到Google Earth之后的结果概览。
通过调整KML文件中干涉条纹图片的透明度,将干涉的结果和Google Earth的光学图像进行叠加和结果。
选择合适的透明度,放大到局部。
到这里就完成了差分干涉的操作。整体看下来,在Skysense软件中,看不出TOPS数据和普通条带数据的处理流程上有太大的区别,软件封装的比较好。唯一需要注意的就是准备足够的磁盘容量,硬盘的剩余空间要足够。
当然,严肃来说,我只是抱着学习的态度使用了Skysense软件,相比于很多老师教授和专业团队的处理结果还存在很大的差距。要想获得更好的处理结果还要对InSAR的原理和参数设置有更深的了解和尝试,但是Skysense很好的一点就是每一步骤处理完都能直接看到图像,在以后的学习过程中还要多多调整参数,观察比较结果,促进学习,希望有一天自己也能处理出漂亮、真实、准确的结果,也希望InSAR能够引起越来越多人的关注,大家一起学习使用这门技术,为地震后的救灾工作贡献一份力量。
最后说一下处理时间,整个处理下来还是很快的,在个人电脑上,大概用了大半天的时间。我还做了一个8*2多视的Subswath 3数据,从导入数据到最后的出KML结果,差不多两个多小时,对于灾害后短时间获取地表形变还是非常有帮助的。