青藏高原东北缘岩石圈有效弹性厚度与地震分布关系

苏子旺,鲁宝亮,2,3*,李柏森,朱武

1 长安大学地质工程与测绘学院,西安 710054 2 海洋油气勘探国家工程研究中心,北京 100028 3 西部矿产资源与地质工程教育部重点实验室,西安 710054

0 引言

青藏高原东北缘与阿拉善地块和鄂尔多斯地块交汇区是现今高原向外扩展的最前缘,印度和欧亚板块碰撞的远程响应,造就了该区域复杂的地形地貌,独特的地质构造,导致了该区域大型断裂发育、中强地震频发(邓起东等,2002).青藏高原东北缘按照活动地块可划分为西南的青藏地块区、西北的阿拉善地块和东部的鄂尔多斯地块(张培震等,2003).受印度板块运动影响,研究区(30°N—42°N,95°E—112°E,图1)中存在大量深大断裂带,如海原断裂、昆仑断裂、鲜水河断裂等(Tapponnier et al.,2001),青藏高原东北缘由昆仑断裂、鲜水河断裂分为祁连地块、柴达木地块、巴颜喀拉地块、羌塘地块,由于地块边界的断裂带切割地壳深度大,差异运动强烈而非连续性更强,有利于应力积累而孕育大地震; 鄂尔多斯块体内部构造活动性微弱,不发育大规模的活动断层,周边的地震活动却十分强烈(张培震等,2003).由于大型活动断裂彼此的共同作用,不同活动地块运动方式和速度的不同,新构造运动强烈,该区域地震活动频发,具有空间分布广、震源深度浅、成因类型杂、灾害损失重等的特点(张培震等,2013).

图1 青藏高原东北缘地形、地震和地质构造

岩石圈有效弹性厚度Te反映了岩石圈在地质时间尺度(>105年)荷载下的力学强度(Watts and Burov,2003),定义为与岩石圈板块中实际应力分布所产生的弯矩相等的理论弯曲弹性薄板的厚度(Forsyth,1985),与地幔物质组成、热状态、壳幔耦合程度等相关(Burov and Diament,1995).较高的Te值表示岩石圈强度较大,较低Te值的地区岩石圈强度较小,岩石圈的强度差异控制着岩石圈在构造力作用下的局部变形,与地震分布有一定的对应关系(Lowry and Smith,1995).Te变化剧烈的区域即Te梯度带有利于应力释放和地震孕育; 同时,地震发生频率较高的地区通常有较低的Te值(10～40 km),频率较低的地区Te值较高(大于80 km)(Chen et al.,2013).因此,选择Te为分析目标,研究其与地震分布的耦合关系,对于认识地震的发震条件具有重要意义.一般利用地形数据、重力数据,采用谱分析法进行Te估计.在谱分析法中需要计算地形和重力异常数据的自功率谱和互功率谱.为了准确地计算谱,学者提出了周期图法(Percival and Walden,1993)、最大熵法(Lowry and Smith,1994)、多窗谱法(Thomson,1982; Simons et al.,2000)、Fan小波谱法(Kirby and Swain,2004)等方法.其中,Fan小波谱法不需要加窗,因此不受窗口大小影响,故采用Fan小波谱法.

本文使用布格重力异常和地形数据,利用基于Fan小波谱法的导纳和相关联合反演,计算青藏高原东北缘的岩石圈有效弹性厚度,分析其分布特征与地质意义; 研究天然地震与岩石圈力学性质间的耦合关系,进一步认识该区域地震的孕育环境,为防震减灾提供科学依据.

1 研究方法

假定岩石圈由均一密度为ρc的地壳组成,可以近似为理想的水平无限延伸的薄弹性板,上覆于软流圈流体上.对于一个漂浮于流体之上的弹性薄板或梁,一维挠曲响应方程一般用四阶偏微分方程表示(Gunn,1947):

(1)

其中,v(x)为弹性薄板的变形,是水平距离x的函数;ρm和ρc分别为下伏流体层的密度(设为地幔密度)和挠曲凹地的填充物质密度(设为地壳密度);g为重力加速度;D为弹性板的挠曲刚度,与有效弹性厚度Te的关系为(Walcott,1970):

(2)

其中,Te为岩石圈有效弹性厚度;E为杨氏模量,取1×1011N·m-2;σ为泊松比,取0.25.

把地球岩石圈视为滤波器,输入为地形数据,经过岩石圈的滤波后,输出为重力异常,定义重力导纳来描述这种滤波作用(Dorman and Lewis,1970):

(3)

为减少实测地形和重力异常数据中误差的影响,一般采用统计的方法(例如分波段平均)计算实测导纳:

(4)

假设地表荷载和地下荷载不相关(即具有随机相位差),经过波数段统计平均,交叉项将被消除,实测导纳展开获得模型导纳:

(5)

其中,HT为地表荷载对地表地形的贡献,HB为地下荷载对地表地形的贡献,ΔgT为地表荷载对重力异常的贡献,ΔgB为地下荷载对重力异常的贡献.

(6)

做相同假设,实测相关展开得到模型相关(Forsyth,1985):

(7)

如图2所示,假定地下荷载加载和均衡面均为Moho面,利用实测地形和重力异常获得地表初始荷载Hi和地下初始荷载Wi,进而得到岩石圈对地表初始荷载的挠曲响应WT、岩石圈挠曲后的地表初始荷载剩余分量HT、岩石圈对地下初始荷载的挠曲响应HB、岩石圈挠曲后的地下初始荷载剩余分量WB,即可计算预测导纳和预测相关.

图2 岩石圈初始荷载与挠曲响应(Watts,2001)

四个变形分量(HT,HB,WT,WB)和地表初始荷载Hi、地下初始荷载Wi在波数域下的关系为(Kirby and Swain,2011):

(8)

其中,φ=Dk4/g+ρm-ρf,Δρ1=ρc-ρf,Δρ2=ρm-ρc;ρf是上覆流体的密度(在海洋中ρf等于水的密度,在陆地上ρf=0).实测地形H和Moho面起伏W是上述分量的总和:H=HT+HB,W=WT+WB.

地表地形谱H和Moho面地形谱W与二者引起的重力异常Δg之间的一级近似关系式(Parker,1973)为:Δg=2πAGΔρ1H+2πGΔρ2e-|k|zmW,其中,A=0则Δg为布格重力异常,A=1则Δg为自由空气重力异常;G是万有引力常数;zm为地壳厚度.利用地形分量和初始荷载的关系式、重力异常谱与地形谱的一级近似式,可得地表荷载和地下荷载产生的重力异常分量(ΔgT和ΔgB)与初始荷载Hi和Wi的关系为

(9)

可以解出地表和地下初始荷载分别为

(10)

进而计算得到预测导纳:

(11)

(12)

使用“均匀f”法(Kirby and Swain,2008),f不随波数变化,全波数下均一.故给定两个独立参数Te和F,即可得到预测的导纳和相关.

为了反演计算Te和F,采用导纳和相关联合反演,使用实测导纳和相关与预测导纳和相关之间的简化卡方准则(reduced chi-squared criterion)来衡量两者的拟合差(Audet,2014),使用非线性最小二乘法搜索满足卡方最小的Te和F,即为最优解.导纳和相关联合反演Te和F的简化卡方准则为

(13)

其中j表示导纳(j=1)和相关(j=2);dij表示方差为εij的实测值,sij表示预测值;N为采样网格的总数量.

2 数据及计算结果

2.1 数据

岩石圈有效弹性厚度的计算需要地形、重力异常、地壳厚度、地壳密度与岩石圈地幔密度.布格重力异常数据来自WGM2012(图3a,Bonvalot et al.,2012),地形数据来自ETOPO1(图1,Amante and Eakins,2009),地壳厚度(图3b)和地壳密度(图3c—e)数据来自CRUST 1.0(Laske et al.,2013),数据分辨率依次为2′×2′、1′×1′、1°×1°,岩石圈地幔密度取地幔平均密度3.2g/cm3.为避免边界效应,实际处理的数据范围大于研究区域.岩石圈被假设为均一密度地壳组成的弹性板,此均一密度取为地壳平均密度(图3f),由上、中、下地壳密度对各层厚度加权平均得到.对数据以107°E为中心子午线做Mercator投影变换,并以10 km为间隔进行网格化,对Te和荷载比F进行逐点反演计算.反演使用PlateFlex(Audet,2019)进行.

图3 研究区布格重力异常、地壳厚度及地壳密度数据

2.2 计算结果

图4 使用不同中心波数得到的反演结果

在中心波数|k0|取5.336的结果中,选择不同构造单元地区,给出曲线拟合状况(图5),在岩石圈Te及F反演过程中,导纳和相关的实测值与预测值总体上吻合程度较好,表明本文利用导纳和相关函数联合反演获得的结果较为可靠.

图5 不同构造单元地区的导纳和相关曲线拟合

3 讨论

3.1 有效弹性厚度特征

研究区内的岩石圈有效弹性厚度Te整体呈明显的东高西低分布,西南部存在局部高值区,中部为Te值高低变化的分界区域,是Te值自东向西变低的过渡带(图6).

图6 青藏高原东北缘Te

东部鄂尔多斯地块的Te分布大体走向为北北西,Te值普遍大于34 km,地块西北部与东南部Te值更是达到40 km以上.前人使用多窗谱法得到的Te值高达110 km(Chen et al.,2013)和80 km(Chen et al.,2022),与本文差别较大,可能是使用方法不同导致的; 而使用Fan小波谱法,郑勇等(2012)发现鄂尔多斯盆地Te较高(30～65 km),南北两端较高、中南略低,Lu等(2020)的Te值大于40 km,均与本文结果相近,形态相似.在鄂尔多斯地块西缘,较强的鄂尔多斯地块插入了较弱的青藏块体,形成了六盘山区域(Tian et al.,2021).秦岭可能存在物质向东流动的通道(Royden et al.,2008),该区域存在Te高值区(>30 km),所以若存在此通道,应该位于岩石圈底部或软流圈深度上(郑勇等,2012).总体而言,本文的计算结果表明鄂尔多斯地块岩石圈力学强度较大,印证了鄂尔多斯地块具有内部变形小、相对稳定的特征(张先康等,2003),对应了鄂尔多斯地块内部构造活动性微弱的性质(张培震等,2003).鄂尔多斯地块存在有相对冷的岩石圈地幔(Guo et al.,2016)也支撑了这一特征.

西部的青藏块体Te呈现8～48 km的剧烈变化,跨度包含了整个研究区域Te的变化范围,块体东缘Te值较低,小于25 km,块体内部如巴颜喀拉地块和羌塘地块同时存在Te的高值区(>40 km)和低值区(<20 km),呈现出岩石圈力学强度强烈的横向变化特征.巴颜喀拉地块岩石圈强度整体较弱,Te数值与李永东等(2013)接近.昆仑断裂的Te低值带可能是较厚的青藏高原地壳向北插入或局部熔融侵入柴达木盆地南部地壳之下(Karplus et al.,2011; Le Pape et al.,2012)的均衡信号(李永东等,2013).巴颜喀拉地块东部与龙门山断裂带西部Te值较低,与李永东等(2013)、杨光亮等(2020)、胡敏章等(2020)结果相近.鲜水河断裂具有明显的低电阻结构,中下地壳存在熔融物质(何梅兴等,2017),岩石圈强度较弱,Te值较低.与东部稳定的鄂尔多斯地块和扬子地台接触区域岩石圈强度明显较弱,前人使用接收函数(王椿镛等,2008; 叶卓等,2018)、Lg波Q值成像(Zhao et al.,2013)等方法发现该区域下地壳物质处于热和软弱状态,可能发生有下地壳流,受四川盆地影响,物质流动的通道可能是向东和东南的.岩石圈强度的剧烈变化与印度板块向欧亚板块的深部俯冲-碰撞-拼贴过程、下地壳存在的局部熔融(李宝春等,2020)以及地壳物质在印度板块推挤和周边刚性地块的阻挡下围绕东构造结发生顺时针旋转并向东挤出逃逸(张培震等,2002; 甘卫军等,2004)有密切关系,形成了青藏高原东缘及东北缘东向与东南向的Te低值带与梯度带.

研究区的荷载比F分布同样东西差异明显,但与Te不同,荷载比F整体是西南角为明显高值,其他区域为低值(图7).青藏块体内部的巴颜喀拉地块内部和鲜水河断裂带存在荷载比F大于0.5的区域、岩石圈受到荷载以地下为主,岩石圈深部物质密度存在较大的横向变化,主导了岩石圈的挠曲变形,可能与青藏高原东部物质向东流动有关.研究区的其他地区均低于0.2,且其中约一半地区的荷载比F小于0.1,整体与Lu等(2020)的结果相近,说明以地表荷载为主,地表物质密度横向变化较大,是该区域岩石圈挠曲变形的主要影响因素.六盘山荷载比F较小,与佘雅文等(2016)结果相近.龙门山断裂西北部的荷载比F较小(约0.15),与杨光亮等(2020)的结果(约0.4)存在差别,但均表明该地区荷载以地表为主.青藏块体与周边块体交汇处的Te与荷载比F的复杂变化反映了该区域复杂的地质构造应力作用.

图7 青藏高原东北缘荷载比F

3.2 与地震分布关系

地震实际上是在区域构造应力作用下,应变在活动断裂带上不断积累并达到极限状态后而突发失稳破裂的结果(Scholz,1998).即地震的发生主要受到内部应力与活动断裂的影响.Te的空间变化反映了岩石圈内部应力的差异,所以弹性厚度的高分辨率地图可能是分析新构造和地震成因的有用工具(Tassara et al.,2007).中国大陆西部由于活动断裂和褶皱广泛发育,地震强度和频度都大于东部(张国民和张培震,2000).所以对于青藏高原东北缘,研究Te、荷载比F与地震分布之间的对应关系具有重要意义.

研究区内60%以上的地震震源深度在10～20 km之间(图8),其所表征的多震层深度与Te关联性较弱,震源深度与Te基本没有明显的规律可循(陈波,2013); 又因不同震级的地震成因各异,较难展开详细分析,所以后续讨论主要针对震级、发生占比与Te、荷载比F及二者梯度的关系.为减少占比较大的5级以下地震对结论的影响,添加了各参数与震级的散点图.将地震震级与Te、荷载比F、Te梯度、荷载比F梯度叠加(图9—12),得到如下结果:

图8 青藏高原东北缘地震震源深度

图9 青藏高原东北缘Te与地震分布

图10 青藏高原东北缘荷载比F与地震分布

图11 青藏高原东北缘Te梯度与地震分布

图12 青藏高原东北缘荷载比F梯度与地震分布

在鄂尔多斯地块,地震主要出现在地块边缘.地块西缘、北缘与东南缘地震分布较多,以6级以下地震为主,该部分岩石圈Te值较高(>30 km),地块西北缘Te值甚至达到了40 km以上,荷载比F和荷载比F梯度均较小,Te梯度相对较大,呈现为明显的梯度带,地震分布与Te梯度带耦合度相对较高,但不明显.地块内部地震分布较少,集中在地块东北部,且绝大多数为5级以下地震,发震区域Te虽大于32 km,但与地块内部其他地区相比为相对低值区,Te梯度、荷载比F、荷载比F梯度很小,地震分布与Te相对低值区耦合度较高,说明地块内部的地震易于发生在Te相对低值区.

研究区内的青藏块体包含祁连地块、柴达木地块、巴颜喀拉地块和羌塘地块,该区域地震活动频发,地震成因复杂,震级跨度大,6级以下地震在整个区域的分布没有明显倾向,6级以上地震则相对集中在块体内部.在块体边缘区域,Te变化明显,但基本都大于20 km;Te梯度普遍较小,块体东缘存在明显梯度带,与海原断裂、昆仑断裂相交; 荷载比F均小于0.5,岩石圈荷载以地表为主; 昆仑断裂与龙门山断裂之间存在荷载比F梯度带; 均与地震分布耦合度较差,说明块体边缘的地震分布受四个参数影响较小.块体内部的地震分布较为发散,块体西部地震稍多于块体东部; 强度较小(Te为20～36 km)的祁连地块西部和柴达木地块5级以上地震较多,小震较少,和地震分布无明显对应关系; 柴达木地块东侧地区Te值明显变低(<20 km),5级以上地震占比减少,地震与Te梯度、荷载比F梯度有一定的对应性,说明在岩石圈弹性很弱的地区,地震易发生在Te梯度、荷载比F梯度相对较大区域; 巴颜喀拉地块边缘小震居多,地块内部6级以上强震增加,与Te梯度、荷载比F梯度对应性较好; 羌塘地块内部6级以下地震居多,与Te梯度、荷载比F梯度有一定对应性,地块东部岩石圈强度较差区域存在多个6级以上强震,但其与Te梯度、荷载比F梯度并无明显对应性.

阿拉善地块较为稳定,所受荷载以地表为主.地块内部地震活动少,地块边缘地震活动相对内部较多,基本都是6级以下地震,该区域地震的分布与Te、荷载比F及二者梯度耦合性较差.

走滑断裂中,昆仑断裂与鲜水河断裂的Te和荷载比F变化大,昆仑断裂和鲜水河断裂分别与荷载比F梯度、Te梯度有一定的对应性,但对应性一般; 海原断裂西部的Te梯度与地震对应性较好,其他区域较差.位于块体边缘的阿尔金断裂和逆冲断裂带龙门山断裂处的Te和荷载比F及其梯度与地震分布则无明显对应性.故走滑断裂带和逆冲断裂带的Te、荷载比F及其梯度与地震分布并无明显关系,对应性总体较差.

图13清晰地表现了岩石圈有效弹性厚度、荷载比、梯度、梯度与地震分布的关系,其中,考虑到面积占比的影响,将地震发生占比与不同参数的区间面积占比之比作为表征地震发生占比与不同参数对应关系的指数(后称指数),绘制图13c—l,对于同一参数,某一区间的指数大说明地震易发生在该区间.从Te与地震发生占比(图13a)看,研究区内1.62%的地震发生在10～14 km的Te区域,7.56%在14～18 kmTe区域,12.39%在18～22 kmTe区域,20.09%在22～26 kmTe区域,33.96%在26～30 kmTe区域,13.29%在30～34 kmTe区域,3.31%在34～38 kmTe区域,4.21%在38～42 kmTe区域,3.35%在42～46 kmTe区域,0.22%在Te46 km以上区域; 75.62%的地震发生在Te小于30 km的区域.结合Te与震级分布(图13b),不同震级地震的发生均偏向于Te值小于30 km的区域,4级以下地震与Te对应性较弱,4级以上地震与Te对应性较强,6级以上强震发生在Te较低区域的倾向更加明显.在地震发生占比/Te区间面积占比(图13c)中,Te值大于26 km的指数明显低于Te值小于26 km.可以得出4级以上地震均易于发生在Te值较小即岩石圈综合力学强度较差区域的结论.

图13 Te、荷载比F、Te梯度、荷载比F梯度与地震震级、发生占比关系

从F与地震发生占比(图13d)看,研究区内76.41%的地震发生在0.2以下的荷载比F区域,7.09%在0.2～0.4荷载比F区域,4.61%在0.4～0.6荷载比F区域,5.19%在0.6～0.8荷载比F区域,看似地震易于发生在荷载比F小于0.2的区域.但结合荷载比F与震级分布(图13e)发现,6级以下地震占比巨大且在0～0.2荷载比F范围内集中分布,在0.2～1荷载比F范围内分布很少且较为均匀; 6级以上地震占比很小且分布存在明显两极分化现象,分别在荷载比F值约0.1、约0.8处集中,说明在极不均衡的荷载作用下,岩石圈容易产生较大的挠曲形变,易于发生强度较大的地震.在地震发生占比/F区间面积占比(图13f)中,荷载比F值较大和较小区间的指数略大于荷载比F值居中的指数,但该对应性不显著,且荷载比F值小于0.1的区间的指数还存在明显降低.表明6级以下地震倾向于发生在地表荷载为主的区域,6级以上强震则倾向于发生在地表、地下荷载不均衡的区域.

从Te梯度与地震发生占比(图13g)看,研究区内71.95%的地震发生在0.0～0.05的Te梯度区域,14.04%在0.05～0.1Te梯度区域,8.35%在0.1～0.15Te梯度区域,2.52%在0.15～0.2Te梯度区域,1.84%在0.2～0.25Te梯度区域,1.22%在0.25～0.3Te梯度区域,0.07%在0.3～0.35Te梯度区域.Te梯度与震级分布(图13h)同样表明不同震级地震的分布均集中于Te梯度较低的区域.但在去除面积占比的影响后(图13i),Te梯度较大区域的指数整体大于较小区域.表明地震易于发生在Te梯度较大的区域.

从F梯度与地震发生占比(图13j)看,研究区内84.91%的地震发生在0.0～0.001的F梯梯度区域,6.66%在0.001～0.002F梯度区域,3.82%在0.002～0.003F梯度区域,3.31%在0.003～0.004F梯度区域,1.04%在0.004～0.005F梯度区域,0.18%在0.005～0.006F梯度区域,0.07%在0.006～0.007F梯度区域.结合F梯度与震级分布(图13k),不同震级地震的分布均集中于F梯度较低的区域.同样地,在去除面积占比的影响后(图13l)可以看出指数随F梯度增加而增加.说明地震易于发生在F梯度较大的区域.

总体上可以认为地震更可能发生在岩石圈Te较低、Te变化剧烈的地区; 震级较小的地震倾向于发生在以地表荷载为主的地区,6级以上地震倾向于发生在荷载以地下为主、荷载比F变化剧烈的地区.

4 结论与认识

本文使用WGM2012布格重力异常数据、ETOPO1地形数据和CRUST1.0地壳模型数据,利用导纳和相关联合反演了青藏高原东北缘的岩石圈有效弹性厚度Te和荷载比F.通过对比分析,研究了青藏高原东北缘岩石圈Te、F及二者梯度与地震分布的耦合关系,得到如下结论和认识:

(1) 研究区Te整体呈明显的东高西低分布,青藏块体内部Te变化剧烈,西部高值(>40 km)和东部低值(<20 km)共存; 鄂尔多斯地块Te值较高(>30 km),变化相对平缓; 荷载比F存在局部西南高,巴颜喀拉和羌塘地块荷载比F较高(>0.5),说明以地下荷载为主,其他地块荷载比F较低(<0.2),说明以地表荷载为主.

(2) 鄂尔多斯地块结构稳定,岩石圈强度较大,Te值较高; 内部构造活动性微弱,深部物质密度横向变化较小,岩石圈所受荷载以地表为主,荷载比F值较低.柴达木地块东部与巴颜喀拉地块东部中下地壳可能发生有下地壳流,岩石圈强度较小,Te值较低; 巴颜喀拉地块内部与鲜水河断裂的荷载比F值较高,岩石圈所受荷载以地下为主,可能是地壳物质向东流动导致岩石圈深部物质密度横向变化较大引起的.

(3) 岩石圈高Te值的地区,地震出现的频次较少;Te值较低、Te梯度较大的区域,地震出现的频次较高.荷载比F值较小的地区,震级较小的地震出现频次高,6级以上强震则较多发生在荷载比F约为0.1和0.8的区域.荷载比F梯度较大的区域易于发生地震尤其是强震.

致谢感谢中国地震台网中心国家地震科学数据中心(http:∥data.earthquake.cn/)提供数据支撑.感谢Audet(2019)提供的PlateFlex程序对本文计算部分的帮助.感谢论文匿名评审专家和编辑部提出的宝贵修改意见.