基于D-InSAR技术的当雄地震形变场提取研究

杨　珍，张永志，吴　然，张文军，叶　凯

(1.长安大学 地质工程与测绘学院，陕西 西安 710000；2.甘肃省测绘地理信息局 地图院,甘肃 兰州 730000)



基于D-InSAR技术的当雄地震形变场提取研究

杨珍1，张永志1，吴然1，张文军2，叶凯1

(1.长安大学 地质工程与测绘学院，陕西 西安 710000；2.甘肃省测绘地理信息局 地图院,甘肃 兰州 730000)

摘要：随着D-InSAR技术在同震形变量测方面的优势日益突出，文中以2008年10月6日当雄Mw6.3级地震为对象，采用ENVISAT搭载的ASAR传感器c波段在图像模式下的5景level 0级影像，以90 m分辨率的SRTM DEM作为外部DEM，使用GAMMA软件采用二轨差分干涉测量的方法获取当雄地震同震形变场，并对其进行分析，确定震中位置为90.372°E，29.734°N，得到同震垂直形变约为33 cm，形变量精度达到厘米级。

关键词：D-InSAR；地震；ENVISAT；二轨差分；地震形变场

2008-10-06T16 p0 min，在西藏自治区当雄县发生了Mw6.3级地震，监测震间微小形变场，研究断裂活动引起地表位移的时空演化特征，捕获可能的震前形变异常，为地震的监测预测提供可靠依据，推进地震预报工作是地震研究迫切需要解决的问题。

合成孔径雷达差分干涉测量(D-InSAR)技术是雷达干涉测量应用的一个拓展，雷达干涉图的差分可以监测雷达视线方向厘米级或更微小的地球表面形变，以揭示许多地球物理现象，如地震形变、火山活动、冰川漂移、地面沉降以及山体滑坡等[1]。

近年来，国内外很多学者利用D-InSAR技术研究地震形变，许才军等利用InSAR数据提取及分析汶川地震形变场[2]，单新建等研究了汶川地震前的垂直形变场变化特征[3]。但是以往的研究大多针对的是大地震，本文选取强震级别的当雄地震为研究对象，基于ENVISAT ASAR数据，采用D-InSAR方法，提取及分析当雄Mw6.3级地震形变场。

1D-InSAR基本原理

InSAR测量的相位由5部分组成[4]：

(1)

等式右端的5个分量分别为地形相位、所求形变相位、平地相位(参考相位)、大气相位、噪声引起相位。地表变化前后的2幅影像的干涉图中同时包含地形信息和形变信息，由式(1)知差分干涉的目的是使用一定的方法去掉一些相位信息，最终获得形变相位的过程[5]。

以二轨法为例，将2幅SAR影像进行干涉生成干涉图，生成的干涉图包含地形变化信息，随后再利用DEM数据进行相位模拟，模拟得到地形信息的条纹图，再将SAR影像干涉图与DEM模拟的条纹图叠加去除地表形变，得到地表形变信息。

2基于GAMMA软件的二轨法数据处理流程

本文采用GAMMA软件来实现二轨法数据处理，具体的实现流程如图1所示，首先在聚焦模块下，利用仪器的参数文件及精密轨道文件对数据进行聚焦得到SLC文件，然后与常规差分处理流程一样，包括影像的配准、滤波、干涉图生成、基线估算、外部相位模拟、差分处理、相位解缠、形变量计算、地理编码等处理步骤。

图1　GAMMA软件二轨法数据处理流程

3数据处理

3.1数据源

研究区域当雄地形起伏在4 000～6 000 m，地形高低错落，沟谷纵横。本文获得了5景ENVISAT ASAR level 0 级数据，天线文件和ASAR仪器参数文件(ESA提供)及精密轨道数据文件(DELFT提供)，数据情况见表1。

表1　ENVISAT ASAR数据列表

DEM数据为SRTM DEM，分辨率30 m，平面基准WGS-84，高程基准EGM96，标称精度±10 m。

3.2数据选择

选择适合干涉的SAR影像是进行任何干涉处理的首要步骤，也是非常关键的步骤。数据源的好坏直接影响到后续步骤，为了得到最好的结果，数据选择需要考虑的参数包括传感器类型、视角、几何基线/空间基线、时间基线、影像获取时刻、相干性、气象条件。

将5景数据处理成SLC数据后，对数据进行滤波处理(前置滤波)，配准以与地震发生时间最近的20080921.slc影像为配准主影像，将其余四景影像配准至20080921，配准方法采用精密轨道数据和强度互相关及偏移量跟踪法，配准精度可达到1/10像元。再根据表2，干涉成像。基线估计、差分处理，得到10幅干涉图，时空基线的结果见表2。

表2　当雄Envisat ASAR干涉组合

根据卫星影像的覆盖范围，InSAR数据处理对空间基线、时间基线的要求，差分干涉的效果，故选取20080921～20081026和20080921～20090104，如图2所示，进行后续操作，最后选取相干性好的一幅来提取同震形变场。

图2　差分处理后的干涉图

3.3滤波处理

在影像进行精确配准前，需要进行频率域的滤波处理，即前置滤波，前置滤波是在方位向和距离向分别进行的。一般而言，如果频谱偏移过大，方位向滤波应该在配准前进行，可提高配准精度，而距离向滤波则必须在配准后的主从影像之间才能进行。方位向和距离向前置滤波的具体流程相类似，首先都是对主从影像在方位向和距离向进行频率域分析，然后通过带通滤波器将两者的频谱非重叠部分滤除。

图2得到差分干涉图去除地形相位及平地效应，但其他因素导致的随机噪声也会影响干涉图的质量，对影像也要进行滤波(后置滤波)。Adaptive滤波运算速度快、效率高、效果好、占用内存小，对噪声较小的区域能够得到正确的解缠值，且能够较多保存边缘信息，对于高分辨率数据处理效果非常好，如TerraSAR-X或COSMO-SkyMed；Boxcar滤波是使用局部干涉条纹的频率来优化滤波器，该方法尽可能的保留微小的干涉条纹；Goldstein自适用滤波方法的滤波器是可变的，提高干涉条纹的清晰度，减少由空间基线或时间基线引起的失相干的噪声，是目前最流行的方法。

本文采取了Goldstein自适用滤波的方法，且进行两次，第一次滤波窗口64×64，第二次滤波窗口32×32。

3.4相位解缠与去除基线误差

相位解缠方法很多，大致分为两类：①基于路径的相位解缠算法；②基于最小范数的相位解缠算法。最常用效果较好的方法有区域增长法、最小费用流法、枝切法、最小二乘法。本文通过比较，最后采用最小费用流法。在进行相位解缠时避开低相干性区域，生成有效掩膜文件，提取相干系数较高的点，建立Delaunay网，使缠绕相位与解缠相位的差异最小化，然后将最小化问题转换为最小费用流问题。

据有关实验研究证明在海拔5 000 m的高原水汽变化引起的误差小于1 mm，当雄海拔高，属于干燥性气候，水汽总含量相对较低，所以大气误差造成的影响本文不考虑。

采用Delft精密轨道数据，轨道精度为5～10 cm，基线误差属于系统误差，本文采用远场最佳二次多项式拟合的方法去除非线性轨道残余相位信息，最后得到去除基线轨道误差的干涉图。

3.5形变量计算

形变量计算是将差分解缠相位转换为视线(LOS)向形变量，LOS向形变同时包含了东西向、南北向、及垂直向3个方向的分量，经过计算分析发现本次地震LOS向形变在东西向和南北向的形变量非常小，经过处理转换最终得到地表垂直方向形变。

3.6地理编码

计算出每个像素对应的地表几何信息后，得到一个在影像坐标系下的点阵图，还需要把各种数据从影像坐标系转换到统一的地理坐标系下。将雷达坐标系下的原始数据通过一定的几何校正方法消除由轨道、传感器、地球模型引起的扭曲和畸变，然后变换到某种制图参考系中，地理编码之后影像坐标系与提供的外部DEM坐标系一致，故形变监测结果为WGS-84坐标系统。图3为地理编码后的形变图，总体上两幅图的同震形变场保存高度一致，图3(a)相干性较好，图3(b)失相干现象比较严重，故本文选取地理编码后的20080921～20081026形变图进行结果分析。

图3　地理编码后的形变图

4当雄地震形变结果分析

根据图4地震同震形变图，从宏观上分析，沉降区的中心位置就是震中位置，形变条纹左上侧圆环沉降，区域面积大于右侧的隆升区，越靠近地表破裂带，干涉条纹越密集，地表形变量越大，越远离地表破裂带，干涉条纹越疏松，地表形变量越小，地表破裂带沉降区的最大变形量约为33 cm，隆升区最大变形量约为5 cm，圆环的中心位置即为震中位置。可判断出震中位置在90.372°E，29.734°N。与乔学军推算的震中位置90.374°E，29.745°N非常接近[6]，相差约640 m。

4.1剖面分析

取两条测线进行剖面分析，两条测线的位置均过震，一条对应雷达的视线方向且近似垂直地表破裂带[7]，如图4中的WE，一条近似平行于地表破裂带，如图4中的NS，剖面图结果如图5所示。

图4　当雄地震同震形变

从图5(a)中可以看出，从起点至大约12.5 km处，形变值约从-0.03变化到-0.32，在大约12.5～19 km处，形变值又由约-0.32变化到0，在约19～31 km处，形变值在0.03左右波动，这说明地表变形从沉降逐渐变为隆升，最大沉降量约为32 cm。从图5(b)中可以看出，剖面图近似呈对称状，在距离起点10.5 km处，沉降量达最大值约33 cm，说明离震中越近，地表位错变换越剧烈。

4.2等直线分析

对同震形变做等值线图[8]，如图6所示，可以看出最大形变量在32 cm左右，根据图6的形变场特征可以看出，与图4、图5中的形变场保持高度一致，等值线中因部分像元无变形值出现“孤岛”现象，导致这种情况的原因有多种，大气改正失败，相位解缠失败，山体滑坡等，另一方面来说，对等值线进行分析，可评判InSAR结果的质量。

图5　同震形变场剖面图

图6　同震形变场等值线

5结束语

本文通过D-InSAR的两轨差分方法提取当雄地震的同震形变场，得到的同震垂直形变约为33 cm，确定的震中位置：90.372°E，29.734°N。采用Goldstein自适应滤波2次，减少斑点噪声，相位数据连续，周期性更加明显。但是本文没有考虑大气误差对形变结果的影响，还有待进一步研究。D-InSAR技术对于跨断层的震间形变监测具有较大潜力，特别是对GPS测站分布稀疏的地区，利用D-InSAR获取的震间形变场对于理解活动断层的地壳形变特征、地球动力学问题具有重要的应用价值。

参考文献：

[1]廖明生，林珲.雷达干涉测量：原理与信号处理基础[M].测绘出版社，2003.

[2]许才军，江国焰，王浩，等.基于GIS的InSAR结果分析方法及在汶川Mw7.9级地震同震解释中的应用[J].武汉大学学报(信息科学版)，2011，36(4):379-383.

[3]单新建，宋小刚，韩宇飞，等.汶川Ms8.0地震前InSAR垂直形变场变化特征研究[J].地球物理学报，2009(11):2739-2745.

[4]吴中海，叶培盛，吴珍汉.2008年10月6日西藏当雄Ms6.6级强震的地震烈度控震构造和发震机理[J].地质通报，2009，28(6):713-725.

[5]张军，尼玛，尚荣波.2008年西藏当雄6.6级地震浅[J].高原地震，2009，20(4):40-44.

[6]乔学军，王琪，杨少敏，等.2008年新疆乌恰Mw6.7地震震源机制与形变特征的InSAR研究[J].地球物理学报，2014(6):1805-1813.

[7]雷广渊，周辉.基于SBAS技术的金属矿山沉陷规律研究[J].测绘工程，2015，24(3):40-46.

[8]谷延超，范菇宇，范东明，等.基于不同滤波算法差异的机载LiDAR数据桥梁提取[J].测绘工程，2014，23(11):67-70.

[责任编辑：李铭娜]

DOI:10.19349/j.cnki.issn1006-7949.2016.08.002

收稿日期：2015-09-06

基金项目：国家自然科学基金资助项目(41374028)

作者简介：杨珍(1984-)，女，博士研究生.

中图分类号：TP751

文献标识码：A

文章编号：1006-7949(2016)08-0005-06

On the extraction of surface deformation caused by Damxung earthquake based on D-InSAR technology

Yang Zhen1,Zhang Yongzhi1,Wu Ran1,Zhang Wenjun2,Ye Kai1

(1.School of Geology Engineering and Geomatics,Chang’an University,Xi’an 710000,China;2.Map Institute,Gansu Administration of Surveying and Mapping Geoinformation,Lanzhou 730000,China)

Abstract：As advantages of D-InSAR technology in coseismic deformation and measurement are becoming more and more prominent,this paper takes Damxung Mw 6.3 earthquake on October 6,2008 as the case,adopting the image of 5 scenes and level 0 of C band in the image mode of ASAR sensor being equipped on ENVISAT,using SRTM DEM with 90 m resolution ratio as the external DEM,and with GAMMA software and the approach of two-pass differential InSAR technology to obtain the coseismic displacement of Damxung earthquake.And analysis is made to determine the location of the epicenter as 90.372°E,29.734°N,gaining about 35 cm of the coseismic vertical deformation,of which the deformation precision reaches the cm-level.

Key words：D-InSAR;earthquake;ENVISAT;two-pass differential;coseismic displacement