地球物理学报  2017, Vol. 60 Issue (6): 2480-2492   PDF    
重力与地形数据揭示的巴颜喀拉块体东缘垂向构造应力场
佘雅文1,3, 付广裕2 , 王灼华2, 高原1,2     
1. 中国地震局地球物理研究所, 北京 100081;
2. 中国地震局地震预测研究所, 北京 100036;
3. 河北省地震局, 石家庄 050021
摘要: 利用基于消去-恢复原理的最小二乘配置方法,对2009—2013年相对重力/GPS联合观测数据与EGM2008模型数据进行融合,更新了巴颜喀拉块体东缘地区的自由空气与布格重力异常场.基于该布格重力异常数据,以CRUST1.0地壳密度模型为初始条件,使用二维多边形棱柱体正演与非线性最小二乘反演方法,获取了巴颜喀拉块体东缘地壳分层密度结构.基于地壳不可压缩和均衡调整原理提出了计算垂向构造应力新方法,并结合上述地壳分层密度结构和地形数据计算了巴颜喀拉块体东缘垂向构造应力分布.结果表明,龙门山断裂带中南段蓄积了较高的正向构造应力(约40 MPa),马尔康周边地区蓄积了较高的负向构造应力(约-30 MPa).对研究区域1970年以来5级(MS)以上地震进行统计发现,地震多发生在垂向构造应力梯度带上,垂向构造应力为正的地区易触发浅源地震,为负的地区易触发深源地震.在地壳横向变形强烈的区域,垂向构造应力与地震深度的对应关系减弱.
关键词: 巴颜喀拉块体      布格重力异常      地壳密度结构      垂向构造应力     
Vertical tectonic stress in eastern margin of Bayan Har block revealed by gravity and terrain data
SHE Ya-Wen1,3, FU Guang-Yu2, WANG Zhuo-Hua2, GAO Yuan1,2     
1. Institute of Geophysics, China Earthquake Administration, Beijing 100081, China;
2. Institute of Earthquake Science, China Earthquake Administration, Beijing 100036, China;
3. Hebei Seismological Bureau, Shijiangzhuang 050021, China
Abstract: In the eastern margin of Bayan Har block, the Gravity/GPS measurement data from 2009 to 2013 and the EGM2008 model gravity data are merged by the least square collocation based on the remove and restore method. Simultaneously, the free-air gravity anomalies (FGAs) and the Bouguer gravity anomalies (BGAs) in the study area are updated. According to the new BGA field above and the Crust 1.0 data, the density structure of the crust is inverted by using the two-dimension polygon method and the nonlinear least square method. Based on the incompressible crust and the isostatic adjustment method, we propose a new method to calculate the vertical tectonic stresses. Furthermore, the distribution of the vertical tectonic stresses around the eastern margin of the Bayan Har block is calculated by using the density structure above. Through the vertical tectonic stress field, we find that the high positive vertical stresses (~40 MPa) are around the southern segment of the Longmen Shan faults, and the high negative vertical stresses (~-30 MPa) are around Barkam. According to the distribution of 99 earthquakes (MS > 5) since 1970 in the study area, we find that the gradient belts of the vertical stresses are corresponding to the lateral distribution of the earthquakes, and the area of the negative and positive vertical stresses corresponds to the deep and shallow earthquakes, respectively. However, when the horizontal deformation of the crust is large, the corresponding relation between vertical tectonic stresses and the focal depth of earthquakes is weak.
Key words: Bayan Har block      Bouguer gravity anomaly      Structure of the crustal density      Vertical tectonic stress     
1 引言

地壳构造应力的分布对于研究地震的孕育和发生具有重要意义.随着现代大地测量技术的进步,观测数据的丰富使得获取区域应力场分布成为可能.巴颜喀拉块体东缘的龙门山断裂带近年来相继发生了2008年汶川MW7.9地震与2013年芦山MW6.6地震,造成巨大人员伤亡和财产损失,同时也使得该区域成为地学研究的热点地区.

近年来不少学者在巴颜喀拉块体东缘展开地壳构造应力研究.Wu等(2009)利用汶川MW7.9地震发生后4个月内观测的3个钻孔的水压致裂应力数据,与之前的数据对比分析,发现龙门山断裂带西南段的应力水平依然维持在较高水平.Zhang等(2009)依据汶川MW7.9地震发生前的地震波形数据,分析了该区域的应力分布状态,发现汶川MW7.9地震前地震活动处于稳定状态,并且小震的微破裂应力场分布与主震一致.汶川MW7.9地震之后,陈群策等(2012)利用水压致裂方法测量了6个孔深在200~500 m之间的钻孔应力,发现龙门山断裂带上盘的应力作用强度高于下盘,其印模定向实验数据表明龙门山断裂带东北段应力方向与西南端不同.刘峡等(2014)利用1997—2007年期间的GPS数据,结合有限元方法研究了汶川MW7.9地震与芦山MW6.6地震发生前的应力积累状态,发现地震发生前龙门山断裂的挤压应力积累为全区最高.高尚华等(2016)提出了一种基于均衡理论的计算垂向构造应力的新方法,并在均匀地壳假设前提下,利用巴颜喀拉块体东缘的重力/GPS联合观测数据,计算了观测路线上的垂向构造应力的分布,并指出地震空间分布与垂向构造应力的梯度带具有对应关系.

上述构造应力领域的研究大多数集中于讨论横向构造应力与地震孕育和发生的关系,而垂向构造应力对于地震的影响鲜有讨论,虽然Wu等(2009)陈群策等(2012)所采用的钻孔应力观测中包含垂向构造应力的数据,但是其钻孔深度较浅,且观测的仅为观测站所在地的应力状态.高尚华等(2016)计算的垂向构造应力为整个岩石圈所承受垂向构造应力分布,然而其采用的地壳模型为均匀地壳,并未考虑地壳密度分层结构的影响,并且数据范围仅局限在观测路线附近.本文首先使用基于消去-恢复原理的最小二乘配置方法,融合巴颜喀拉块体东缘的重力/GPS观测数据与EGM2008重力模型数据,更新该地区布格重力异常场.在此之后,以Crust1.0地壳模型为初始值反演区域地壳密度分层结构.最后,基于地壳不可压缩和均衡调整原理,并结合以上地壳密度结构精细计算了巴颜喀拉块体东缘垂向构造应力分布,同时讨论了该区域垂向构造应力分布与地震空间分布的关系.

2 巴颜喀拉块体东缘地表观测与模型重力数据融合

图 1为收集整理的巴颜喀拉块体东缘重力/GPS联合观测数据的测线位置,其中松潘至武胜测线和马尔康至内江测线的观测时间为2009年(Zhang et al., 2014),测线横跨龙门山东部地区,与龙门山走向垂直,通过映秀—北川断裂、灌县—安县断裂和汶川—茂汶断裂,且马尔康至内江断裂通过汶川MW7.9地震震中.四川盆地测网的观测时间为2012年(Fu et al., 2014; Fu and Zhang, 2014),测网覆盖四川盆地西部以成都为中心包括龙泉山断裂的大部分地区.小金至雅安和康定至雅安测线的观测时间为2013年(杨光亮等, 2015),其中小金至雅安测线与Zhang等(2014)的剖面走向相似,同样跨越龙门山断裂带,并且通过芦山MW6.6地震震中,康定至雅安测线跨越鲜水河断裂、汶川—茂汶断裂和灌县—安县断裂.通过以上描述可知,巴颜喀拉块体东缘的重力/GPS联合观测数据对该地区的断裂带具有较好的覆盖,同时也覆盖了龙门山大部分地区.以上观测数据在经过正常重力改正和高度改正后得到测线的自由空气重力异常数据.

图 1 巴颜喀拉块体东缘重力与GPS联合观测数据分布图 其中,红线标示断层,黑色三角标示城镇,黑色五角星与沙滩球标示2008年汶川MW7.9地震与2013年芦山MW6.6地震以及其震源机制. Fig. 1 Distribution of measurement profiles in eastern margin of Bayan Har block The black stars and beach balls show the locations and focal mechanisms of the 2008 MW7.9 Wenchuan earthquakeand the 2013 MW6.6 Lushan earthquake. Black triangles denote cities. Red lines are faults.

实测数据中跨越龙门山及其以西地区的重力/GPS联合观测受限于当地陡峭的地形和交通状况,基本沿峡谷进行观测所得.使用观测数据计算的自由空气重力异常的精度是可以保证的,然而由于自由空气异常与地形的相关性,对于地形起伏较大的地区,重力异常变化同样较为剧烈.具体来说,高程的变化会产生较大的重力变化,其量值约为0.3086 mGal·m-1,这时使用普通的插值方法获取重力异常分布会出现较大的误差,甚至谬误.因此,在地形起伏较为剧烈的地区,直接采用观测数据插值计算自由空气异常分布一般精度不足(She et al., 2016).

随着CHAMP、GRACE和GOCE等重力卫星的发射,地球低阶重力场精度被极大地提高,基于以上卫星重力观测数据,并结合测高卫星等数据构建的重力模型也处于不断的更新中.EGM2008(Pavlis et al., 2012)地球重力模型是美国国家地理空间局NGA于2008年4月发布的全球超高阶地球重力场球谐模型,该模型融合地表重力观测、卫星测高、重力卫星以及航空观测等数据,其球谐系数展开至2190阶和2159次.该模型自发布以来一直处于更新中,到目前为止发布的最新版本为V23.1(http://topex.ucsd.edu/).为了比较不同重力模型的差异,分别从EGM2008和EIGEN-6C4(Förste et al., 2014)重力模型中提取图 1所示测线的自由空气重力异常值,并与观测值进行对比分析.发现EGM2008和EIGEN-6C4与实测值差异的标准差分别为47.98 mGal和57.81 mGal,可见在研究区域EGM2008相对于EIGEN-6C4有较高的精度,因此本文后续研究工作选取EGM2008重力模型进行数据融合.

在EGM2008重力模型的构建中,对于重力观测数据较为稀少的地区(如两极、俄罗斯、中国和非洲等地区),重力值为地形数据推算得到(Pavlis et al., 2012).由于中国地区为重力专利区,所以EGM2008包含的研究区域高精度重力值分辨率较低,大部分地区的重力值为地形推算结果.杨金玉等(2012)付广裕等(2013)分别对中国地区的EGM2008模型数据和实测数据进行了对比分析,结果显示川西地区模型值和实测值差异为10 mGal (付广裕等, 2013);中国东部地区两者差异为10 mGal左右,而在西部地区差异为50 mGal左右(杨金玉等,2012).由以上分析可知,EGM2008在中国地区存在系统性的误差,且高程起伏对于该系统性误差具有较大影响.

据以上对实测数据和EGM2008模型数据的分析可知,由于自由空气重力异常与地形正相关,直接插值获取该区域自由空气重力异常会产生较大的误差,甚至谬误,而EGM2008重力模型数据在研究区域的精度又低,且其精度受地形起伏的影响.因此,单一使用以上任何一种数据都很难获取高精度、高分辨率的自由空气重力异常场.鉴于此,本文引入基于消去-恢复原理的最小二乘配置方法,融合地表实测和EGM2008自由空气重力异常数据,以获取巴颜喀拉块体东缘的自由空气重力异常分布.

基于消去-恢复的最小二乘配置方法被广泛用于将地表观测数据融入EMG2008重力模型的过程中,同时此方法被普遍地用于不同观测类型重力数据的融合,如船载重力观测数据与卫星观测数据的融合(Hwang and Parsons, 1995; Zhao et al., 2017).在EGM2008的构建过程中,基于高程模型数据DTM2006.0(Pavlis et al., 2012)导出的重力异常场被作为消去-恢复法的参考场,DTM2006.0是一个全球数字地形模型,包含陆地和海洋的高程、深度和冰川的厚度等信息.参考EGM2008的构建过程,本文将EGM2008数据作为消去-恢复法的参考场来处理.观测数据与EGM2008数据的融合主要包括以下三个步骤.

首先,从EGM2008重力模型中提取对应观测点位置的自由空气重力异常数据gegm,并从实测自由空气重力异常数据gm中扣除以上提取的数据,以获取实测与模型数据的残插值Δgres.

(1)

其次,获取残差与距离的统计关系.计算残差值Δgres之间的自协方差值Cres

(2)

公式(2) 中N为全部观测点的数量,d代表观测点之间的距离,本文取d为10 km,因此距离范围可以划分为0~5 km、5~15 km、15~25 km,依次类推.Δgresl(d)和Δgresm(d)代表落入相同距离范围的残差值(lm),n为落入相同距离范围残差值的数量.Cres随相应的空间距离变化的统计信息如图 2a所示,利用此统计信息的拟合结果(图 2b中实线)构建Δgres与推估值之间的协方差矩阵.需要注意的是协方差的负值在残差值与距离的关系中没有实际意义,故在此仅对正值进行拟合.

图 2 实测与模型自由空气重力异常残差值协方差与拟合结果 (a)残差值协方差Cres; (b)残差值协方差拟合结果. Fig. 2 The covariance of the FGAs residual between measurements and model data, and the fitting result of the covariance (a) The covariance of the residual (Cres); (b) The fitting result of Cres .

最后,使用最小二乘配置方法(Moritz, 1980) 参考以上构建的协方差矩阵推估研究区域0.05°×0.05°分辨率的残差数据:

(3)

公式(3) 中LDΔ分别为与其平均值的差值和误差矩阵(观测误差对角阵),将加回研究区域相同分辨率的EGM2008自由空气重力异常数据gegm:

(4)

从而得到了两种数据融合的结果gmerge.EGM2008和使用最小二乘配置方法融合后的自由空气重力异常分布如图 3a图 3b所示.融合后的自由空气重力异常数据,在龙门山及其以西地区明显小于重力模型的结果,表明最小二乘配置方法使用观测数据的先验信息有效地压制了在地形起伏剧烈地区的重力模型数据.

图 3 巴颜喀拉块体东缘自由空气重力异常分布 (a) EGM2008自由空气重力异常分布;(b)实测和EGM2008融合的自由空气重力异常分布.黑色三角标示城镇. Fig. 3 Free-Air gravity anomalies around the margin of Bayan Har block of EGM2008 data and in situ observations. Black triangles denote cities. (a) Free-Air gravity anomalies of EGM2008; (b) Free-Air gravity anomalies derived from the combination
3 巴颜喀拉块体东缘地壳密度结构

图 3b所示的融合后的自由空气重力数据进行层间改正,并使用ASTER GDEM 2011数字高程数据分别以1″×1″、5″×5″、10″×10″网格数据对近场(0~2 km)、中场(2~20 km)和远场(20~167 km)进行地形改正,最后获取布格重力异常数据(图 4).布格重力异常分布整体呈西北低东南高的趋势,且梯度变化较大,研究区域的地形与布格重力异常呈负相关分布.

地震层析成像研究表明,巴颜喀拉块体东缘具有明显的层状构造(王有学等,2005邓文泽等,2014赵盼盼等,2015; Liu et al., 2014),因此不应忽略地壳层状结构对于后续研究的影响.本文选取CRUST1.0(Laske et al., 2013)全球地壳模型作为初始地壳密度结构,并结合布格重力异常数据,反演获取巴颜喀拉块体东缘地壳密度结构.首先将研究区域划为9个东西走向的剖面,经度101°至107°,纬度29°至33°,剖面间隔为0.5°.其次,从CRUST1.0中提取各剖面的地壳密度结构数据,同时从布格重力异常图像中提取各剖面对应的重力异常值.需要注意的是提取的布格重力异常数据的间隔为0.1°,并对布格重力异常进行了50 km的低通滤波处理,以消除高频信号的影响.最后,采用二维多边棱柱体模型(Talwani et al., 1959)作为正演模型,反演方法使用非线性最小二乘方法(Marquardt, 1963),以获取每条剖面对应的地壳密度结构,同时采用三次样条插值对该剖面数据进行插值以获取该地区的三维地壳密度结构(图 5).

图 4 巴颜喀拉块体东缘布格重力异常分布图 黑色三角标示城镇. Fig. 4 The distribution of Bouguer gravity anomalies in eastern margin of Bayan Har block Black triangles denote cities.
图 5 巴颜喀拉块体东缘地区2-D和3-D地壳密度结构 Fig. 5 The 2-D and 3-D structure of the crustal density in eastern margin of Bayan Har block

图 5给出了巴颜喀拉块体东缘地壳密度构造图.该结果中莫霍面的深度与Liu等(2014)Zhang等(2014)基本一致,研究区域东端莫霍面深度为40 km左右,西端莫霍面深度为60 km左右,龙门山地区为显著的莫霍面梯度带.北部地区莫霍面在经度102°—104°地区存在明显的起伏变化,而在南部地区莫霍面变化趋势较为一致,整体趋势呈现为自东向西逐渐变深.上地壳、中地壳和下地壳的变化趋势与莫霍面基本一致,沉积层主要分布于四川盆地内部,并且深度较浅,最深为10 km左右.

4 巴颜喀拉块体东缘垂向构造应力场

在均衡理论中假定地壳处于静水力平衡状态,即其所受重力与浮力相等,所受重力为地壳自重,而浮力是由于地壳排开地幔物质而引起的.为了保持这种平衡状态,密度界面的扰动就需要其他密度界面产生相应的挠曲进行补偿.因此,均衡问题实际上是一个平衡密度界面扰动的加载问题.

由以上介绍可知,在地壳不可压缩、为弹性体并且符合浮力原理的假设下,不同密度界面的扰动都会造成其他界面的挠曲来进行补偿.需要注意的是下文中对于密度界面的加载和补偿的研究均在频域下进行,其为k波数(k=2π/λ, λ为波长)的函数.Forsyth(1985)首先对该问题进行了研究,他将岩石圈分为地壳地幔两层进行处理,密度界面为地表与莫霍面.地壳起伏HI造成的加载由莫霍面挠曲WT来补偿,莫霍面起伏WI则由地表挠曲HM进行补偿,且假设这两个过程同时进行(图 6).根据以上设定,均衡调整后的界面HW即为力平衡状态.

图 6 均衡调整过程示意图 Fig. 6 Sketch of the isostatic adjustment

Simons等(2003)利用Turcotte和Schubert(1982)的平衡方程以及Forsyth(1985)的方法在频域给出了多层密度结构地壳均衡调整公式:

(5)

其中,HTWT是起始地表加载HI的均衡状态,相应地其他各层密度界面的起始加载WIi所形成的均衡状态为HMiWMi.Δi为相邻两层介质的密度差,ij均为地壳密度界面层数,D为挠曲刚度,E为杨氏模量,Te为有效弹性厚度,υ为泊松比.依据Forsyth的思路,起始加载与均衡状态的关系式:

(6)

将公式(5) 代入公式(6) 可得公式(7):

(7)

利用该式可计算频域下不同密度界面起始加载所对应的地表均衡状态HTHMi.

为了计算均衡调整后界面分布情况,延续Forsyth(1985)的思路,给出均衡状态与均衡调整结果的关系式:

(8)

其中H1Wi为均衡调整后的频率域地形以及其他各层密度界面.需要注意的是,由于地壳的不可压缩性,地壳内部密度层扰动造成的其他密度界面的补偿都是相等的,即其他界面对应HT的补偿均为WTWMi对应的其他各密度界面补偿应为HMi.将公式(5) 和(7) 代入公式(8) 得到求解均衡调整后的地壳密度结构的公式:

(9)

通过公式(9) 可根据现有的任意层数的地壳密度结构计算均衡调整后的结果.需要注意的是,以上公式都是在频域中计算的,输入为频域下的各层密度界面数据,因此在开始计算时要先将空间域数据通过二维傅里叶变换转到频率域,并且最后需要对输出的均衡调整后的频域数据进行逆傅里叶逆变换转为空间域数据.

为了更直观地了解上述方法,本文模拟了三层密度界面不同加载位置的结果.上地壳、下地壳和地幔密度分别设定为2.6 g·cm-3、2.9 g·cm-3以及3.3 g·cm-3.上下地壳界面深度为15 km,莫霍面深度为35 km,有效弹性厚度取0,泊松比取0.25,杨氏模量为1.78×1011 Pa.首先生成高斯分布的频域数据Z(k),作为模拟随机地形的基础.其次,考虑到真实地形的分布并不是完全的随机分布,地形的特征为在频域中高频域的振幅较低,低频域的振幅较大.因此引入地球物理学中使用较多的各向同性Matérn频谱模型作为生成模拟数据的功率谱密度S(k)(Goff and Jordan, 1988; Stein, 2012; Guttorp and Gneiting, 2006; Paciorek, 2007),

(10)

公式(10) 中,σ2, v, ω分别为Matérn频谱模型的第一,第二和第三参数.三个参数分别控制生成地形的振幅,平滑程度和整体分布.最后,使用公式

(11)

得到模拟地形的频域分布H(k),再进行二维傅里叶逆变换获得空间域的界面数据.在此利用上述方法模拟了三组地壳界面数据,每组数据含有地表、上下地壳三层结构,并且依次仅有一层界面具有起伏,其他界面起伏为零.图 7a7b7c分别是三组数据中具有起伏的地表、上下地壳界面与莫霍面.图 7d7e7f依次对应的是经过均衡调整后的以上三个界面,图 7g7h7i为对应的均衡调整前后各层界面的差异,由于地壳的不可压缩性,所以三层界面调整前后的差异是一致的,故三组模拟数据的其他界面变形不再列举.从图 7g可见,上下地壳界面与莫霍界面对地表加载进行了补偿,地形较高处形成了对应的山根,地形较低处进行了相应的抬升.图 7i反应的是上下地壳界面与地形对莫霍面加载的补偿,莫霍面较深地区对应地表发生抬升,莫霍面较浅地区对应地表发生沉降.同上,图 7h是莫霍面与地表对上下地壳界面加载的补偿.从以上模拟计算可知,该方法可以计算来自任意界面加载所对应的均衡调整结果,适用于多层地壳密度界面的均衡调整计算,为后续进一步计算密度分层地壳的垂向构造应力提供了基础.

图 7 不同界面加载的均衡调整模拟结果 Fig. 7 The simulation of the isostatic adjustment with the different loading position

基于公式(9) 和图 5给出的地壳密度结构,同时设定有效弹性厚度Te为0(即符合艾里均衡理论),泊松比取0.25,杨氏模量为1.78×1011 Pa,可计算均衡调整后的龙门山地壳密度结构,继而得到均衡调整前后的界面变化(图 8).由于以上方法假设地壳是不可压缩的且符合均衡原理,故地壳在均衡前后的受力状态可分别表示为

图 8 巴颜喀拉块体东缘均衡调整前后差异分布 白色虚线为块体边界,黑色三角标示城镇. Fig. 8 The distribution of the differences between the before and after isostatic adjustment in the eastern margin of Bayan Har block The white dash lines denote the boundary of the blocks, and black triangles denote cities.

(12)

(13)

其中ρM为地幔密度,取3340 kg·m-3g取9.81 m·s-2hmohohiso分别为均衡调整前后的莫霍面厚度,ρihi分别为地壳各层的密度和厚度.由公式(12) 和(13) 易得垂向构造应力为

(14)

其中Δh即为均衡调整前后界面变化(图 8).利用公式(14) 得到巴颜喀拉块体东缘的垂向构造应力分布,如图 9所示.从计算结果可见,松潘、汶川、马尔康、宁强、内江和康定附近地区垂向构造应力向下,其余地区垂向构造应力向上,其中大邑至雅安地区的垂向构造应力值较大.汶川和芦山地震震中与该局部垂向构造应力极大值较为靠近,并且位于垂向构造应力的梯度带上.龙门山断裂带中南段蓄积了较高的正向构造应力(约40 MPa),马尔康周边地区蓄积了较高的负向构造应力(约-30 MPa).

图 9 巴颜喀拉块体东缘垂向构造应力与地震深度分布 红线为断层,白色虚线为块体边界,黑色三角标示城镇.圆形代表地震位置,其大小代表震级,其颜色代表地震深度. Fig. 9 The distribution of the vertical stress and earthquakes in the eastern margin of Bayan Har block The red lines denote faults, white dash lines denote the boundary between the blocks, and black triangles denote cities. The circles denote the location of the earthquakes (MS > 5), sizes of the circles denote the earthquake magnitude, and the color denotes the depth of the earthquakes.

需要注意的是,在忽略地壳横向构造活动的前提下,垂向构造应力为正表示该区域地壳运动趋势为抬升;反之,垂向构造应力为负表示地壳运动趋势为下沉.上述计算垂向构造应力的方法是基于地壳不可压缩和均衡调整原理,利用均衡调整前后的地壳界面变化获取垂向构造应力;高尚华等(2016)提出的计算垂向构造应力的方法基于均衡理论,利用高程以及重力数据计算均衡面和莫霍面的差异,进而得到垂向构造应力.经过计算发现两种方法的结果是一致的,但是前者考虑了地壳的不可压缩并可用于计算均衡调整前后的地壳界面变化,显然该方法更为完备.

5 巴颜喀拉块体东缘垂向构造应力场与地震分布的对应关系

地壳的应力分布状态应与地壳内部地震分布有直接关联,本文以巴颜喀拉块体东缘地区垂向构造应力分布为基础定性地讨论其与地震空间分布的关系.图 9绘制了研究区域1970年以来5级(MS)以上的历史地震以及垂向构造应力的分布,其中地震标示的颜色代表震源深度.

在艾里均衡模型中,地壳为塑性,此时地壳有效弹性厚度Te为0,在此前提下地震更容易发生在莫霍面与均衡面差异的梯度带上(张永谦等,2010).根据本文计算垂向构造应力的方法,两个界面差异的梯度带即为垂向构造应力的梯度带(图 9).从图 9的结果来看,一方面地震的水平方向上的空间分布与垂向构造应力的梯度带分布具有对应关系,这与前人的研究一致(高尚华等,2016).另一方面,研究区域部分地震的深度分布与垂向构造应力的方向分布有一定的相关性,即负的构造应力对应深源地震,正的构造应力对应浅源地震.具体来说,研究区域西北部、中北部和西南地区的地震深度较深,垂向构造应力为负,其他垂向构造应力为正的区域地震深度较浅.

若定义大于研究区域平均地震深度(15.6 km)为深震,小于为浅震,符合上述地震深度与垂向构造应力对应关系的地震占比为53%.符合对应关系的地震事件占比较低的原因可能是本文讨论的垂向构造应力未考虑横向构造活动.由于地壳的横向构造活动是影响地震活动的重要因素,因此将地震以构造区域分为三组:龙门山逆冲断裂带、鲜水河走滑断裂带以及位于巴颜喀拉块体内部的松潘—马尔康及其以西区域进行讨论.

分别统计以上三个区域内符合上文对应关系的地震事件的占比,其中龙门山断裂带地区占比为50%,鲜水河断裂带占比为54%,松潘—马尔康地区占比为69%.从以上结果和图 10中近十年来相对华南地块的GPS水平速度场(Zhang et al., 2013)可见,龙门山和鲜水河断裂带符合比例较低,可能是由于断层两侧地壳处于强烈的横向变形造成的,即横向构造活动对该地区地震活动的影响较大;而在松潘—马尔康及其以西地区地壳水平运动速率和方向较为一致,说明横向构造活动对该区域影响较小,故符合地震深度与垂向构造应力的对应关系的地震事件占比较高.

图 10 巴颜喀拉块体东缘GPS水平速度场,相对华南块体 红线为断层,黑色三角位地名,黑色虚线为块体边界. Fig. 10 Horizontal GPS velocities of the eastern margin of the Bayan Har block, respected to South China block Red lines denote the faults, black triangle denote cities.

为了解释垂向构造应力与地震深度的关系,需要了解地壳内部的应力分布状态.一方面,地壳受到梯度较大的垂向构造应力的影响时,地壳内的不连续面即断层必然是容易发生相对运动的,即触发地震.但是以上机制仅能解释地震在水平面上的位置分布,而不能解释垂向构造应力与地震深度的关系.另一方面,当地壳受到垂向构造应力的影响时,地壳内部的应力分布会发生变化.McNutt(1980)Turcotte和Schubert(1982)对该问题进行了详细的研究,他们给出了弹性板在受到垂向应力的作用时,其内部的正应力分布(图 11).从图 11可见,当垂向构造应力为负时,地壳浅部呈压缩状态,地壳深部呈拉长状态,根据库伦破裂准则,地壳深部库仑应力增大,更易触发地震;同理,当垂向构造应力为正时,地壳浅部库伦应力增大,更易触发地震.由于该研究区域范围有限,对于垂向应力与地震深度分布的关系尚需进一步的研究进行解释验证.

图 11 垂向构造应力影响下的地壳内部应力状态示意图 Fig. 11 Sketch of the stresses in the crust which is controlled by the vertical tectonic stress
6 结论

首先,本文基于消去-恢复原理的最小二乘配置方法,对地表重力/GPS联合观测数据与EGM2008重力模型数据进行融合,更新了巴颜喀拉块体东缘地区的自由空气重力异常场.对自由空气重力异常进行重力归算,获取了巴颜喀拉块体东缘布格重力异常场.然后基于上述布格重力异常数据,采用二维多边棱柱体正演公式和非线性最小二乘反演方法,以CRUST1.0地壳模型为初始值,反演得到了巴颜喀拉块体东缘的地壳密度结构.

其次,依据Forsyth(1985)Simons等(2003)关于均衡调整的公式,给出了计算任意密度分层结构地壳的均衡调整计算公式.同时,提出了基于地壳不可压缩和均衡调整原理的计算垂向构造应力的新方法,并得到了巴颜喀拉块体东缘的垂向构造应力分布.结果表明,龙门山断裂带中南段蓄积了较高的正向构造应力(约40 MPa),马尔康周边地区蓄积了较高的负向构造应力(约-30 MPa).

经比较发现,巴颜喀拉块体东缘垂向构造应力分布与5级以上地震空间分布存在对应关系.具体地说,在水平方向上,垂向构造应力的梯度带与地震分布有对应关系;在垂直方向上,垂向构造应力为正的地区更易触发浅源地震,垂向构造应力为负的地区更易触发深源地震.需要注意的是,将地震以构造区域分组并结合GPS水平速度场的分析结果显示,横向构造运动对以上结论有较大影响,在地壳横向变形强烈的区域,垂向构造应力与地震深度的对应关系减弱.依据弹性板模型与库伦破裂准则对以上对应关系进行了机理分析,发现垂向构造应力方向的不同造成地壳内部正应力随深度的变化,进而影响了地震的触发深度.由于该研究区域范围有限,对于垂向应力与地震深度分布的关系尚需进一步的研究进行解释验证.

参考文献
Chen Q C, Feng C J, Meng W, et al. 2012. Analysis of in situ stress measurements at the northeastern section of the Longmenshan fault zone after the 5.12 Wenchuan earthquake. Chinese J. Geophys. , 55(12): 3923-3932. DOI:10.6038/j.issn.0001-5733.2012.12.005
Deng W Z, Chen J H, Guo B, et al. 2014. Fine velocity structure of the Longmenshan fault zone by double-difference tomography. Chinese J. Geophys. , 57(4): 1101-1110. DOI:10.6038/cjg20140408
Förste Christoph, Bruinsma Sean L, Abrikosov Oleg, et al. 2014. EIGEN-6C4 The latest combined global gravity field model including GOCE data up to degree and order 2190 of GFZ Potsdam and GRGS Toulouse. GFZ Data Services. http://doi.org/10.5880/icgem.2015.1.
Forsyth D W. 1985. Subsurface loading and estimates of the flexural rigidity of continental lithosphere. Journal of Geophysical Research: Solid Earth, 90(B14): 12623-12632. DOI:10.1029/JB090iB14p12623
Fu G Y, Gao S H, Freymueller J T, et al. 2014. Bouguer gravity anomaly and isostasy at western Sichuan Basin revealed by new gravity surveys. Journal of Geophysical Research: Solid Earth, 119(4): 3925-3938. DOI:10.1002/2014JB011033
Fu G Y, Zhu Y Q, Gao S H, et al. 2013. Discrepancies between free air gravity anomalies from EGM2008 and the ones from dense gravity/GPS observations at west Sichuan Basin. Chinese J. Geophys. , 56(11): 3761-3769. DOI:10.6038/cjg20131117
Fu G Y, Zhang G Q. 2014. Significant isostatic imbalance near the seismic gap between the M8.0 Wenchuan and M7.0 Lushan earthquakes. Chinese Science Bulletin, 59(34): 4774-4780.
Gao S H, She Y W, Fu G Y. 2016. A new method for computing the vertical tectonic stress of the crust by use of hybrid gravity and GPS data. Chinese J. Geophys. , 59(6): 2006-2013. DOI:10.6038/cjg20160607
Goff J A, Jordan T H. 1988. Stochastic modeling of seafloor morphology: Inversion of sea beam data for second-order statistics. Journal of Geophysical Research: Solid Earth, 93(B11): 13589-13608. DOI:10.1029/JB093iB11p13589
Guttorp P, Gneiting T. 2006. Studies in the history of probability and statistics XLIX on the Matérn correlation family. Biometrika, 93(4): 989-995. DOI:10.1093/biomet/93.4.989
Hwang C, Parsons B. 1995. Gravity anomalies derived from Seasat, Geosat, ERS-1 and TOPEX/POSEIDON altimetry and ship gravity: A case study over the Reykjanes Ridge. Geophysical Journal International, 122(2): 551-568. DOI:10.1111/gji.1995.122.issue-2
Laske G, Masters G, Ma Z T, et al. 2013. Update on CRUST1.0—A 1-degree global model of Earth's crust. Geophys. Res. Abstracts, 15: EGU2013-2658.
Liu Q Y, Van der Hilst R D, Li Y, et al. 2014. Eastward expansion of the Tibetan Plateau by crustal flow and strain partitioning across faults. Nature Geoscience, 7(5): 361-365. DOI:10.1038/ngeo2130
Liu X, Sun D Y, Ma J, et al. 2014. Present-day deformation and stress state of Longmenshan fault from GPS results—comparative research on active faults in Sichuan-Yunnan region. Chinese J. Geophys. , 57(4): 1091-1100. DOI:10.6038/cjg20140407
Marquardt D W. 1963. An algorithm for least-squares estimation of nonlinear parameters. Journal of the Society for Industrial and Applied Mathematics, 11(2): 431-441. DOI:10.1137/0111030
McNutt M. 1980. Implications of regional gravity for state of stress in the Earth's crust and upper mantle. Journal of Geophysical Research: Solid Earth, 85(B11): 6377-6396. DOI:10.1029/JB085iB11p06377
Moritz H. 1980. Advanced Physical Geodesy. Karlsruhe: Wichmann, 1.
Paciorek C J. 2007. Bayesian smoothing with Gaussian processes using Fourier basis functions in the spectralGP package. Journal of Statistical Software, 19(2): nihpa22751.
Pavlis N K, Holmes S A, Kenyon S C, et al. 2012. The development and evaluation of the Earth Gravitational Model 2008 (EGM2008). Journal of Geophysical Research: Solid Earth, 117(B4): B04406.
She Y W, Fu G Y, Wang Z H, et al. 2016. Gravity anomalies and lithospheric flexure around the Longmen Shan deduced from combinations of in situ observations and EGM2008 data. Earth, Planets and Space, 68(1): 163. DOI:10.1186/s40623-016-0537-7
Simons F J, Van der Hilst R D, Zuber M T. 2003. Spatiospectral localization of isostatic coherence anisotropy in Australia and its relation to seismic anisotropy: Implications for lithospheric deformation. Journal of Geophysical Research: Solid Earth, 108(B5): 2250.
Stein M L. 2012. Interpolation of Spatial Data: Some Theory for Kriging. New York: Springer.
Talwani M, Worzel J L, Landisman M. 1959. Rapid gravity computations for two-dimensional bodies with application to the Mendocino submarine fracture zone. Journal of Geophysical Research, 64(1): 49-59. DOI:10.1029/JZ064i001p00049
Turcotte D L, Schubert G. 1982. Geodynamics: Applications of Continuum Physics to Geological Problems. New York: Wiley.
Wang Y X, Han G H, Yuan X C, et al. 2005. Crustal p-wave velocity structure from Altyn Tagh to Longmen mountains along the Taiwan-Altay geoscience transect. Chinese J. Geophys. , 48(1): 98-106.
Wu M L, Zhang Y Q, Liao C T, et al. 2009. Preliminary results of in-situ stress measurements along the Longmenshan Fault Zone after the Wenchuan MS8.0 Earthquake. Acta Geologica Sinica, 83(4): 746-753. DOI:10.1111/j.1755-6724.2009.00098.x
Yang G L, Shen C Y, Wu G J, et al. 2015. Bouguer gravity anomaly and crustal density structure in Jinchuan-Lushan-Qianwei profile. Chinese J. Geophys. , 58(7): 2424-2435. DOI:10.6038/cjg20150719
Yang J Y, Zhang X H, Zhang F F, et al. 2012. On the accuracy of EGM2008 earth gravitational model in Chinese Mainland. Progress in Geophysics , 27(4): 1298-1306. DOI:10.6038/j.issn.1004-2903.2012.04.003
Zhang Y Q, Wang Q S, Teng J W. 2010. The crustal isostatic anomaly beneath eastern Tibet and western Sichuan and its relationship with the distribution of earthquakes. Chinese J. Geophys. , 53(11): 2631-2638. DOI:10.3969/j.issn.0001-5733.2010.11.011
Zhang Y Q, Teng J W, Wang Q S, et al. 2014. Density structure and isostatic state of the crust in the Longmenshan and adjacent areas. Tectonophysics, 619-620: 51-57. DOI:10.1016/j.tecto.2013.08.018
Zhang Z Q, McCaffrey R, Zhang P Z. 2013. Relative motion across the eastern Tibetan Plateau: Contributions from faulting, internal strain and rotation rates. Tectonophysics, 584: 240-256. DOI:10.1016/j.tecto.2012.08.006
Zhang Z W, Cheng W Z, Ruan X, et al. 2009. The seismicity and tectonic stress field characteristics of the Longmenshan fault zone before the Wenchuan MS8.0 earthquake. Earthquake Science, 22(2): 119-128. DOI:10.1007/s11589-009-0119-x
Zhao P P, Chen J H, Liu Q Y, et al. 2015. Fine structure of middle and upper crust of the Longmenshan Fault zone from short period seismic ambient noise. Chinese J. Geophys. , 58(11): 4018-4030. DOI:10.6038/cjg20151111
Zhao Q, Wu W W, Wu Y L. 2017. Using combined grace and GPS data to investigate the vertical crustal deformation at the northeastern margin of the Tibetan Plateau. J. Asian Earth Sci., 134: 122-129. DOI:10.1016/j.jseaes.2016.11.010
陈群策, 丰成君, 孟文, 等. 2012. 5.12汶川地震后龙门山断裂带东北段现今地应力测量结果分析. 地球物理学报, 55(12): 3923–3932. DOI:10.6038/j.issn.0001-5733.2012.12.005
邓文泽, 陈九辉, 郭飚, 等. 2014. 龙门山断裂带精细速度结构的双差层析成像研究. 地球物理学报, 57(4): 1101–1110. DOI:10.6038/cjg20140408
付广裕, 祝意青, 高尚华, 等. 2013. 川西地区实测自由空气重力异常与EGM2008模型结果的差异. 地球物理学报, 56(11): 3761–3769. DOI:10.6038/cjg20131117
高尚华, 佘雅文, 付广裕. 2016. 利用重力/GPS联合观测数据计算地壳垂向构造应力的新方法. 地球物理学报, 59(6): 2006–2013. DOI:10.6038/cjg20160607
刘峡, 孙东颖, 马瑾, 等. 2014. GPS结果揭示的龙门山断裂带现今形变与受力——与川滇地区其他断裂带的对比研究. 地球物理学报, 57(4): 1091–1100. DOI:10.6038/cjg20140407
王有学, 韩果花, 袁学诚, 等. 2005. 台湾—阿尔泰地学断面阿尔金—龙门山剖面的地壳纵波速度结构. 地球物理学报, 48(1): 98–106.
杨光亮, 申重阳, 吴桂桔, 等. 2015. 金川—芦山—犍为剖面重力异常和地壳密度结构特征. 地球物理学报, 58(7): 2424–2435. DOI:10.6038/cjg20150719
杨金玉, 张训华, 张菲菲, 等. 2012. EGM2008地球重力模型数据在中国大陆地区的精度分析. 地球物理学进展, 27(4): 1298–1306. DOI:10.6038/j.issn.1004-2903.2012.04.003
张永谦, 王谦身, 滕吉文. 2010. 川西藏东地区的地壳均衡异常及其与地震分布的关系. 地球物理学报, 53(11): 2631–2638. DOI:10.3969/j.issn.0001-5733.2010.11.011
赵盼盼, 陈九辉, 刘启元, 等. 2015. 龙门山断裂带中上地壳速度结构的短周期环境噪声成像. 地球物理学报, 58(11): 4018–4030. DOI:10.6038/cjg20151111