宜宾地震同震形变场双轨法D-InSAR技术的研究

第44卷第6期
测绘与空间地理信息Vol.44，No.6
Jun.，2021 2021年6月GEOMATICS&SPATIAL INFORMATION TECHNOLOGY
宜宾地震同震形变场双轨法D-InSAR技术的研究
李铭
(河南理工大学测绘与国土信息工程学院，河南焦作454000)
摘要：为了研究2019年6月17日宜宾Ms6.0级地震引起的地表形变量，在欧空局下载覆盖研究范围内的Sen-tinel-1A卫星C波段雷达影像数据，使用SRTM(3") DEM，采用双轨法差分干涉测量技术进行地表形变信息的提取，对比不同相位解缠方法对研究结果精度的影响，最终得到宜宾地震在视线方向上的形变图、形变量。

结果表明：双轨法D-InSAR技术在地表形变信息提取的应用中，使用合适的解缠方法能够提高成果的精度。

关键词：D-InSAR；地表形变；相位解缠；形变场
中图分类号：P225.1文献标识码：A文章编号：1672-5867(2021)06-0214-04
Research on D-InSAR Technology of Double Track Method for Coseismic Deformation Field of Yibin Earthquake
LI Ming
(School of Surveying and Landing Information Engineering,He'nan Polytechnic University,Jiaozuo454000,China)
Abstract:In order to study the surface deformation caused by the Yibin Ms6.0earthquake on June17,2019,the C-band radar image data of sentinel-1a satellite in the study area were downloaded from ESA.The SRTM(3")DEM was used to extract the surface de­formation information by using the dual orbit differential interferon technology.The influence of different phase unwrapping methods on the accuracy of the research results was compared,and the Yibin ground deformation information was finally obtained The deformation map and deformation variable of the earthquake in the direction of line of sight.The results show that：in the application of dual orbit D -InSAR technology in surface deformation information extraction,using appropriate unwrapping method can improve the accuracy of the results.
Key words：D-InSAR；surface deformation；phase unwrapping；deformation field
0引言
2019年6月17日22时55分在四川省宜宾市长宁县发生Ms6.0级地震，震源深度16km,属于“浅源”地震。

据不完全统计，此次地震共造成长宁县、琪县、高县等区县受灾，因灾死亡13人,220人受伤，受灾人数243880人,倒塌和严重破坏房屋9532户20185间。

虽然损失不是特别严重，但迅速、准确地对震后引起的地表形变进行研究是非常重要的，可以为地震多发区提供技术参考和理论依据。

传统的地表形变监测方法，并不适合在城市、山地等进行大面积区域的形变监测。

合成孔径雷达(Synthetic Aperture Radar,SAR)由于其全天候、大面积、低成本等优势［1］,已成为大范围滑坡监测［2］、积雪深度反演⑷以及地表形变等重要监测手段。

合成孔径雷达差分干涉测量(D-InSAR)技术是在In­SAR的基础上发展而来的,其原理是对覆盖同一区域形变前后的两幅SAR影像，经过差分干涉处理获得有效的地表形变信息。

D-InSAR技术可用于监测厘米级甚至毫米级的地球表面形变［4］,在地壳运动、震源机理、火山喷发、地面沉降、冰川移动等研究领域越来越受关注，其技术应用也越来越成熟。

早在1993年Massonnet等［5］利用 ERS-1的SAR数据对1992年加利福尼亚Landers地震同震位移场数据做出了提取，并发表在《Nature》杂志上，引起了国际上的震惊。

在国内，徐小波、单新建等［6］对柴达木盆地地震群进行了同震形变场的反演;臧妻斌［7］利用 不同相位解缠方法对意大利拉奎拉6.3级地震同震形变场做了研究;杨亚莉［8］利用D-InSAR技术分析了神东矿区地表形变规律;王钰⑼采用D-InSAR技术结合数值模拟对张双楼开采沉陷规律进行了研究；王燕燕［10］对九寨
收稿日期：2020-01-13
作者简介:李铭(1991-)，男，河南濮阳人，测绘工程专业硕士研究生，主要研究方向为InSAR 技术
第6期李铭:宜宾地震同震形变场双轨法D-InSAR技术的研究215
沟同震形变场进行了三维反演等。

本文以欧空局发布的Sentinel-1A卫星影像为数据来源，选取震前、震后两幅SAR影像，采用双轨D-InSAR技术对宜宾地震引起的地表形变进行研究,并对不同相位解缠方法得到的形变结果进行精度分析，最终得到该研究区域在雷达视线方向上的形变图、形变量。

1研究区概况及数据选取
1.1研究区概况
2019年6月17日发生的宜宾Ms6.0级地震，主震过后,一天之内共记录到余震68次,其中包括一次5.1级、一次5.3级。

从地理位置上看，宜宾属于华蓥山地震带系统，该断裂带历史上最大地震等级为Ms6.3。

近年来，此地震带上的宜宾、长宁、荣昌等地区都发生过4—6级地震，因此，对该地区进行震后地表形变分析是必要的，可为此方面的研究提供技术支持和参考。

本文以宜宾震中位置为中心进行区域选择,地理范围如图1所示。

黑色方框为哨兵数据单景SLC影像覆盖范围，白色线框表示裁剪后研究区域地理范围。

图1研究区域地理范围图
Fig.1Geographical location map of research region
1.2数据选取
为了对宜宾地震同震形变场进行研究，本文选取时间基线较短的20190614（震前）、20190708（震后）两幅Sentinel-1A卫星影像进行差分干涉处理,SAR影像数据各参数信息见表1。

该卫星载有C波段成像系统，重访周期12d,4种采集模式（SM、IW、EM、WAVE）。

本文选取成像模式为IW—干涉宽幅模式，该模式有4种极化方式，采用中等分辨率（5mx20m），数据覆盖范围幅宽250km,通过递进的地形扫描方式（TOPSAR,Terrain Observation with Progressive Scans）获取3个子条带［11］。

本文采用双轨法差分干涉测量（D-InSAR）进行地表形变信息的提取，使用3〃即90m空间分辨率的SRTM DEM（数字高程模型）模拟地形相位，从而消除干涉相位所包含的地形相位信息，最终得到雷达视线方向（LOS）上的地表形变信息。

1.3方法原理和技术流程
实验借助Envi SARscape数据处理平台对覆盖研究区域形变前后的两幅SAR影像进行干涉处理，通过3种相位解缠方法得到研究区域雷达视线上（LOS）的地表形变信息，然后对其精度进行分析对比,使用精度最高、对形变信息误差影响最小的解缠方法得到地表形变图，最后把形变图导入ArcGIS中进行分后续处理。

双轨法D-InSAR重复轨道干涉测量模式的干涉相位组成为［12］:
①=①切+①de/+①at”+①”。

ise+①曲（1）
式中，①切为地形相位；O de/为视线方向变形引起的相位,是有效相位；叽大气延迟相位；O nQise为各种斑点噪声、热噪声引起的噪声相位；％为平地相位。

双轨法差分干涉测量技术采用外部参考DEM进行地形相位反演，通过干涉相位和反演的相位差分得到去除地形相位的地表形变图。

具体流程如图2所示。

图2双轨法D-InSAR流程图
Fig・2D-InSAR processing flow chat
2数据预处理
Sentinel-1A雷达影像原始SAR数据首先需进行聚焦处理和基线估算，得到可以进行干涉处理的单视复数影像一一SLC数据。

通过下载和原始数据时间相对应的POD精密定轨星历数据进行聚焦处理，可消除轨道误差引起的系统性误差。

根据试验区影像范围，对SLC数据按照地理坐标进行裁剪,再进行基线估算,估算结果表明空间垂直基线远小于临界基线长度（见表2），可以进行差分干涉处理。

表1实验数据主要参数
Tab・1Main parameters about experiment data
产品级别主影像辅影像极化方式成像模式
入射角度
（°）SLC_12019.6.142019.7.8VV IW39.38
表2实验数据的基线特征
Tab・2Baseline characteristics of experimental data 主影像从影像
时间基线空间基线临界基线
(d)(m)(m) 2019年6月14日2019年7月8日24-55.5076
091.829
216测绘与空间地理信息2021年
3地表形变信息提取
1)干涉图滤波和相干性计算
根据SLC数据对的像素间距及成像入射角大小，设
置制图分辨率为20,对距离向和方位向分别进行1倍、5
倍的多视处理，主、副影像经过对应像元之间复数的共轭
相乘得到干涉条纹图(如图3所示)o此时的干涉图中含
有大量噪声误差，对相位解缠的效率和精度有很大影响，
因此必须进行去噪处理。

本文选取软件自带的Goldstein
滤波方法去除噪声误差,相对于Adaptive和Boxcar滤波，
此方法既能提高干涉图信噪比，又可以减少基线引起的失相干噪声。

滤波后的干涉图和相干性系数图如图4、图5所示。

图3干涉图
Fig・3Interferogram
图4滤波后干涉图
Fig・4Interferogram after filtering
2)相位解缠
相位解缠是双轨法D-1nSAR技术处理的一个重要环节,解缠的效果是成像质量好坏的关键之一。

实验选取3种相位解缠方法进行地表形变信息的提取，设置相干性阈值为0.2,只对大于该相干性的像元进行相位解缠,小于 0.2的做掩膜处理。

论文对比了区域增长法(Region Growing)、最小费用流(minimum cost flow MCF)和Delaunay MCF共3种常见解缠方法的结果。

加载滤波后三种相位解缠方法对应的解缠图，如图6所示。

图5相干性系数图
Fig・5Coherence coefficient diagram
(a)Del a unay MCF(b)最小费用流(c)区域增长法
图6相位解缠图
Fig・6Phase unwrapping map of different methods
3)轨道精炼与重去平
轨道精炼之前需要通过选择控制点对轨道精炼的参数进行修正,为了提高控制点的质量,人工选取远离形变区、干涉效果良好的区域,多次检查精炼之后生成的文件并剔除质量低的控制点。

轨道精炼和相位偏移的计算可以消除斜坡相位,矫正卫星轨道及相位偏移，是解缠后的相位能否转化为形变值的关键步骤。

轨道精炼与重去平后，对3种相位解缠方法形变结果的精度分析对比，见表3,Delaunay MCF解缠方法精度最高,像素DN值在轨道精炼与重去平前后一致性最好,结果最可靠。

表3不同相位解缠方法精度指标对比
Tab・3Comparison of precision indexes of
different methods
解缠方法
精度
平均值
最小值最大值
Delaunay MCF-0.149169 2.0654020.119765
最小费用流-0.861952 4.7123320.240649
区域增长法-0.149462 2.8695640.170006
4)地理编码与形变图生成
地理编码是差分干涉处理的最后一个流程，差分干
涉测量经过一系列技术流程最终得到的产品是在雷达坐
第6期李铭:宜宾地震同震形变场双轨法D-InSAR技术的研究217
标系下完成的,不能够直接使用,为了方便解译和理解,需将产品经过地理编码转换至地理坐标系下。

地理编码可以纠正由地形起伏造成的SAR影像几何形变。

经过地理编码得到雷达视线方向(LOS)地表的形变图和DEM (如图7所示)o为了形象描述本次地震引起的地表形变,以震中为中心，截取在东南、西北方向上的剖面图，如图8所示。

编码后LOS方向形变图生成的DEM图
图7D-InSAR视线方向形变结果图
Fig・7D-InSAR line of sight deformation result graph
图8视线向形变剖面图
Fig・8Line-of-sight deformation profile
形变图的像素值表示地表的起伏变化，正号代表地面上升，负号代表地面下沉。

图7可知，震中位置并没有导致地面明显下沉,反而有所抬升，震中地表形变量约3.5 cm；震中东南方向约6km处引起地表最大抬升7.8cm，位于震中西北方向23km处地表沉降量为6.7cm,原因是丰富的余震引起的地面大面积下沉。

研究表明，此次地震引起地表最大沉降量8.4cm,最大抬升量9.3cm,连续的余震造成地表抬升与沉降区形成一条长约38km东南一西北走向的断层。

4结束语
1)本次研究通过双轨法差分干涉测量技术(D-InSAR)获得了2019年6月17日宜宾Ms6.0级地震引起的地表形变量,得到了该区域地震同震形变场,分析了以震中为中心在各个方向上的形变量，为抗震决策和地震风险评估提供了参考依据。

研究结果表明，Sentinel-1卫星SAR数据结合90m分辨率的SRTM DEM可以很好地进行大面积区域的D-InSAR地表形变监测，监测精度能够达到毫米级，在地表形变监测方面有非常大的应用潜力。

2)相位解缠是双轨法D-InSAR技术的重要环节，对3种相位解缠方法形变结果的精度进行了分析对比，实验证明,选取合适的解缠方法能够提高实验精度并影响形变监测结果。

3)本次研究获得的只是雷达视线方向上的形变量，而震区实际地表形变可能是任意方向，还需与GPS、高精度水准测量等真实观测数据进行对比分析，建立模型模拟反演,才能精确地了解此次地震形变特征,确定形变机制。

另外,滤波并不能完全消除噪声相位残差，大气延迟误差的改进需进一步处理，外部DEM误差对形变结果的精度也有影响，因此对于D-InSAR技术的研究还有很长的路要走,关键技术还需要更深层次的探索。

参考文献：
[1]刘云华，屈春燕,单新建，等・SAR遥感图像在汶川地震
灾害识别中的应用[J].地震学报，2010,32(2)：
214-223.
[2]王桂杰,谢谟文,邱骋,等.D-INSAR技术在大范围滑坡
监测中的应用[J].岩土力学，2010,31(4)：1337-
1344.
[3]李晖，肖鹏峰,冯学智,等.基于重轨InSAR的积雪深度
反演方法[J].冰川冻土，2014,36(3)：517-526.
[4]HAO K X,SI H,FUJIWARA H,et al.Coseismic surface
-ruptures and crustal deformations of the2008Wenchuan
earthquake Mw7.9，China[J].Geophysical Research
Letters,2009,36(11)：192-200.
[5]MASSONNET D,ROSSI M,CARMONA C,et al.The dis­
placement field of the landers earthquake mapped by radar
interferometry[J].Nature,364(6433):138-142.
[6]徐小波，单新建,屈春燕，等.2008、2009年柴达木盆地
地震群InSAR分析及同震形变场反演[J].地震地质，
2017,39(3)：485-496.
[7]臧妻斌，黄腾，郭献涛.基于DInSAR技术的意大利拉奎
拉同震形变研究[J].河南理工大学学报(自然科学
版),2017,36(1)：52-57.
[8]杨亚莉.神东矿区开采沉陷干涉测量与时空变化规律
分析[D].焦作:河南理工大学,2016.
[9]王钰.基于D-InSAR及数值模拟的张双楼矿开采沉陷
规律研究[D].焦作：中国矿业大学,2019.
[10]王燕燕，于海洋.九寨沟M_S7.0级地震三维同震形变
场提取[J].测绘科学,2019,44(5)：8-13.
[11]杨金明，刘志辉.Sentinel-1卫星数据产品应用探讨
[J].地理空间信息,2016,14(12):18-20.
[12]杨俊凯,范洪冬,赵伟颖,等.基于D-InSAR技术和灰色
Verhulst模型的矿区沉降监测与预计[J].金属矿山，
2015(3):143-147.
[编辑:任亚茹
]。