利用GPS数据反演柴北缘断裂带现今闭锁程度与滑动亏损

刘 洋,邱雨轩,王君毅,李航昊,张 宇,温扬茂,许才军

1. 武汉大学测绘学院,湖北 武汉 430079; 2. 武汉大学地球空间环境与大地测量教育部重点实验室,湖北 武汉 430079; 3. 自然资源部地球物理大地测量重点实验室,湖北 武汉 430079

青藏高原东北缘作为高原向外扩张的前缘部位,是研究高原地壳运动学模式、探索地震孕育和发生机制的天然实验室[1-2]。青藏高原东北缘的断裂闭锁程度、滑动亏损及地震危险性的研究主要集中于块体边界的大型断裂,如海原断裂[3-5]、东昆仑断裂[4,6-7]。针对其他断裂闭锁程度、滑动亏损及地震危险性的研究较少[8]。柴达木盆地北缘断裂带(简称“柴北缘断裂带”)是青藏高原东北缘的重要地质边界,分隔了柴达木盆地和祁连山地块,具有几何形态复杂、断裂数量多、活动性强等特点[2,9-10](图1)。柴北缘断裂带全长约700 km,总体呈NW向,由赛什腾山断裂、绿梁山断裂、大柴旦断裂、锡铁山断裂、阿木尼克山断裂、牦牛山断裂和鄂拉山断裂等组合而成[11-14]。20世纪60年代以来,柴达木盆地北缘已记录多个中强震事件(中国地震台网中心,https:∥news.ceic.ac.cn/):1962年Mw6.8阿木尼克山地震,1977年Mw6.3锡铁山地震,2003年Mw6.6德令哈地震,2008年、2009年Mw6.3大柴旦地震群,2022年Mw6.0德令哈地震。

图1 柴北缘断裂带区域构造Fig.1 Regional tectonic settings of the fault zone in northern Qaidam Basin

目前,基于航卫片解译、野外考察等地质学方法对柴北缘断裂带的活动性已开展了大量研究。文献[15—17]研究得到鄂拉山断裂水平、垂直滑动速率分别为1～4 mm/a、0.15 mm/a。文献[12—13]研究得到锡铁山断裂全新世晚期以来的水平滑动速率为1.81～2.1 mm/a,平均逆冲速率为0.33～0.38 mm/a;阿木尼克山断裂全新世平均垂直速率为0.43±0.02 mm/a。文献[18]研究得到牦牛山断裂晚更新世以来的垂直滑动速率为0.19～0.51 mm/a。然而地质学方法得到的断裂运动特征是百万年尺度的平均值,难以准确刻画其现今运动特征。

现代大地测量数据(如GPS、InSAR)能有效反演断层的现今闭锁程度及滑动亏损,是评价断裂孕震潜能的重要手段[3,19-28]。文献[4,29—30]利用GPS数据反演得到鄂拉山断裂的右旋走滑速率为1～4 mm/a,但文献[29]未考虑块体内部应变对震间速度场的影响,文献[4,30]对青藏高原整体进行运动学块体建模,忽略了区域内次级块体和次级断裂。文献[31—32]利用GPS数据反演得到青藏高原东北缘断裂平均闭锁深度为17～22 km,以及主要断裂的滑动特征,但未研究滑动速率沿走向的变化特征。文献[8]利用GPS数据,基于贝叶斯方法反演得到鄂拉山断裂的闭锁深度(约15 km)及其滑动特征,但未考虑周边断层闭锁效应的影响。

综上所述,不同研究采用不同方法得到的柴北缘断裂带活动特征存在差异。综合分析柴北缘断裂带的现今闭锁程度、滑动亏损及评估其地震危险性具有重要科学意义。本文收集并融合得到较密集的柴北缘GPS数据,基于GPS速度剖面法分析断层运动特征,进而利用TDEFNODE软件建立负位错模型,采用网格搜索和模拟退火方法定量研究柴北缘断裂带的现今闭锁程度、滑动亏损等变形特征,以全面了解该地区块体运动和构造变形模式,并结合历史地震分布分析区域地震危险性。

1 GPS数据及跨断层剖面

1.1 GPS数据

本文采用的GPS数据来自文献[33—35],参考框架均为欧亚参考框架,时间跨度分别为1999—2016年、1991—2015年和1998—2018年。其中,文献[33]通过GAMIT软件处理得到了中国大陆及周边地区共2403个GPS站点的速度场;文献[34]通过PANDA软件处理得到了印度-欧亚碰撞带内2576个GPS站点的速度场;文献[35]采用Bernese软件处理得到了中国大陆及周边地区2260多个GPS站点的速度场。

本文对上述3个数据集进行融合,以得到较密集的柴北缘GPS数据。数据融合中,根据相同GPS测站之间的速度差求解欧拉旋转参数,以文献[33]为参考,将文献[34—35]的GPS数据进行转换。进一步地,本文剔除与区域速度特征差异显著的个别GPS站点,最终选取63个站点的数据参与后续研究(图2)。

图2 GPS数据及跨断层速度剖面的位置Fig.2 GPS data and location of cross-fault velocity profile

1.2 GPS跨断层剖面

为了分析柴北缘断裂带各段的运动特征,本文对GPS速度场数据进行坐标系转换后,采用速度剖面投影法对不同断裂段进行分析。速度剖面投影法是一种简单有效的活动断裂研究方法,通过投影得到GPS站点平行和垂直于断层走向的速度分量,反映了断裂两侧站点的速率随着站点与断裂之间距离的变化情况[36-37]。其原理可表示为

Vv(X)=Ve(X)×cosα-Vn(X)×sinα

(1)

Vp(X)=Ve(X)×sinα+Vn(X)×cosα

(2)

式中,Vv(X)和Vp(X)分别为测站X垂直和平行于断裂走向的速度分量;α为断裂走向;Ve(X)和Vn(X)分别为测站X的东向和北向速度分量。本文横跨柴北缘断裂带共绘制了9条剖面(图2),剖面长为300 km,其中,剖面1横跨赛什腾山断裂,剖面2横跨绿梁山断裂,剖面3横跨锡铁山断裂,剖面4、5横跨阿木尼克山-牦牛山断裂,剖面6、7、8、9横跨鄂拉山断裂。将断裂带两侧站点的速度分别沿平行和垂直断层走向进行投影,利用螺旋位错模型估算断层的走滑速率,根据剖面两侧站点速度平均值之差估算断层的张/压速率[36](图3、图4)。

图4 沿剖面1—9的速度分量(垂直于断层方向)Fig.4 Velocity component along profiles 1—9 (perpendicular-fault direction)

结果表明:剖面1、2、3区域显示为右旋挤压运动,滑动速率由西北向东南逐渐增大,依次约为0.8、1、1.2 mm/a,挤压速率约为0.5、3、1 mm/a;剖面4、5区域显示为左旋挤压运动,滑动速率约为1.5、1 mm/a,挤压速率约为2、1.2 mm/a;剖面6区域显示为右旋拉张运动,滑动速率约为3 mm/a,拉张速率约为0.5 mm/a;剖面7、8、9区域显示为右旋挤压运动,滑动速率为2～3 mm/a,挤压速率从北向南由1 mm/a逐渐增大至2 mm/a。需要特别指出的是,剖面5、6、7同时受阿木尼克山-牦牛山断裂、北鄂拉山断裂等周围断裂闭锁效应的影响,单一剖面无法准确刻画相应断层的运动特征,需要开展后续块体模型反演研究。

根据图3、图4中的剖面结果,结合各段走向、倾向及断层上下盘的差异[11,38-39],本文将柴北缘断裂带分为4段:剖面1、2位于赛什腾-绿梁山断裂段,由茫崖市和平乡至大柴旦西北部;剖面3、4、5位于锡铁山-阿木尼克山-牦牛山断裂段,由大柴旦西北部至都兰县东北部;剖面6、7位于北鄂拉山断裂段,由天峻县阳康乡至都兰县东北部;剖面8、9位于南鄂拉山断裂段,由都兰县东北部至玛沁县西北部。

2 研究方法

2.1 反演理论

假定块体内部点的运动为块体刚性旋转、断层闭锁和块体内部永久应变对应的运动分量之和[40-41],可表示为

(3)

反演结果的优劣采用χ2来评估,其定义可表示为

(4)

式中,dof是自由度;r是观测值的残差;s是观测值中误差;F是观测值权重。χ2约等于1时表示反演结果较好,既没有过度拟合,也没有遗漏必要信号[40,44]。

2.2 块体模型、断层几何及反演参数

基于研究区域(34°N—40°N,92°E—103°E)内山脉、断层的分布及已有研究[30,45],以柴北缘断裂带作为主要边界,本文将研究区域划分成3个块体:阿克塞-德令哈块体、柴达木块体和鄂拉山块体(图5)。本文试验、调整不同的初始条件(各闭合块体边界坐标、边界断层几何形状、倾角和闭锁深度初值等),最终获得较优模型(χ2=1.003)。

注:1-1、1-2、1-3依次表示赛什腾-绿梁山断裂段的西北端、中部和东南端;2-1、2-2、2-3依次表示锡铁山-阿木尼克山-牦牛山断裂段的西北端、中部和东南端;3-1、3-2、3-3依次表示北鄂拉山断裂段的北端、中部和南端;4-1、4-2、4-3依次表示南鄂拉山断裂段的北端、中部和南端。图5 块体及断层划分Fig.5 Block and fault division

反演计算以断层几何形状作为断层参数的输入值,通过分布在断层面上沿等深线分布的节点表示断层的三维结构。模型沿断裂带共有4条等深线,每条等深线上有21个节点,深度依次为0.1、6、15、22 km,各深度的初始闭锁系数依次设置为0.8、0.6、0.4、0。根据1.2节柴北缘断裂带的划分,参考已有研究并通过迭代寻优得到最优模型,其中断层倾向设置如下:赛什腾-绿梁山断裂段倾向NE、倾角70°,锡铁山-阿木尼克山-牦牛山断裂段倾向NE、倾角55°,北鄂拉山断裂段和南鄂拉山断裂段倾向均为SW、倾角60°[11,38-39]。断层节点φ值约束如下:节点均为独立节点,闭锁系数从地表向下单调递减,22 km深度处的节点设置为自由滑动(φ=0)。

反演过程中,采用网格搜索和模拟退火方法计算各节点的闭锁系数及各块体的欧拉矢量;根据双线性插值方法获得相邻节点之间断层网格(网格大小1 km×1 km)区域的闭锁程度;将闭锁系数与根据欧拉矢量计算得到的断层滑动速率相乘,从而获得滑动亏损速率的空间分布。进一步地,使用蒙特卡洛(Monte Carlo)方法从统计学的角度估计了模型反演结果的不确定度。对GPS数据加入以观测误差作为标准差生成的随机噪声,得到100组带有扰动的速度场数据。基于扰动速度场数据重复反演,统计分析结果,得到反演结果的不确定度。

图6给出了GPS速度拟合残差结果。可以看出:拟合残差方向具有随机性,约83%的GPS站点的残差小于1.5 mm/a,且残差在东西、南北方向上的分量均符合正态分布,表明模型未遗漏明显信号,拟合效果较好[46-47]。残差值较大的站点主要集中在:①柴达木块体北边界以南的近场区域,可能是由于断裂带倾角沿走向有一定的变化,断层近场区域应力应变复杂从而产生局部运动分量[48],而本文采用的断层几何不足以刻画所有局部形变特征;②东昆仑断裂的近场区域,可能是由于模型反演中未考虑东昆仑断裂的闭锁效应对地表形变的影响[6-7]。

图6 GPS数据拟合残差Fig.6 Fitting residual of GPS data

3 反演结果

图7给出了柴北缘断裂带的闭锁程度分布。结果表明:赛什腾-绿梁山断裂段中部处于强闭锁状态(15 km深度内平均闭锁系数约0.7),闭锁深度达到20 km,从中部向两端闭锁程度逐渐减弱。锡铁山-阿木尼克山-牦牛山断裂段从西北向东南闭锁深度由18 km增大到20 km;西北端0～6 km深度平均闭锁系数为0.9,6～15 km深度平均闭锁系数减小为0.6,15～22 km深度基本为蠕滑状态;中部和东南端0～15 km深度为强闭锁状态(平均闭锁系数约0.99),从15～22 km深度闭锁程度逐渐降低,由强闭锁转变为蠕滑。北鄂拉山断裂段北端闭锁程度较弱,从5～22 km深度基本处于蠕滑状态;中部和南端处于强闭锁状态,闭锁深度达到19 km,0～15 km深度平均闭锁系数分别约0.99和0.8。南鄂拉山断裂段北端和南端处于强闭锁状态,闭锁深度达到20 km,0～15 km深度平均闭锁系数分别约0.9和0.99;中部闭锁程度较弱,从5～22 km深度基本处于蠕滑状态。本文得到的柴北缘断裂带闭锁深度与文献[31—32]通过GPS数据确定的青藏高原东北缘断层平均闭锁深度17～22 km较为吻合,基本在文献[49—50]得到的区域震源深度18±8 km范围内。

图7 柴北缘断裂带闭锁程度Fig.7 Locking degree of the fault zone in northern Qaidam Basin

图8给出了滑动亏损分布。结果显示:赛什腾-绿梁山断裂段的滑动亏损较低,滑动亏损从中部(0～15 km深度的平均大小约为1.26 mm/a)向两端逐渐减小,滑动亏损积累较少。锡铁山-阿木尼克山-牦牛山断裂段的滑动亏损从西北向东南逐渐增大,西北端从地表到22 km深度滑动亏损逐渐减小,0～15 km深度的平均大小约为2.30 mm/a,15～22 km深度的平均大小约为1.08 mm/a;中部和东南端0～15 km深度的滑动亏损较西北端有所提高,平均大小分别约为3.34 mm/a和3.80 mm/a,15～22 km深度滑动亏损较小。北鄂拉山断裂段最北端0～5 km深度平均滑动亏损大小约为3.12 mm/a,5～22 km深度滑动亏损逐渐减小,中部15 km深度内平均滑动亏损大小约为2.23 mm/a, 其他部分平均滑动亏损较小。南鄂拉山断裂段北端在地下15 km深度平均滑动亏损大小约为2.95 mm/a,15 km深度以下滑动亏损逐渐减小,20～22 km深度处几乎无滑动亏损积累,其他区域滑动亏损较小。

图8 柴北缘断裂带滑动亏损Fig.8 Slip deficit of the fault zone in northern Qaidam Basin

表1给出了基于较优模型得到的阿克塞-德令哈块体、柴达木块体、鄂拉山块体的刚性旋转参数、内部应变参数。结果表明:3个块体的运动趋势相一致,均为顺时针旋转;鄂拉山块体的内部应变表现为北东-南西方向的挤压和北西-南东方向的拉张,主压应变率为22.49×10-9/a、主张应变率为12.70×10-9/a;柴达木块体、阿克塞-德令哈块体的内部应变表现为北北东-南南西方向的挤压和北西西-南东东方向的拉张,柴达木块体的主压应变率为15.08×10-9/a、主张应变率为10.04×10-9/a,阿克塞-德令哈块体的主压应变率为12.71×10-9/a、主张应变率为7.53×10-9/a。

表1 欧拉矢量与块体内部应变率反演结果

本文使用的负位错模型假设块体运动由块体刚性旋转、断层闭锁和块体内部永久应变3部分组成。表1结果显示,3个块体均呈顺时针方向旋转,块体内部呈北北东-南南西或北东-南西方向的挤压应变特征,揭示了在印度板块与欧亚板块汇聚作用和鄂尔多斯、阿拉善刚性块体阻挡的影响下,区域发生了顺时针方向的旋转、北北东-南南西或北东-南西方向的缩短,区域断裂共同调节并吸收了相邻块体间的差异运动,区域变形具有局部旋转和整体弥散变形的特征[45,51-52]。但由于模型假设块体内部应变为均匀应变,未顾及其非均匀性,反演结果可能受相关因素的影响[53]。

4 讨 论

4.1 滑动速率特征对比分析

将模型反演结果与GPS剖面结果进行对比。模型反演结果表明,锡铁山-阿木尼克山-牦牛山断裂段自西北向东南、北鄂拉山断裂段自北向南均呈右旋走滑兼挤压运动特征,锡铁山-阿木尼克山-牦牛山断裂段右旋走滑、挤压速率分别约2.67～3.13 mm/a、0.70～2.32 mm/a,北鄂拉山断裂段右旋走滑、挤压速率分别约1.27～3.91 mm/a、0.35～0.67 mm/a。GPS剖面结果表明,锡铁山-阿木尼克山-牦牛山断裂段自西北向东南由右旋挤压转为左旋挤压,其中,右旋走滑速率约为1.2 mm/a,左旋走滑速率约为1～1.5 mm/a;北鄂拉山断裂段自北向南由右旋拉张转为右旋挤压,其中,拉张速率约为0.5 mm/a,挤压速率约为1 mm/a。二者结果存在差异,其原因可能为:基于GPS剖面确定断层滑动速率时,通常假定断层倾角垂直且地表形变仅受单一断层影响;而基于负位错模型反演时能顾及较为真实的断层几何及多断层运动对地表形变的影响,从而得到较为准确的反演结果。

将本文与先前研究结果进行对比(表2)。赛什腾-绿梁山断裂段的右旋走滑速率为1.55～1.96 mm/a,介于文献[31—32]给出的1～2.9 mm/a之间;挤压速率为1.59～2.50 mm/a,显著小于文献[31]给出的7.3 mm/a,且与文献[32]给出的拉张速率1 mm/a有较大差异。锡铁山-阿木尼克山-牦牛山断裂段的右旋走滑速率为2.67～3.13 mm/a,略大于1.81～2.1 mm/a[12-13];挤压速率为0.70～2.32 mm/a,大于0.19～0.51 mm/a[12-13,18]。北鄂拉山断裂段的右旋走滑速率为1.27～3.91 mm/a,南鄂拉山断裂段的右旋走滑速率为1.48～2.24 mm/a,介于前人给出的鄂拉山断裂全段右旋走滑速率0.6～5 mm/a[4,8,15,17,29-32]之间。北鄂拉山断裂段的挤压速率为0.35～0.67 mm/a,介于前人给出的0～1.1 mm/a[15-16,31-32]之间。南鄂拉山断裂段由挤压速率0.27～0.71 mm/a转为拉张速率0.33 mm/a,该特征与前人给出的全段挤压速率0.15～1.7 mm/a存在差异[15-16,31-32]。本文研究结果与先前研究结果整体上一致,但存在局部差异,其原因可能为:①研究区域、断层几何及块体模型不同;②采用的数据源不同;③反演存在不唯一性。

表2 不同研究给出的柴北缘断裂带滑动速率

4.2 断层闭锁分辨率试验

为了研究GPS数据对闭锁程度反演结果的分辨能力,本文进行了分辨率试验。该方法主要分为正演和反演两步:①正演:给定每个节点处的闭锁系数φ,正演得到各GPS站点的模拟速度场;②反演:将各GPS站点的中误差作为白噪声加入到模拟速度场中,进而反演断层闭锁程度[46-47]。与正演模型(图9(a)、(b))相比,由反演结果(图9(c)、(d))可知:①在赛什腾-绿梁山断裂段,反演结果较好地识别出高闭锁区域集中在中部0～15 km深度范围的特征;②在锡铁山-阿木尼克山-牦牛山断裂段,反演结果识别出了中部的高闭锁区域,但分布相对分散,在东南端闭锁程度略高于正演模型;③在北鄂拉山断裂段,反演结果较好地刻画了中部的高闭锁区域;④在南鄂拉山断裂段,反演结果显示高闭锁区域集中在南端0～5 km、中部0～15 km深度范围内,与正演模型设置的集中在中部0～15 km深度范围略微存在差异。综上,本文使用的GPS数据能较好约束块体模型反演,能够有效识别出断层闭锁的分布特征。

图9 柴北缘断裂带闭锁程度分辨率试验Fig.9 Resolution test of the locking degree of the fault zone in northern Qaidam Basin

4.3 区域地震危险性分析

断层现今闭锁程度与滑动亏损、地震危险性密切相关[1-2,54-56],如2008年汶川地震、2010年玉树地震、2013年芦山地震、2015年廓尔喀地震等发震断层在震间阶段的闭锁系数接近于1[42,56-57]。基于中国地震台网中心(https:∥news.ceic.ac.cn/)地震目录,本文给出了该区域1970—2016年的历史地震分布(图10)。根据历史地震统计,柴北缘断裂带近场共发生M≥5级地震104次,其中M≥6级地震17次。综合反演结果与历史地震分布分析,可以看出:

图10 研究区域历史地震分布Fig.10 Historical earthquakes distribution in the study area

赛什腾-绿梁山断裂段闭锁程度从中部向两端逐渐减小,全段滑动亏损较小,应变难以积累。锡铁山-阿木尼克山-牦牛山断裂段闭锁程度在中部和东南端较高,闭锁深度达到20 km,应变在此积累;中部和东南端现今小震活动相对西北端较弱,应变释放程度低,地震危险性值得关注[56,58]。北鄂拉山断裂段闭锁程度在中部和南端较高,闭锁深度达19 km,滑动亏损在中部较大,应变在此积累,且中部现今小震活动不活跃,为历史地震的未破裂段,应变释放较弱,地震危险性值得关注。南鄂拉山断裂段闭锁程度在北端和南端较高,闭锁深度达20 km,滑动亏损在北端较大,应变在此积累,且北端是历史地震的未破裂段,现今小震活动活跃程度较低,释放的应变较少,地震危险性值得关注。

综上,在闭锁程度较高、滑动亏损较大、小震活动较弱的锡铁山-阿木尼克山-牦牛山断裂段中部和东南端、北鄂拉山断裂段中部、南鄂拉山断裂段北端应变积累且释放较弱,其地震危险性可能较高,须加以关注。

5 结 论

本文利用融合得到的较密集的柴北缘GPS数据,基于GPS速度剖面法初步分析断层运动特征,进一步通过负位错模型定量研究了断裂带的滑动速率、闭锁程度及滑动亏损等,结合历史地震对区域地震危险性进行分析,主要结论如下:

(1) 滑动速率反演结果表明,赛什腾-绿梁山断裂段、锡铁山-阿木尼克山-牦牛山断裂段、北鄂拉山断裂段呈右旋走滑兼挤压运动特征,南鄂拉山断裂段的运动特征自北向南由右旋挤压转为右旋拉张。

(2) 闭锁程度反演结果显示,锡铁山-阿木尼克山-牦牛山断裂段中部和东南端、北鄂拉山断裂段中部和南鄂拉山断裂段南端闭锁显著,15 km深度内的平均闭锁系数约0.99。赛什腾-绿梁山断裂段中部、锡铁山-阿木尼克山-牦牛山断裂段西北端、北鄂拉山断裂段南端、南鄂拉山断裂段北端的闭锁较强,15 km深度内的平均闭锁系数约为0.7～0.9。

(3) 滑动亏损反演结果显示,赛什腾-绿梁山断裂段中部的滑动亏损相对较大,15 km深度内的平均大小约为1.26 mm/a,向两端逐渐减小。锡铁山-阿木尼克山-牦牛山断裂段的滑动亏损自西北向东南均逐渐增大,中部和东南端较大,15 km深度内平均大小约为3.34 mm/a和3.80 mm/a。北鄂拉山断裂段中部滑动亏损较大,15 km深度内平均大小约为2.23 mm/a。南鄂拉山断裂段北端滑动亏损较大,15 km深度内平均大小约为2.95 mm/a。

(4) 闭锁程度、滑动亏损反演结果及区域历史地震分布表明,锡铁山-阿木尼克山-牦牛山断裂段中部和东南端、北鄂拉山断裂段中部、南鄂拉山断裂段北端具有闭锁程度较高、滑动亏损较大、小震活动相对较弱的特点,地震危险性可能较高,是值得关注的区域。