InSAR数据约束的2022年1月8日青海门源MS6.9地震发震构造研究

颜丙囤, 季灵运, 蒋锋云, 殷海涛, 陈其峰, 连凯旋

(1. 山东省地震局 聊城地震监测中心站, 山东 聊城 252000; 2. 中国地震局第二监测中心, 陕西 西安 710054;3. 山东省地震局, 山东 济南 250014)

0 引言

根据中国地震台网官方公布,2022年1月8日1时45分,青海省海北州门源县(101.26° E,37.77° N)发生MS6.9地震,震源深度10 km。此次地震造成约22 km的地表破裂,等震线长轴呈NWW走向,长轴约200 km,短轴约153 km。震中位于青藏高原东北边缘的冷龙岭断裂带南部,是继1986年门源MS6.5和2016年门源MS6.4地震后该地区的又一次破坏性地震。冷龙岭断裂带由一组近平行的NWW向左旋逆走滑断裂组成,长度约120 km[1-2]。自21世纪以来,在此次地震震中100 km范围内发生过13次5级及以上的中强地震(https://earthquake.usgs.gov/)。门源地震的区域构造背景如图1所示。详细震源参数信息列于表1,表中全球矩心矩张量(GCMT)计划、美国地质调查局(USGS)测定的断层几何及震源机制由于采用数据不同而存在差异。

图1 2022年门源地震的区域构造背景Fig.1 Regional tectonic background of the Menyuan earthquake in 2022

表1 2022年1月8日青海门源MS6.9地震震源参数

冷龙岭断裂位于青藏高原隆起区的东北边缘,全新世断裂活动强烈,主要表现为左旋走滑运动,伴有正倾滑,断裂地貌特征明显[3-5]。前人在青藏高原东北缘已经开展了大量研究地壳运动的工作。葛伟鹏等[6]认为青藏高原东北缘块体旋转运动是阿尔金断裂左旋走滑作用与沿鄂拉山到阿拉善地块的地壳缩短及海原断裂左旋走滑共同作用的结果。李煜航等[7]研究了青藏高原东北缘的主要活动断裂和广义海原断裂带的空间亏损滑动分布及其耦合变化。郑文俊等[8]认为青藏高原东北缘新生代中晚期以来发生了有序的向外扩展。近几年对海原断裂左旋滑动速率的研究成果相继发表:Zhang等[9]和Burchfiel等[10]认为海原断裂中卫段左旋滑动速率约为8～10 mm/a;Lasserre等[11]认为海原断裂冷龙岭段滑动速率为(19±5) mm/a,老虎山—毛毛山一带为(12±4) mm/a;Gaudemer等[12]认为海原断裂中段滑动速率为(11±4) mm/a;何文贵等[3]、袁道阳等[13]认为冷龙岭一带左旋滑动速率约为4～5 mm/a;Yuan等[14]认为鄂拉山断裂的晚第四纪右旋速率为(1.1±0.3) mm/a、日月山右旋速率为(1.2±0.4) mm/a;Li等[15]通过阶地地貌详细测图和年代样品测试确定海原断裂带左旋滑动速率为(4.5±1.1) mm/a;Zhang[9]等认为海原断裂带滑动速率向东衰减到4～6 mm/a。

此次门源MS6.9地震的孕震区位于高原山区,常规大地测量方法难以获得连续、全覆盖的地表形变场。近年来发展起来的差分干涉合成孔径雷达(Differential Interferometric Synthetic Aperture Radar,D-InSAR)技术是基于InSAR原理,采用差分方法获得的合成孔径雷达差分干涉测量技术。它具有全天候、全天时、覆盖范围广、分辨率高、成本低、精度高(可达厘米或毫米级)的优点,广泛应用于火山运动、冰川运动、地面沉降、地震变形和滑坡的监测[16-18]。D-InSAR技术不仅能减少现场测点的工作量和震后余震现场观测的危害,而且能够反映几何参数、滑动参数和动力学参数。通过分析所反映出的数据参数,反演震源破裂过程,是理解地震成因机制的有效途径之一[19-22]。

基于此,本文利用D-InSAR技术获得2022年1月8日青海门源地震的同震形变场。在此基础上,利用弹性半空间位错模型反演发震断层参数和同震滑动分布,并对形变场、断层活动和发震构造特征进行分析和讨论。

1 InSAR数据处理与同震形变场

1.1 InSAR 数据处理

下载Sentinel 1A干涉宽幅模式SLC IW L1 .1 C波段合成孔径雷达(Synthetic Aperture Radar,SAR)图像作为数据源,数据极化方式为VV极化,距离和方位分辨率为5 m×20 m。地震前后升降轨道的图像参数如表2所列。升轨观测时间为2022年1月5日和1月17日,升轨干涉相对的空间基线为38.3 m,时间基线为12 d;降轨观测时间为2021年12月29日和2022年1月10日,降轨干涉相对的空间基线为55.1 m,时间基线为12 d。短时间基线有助于保持干涉图的相干性,获得高时空采样率的同震形变场。D-InSAR技术采用SARscape软件处理模块,使用的外部参考数字高程模型(Digital Elevation Model,DEM)数据为NASA SRTM DEM(数据精度为30 m),用于消除地形相位的影响。具体步骤是:首先,采用Gold-Stein滤波方法对去地形后的干涉图进行滤波,以提高干涉条纹的信噪比;然后,采用最小费流量相位解缠法实现绝对相位恢复,经过轨道精炼和轨道误差平面去除;最后,进行地理编码将SAR坐标系统结果投影到WGS-84地理坐标系下的地震InSAR同震形变场[23-24]。

表2 升降轨差分干涉影像参数

1.2 InSAR 同震形变场

图2所示为利用双轨D-InSAR技术获得的2022年1月8日MS6.9地震升降轨InSAR干涉图。从图中可以看出,同震特征明显,升降轨InSAR形变干涉条纹具有良好的相干性。InSAR同震形变场(LOS方向)由升降轨数据获得(图3)。如图3所示,形变场整体上基本沿冷龙岭断裂带呈对称分布,覆盖本次地震震中,中间长轴大体呈NWW-SEE向迹线。结合研究区的地质构造背景,推测冷龙岭断裂可能为本次地震的发震断层,断层迹线两侧的形变场存在明显差异。升轨形变场范围南北约20 km,东西约30 km,正值表示它向卫星方向移动,地面显示隆起,最大视距形变量约为43 cm;负值表示它远离卫星,地面显示沉降,最大视距形变量约为61 cm。降轨形变场范围南北长约20 km,东西长约30 km,南北两盘最大视距形变量分别约56 cm和62 cm。升轨观测的沉降值略小于降轨,这种差异主要是升、降轨观测角度不同造成的。根据本次地震的发震构造与升降轨LOS沉降和抬升的关系,同一轨道形变场的上、下盘呈现相反的运动状态。这一现象表明地震引起的地表形变主要是水平运动,符合走滑断层地震的运动特征。同震视线向形变场的断层附近出现了较为严重的失相干现象,主要原因是发震断层向地表断裂,形变梯度大,形变相位不连续,这一现象与现场调查结果一致。因此认为2022年1月8日青海门源MS6.9地震具有较大的走滑分量,且发震断层呈NWW-SEE走向的左旋走滑特征。

图2 2022年门源MS6.9地震InSAR升降轨干涉条纹图Fig.2 InSAR ascending and descending interferogram for the MenyuanMS6.9 earthquake

图3 2022年门源MS6.9地震InSAR升降轨同震形变场(LOS方向)Fig.3 InSAR ascending and descending coseismic deformation field for the Menyuan MS6.9 earthquake (LOS direction)

2 InSAR震源机制反演及发震构造分析

2.1 震源机制反演方法

本文在InSAR同震形变分析结果的基础上,以青海门源MS6.9同震形变场为约束条件,利用Okada弹性半空间位错模型反演研究发震断层的几何参数和同震滑动分布特征[25-26]。主要技术路线包括:第一步,以升降轨LOS形变场降采样结果为约束,采用非线性方法反演均匀滑动断层的几何参数;第二步,基于非线性反演方法得到的均匀滑动断层的几何参数,采用线性反演获得青海门源2022年MS6.9地震断层面的精细滑动分布。

为了获得此次门源地震的同震滑动分布特征,采用ENVI遥感图像处理软件SARscape雷达图像处理模块的Displacement Modeling处理功能,利用Okada弹性半空间位错模型,反演门源地震发震断层的位置、走向、倾角、倾向等几何参数以及滑动分布特征,分析此次地震的断层参数和滑动机制,以深入了解断层活动和深部断裂的特征。由于InSAR同震形变场数据量庞大,且形变结果在空间上是高度相关的,因此在断层几何参数反演之前,首先对干涉图进行降采样处理来获取适当大小的InSAR形变数据集。在实际采样过程中,对近场区域采样点选取相对密集,对于远场区域采样点选取相对稀疏,这样能够最大程度地保留原始形变场的空间特征。采用矢量图层圈定形变场(小区域,采样间隔0.5 km)和模型反演范围(大面积,采样间隔2 km),对升降轨形变场采用均匀矢量网格降采样处理方法,得到1 366个升轨形变数据点和2 360个降轨形变数据点。根据形变场的模拟特征,确定了走向为“NWW-SEE”的发震断层,断层上下两盘分布所示升降轨的反演、模拟结果分别如图4(a)、(b)所示。将反演结果和模拟结果进行拟合,得出残差分布,拟合效果如图4(c)所示。残差绝对值越小的模型,拟合度越好,由此可见二者的拟合度较好。

图4 升降轨形变场采样及反演残差Fig.4 Deformation field sampling and inversion residuals

基于降采样InSAR视线形变数据的点位形变场结果,通过Displacement Modeling模块的非线性反演Non-Linear Inversion工具,结合GCMT推荐的断层几何参数,设置搜索间隔,并利用Okada模型搜索最优解,得到搜索结果。最后,将非线性反演得到的均匀滑动断层的几何参数(包括经纬度、走向、倾角、滑动角、深度等)引入到线性反演模型中,设置线性反演参数,使断层滑动信息得以完整显示。在反演时,InSAR形变场已经破裂到地表,断层的上边界设置为0,走向变化设置为80°～130°,长度变化设置为3～10 km,宽度变化设置为2～8 km,断层倾角变化设置为68°～90°。均匀反演得到的断层参数列于表3。

表3 均匀滑动反演所得断层参数

为了获取断层面上的精细滑动分布,将发震断层沿着走向扩展为30 km,沿着倾向扩展为20 km。通过将反演设置为固定滑动角,将走向和倾角细分为15×15来划分断层;设置阻尼系数为0.05,断层精细反演相对平稳。从图5所示的精细滑动分布反演结果可以看出,本次地震宏观震中为101.29°E、37.78°N,矩震级为MW6.6;断裂规模模型断层长30 km,宽18 km,断裂深度主要集中在0～8 km范围内;断层走向110°,倾角80.5°,平均滑动角12.7°。这次地震是左旋走滑地震,属于高倾角走滑型地震。断层最大滑动量为3.5 m,出现在地下约4 km处,滑动主要集中在地下1～7 km的区域。图5显示的同震滑动分布结果表明,InSAR观测到的升降轨道形变场与相应的拟合结果基本一致,在近场中断层两侧的残余误差明显,主要是由于内插值错误以及无效形变测量值造成的。拟合残差有正负值,在残差图中,残差绝对值越小的模型拟合效果越好,残差点所在的带状区域的宽度越窄拟合效果越好。

图5 同震滑动分布Fig.5 Co-seismic slip distribution

2.2 发震构造分析

上文利用InSAR技术获得的2022年门源地震的同震形变场(图3)表明,同震形变场出现在冷龙岭断裂附近,形变场的长轴方向为NWW,反演确定的震源机制与冷龙岭断层的运动性质一致[27],初步判定此次地震的发震断层主要为冷龙岭断裂。近期研究发现,冷龙岭断层伴随着海原断裂带的运动特征,在101°45′E 处衍生出一条伴生断层——冷龙岭北断层。由于该断层发育至地表,故一直未被发现,直至1986年门源发生地震,才被地质专家研究发现。该断层总长47 km,整体展示出逆冲左旋的运动趋势,其与冷龙岭断层总体呈环形结构,于101°20′E处汇入[28]。滑动分布建模表明该断层为左旋走滑型,与本文的InSAR反演结果一致,表明本文得到的同震形变场反演结果是科学、合理的。

距离此次MS6.9地震最近的一次地震是2016年MS6.4门源逆冲型地震,其震源机制完全不同。距离2016年地震震中位置最近的已知断裂为冷龙岭断裂,该断裂晚第四纪以来主要表现为左旋走滑运动,局部兼具倾滑分量。但2016年地震的震源机制解显示为逆冲型地震,与人们认知的冷龙岭断裂走滑运动性质有所差异。刘云华等[29]认为2016年地震由冷龙岭北侧的弧形次生断层引起,该次生断层与冷龙岭断裂一起构成正花状构造,冷龙岭主干断层近直立插入基底,夹持部分形成隆起断块,推测该MS6.4地震是青藏高原向NW推挤生长,在压扭性作用下隆起断块的一种表现。Wang等[30]利用InSAR技术反演了2016年MS6.4地震发震断层的几何参数和滑动分布,指出冷龙岭断层复杂的聚集形式导致发震断层的滑动集中在小范围内,导致应力降水平较高,加剧了局部应力积累,增加了区域地震灾害的风险,2016年门源地震对2022年门源地震的发生有一定的促进作用。

3 结语

本文应用Sentinel-1A资料,基于2022年1月8日门源MS6.9地震LOS方向的同震形变场,利用D-InSAR技术对此次地震进行了重建,并利用Okada模型反演了断层滑动的分布特征,得出以下结论:

(1) 使用D-InSAR技术获取2022年1月8日门源MS6.9地震同震形变场,发现升降轨数据干涉图的相干性较好,地表破裂沿NWW-SEE方向延伸,地震引起的地表形变范围为30 km×20 km,最大变形量为62 cm。

(2) 基于同震形变场的线性反演,断层滑动分布模型进一步表明本次地震的发震断裂主要为冷龙岭断裂西段,断裂沿NWW-SEE向展布,地震破裂到达地表,与现场地质调查结果一致。

(3) 本次地震宏观震中为101.29°E,37.78°N,矩震级为MW6.6,断裂深度主要集中在0～8 km范围内,断层最大滑动量为3.5 m,大约出现在地下4 km处,滑动主要集中在地下1～7 km的区域。

致谢:感谢欧洲航天局(ESA)免费提供的Sentinel-1A雷达卫星图像！