深地震宽角反射折射剖面联合迭代法的速度反演

深地震宽角反射折射剖面联合迭代法的速度反演

龙门山断裂带地震地质解释龙门山断裂带位于四川盆地西部，青藏高原东部，南起四川天全，北至东南40-50延伸至地壳，全长约500公里，宽约30-40公里（李勇等，2009a）（图1）。龙门山断裂带主要由三部分组成：四川-茂县断裂（后山断裂）、四川-映秀断裂（中央断裂）和江油-灌木断裂（前山断裂）。这三件古事几乎平行地展开。经过长期和多阶段的构造发育，该地区的地质构造变得越来越复杂（徐杰等人，2010；李志伟等人，2011）。龙门山断裂带作为一个稳定的扬子地块与活动的松潘-甘孜地块之间的过渡带，不仅显示了地形的严重变化，而且两侧地壳厚度（包括沉积层，下同）也有类似的变化模式。换句话说，地壳厚度从松潘-甘孜地块的57-58km减少到四川盆地的43-44km。地形和地壳厚度之间的巨大差异表明该区域发生了强烈的构造变形（李志伟等人，2011）。据历史地震资料记载,中等强度及以上的地震主要集中在龙门山断裂带的南段及中南段,且龙门山断裂带的3条断裂上均已发生过M6以上的中等强度地震.李勇等(2009b)根据历史地震记录分析推测,龙门山断裂带近几十年来已进入地震高发期,3条主断裂均具备发生更大强度地震的潜在危险.龙门山断裂带长期以来一直是地球科学家们关注的焦点.2008年5月12日,位于龙门山中央的映秀—北川断裂上的汶川县映秀镇发生M自20世纪80年代中期起,龙门山地区实施了一系列壳幔构造的探测计划,布设了多条人工地震剖面,仅横穿龙门山断裂带中段的测线就有两条.但是,由于当时观测系统不完善,地震仪器落后,一直未能得到翔实可靠的地下深部构造数据资料,难以对龙门山地区的地壳结构进行细致的研究.此外,李志伟等(2011)关于该区天然地震的观测研究显示,低速层分布和形态分布等方面存在着不一致.而汶川地震的发生,推翻了之前对于龙门山断裂带的一些认知,提出了一系列新的科学问题.为了更加深入细致地研究该区的壳幔精细速度结构,中国地震局地球物理勘探中心布设了一条横穿龙门山断裂带的人工深地震测深剖面.本文拟对横穿龙门山断裂带的这条人工深地震测深剖面数据进行反演计算和解释.首先用每一炮的多个震相反演一维速度模型,利用多炮一维模型插值出一个较为粗略的二维初始速度模型;然后使用射线追踪法计算反射、折射波的理论到时,并且根据理论走时与观测走时的拟合程度(试错法),调整不同层位的速度和深度,得到一个较为接近真实地壳上地幔顶部的二维速度模型.在试错法建立的速度模型基础上,本文提出一种对深地震测深剖面数据进行联合迭代反演(simultaneousiterativereconstructiontomography,简写为SIRT)的方法,并使用该方法对横跨龙门山断裂带中段的这条深地震宽角反射/折射剖面进行反演,以期获取该剖面的二维地壳速度结构.1反射/折射观测系统2008年汶川M该测线横穿龙门山断裂带中段,途经绵竹、茂县等极震区,向西北经过黑水,最终到达阿坝州.整条测线布设12个炮点,炸药量为600—2800kg(嘉世旭等,2014),采用井下爆破的方式,井深范围为15—61m.龙门山断裂带及其附近的炮点较密集,最小炮检距约为10km,而在远离龙门山断裂带的四川盆地和高原地区炮检距最大约为110km.图2给出了测线所采用的宽角反射/折射观测系统,通过该系统可以实现对地下介质进行多次覆盖和交叉采样,保证壳内波组可以连续对比追踪(张恩会等,2013),从而压制部分干扰波,提高地震记录的信噪比.本文选取12炮中震相清晰、容易识别且较为典型的6炮数据进行处理分析,具体列于表1.2震相的识别与识别整条测线的原始数据分为12个组段,首先,从其中挑选出同一炮点的数据,重新分组形成共炮点数据文件;然后,对每一炮的地震数据进行分析,识别出清晰可靠的震相;最后,进行走时拾取,得到观测走时数据.识别震相时,不仅要对多炮的对应波组进行对比,而且要分析震相的运动学和动力学特征.震相识别是所有工作的基础,它不仅是构建初始速度模型的依据,同时还对反演的可靠性有决定作用.因此,对于不太清晰的或者无法确定的震相应采取宁缺毋滥的原则予以剔除.最终分析识别得到的主要震相为:结晶基底折射波Pg,莫霍面反射波P2.1沉积相与震相波组的构造特征Pg震相是地壳内结晶基底的折射波.测线东侧(100—270km桩号段)追踪区间为10—100km,以清晰可靠的初至波出现(图3).在炮点附近视速度迅速增加,远离炮点时视速度逐渐趋于稳定,大约为5.8—6.2km/s,折合到时约为0—1.8s,且其能量随着炮检距的增大而逐渐减弱.这些特征体现了扬子地块稳定的结构和较厚的沉积层.测线中部(270—460km桩号段)位于龙门山断裂带,复杂的地下地质构造使得Pg波组的追踪范围大大缩小,仅在炮点附近可以追踪到,折合到时降至约0.1—0.4s.测线西侧(460—600km桩号段)的Pg波折合到时约为0.3—1.0s,显示了较浅的沉积层厚度和高速的地下速度结构.Pg震相波组的上述特征是由地表沉积盖层厚度和地壳结晶层顶部介质速度结构所决定的,而局部折合到时的超前、滞后一般是由地表局部的隆起和凹陷构造所引起.2.2莫霍界面的反射波p作为壳幔分界面莫霍面的反射波,由于上下岩性差异较大,速度差异也较大,P2.3上地禾地震追踪Pn波是来自上地幔的折射波,能够反映上地幔顶部的速度结构.由于该波组的能量较弱,连续对比追踪比较困难,除了SP8可以在炮点两侧追踪到外,其余炮仅能在单侧追踪.Pn波出现的最小距离为140km,最远可达300km,追踪长度最小仅有十几千米,显示了复杂的上地幔地质构造.以龙门山断裂带为界,测线西段川西北高原的折合到时在约190km处为0;而在东段的四川盆地,其折合到时在约260km处为0,显示了高速稳定的上地幔速度结构,视速度明显大于川西北高原.2.4地壳内反射波pP位于四川盆地的炮点SP1西侧和SP5东侧的震相显示,有3组反射波震相较明显,为P3速度模型的构造构建地下速度模型需要自上而下确定地壳各层的埋深和速度值.地表的速度值利用炮点附近的Pg波走时得到,界面深度及其上覆层介质的平均速度由x龙门山断裂带构造变形剧烈.横跨该断裂带的地壳和上地幔速度结构存在强烈的横向不均匀性,因此使用简单的一维地壳模型不可能很好地描述真实的地下速度结构.从图4可以看出,位于该区的SP5的右支(图4b),SP8的两支(图4c)和SP9的左支(图4d)均较难拟合,比较突出的是SP8左支和SP9左支的P首先利用每一炮的一维速度模型,使用插值法构造出一个较为粗略的二维初始速度模型;然后应用“剥皮法”,由浅入深地对沉积层、上地壳、中地壳、下地壳和上地幔顶部,采用射线追踪法计算反射、折射波的理论到时(ˇCervenue7b3,P2enˇcík,1984);最后根据理论走时与观测走时的拟合程度,调整不同层位的速度和深度,即可得到一个较为接近真实地壳上地幔的二维速度模型.4维地壳速度结构的反演.由于试错法存在耗时和不能给出定量误差的缺点,仅仅借助试错法,很难得到高精度的地壳速度结构.为此,本文将使用SIRT方法对观测走时进行反演.首先对二维初始模型进行网格化并建立射线路径上任意节点处的速度与周边网格节点速度之间的关系,然后推导出走时相对于网格节点速度的偏导数,最后建立走时反演方程组.4.1维网格化速度模型首先使用矩形网格对二维速度模型进行网格化.每一层介质为一个独立的网格单元,单元的左、右边界与研究区域的左、右边界重合,网格的上、下边界将所在层位的上、下两个界面包含在内.每个层位的矩形网格尺度可以不一致,网格的大小取决于二维初始模型速度的纵向和横向变化,速度梯度较大时采用较密的网格,反之,则采用较稀疏的网格.对之前得到的二维初始速度模型进行横向和纵向线性插值,可以得到用于反演的二维网格化速度模型.每个网格有4个节点,节点速度用v对网格内的射线节点采用反距离平方法进行线性插值.根据射线追踪结果可知,射线2在图5所示的矩形网格内没有节点,它在网格外的两个节点分别为V射线路径上任意相邻两点A和B的走时差t式中其中s对于每一条射线,沿射线路径周围节点速度的改变会导致射线终点走时的改变.用Δt上式中每一条射线的走时是该条射线上所有射线段走时的总和.将方程组(4)等号左右两边互换,可用矩阵形式表示为此即为用于反演的方程组,其右边为观测走时与理论走时的差,左边系数矩阵的每个元素是理论走时对模型节点速度的偏导数.在完成射线追踪后,系数矩阵可由式(1),(2)和(3)求出;等号左边的矢量(Δv4.2维速度模型反演大型稀疏矩阵方程组(5)可通过许多方法求解,例如,直接解法,雅克比共轭梯度法,不完全柯列斯基(Cholesky)共轭梯度法,预条件共轭梯度法等.这些方法有的要求高内存,有的要求计算速度快,有的可能存在收敛和稳定性问题.本文将采用层析成像理论中的代数迭代法和联合迭代法(SIRT).利用代数迭代法计算一条射线之后就可以将走时残差所对应的速度变化分配到射线路径上的相关节点;利用SIRT法则是计算完所有射线后将走时残差所对应的速度变化量分配到相关的节点上.由于SIRT法要求模型更新的次数较少,所以比较稳健.本文主要使用SIRT进行迭代反演.将走时拟合后得到的二维速度模型作为初始反演模型,使用上述反演方法对2668个走时数据(见表2)进行反演,可以进一步缩小理论走时与观测走时之间的误差.Pg波组的到时记录用来约束地壳顶部速度结构,P图7给出了反演前与反演迭代第11次的对比结果,可以看出,速度模型有了较大的改善.SP1(图7a)和SP13(图7f)炮点模型反演前实测走时与理论走时的残差较小,拟合效果比较好;反演后速度模型变化不大.SP5(图7b)反演前左支走时残差较小,右支P5地壳厚度和速度结构通过反演得到了二维P波速度结构剖面图,如图8所示.可以看出,整条测线下方的地壳结构以龙门山为界,两侧的厚度和速度结构均存在较大的差异.测线东段(100—270km桩号段)位于四川盆地,地壳厚度约为43—44km;测线西段(460—600km桩号段)位于川西北高原,地壳厚度约为57—58km;测线中段(270—460km桩号段)位于松潘—甘孜构造转换带,强烈的构造活动使得该区界面起伏明显并伴有错断,地壳厚度介于东、西两段之间.5.1龙日坝内盖层沿着测线自西向东,东、西两段沉积盖层的特征存在巨大差异.川西北高原沉积盖层的厚度为2.5—3.0km,介质速度为3.60—5.25km/s;沿测线向东经过龙日坝后,沉积盖层逐渐变薄;而从桩号280km左右开始,盖层厚度急剧增大;当测线到达德阳市时,盖层厚度达到5.5—6.5km,并延伸至遂宁,介质速度为3.55—5.34km/s.龙门山断裂带两侧沉积盖层厚度的差异与扬子板块和松潘—甘孜板块的地质年龄相吻合.5.2上地壳结构上地壳位于沉积盖层G与C5.3中地壳结构中地壳位于C5.4下地壳结构下地壳位于C5.5上地潭顶部地下水震莫霍面作为壳幔边界,其深度沿测线变化很大,由川西北高原的53—55km,自西向东经龙门山褶皱带,逐渐抬升至四川盆地的43—44km.上地幔顶部的速度Pn在横向上存在比较剧烈的变化,川西北高原速度约为7.95km/s,而四川盆地高达8.10km/s.将汶川地震的震源垂直投影至测量剖面上,如图8红色星形所示,其在地面的投影坐标为(103.7830°E,31.6879°N),震源深度为(15.5±0.3)km,桩号为332.9km,可见其位于上地壳底部的地层抬升区.6下地壳厚度及厚度本文研究的测线横跨西部松潘—甘孜地块、中部龙门山断裂带和东部扬子地块,在地壳厚度、壳内速度结构方面存在着巨大的差异,低速层、高速体的分布特征也有其独特的地质成因.测线剖面显示:沉积盖层在东部四川盆地的分布明显厚于中部褶皱带和西部高原,而中部褶皱带部分地区出现基岩裸露,表明扬子板块的地质年龄要老于松潘—甘孜板块,中部褶皱带处于活跃的地质构造期;地壳内构造转换带两侧的地层分界面近于水平层状分布,转换带西侧的中、下地壳内各存在一个层间速度间断面,而东侧则不存在类似的间断面.对应上、中、下地壳的分界面,西侧比东侧的速度要小,而且自上而下逐渐增大,下地壳底部莫霍面的速度西侧比东侧要低0.3km/s.在厚度方面,西侧比东侧的上、中、下地壳厚度分别厚3—7km,6—8km,6—9km,可见地壳的增厚主要由被改造的中下地壳介质形变增厚所致