华北地区是中国大陆地震活动强烈的地区之一, 自新生代至今, 华北裂陷盆地的活动依旧强烈, 近期显示出频繁地震活动的特征, 华北地区由此成为国际地球科学界关注的热点区域(王椿镛等, 2016)(图 1).张家口—渤海地震带(以下简称为“张渤地震带”)位于华北盆地、燕山隆起和太行山隆起三个构造单元的结合处(张四昌, 1993; 徐杰等, 1998), 区内断裂相交和切割, 构成复杂的孕震构造格局.张渤地震带西起太行山北缘, 沿横贯燕山与华北盆地结合区域分布, 东入渤海, 是我国东部地区一条重要的NW向地震活动带, 该地震带上发生过26次M6以上的地震(顾功叙, 1983; 赵博等, 2011).从构造背景上, 张渤地震带位于显生宙以来遭受破坏的华北克拉通中东部(Wu et al., 2005; 赵越等, 2010; 朱日祥等, 2012).华北克拉通破坏地区主要发生在太行山以东地区, 既包括岩石圈地幔的整体性破坏以及地壳的强烈改造和减薄作用, 也包括岩石圈地幔物理化学性质的明显改变(朱日祥等, 2011, 2012).华北克拉通破坏的机制是热-化学侵蚀还是拆沉作用还存在一定的争论; 目前一种倾向性的观点认为, 在华北的周边地区华北克拉通减薄的主要机制是板块俯冲导致的拆沉作用; 在华北内部更为重要的作用可能热侵蚀(吴福元等, 2008).但考虑该区中新生代岩浆活动极为强烈, 也不排除该区发生拆沉后, 地壳底部发生大规模岩浆底侵的可能性(朱日祥和郑天愉, 2009).

图 1

Fig. 1

图 1 华北地区构造背景和沿张渤地震带位置(白色虚框)示意图(王椿镛等, 2016) 图中所示地震为区域历史强震(≥6.0级)红色线AB为剖面线, 黑色虚线为主要地质构造线. Fig. 1 Geological settings of North China and section position of Zhangbo seismic belt(White Dashed Box) The hollow circles denote the historical strong earthquakes(≥6.0M), red line AB denotes vertical section,main geological structures represented by black dashed lines.

华北地区地震层析成像研究揭示了大尺度的地壳上地幔横向不均匀性的主要特征, 壳幔边界结构的地震观测结果(Zheng et al., 2008; 陈凌等, 2010a, 2010b), 揭示了华北克拉通破坏过程中不同地区经历了不同的岩浆作用和地壳改造过程(朱日祥等, 2012).可是这些图像也遗留了些疑问, 如P波速度异常是由温度的变化、物质组成的成分引起的, 或两者兼而有之, 还是孔隙流体的存在造成的？地球物理观测手段提供的下地壳、上地幔的物质成分和温度的变化也是相当有限的(孙若昧和刘福田, 1995).所获得的地壳厚度和泊松比结果(许卫卫和郑天愉, 2005; 罗艳等, 2008; 王峻等, 2009; 嵇少丞等, 2009; 葛粲等, 2011; 齐刚和陈棋福, 2015), 也难以精细探讨张渤地震带地壳结构所携带的构造演化信息.考虑到V_P/V_S波速比是反映岩石物理性质和流变学特征的重要参数, 比单纯的P波或S波速度参数提供了更为丰富的信息, 对提高地壳和上地幔物质的固有特征、动力学性质以及地质演化的认识有着重要的意义(张学民等, 2004).本文充分收集了华北地区区域数字地震台网176个固定台站所记录到的观测数据, 应用近震体波层析成像方法, 获得了华北地区地壳三维P波速度和V_P/V_S波速比结构, 在此基础上, 分析和讨论了张渤地震带深部孕震构造环境和地幔岩浆对地壳底侵作用等问题.

本文研究区域为110—120°E和36—43°N; 资料主要取自中国地震台网中心和河北省地震局台网中心的地震到时资料.研究地区的观测台站和所选用地震的震中分布图如图 2.

图 2

Fig. 2

图 2 本研究所用的地震(灰色圆形)和台站分布图(黑三角) (灰色线条为主要地质构造线) Fig. 2 The earthquake epicenters (blue circles) and the distribution of stations (red triangles) used in the study (The main geological structures represented by grey lines)

为了提高数据的质量, 依据走时残差分布检验数据可靠性的原则挑选数据.近震的走时残差可以表示为

其中, t_{i, j}^obs表示第i个台站接收的第j个事件的观测走时, t_{i, j}^cal表示射线追踪计算的理论走时.通过正演部分即伪弯曲射线追踪方法(Thurber, 1983; Haslinger and Kissling, 2001)获得走时残差.如果其残差值大于3 s, 则基本可以认为是数据存在错误, 需要查找原因或剔除, 尽可能消除数据对反演结果的影响.经过严格筛选, 最终获得了176个地震台站, 677个地震(其分布见图 2); 地震事件的震级分布如图 3, 震级在M2.5~4.7之间分布.筛选所得的P波和S波走时数据, 如图 4所示:P波和S波的震相分别为Pg、Pn、Sg和Sn(如图 4a).在大约270 km以后, 莫霍界面折射Pn、Sn震相首先到达, 成为首波, 所以在270 km以后只保留Pn、Sn震相; 本文将大于270 km的Pg、Sg震相数据作进一步筛选并剔除, 最终获得反演所用的P波、S波走时数据如图 4b, 其中P波走时43354条:Pg波32500条, Pn波10854条; S波走时30407条:Sg波23154条, Sn波7253条.需要指出的是震中距150~300 km之间存在一个非首波震相, 但并非反射震相, 而是Pn首波在成为初至波以后, 紧随其后的Pg震相.在数据筛选时之所以将其保留, 一是因为突起部分的台站没有记录到明显的Pn震相, 二是增加一部分中地壳的射线数, 以此来提高这一深度范围内的分辨率.

图 3

Fig. 3

图 3 所选地震震级分布图 Fig. 3 Magnitude distribution

图 4

Fig. 4

图 4 地震走时数据震相分布(a)图和参与反演数据站P波和S波走时曲线(b) Fig. 4 Travel time curves of P and S waves

本研究采用的网格节点对研究区进行了模型化分; 水平方向采用0.5°×0.5°均匀划分网格, 模型边缘间隔为1°; 在深度方向从地表至45 km之间, 每层的深度分别设为1、5、10、20、35、45 km, 地壳中初始速度模型参考人工地震测深的结果(嘉世旭和张先康, 2005; 嘉世旭等, 2005).本研究所用的一维P波初始速度模型如图 5; 为获取本区的V_P/V_S波速比, 利用P波, S波走时数据, 用和达法进行拟合得到图 6, 本区的V_P/V_S波速比在1.63~1.81之间, 平均值为1.73.

图 5

Fig. 5

图 5 P波一维初始速度模型 Fig. 5 1D P-wave velocity model

图 6

Fig. 6

图 6 和达法拟合波速比 Fig. 6 Wave velocity ratio fitting by Multi-station method

为了选取最为合适的阻尼系数, 分别用不同的阻尼值进行了多次反演试算, 权衡了走时残差的大小和速度模型的光滑程度, 在反演过程中阻尼满足了模型和数据方差均为最小时确定P波反演的阻尼系数, 通过绘制由不同阻尼因子产生的数据方差和模型方差的“折衷曲线”来观察确定阻尼因子的最佳取值点(Eberhart-Phillips, 2006)得到图 7中所示的trade-off曲线, 选取P波damp=100, V_P/V_Sdamp=100作为最终反演计算的阻尼系数.

图 7

Fig. 7

图 7 (a) P波trade-off曲线(damp=100);(b) VP/VS trade-off曲线(damp=100). Fig. 7 The tradeoff curves for P wave (a) and VP/VS (b)

地球物理反演的结果不是唯一的, 需要对反演结果的质量进行评价.每个节点附近通过的射线数也可以作为衡量解的可靠性的一个估计, 在地震射线密集并且射线穿越方向分布好的地方, 反演结果的分辨率就好, 通过地震射线分布图 8可知, 射线分布主要分布于张家口—渤海湾一带地震和台站均密集的地带以及太行山一带, 华北平原内部, 由于少地震发生, 射线分布明显稀疏.

图 8

Fig. 8

图 8 地震射线和地震台站(灰色三角)分布 Fig. 8 Ray distributionand distribution of stations (grey triangles) used in the study

图 9(a, b, c)是近震层析成像所用数据在不同深度上地震射线的分布.可以看到, 地壳中、上部的射线密度较大, 地壳下部的射线主要集中在京津唐和太原石家庄一带上.在地震层析成像过程中仅对穿过射线条数大于10的模型单元进行反演计算.图 9(d, e, f)显示了反演结果的分辨情况.图中给出的分辨率是与各个速度模型相对应的分辨率矩阵中的对角元素值.分辨率取值0~0.55, 数值越接近0.55, 分辨能力越好, 可以看出地壳中、上部的分辨率明显优于地壳下部的分辨率, 这和射线分布类似.比较图 9(a, b, c)和图 9(d, e, f), 可以看出地震射线数的分布和分辨率的分布有很大的相关性, 射线密集的地方具有较高的分辨能力.在射线方向分布较好的区域进行层析成像时, 用射线分布和分辨率分布来评价反演结果具有同样的效果(丁志峰, 1999).

图 9

Fig. 9

图 9 P波5 km, 15 km和30 km深度地震射线数分布(a, b, c)和反演结果的分辨率(d, e, f) Fig. 9 Ray distribution (a, b, c) and resolution test (d, e, f) at 5 km, 15 km and 30 km depth

图 10给出了5 km、15 km和30 km等不同深度上的P波速度(图a, b, c)和V_P/V_S波速比(图d, e, f)的水平切片图像.图中不同的P波速度值用不同颜色表示, 其中红色表示低速, 蓝色表示高速.图中不同的V_P/V_S波速比值则用红色表示高速度比, 蓝色表示低速度比.

图 10

Fig. 10

图 10 P波速度(a, b, c)和VP/VS波速比分布图(d, e, f) Fig. 10 The horizontal slices of the tomography results:(a, b, c) and (d, e, f) is for P wave and VP/VS, respectively

图 10a显示5 km深度处的P波速度异常分布与地表的地形密切相关, 华北盆地大部分区域P波速度特征表现为低速异常, 而华北盆地北西缘和燕山隆起带则显示出高速异常分布, 研究区西部太行山隆起则表现为低速异常.图 10d显示该深度处V_P/V_S波速比展布特征, 华北盆地大部分区域显示出较高的波速比, 而华北盆地北西缘及燕山隆起带则表现为较低的波速比.西部太行山和北部的燕山隆起P波正异常和偏低的V_P/V_S波速比相对应, 与地表广泛出露的基岩分布一致.由于本文采用的反演算法涉及到了模型协方差和拉格朗日因子的选择, 不同方法在这两个参数的选择不一样, 可能会导致模型的平滑或者光滑程度不一样, 同时不同的拉格朗日因子的选择模型的权重就不一样, 因此本文得到的异常体在速度值上与DSS结果会有一定的差异, 是反演方法和反演参数选择引起的, 但是无论什么方法, 得到的高速或低速整体的异常特征类似.

在15 km深度上, 华北盆地大部分区域P波速度特征表现为低速异常, 沿燕山山脉为高速异常分布, 太行山山脉则表现为高低速相间的特征(图 10b), 表明研究区上地壳介质明显的横向不均匀性.山间的盆地在周围岩石的高速背景下呈现出小范围的低速体, 虽然唐山地区位于华北盆地, 但它是张渤地震带上一个隆起的块体, 在速度图的上地壳层位上(图 10a和b)都显示为高速异常.这些说明在该区内的地壳速度结构存在明显的横向不均匀性.华北盆地和太行—燕山隆起带波速比分别表现为高、低异常展布特征.华北盆地、燕山和太行山隆起区显示出明显不同的速度结构和不同的波速比结构特征, 以太行山为界, 华北盆地P波速度低速异常显著, 同时V_P/V_S波速比较高, 这可能与松散巨厚的沉积层及其孔隙中含水有关(Eberhart-Phillips et al., 2005; Li et al., 2009).从图中还可以看出, 三河平谷8.0级地震、唐山7.8级地震和张北6.2级地震位于高低速异常分界线或偏向高速体一侧.

在30 km深度上, 研究区范围P波低速异常分布特征明显, 间或夹杂局部的高速异常, 华北盆地和太行—燕山隆起带波速比分别表现为高、低异常展布特征.近些年来发生的M5.0以上地震大多位于低速区或低速至高速的急剧过渡区, 高速区仅有一些震级较小的地震发生.在30 km深度上, 华北地区总体呈现P波低速异常分布, 沿燕山造山带、唐山地震区均表现为低速异常.已有研究结果表明华北地区大部分地震的震源深度小于25 km, 中下地壳低速异常分布被认为是强烈地震发生的深部构造环境.图 10c还揭示了震源体下方均存在广泛的低速异常分布, 如三河平谷8.0级地震、唐山7.8级地震和张北6.2级地震, 据此我们推测这些低速异常的存在与流体有关, 中下地壳中的流体容易引起上地壳发震层的弱化和应力集中, 进而导致大地震的发生.

文中根据P波速度和V_P/V_S波速比的平面分布图像获得了各构造单元速度结构的一些认识, 为了更好地揭示速度和波速比异常在深度上的变化特征, 圈定速度异常的分布范围和展布特征, 研究张渤带地区的孕震环境, 本文沿AB测线(41.5°N, 114°E; 39°N, 119°E)切出剖面图, 剖面位置如图 1所示, 剖面穿过了华北盆地和燕山隆起两大构造单元.

图 11a剖面揭示了张渤地震带下方的P波速度结构垂向特征, 结果显示张渤带地壳内部在20 km下方存在一低速带, 这一低速异常特征在以前的成像结果中(黄金莉和赵大鹏, 2005)也有显示.大地电磁测深资料(刘国栋等, 1983)确定了该唐山地区的地电结构, 揭示在唐山地震震源区下约20多公里深处存在高导异常, 这一结果与本文的层析成像结果相吻合.刘昌铨和嘉世旭(1986本条文献在文后文献中未体现)对穿过唐山地震区NW-SE走向的柏各庄—丰宁—正蓝旗DSS剖面解释所获得的二维速度结构显示:唐山M_S7.8地震位于莫霍界面局部上隆的东南侧; 震中区上地壳下部约20 km深度附近存在低速体, 且延续至下地壳; 大地震发生在高低速交界部位偏高速区的一侧, 并认为震源下方存在的低速异常体可能促进了地震成因; Huang和Zhao(2009)对首都圈地区强震发生的深部构造环境研究发现, 多数大地震都发生在高速块体的边侧, 而在震源区的下方存在明显的低速体分布, 并认为这些低速异常与流体有关; 深反射CDP叠加剖面显示三河—平谷8级大地震在震源区下方存在有低速异常体的存在, 莫霍界面起伏变化和较厚的反射叠层, 以及局部复杂的楔形反射带的存在, 均暗示地壳结构具有强烈挤压、变形的背景, 被认为是地震发生的深部环境(刘保金等, 2011b).结合从1970—2017年地震震中分布, 本文的研究结果也揭示出张渤地震带地震主要分布在上地壳, 从震源深度分布及地壳物性结构分析上地壳具有发震的构造背景, 也有发震的物性基础, 是一个易震层和多震层(段永红等, 2016).

需要说明的是由于该方法所限, 沿剖面见不到Moho界面的变化(图中灰色虚线为该区域的莫霍面深度, 数据参考(王帅军等, 2014; 段永红等, 2015; 刘志等, 2015)), 这主要是由于深纵向分层的大框架设计, 未能对该深度区间的下地壳及上地幔细构造反演计算, 其结果模糊了壳幔边界的深度变化.

刘昌铨和嘉世旭(1986)DSS剖面解释所获得的二维速度结构还显示:唐山M_S7.8地震震中区下方对应上地幔顶部的隆起区, 同时上地壳下部和下地壳的界面均显示局部上隆, 由此推测震源下方存在上地幔的局部活动, 地幔物质上升活动给地壳施加局部应力场, 从而导致地震的发生.接收函数结果显示:在唐山菱形地块下方上地幔顶部存在高达10 km的异常隆起; 唐山震区中上地壳S波速度存在明显的低速异常, 且上地幔顶部有明显异常, 故推测唐山大震区下方应有幔源物质的大规模侵入, 并在中上地壳形成壳内低速体; 唐山大地震能量的积累并非限于震源附近的地壳深度, 而是与地幔物质的运移和变形密切有关.丰南深反射剖面(FN2011)CDP剖面还显示, 下地壳具有良好的反射性质, 并在下地壳和壳幔过渡带顶部出现斜列的楔形或弧形强反射体, 可能反映了岩浆活动对下地壳物质和结构的强烈改造(刘保金等, 2011a; 王椿镛等, 2016).三河—平谷M8.0地震区深地震反射剖面显示的地壳界面强反射带在横向上的中断, 壳内出现的一些局部强反射事件可以认为是上地幔物质沿此深断裂上涌, 以及岩浆侵入从而导致与周围介质有较强的波阻抗差所致, 尤其在震源区下方左侧存在明显的叠层状局部反射事件, 可视为岩浆侵入的异常区域, 莫霍界面起伏变化和较厚的反射叠层, 以及局部复杂的楔形反射带的存在, 均暗示地壳结构具有强烈挤压、变形的背景; 由此推测岩浆活动对下地壳进行了物质和结构的强烈改造(刘保金等, 2011a; 王椿镛等, 2016).华北盆地地震深反射探测显示该区莫霍面有2~3 s近水平的断续强反射带, 这意味着岩浆沿地壳底部的分层底侵, 可能是地幔岩浆底侵地壳底部的结果(刘保金等, 2011b).

下地壳是联系地壳深部与上地幔的纽带, 是壳幔物质交换的重要场所, 也是研究华北克拉通破坏机制的关键所在.在地幔物质上涌区、热点区等地区, 剪切波的传播特性与P波有着很大的差异, P波速度改变的同时, 剪切波有时受到的影响会比P波大, 波速比出现高值; 但有些情况下二者的波速均下降的同时, 其波速比又不会发生很大改变.利用V_P/V_S波速比结果, 结合各种岩石的物理传播性质, 能对研究区域地下介质特性、组成成分以及热物质的物理性质加深认识.图 11b剖面揭示了张渤地震带下方的V_P/V_S波速比结构特征, 结果显示张渤带地壳内部在20 km下方存在偏高的波速比特征, 这些很可能反映了在该区域内中、上地壳幔源物质的侵入和热状态的岩体在横向上的变化(嘉世旭, 1996).地球化学的研究成果也揭示了该区域热状态的岩体在横向上的变化, 地球化学中氦(He)是U系衰变物, 具有很强的化学惰性和物理稳定性, 在水中溶解度低, 能从地下深部沿裂隙运移至地表, 延怀盆地内土壤气He平均含量明显高于大气含量, 前者是后者的3倍, 具有西低东高的特征(李营等, 2009); 龙庆峡水库有来自地壳的氦(谷元珠等, 2001)和延怀地区井泉水中检测到来自地幔源的氦(王广才等, 2003).这些地壳和地幔源的He能逸散到地表, 表明区内地壳深部热流体的供给量是不一致的.

图 11

Fig. 11

图 11 L剖面P波速度结构(a)和VP/VS波速比分布图(b) 图中空心圆由左向右分别代表 1998年张北6.2级地震、1679年三河平谷8.0级地震, 1976年唐山7.8级地震.ZB代表张北、BJ代表北京、SH代表三河和TS代表唐山, 灰色虚线代表莫霍面. Fig. 11 The vertical sections of the tomography results along the L profiles The figureabove represents P wave and the figure below represent VP/VS.The hollow circles from left to right represent the Zhangbei M6.2 earthquake in 1998, the Sanhe-Pinggu M8.0 earthquake in 1679, and the Tangshan M7.8 earthquake in 1976.ZB, BJ, SH, TS represents Zhangbei, Beijing, Sanhe and Tangshan respectively, grey dashed lines represent Moho.

综合P波速度结构、V_P/V_S波速比、大地电磁测深资料和地球化学资料揭示张渤带下方地壳内存在着低速层和高导层, 具有高波速比特征, 反映了幔源物质侵入, 壳内低速、热状态的岩体横向上的变化.虽然唐山地区位于华北平原, 但它是张渤地震带上一个隆起的块体, 在速度图的上地壳层位上(图 10a, b)都显示为高速异常.在深度30 km的水平切片上可以看到唐山地区表现为低速地带; 图 10(c, f)和图 11显示该地区低的P波速度和高波速比, 推测为幔源物质的上涌通道.张渤地震带地壳内部在20 km下方存在偏高的波速比特征, 这些可能反映了在该区域内中、上地壳幔源物质的侵入和热状态的岩体在横向上的变化, 可能是地幔岩浆长期底侵作用的结果; 推测由于幔源物质侵入, 使得地壳深部流体的供给量增加.在地壳发震层下长期存在的流体会影响断裂带的结构, 降低断裂带的强度, 使区域应力场发生变化(Hickman et al., 1995; Sibson, 1992)从而导致断裂带上应力的集中, 进而引发地震发生.

由于本次所使用的近震记录的射线主要分布于中上地壳深度范围, 成像结果未能得到更深处的速度结构, 今后将加入远震事件进行联合成像, 并通过研究对流体作用更为敏感的S波速度结构, 从不同深度来揭示研究区精细速度结构及介质的横向不均匀性, 期待获得的对该研究区的深部介质环境和地震孕育机理更加深入的认识, 这些无疑对于理解华北地区地震孕育的深部地球动力机制以及评估未来地震活动趋势提供重要的科学依据.

V_P/V_S波速比是反映岩石物理性质和流变学特征的重要参数, 本文利用区域地震走时数据, 应用近震体波反演方法, 获得华北地区地壳P波三维速度结构和V_P/V_S波速比, 并结合人工地震测深剖面、大地电磁测深、地球化学等成果, 讨论了张渤地震带的孕震环境和动力学机制等科学问题, 获得如下结果:

(1) 5 km深度处的P波速度异常分布与地表的地形密切相关, 华北盆地大部分区域P波速度特征表现为低速异常, 而燕山隆起带则表现为高速异常分布, 太行山隆起则呈现出低速异常.15 km深度处华北盆地大部分区域P波速度特征表现为低速异常, 而沿燕山山脉为高速异常分布, 太行山山脉则表现为高低速相间的特征, 表明研究区上地壳介质具有明显的横向不均匀性.在30 km深度上, 研究区范围P波低速异常分布特征明显, 间或夹杂局部的高速异常.5 km深度处V_P/V_S波速比展布特征, 华北盆地大部分区域表现为较高的波速比, 而燕山隆起带则表现为较低的波速比.西部太行山和北部的燕山隆起P波高速异常和偏低的V_P/V_S波速比相对应, 与地表广泛出露的基岩分布一致.15 km和30 km深度处华北盆地和太行—燕山隆起带波速比分别表现为高、低异常展布特征.

(2) 研究结果还表明华北地区强震位于高低速异常分界线附近或偏向高速体一侧, 震源体下方均存在广泛的低速异常分布, 据此我们推测这些低速异常的存在可能与流体有关, 本文的研究结果也揭示出张渤地震带地震主要分布在上地壳, 从震源深度分布及地壳物性结构分析上地壳具有发震的构造背景, 也有发震的物性基础, 是一个易震层和多震层.张渤地震带地壳内部在20 km下方存在偏高的波速比特征, 这些可能反映了在该区域内中、上地壳幔源物质的侵入和热状态的岩体在横向上的变化, 可能是地幔岩浆长期底侵作用的结果.由于幔源物质侵入, 使得地壳深部流体的供给量增加, 在地壳发震层下长期存在的流体会影响断裂带的结构, 降低断裂带的强度, 使区域应力场发生变化从而导致断裂带上应力的集中, 进而引发地震发生.

感谢中国地震台网中心和河北省地震局台网中心提供的地震走时数据.在论文的撰写过程中先后与李大虎博士进行了多次有益的交流和探讨, 审稿专家对本文提出了宝贵的修改意见, 在此表示衷心的感谢！