联合slipBERI与D-InSAR技术的青海省玛多县地震形变场提取与模拟

付玉蒂, 薛东剑, 张小轩

(成都理工大学地球科学学院, 成都 610059)

青海省玛多县于2021年5月22日2:04发生里氏7.4级地震,震源深度可达17 km,震中坐标为(34.59°N,98.34°E)。由于地震震源浅、震级大、地面震动强烈,使得震中地区房屋和牲畜棚遭到破坏,道路、桥梁和其他基础设施出现不同程度的隆起或倒塌。对地震形变进行监测能有效地评估震区破坏分布以及滑动断层的特征,用于地震形变监测的传统技术包括GPS技术和传统的水准监测技术等,但其对大面积形变场的监测始终存在空间分辨率不足等问题[1-2]。基于这一问题,应用差分干涉合成孔径雷达(differential interferometric synthetic aperture radar,D-InSAR)进行大面积地表形变监测成为近年来的一大趋势。它具有毫米级的视线精度,特别适用于监测快速而剧烈的地表形变,如地震、火山等。

Massonnet等[3]论证了采用D-InSAR提取地震同震位移场的有效性,并且成功提取了加州Landers地区形变场,让InSAR技术开始受到广泛关注;Budi等[4]采用D-InSAR技术提取龙博地区的地震形变场,分析得到其地震分布特征。余祥伟等[5]采用D-InSAR方法于2020年提取并模拟宜宾市地震形变场,验证了观测结果的可靠性;李成龙等[6]于2021年结合升降轨InSAR形变场,反演新疆伽师县地震断层参数以及滑动分布,并研究了其发震构造;韩炳权等[7]提取了2022年四川泸定地震形变场,并且借鉴Okada弹性位错模型思想,进行滑动断层反演,以验证观测值的有效性;高二涛等[8]采用D-InSAR提取了九寨沟升降轨地震形变场,并对其位置、大小进行分析,为下一步研究断层位置提供了条件。

综上所述,D-InSAR在监测瞬时而剧烈的形变方面(如地震)具有精度高、结果准确等优势,现将slipBERI方法与D-InSAR技术相结合,提取青海省玛多县地震同震形变场,采用slipBERI方法对地震滑动断层进行反演,结合模型和形变参数,提高低相干区域的精度,为研究此次地震的断层滑动、玛多地区地震防灾减灾等方面提供一定依据。

1 研究区概况

1.1 地理位置

玛多县地处青海省果洛藏族自治州西北部,地理坐标为33°50′N～35°40′N,96°50′E～99°20′E,北邻都兰县,东临海南自治州兴海县、玛沁县,南临达日县、四川省石渠县,西部为玉树自治州曲麻莱县,西南面则与巴音喀拉山、玉树藏族自治州称多县接壤,如图1所示。玛多县境内有扎陵湖和鄂陵湖两大主要湖泊,中国第二大河黄河也流经其境内。

图1 玛多县地理位置示意图Fig.1 A map of Maduo County’s geographic location

1.2 地质构造及地震背景

玛多县地属高平原区,大部分地区海拔在4 500 m以上,地势较为平坦。大地构造带属巴颜喀拉褶皱带,构造线自西北向东南延伸。近年来,巴颜喀拉地块先后经历了多次大地震,使其及周边地块成为地震高发区,如2001年青海昆仑山口Ms8.1地震(边界带中部)、2008年新疆和田Ms7.3地震(西段)、2010年玉树Ms7.1地震(区块南部边界)、2013年芦山Ms7.0地震(区块东部边界附近)等[9]。处于青藏高原东北部的玛多地震,位于玛多-甘德断裂带上,距巴颜喀拉块体边缘的东昆仑断裂带70 km[10],是发生于巴颜喀拉块体的又一次大地震,与板块周边断裂带高度相关[11],玛多地震及其周围地质构造如图2所示。

2 数据概况及研究方法

2.1 数据概况

(1)合成孔径雷达(synthetic aperture radar,SAR)数据。本次研究获取覆盖玛多县2景地震前、后Sentinel-1A卫星降轨SAR数据,详细参数如表1所示。

表1 SAR数据参数Table 1 SAR data parameters

(2)精密轨道数据及数字高程模型(digital elevation model,DEM)数据。精密定轨星历数据(precise orbit ephemerides,POD)是sentinel-1卫星最精确的轨道数据,在全球导航卫星系统(global navigation satellite system,GNSS)下行21 d后才能使用。每天生成一个文件,可覆盖26 h,定位精度优于5 cm。POD精密定轨星历数据用于参考校正,能有效地减少形变相位误差,提高监测精度。同时,处理过程中采用DEM数据进行地形校正和地理编码,基于研究区域的地理位置,选取30 m分辨率的航天飞机雷达地形测绘使命-数字高程模型(shuttle radar topography mission-digital elevation model,SRTM-DEM)数据作为参考DEM数据,如图3所示。

图3 玛多县DEM图Fig.3 DEM of Maduo County

2.2 研究方法

2.2.1 双轨差分测量原理

1.2.2 计量分析法 主要使用CiteSpaceV分析系统绘制知识图谱，它是基于JAVA平台的可视化应用软件，通过该软件对文献数据分析处理绘制出直观易懂的科学图谱，分析探测相关研究的重点、热点、前言、变化趋势等[2]，也是一款着眼于分析科学分析中蕴含的潜在知识，并在科学计量学、数据和信息可视化背景下逐渐发展起来的一款引文可视化分析软件[3]。借助知识图谱可视化技术，探究我国体育教学评价领域的研究现状、热点及发展趋势，把握该领域发展动态及演化规律，从而为我国体育教学评价的进一步研究提供理论支撑。

针对地震这种具有大尺度、瞬时等特点的形变,采用双轨差分干涉测量方法对其形变场进行提取。双轨法需要地震前后各一景SAR影像,基本思想是对影像做干涉处理,然后,以外部DEM为参考,对地形相位进行反演,在干涉相位图中将其去除,提取出震区的形变信息。基本原理如图4所示。

A1、A2为传感器在地震前后两次观测P点的位置;R1、R2为传感器到P点之间的距离;A1、A2之间的线性距离是基线L;α为基线L与水平面的夹角;γ为入射角;β⊥为垂直基线;β∥为平行基线;h为观测点P的高程;H为卫星飞行高度;δ为R1、R2之间的距离差图4 D-InSAR基本原理图Fig.4 D-InSAR basic schematic diagram

获取形变信息的关键是从干涉图中去除地形信息,以数字高程模型数据为参考对地形相位进行反演,并在干涉图中进行二次差分,以消除地形相位,获得震区的形变相位。采用双轨差分提取地震形变场,具有数据易获取、运算速率快、精度高等优点。记Δφ为两幅SAR图像的相位差,在SAR影像的获取过程中,大气等各种因素也会影响电磁波的传输,即Δφ表示为

Δφ=φflat+φtopography+φmove+φatomosphere+φnoise

(1)

式(1)中:φflat为干涉相位受平地效应的影响;φtopography为受地形相位影响的部分;φmove为遥感平台移动时引起的干涉相位误差;φatomosphere为大气相位延迟;φnoise为噪声引起的相位误差。其具体流程如图5所示。

图5 D-InSAR技术流程图Fig.5 D-InSAR technical flow chart

2.2.2 数据处理关键步骤

1)滤波

去平后的干涉图(图6)中还有大量的噪声斑点,这难免会对形变信息的提取产生一定影响,导致干涉图的相位值产生偏差或不连续,这不仅会影响后续解缠质量,还会影响最终结果,影响监测精度。因此,应通过滤波去除噪声斑点。3种常用滤波方法包括Adaptive、Goldstein以及Boxcar,其中Adaptive滤波用于处理高分辨率影像;Boxcar滤波采用局部滤波原理,可以有效地保留微小条纹;Goldstein滤波方法使用率最高,它之所以能有效抑制噪声斑点、改善条纹可见度,是因为其具有可变的滤波器。在处理过程中,考虑到图中条纹的清晰度和细条纹的存在,选择了Goldstein方法,将窗口值扩大至4,改善滤波效果,滤波后的差分干涉图(图7)噪声显著减少,干涉条纹清晰,有利于后续相位解缠处理。

图6 玛多地震区差分干涉图Fig.6 Differential interferogram of Maduo earthquake area

图7 Goldstein滤波后的差分干涉图Fig.7 The difference interferogram after Goldstein filtering

相位解缠在整个D-InSAR处理过程中十分关键,相位解缠的结果与形变信息提取的效果密切相关。在未解缠时,干涉图的相位介于[-π,π],此时相位与实际情况相差2nπ,对它解缠是将此时的相位值转换为真实相位值的过程,解缠效果也会直接影响形变信息的质量。根据研究区域的特点,本文研究中采用最小费用流(minimum cost flow,MCF)法进行相位解缠,得到解缠相位(图8)。

图8 相位解缠图Fig.8 Phase unwrapping

3 形变场分析

采用双轨差分干涉方法对玛多震前震后两景SAR数据进行处理后,通过干涉、滤波、相位解缠等各项优化,最终得到玛多地震形变场,如图9所示,其形状近似于一个呈现西北-东南向分布的不规则椭圆,这与玛多县的构造线(西北至东南)走向一致。以AA′为参考界线,AA′下侧为隆升区,剖面上侧为沉降区。从中部往南北方向,形变增大,且北部形变小于南部形变,最大视线向(line of sight,LOS)向上升形变量为0.81 m,最大下沉量为0.65 m。AA′形变剖面(图10)沿发震的断层展布,其形变值在0左右,变化极小;BB′形变剖面(图11)穿过隆升和沉降区,最大沉降量为0.55 m左右,从B开始穿过发震断层,形变量从负值快速增加到该剖面最大隆升值约0.45 m。这与华俊等[12]的InSAR形变场结果(最大形变约0.9 m)较一致。

图9 青海省玛多县地震形变图Fig.9 Earthquake deformation map of Maduo County, Qinghai Province

图10 AA′形变剖面图Fig.10 AA′ deformation profile

图11 BB′形变剖面图Fig.11 BB′ deformation profile

玛多地震有明显的地表破裂痕迹,这与李智敏等[13]的研究一致,破裂带北盘沿视线方向呈拉伸运动(变形约0.65 m),南盘沿视线向呈缩短运动(变形约0.81 m),反映了玛多地震的主要变形场为东西向水平运动,具有明显的左旋走滑特征。断层两侧的相对视线向位移可达1.50 m,同震形变场的长轴方向(断层走向)整体上为北-西-西(NWW),表明地震引起的地表破裂错动位移明显。断层滑动特征与杨君妍等[14]和王守文等[15]研究结果一致。

根据此次地震的形变信息、形变特征以及断层走向,可以判断该破裂带位于巴颜喀拉块体,为昆仑山口-江口断裂。此次地震导致巴颜喀拉地块向东强烈挤压,从而导致北部走滑断层在构造作用下左旋运动,与之前发生在巴颜喀拉地块外围边界的中强地震不同,玛多7.4级地震是近年来该区块内发生的唯一一次强震。

4 玛多地震断层滑动分布反演

4.1 反演方法

将slipBERI[16-17](slip from bayesian regularized inversion)滑动断层反演方法用于玛多地震滑动断层的反演,该方法与常用断层反演方法如梯度下降法、Okada弹性位错模型不同,其使用Okada模型计算表面变形的解析解,再引入冯·卡曼自相关检验与地震的分型特性相结合,利用贝叶斯方法对输入的地震断层数据(InSAR)的滑动进行求解,在滑动模型的建立中,能对光滑度以及反演模型进行选择,并且其自带模拟退火算法用于选择贝叶斯方法起始值,整体操作简便,适合用作断层滑动分布反演。

4.2 断层滑动分布分析

反演原始数据使用D-InSAR处理得到的玛多地震形变数据,在反演过程中采用冯·卡曼正则化结合贝叶斯方法,通过马尔科夫链蒙特卡罗(Markov Chain Monte Carlo,MCMC)链进行数字近似求解,光滑方式选择冯·卡曼正则化配套光滑方式进行反演。断层相关参数如宽度及深度使用Mai等[18]的公式,计算公式为

awid=1.860+0.34alen

(2)

adip=-0.390+0.43awid

(3)

式中:alen为断层长度;awid为断层宽度;adip为断层深度。

断层滑动分布反演结果如图12所示,从图12中可见,左旋走滑是玛多地震同震滑动主要走向,与D-InSAR得到形变场结果相同,沿断层走向破裂为85 km,接近地表滑动量大多在0～1 m,最大滑动量约为5 m处于断裂带东部,位于地下10 km左右,东侧滑移量明显大于西侧滑移量,断层主体破裂带位于地下5～15 km,断层参数与全球质心矩张量目录(global centroid-moment-tensor project,GCMT)给出的断层参数基本一致,反演得出地震矩约为1.108 8×1020NM,根据式(4)和式(5)得出震级约为Mw7.36,与GCMT所给出震级Mw7.4基本相当,如表2所示。

表2 滑动断层参数Table 2 Slip fault parameters

图12 反演得到的玛多地震滑动断层分布图Fig.12 Earthquake slide fault distribution map obtained by inversion

(4)

(5)

4.3 反演结果分析

InSAR数据滑动断层分布反演拟合情况如图13所示。图13(a)为地震原始形变图,图像中红线为此次断层轨迹,沿着这一轨迹附近反演效果良好;图13(b)为slipBERI方法结合InSAR数据模拟得到的形变图,断层上部为沉降区,下部为隆升区,断层附近形变量在-0.2 m和0.2 m左右,与原始数据基本吻合;图13(c)为前两幅图像拟合得到的残差图,大部分区域残差小于10 cm,但图像右下角区域拟合效果不佳,对照地震形变数据,该区域形变复杂,断层可能还存在另一条分支断层,这与余鹏飞等[19]的研究结果一致。

图13 slipBERI模拟结果Fig.13 slipBERI simulation results

5 结论

综合D-InSAR技术监测和slipBERI模拟得到的玛多地震形变场结果,可以得出以下结论。

(1)玛多地震形变场形状近似于一个不规则椭圆,受西北-东南走向的左旋走滑断层控制,断层下侧为隆升区,最大视线向隆升量达0.81 m;断层上侧为沉降区,其最大视线向沉降量为0.65 m。此次形变场断层整体呈北-西-西(NWW)走向,根据形变信息、形变特征和断层走向,可以判断出该断裂带位于巴颜喀拉块体内,属昆仑山口-江口断裂。

(2)slipBERI模拟得到的形变场与观测结果相似,大部分地区的残差小于10 cm,这不仅验证了观测结果的可靠性,而且反演结果连续完整,形变场边界清晰,有助于恢复大气造成的观测结果缺失等问题。