2. 德国GFZ地学研究中心深部探测部, 波茨坦 14473, 德国
2. Geophysical Deep Sounding Sector, GFZ, Telegrafenberg, Potsdam 14473, Germany
龙门山地区是南北地震带的核心部分,龙门山断裂是汶川地震的发震断裂,在青藏高原向东北方向挤压的背景下龙门山地区还发生了2013年芦山地震(易桂喜等,2016)和2017年九寨沟地震(杨宜海等,2017)等7级以上的大地震,目前急需解决龙门山断裂的活动规律和结构特征,包括断裂带两侧和深部的地壳结构特征,如上地壳彭灌杂岩的赋存状态,中地壳滑脱层与龙门山断裂的关系,下地壳的低速体分布范围以及它对上地壳活动的影响等.汶川地震前该地区的地球物理探测剖面屈指可数(Wang et al., 2007),然而汶川地震发生以后,该地区获得了大量的地质和地球物理成果(Zhang et al., 2009, 2010; Robert et al., 2010; Li et al., 2012, 2014a, b;Cook et al., 2013; Guo et al., 2013; Liu et al., 2014; Wang et al., 2015; Xu et al., 2015).不过对于龙门山断裂的深部延展目前还没有统一的认识,对于其深部驱动机制更是争论不休.
密集台站和高密度的近震资料是该地区高精度层析成像的数据基础,近震的优点首先在于分布广,除了沿着龙门山断裂两侧分布,松潘地区常有4级以上的地震活动(徐晶等,2017),四川盆地内部也能记录到3级的地震.其次地震深度分散,从浅部地表到莫霍面以上均有分布,不过由于龙门山断裂的活动,龙门山地区的近震主要集中在5~30 km.还有同一个地震可以有不同的震相记录,如果台站分布比较合理将可以记录到各种不同类型的震相,包括台站正下方的上行直达波,远离台站的下行直达波,以及更远但却能先记录到的折射波,在某些震中距还能记录到莫霍面反射波.S波也存在上述四种震相,由于到时比P波晚,难于识别.不过其频率特征比较明显,通过互相关分析也能获取走时信息.
通过对密集台阵的近震进行分析,可以获取台阵下方密集互相穿插的射线分布,不同于平行入射的远震射线,近震的这种射线分布更有利上地壳部分的速度结构成像, 有更高的水平分辨率.近震层析成像可以很好地拟合近震走时和台站下方的速度结构分布,这些方法一般都自带三维速度结构重定位模块,对地震位置进行四维空间搜索的同时对速度结构进行反演.与三维速度结构反演耦合的近震层析方法中三维定位的一般要求近震位置比较精确,在台站分布没有完全包裹地震或射线分布不均匀的条件下,反演结果往往很不稳定.
由于震源深度的增加、速度的减小以及发震时刻的提前都会导致走时增加,耦合的反演方法没有办法完全区分这些因素的叠加,因此需要在三维耦合近震层析之前就确定高精度的地震位置和时间,从而在三维反演过程中减小地震位置的不确定性或者不对震源进行反演.这样的策略有时能取得更稳定的结果,结果精度的关键还在于震源位置的精确性.常用基于绝对走时的近震层析方法包括Tomo3D(Zhao et al, , 1992)、Lotos(Koulakov, 2009)、FMM3d(Rawlinson and Urvoy, 2006).而双差层析tomoDD(Zhang and Thurber, 2003)的定位方法则基于射线路径差异不大的地震簇.
三维定位之前的地震定位一般采用一维速度模型的初步定位,依据是否对速度模型进行耦合分为两类,耦合的方法包括Velest(Kissling et al., 1994),不耦合的方法则更多如locsat(Bratt and Bache, m 1988)、hypoinverse等.反演方法大多为采用Greg反演,也有采用全局搜索的方法或者走时查找表的方法等.对于一维定位速度模型需要进行展平变换来弥补地球曲率对大震中距的畸变作用的影响,而远震则可能需要采用球坐标.展平变换的优点是走时和距离微分方程简单,但是接近地心的震相深度和速度同时接近无穷大,造成内核震相存在较大的误差是该方法的主要缺点.相对而言,球坐标则只存在球心一个奇点需要特殊考虑,其他位置均采用真实速度和深度,但是其微分方程比较复杂,仅对某些特殊的层位速度分布有解析解.采用解析解的优点是各物理量意义明确,对震源参数,射线参数和速度模型的导数可解析取得,而无需耗时的差分计算.这些导数是反演方法所必须的信息,同时导数优化也是各种线性化反演快速收敛的主要途径.
1 解析定位方法
依据地球速度结构的Bullen模型(Bullen and Bot, 1985), 在(r1, r2)两个深度构成的层位内速度随深度z的变化满足v=ar1-b.其中
某一层位内部的角距离和时间为{x, t}在r1, r2位置的函数值之差,相应层位的导数也是导数公式在r1、r2位置的函数值之差.当该层位包含临界深度时f=0, 临界深度的角距离和时间函数值均为0,导数值也相应全为0,只需要计算层位顶部的函数值和导数值.整体射线的角距离和时间为各层位函数值的求和,导数值也是如此,对于上下层的中间节点导数需要叠加,而顶底节点则没有叠加的问题.当层位内部包含震源或检波点时,需要计算射线对震源深度或检波点速度的导数,供以后反演震源位置和地表静校正使用.p=0时角距离函数可取极限π/2b, 球坐标唯一的奇点是z=0且p=0,此时射线经过球心,距离时间函数均取0, 导数值也相应全为0,相当于整个地球的临界深度.
对不同震相进行编码需要指定开始深度,结束深度和波的类型,同时限制射线所在的深度范围.图 1是检波点和震源深度为0的情况时,不同射线参数对应iasp91模型不同震相的距离时间解析计算结果,并与Taup软件包(Crotwell et al., 1999)的计算结果进行差值的图示.结果表明解析计算的距离误差小于0.035°, 时间误差小于0.6 s.除了散射波不能解析计算外, 其他波均可通过编码表示, 因此对于全球主要震相的走时和距离计算可以用这种方法完成,并可以通过通用的导数优化方法进行射线追踪,由此可以构建相应的定位算法.
一维射线追踪可以定义为解方程X(r, v, p, d, z, w, v0)=X0中的p射线参数,其中r、v是确定速度模型的节点序列,d为震中距,z为震源深度,w和v0为检波点位置和速度, X为各层xi求和的结果.通过导数公式可以看到球坐标条件下各层存在t′i(xi)=p,最终求和也存在T′(X)=p, 由此可以建立射线追踪后的时间对模型和震源位置以及地表速度的导数矩阵,必要时可以附加全部为1的矩阵作为震源时间校正或台站静校正时间量的导数(Qian et al., 2016).
求解射线参数的过程是单变量反演,可以采用最速下降法与黄金分割搜索法结合的方式,能够较快的收敛在给定的取值范围内.而耦合模型的定位反演是多变量优化,由于存在诸多的约束条件,适合采用含导数信息的信赖域优化方法.
2 龙门山近震定位中国地质科学院地质研究所与德国GFZ合作于2012年6月至2013年10月在龙门山地区布设了80个宽频地震流动台站,台站位置如图 2所示.台站范围覆盖了龙门山断裂的主体部分,并记录到了芦山地震前后龙门山地区10~20 km间距网格内的地震活动的情况,为探索龙门山地区的地震孕育机制奠定了一定的数据基础.本次工作采用的近震定位方法是如前所述的球坐标解析追踪附带导数矩阵信息和限制条件的信赖域二次规划优化方法,与Velest方法对比优点在于球坐标可用于远震,处理同一地震的P、S波的多种震相并优化P和S速度模型,并用近地表速度一致性校正代替了台站静校正,结果更符合射线追踪原理.与所有线性化优化方法一样,多解性无法避免,因此反演中选择了更多的地质地球物理约束加入二次规划的限制条件中,保证结果的合理性.
远震由于射线路径较长, 在高精度近震层析中无法运用.而在远近震联合层析中,远震也只是使用相对残差数据.此次工作中,采用了远震近震化处理的方式利用远震绝对走时数据,即按标准地球模型计算地震到达某一深度的参考深度走时和射线穿刺点,然后将台站走时减去参考深度走时获取类似近震的走时数据和穿刺点震源,再通过近震重新定位处理可以获得定位后的位置和震源时间校正量,参考深度外的走时异常通过震源时间校正的方式部分剔除.这种方法的优点是能明显增加深部速度结构的约束,缺点是远震近震射线路径的差异导致外围的走时异常不能完全消除.今后可以发展基于近震三维射线追踪联合球坐标解析追踪的远震定位技术,为远震的纳入提供更为直接的方法,但仍然存在较大的误差,最终的方案将是自适应网格划分的三维射线追踪方法既能保证精度其计算速度也能接受.
台阵观测期间台阵记录到明确P和S走时的2级以上近震814个,射线53618条,通过远震近震化处理获得射线16781条.模型建立的原则是有实际观测走时的地震台站之间存在射线,由Pn和Sn走时曲线的斜率可获取Moho面相应的波速和深度,然后逐步建立初始模型.通过解析射线追踪方法初步定位可以取得误差较小的震源位置,同时还包括震源时间校正,一维速度模型校正和台站近地表速度校正(图 3).
经过地震初步定位获得一维平均速度结构和大致相适应的地震位置和走时,可以作为三维近震层析的输入.采用的近震层析方法是带三维精定位功能的层析方法Lotos,其中网格划分根据射线密度进行,射线追踪采用弯曲法.依据初步定位的数据分布,按照5 km网格进行射线密度统计,采用检测板方法对本次射线分布的分辨能力进行检测(如图 4).可见射线分布对于浅部的大的块体能够有效地分辨,而到30 km以下分辨能力稍差.
经过层析成像反演,获得的走时残差数据方差改进为28.3%,根据分辨较好的切片和剖面位置,获得了P波和S波速度扰动的图像(图 5).整体上看P波速度,20 km深度以上四川盆地以低速物质为主,松潘地区浅表为高速夹部分山间盆地低速物质,而上地壳也表现为低速.二者之间的分界龙门山断裂带浅部表现为高速,并且明确地分割了盆山的结构体系.20 km深度以下,断裂带的分割作用体现不明显,似乎可见彭灌杂岩体的根部,汶川—茂县断裂控制了岩体的边界.此深度上代表扬子板块的高速体逐渐取代四川盆地的低速物质,而松潘地区中下地壳仍广泛发育低速体,在C、D剖面可见低速物质向盆地内部的扩展,构成龙门山山根断裂系的滑脱带.
三维重新定位后浅震主要分布在龙门山断裂和四川盆地的正常地壳的交界处.邛崃附近存在切割龙门山的近东西向31.25°N位置低速分布,此处可以作为汶川地震应力向南释放被吸收的重要标志,浅部的地震还沿着这一通道拐入松潘地区,并在松潘低速体的周围有规律的展布.
S波速度与P波速度相一致的特征有:(1)紧邻龙门山断裂都江堰附近分布的四川盆地浅部低速物质; (2)深度15 km以上走向近东西31.25°N位置切割松潘地体的低速条带;(3)断裂带下方中下地壳普遍分布的低速体.
从速度比看(图 6),浅部断裂带体现为高速度比,盆地和山区均为低速度比;中地壳在汶川一带存在山区向盆地穿插的一条近东西向低速度比条带,20 km以下消失;深部下地壳盆地以高速度比为主,高原内部则为低速度比物质,而汶茂断裂和映秀—北川断裂带的深部是延续盆地的高速度比分布.
对比前人关于龙门山断裂带探测的结果,结合本次层析工作的探讨和体会包括:
(1) 与穿越龙门山断裂带的电性结构对比(Wang et al., 2014)发现电性结构存在盆地浅部位低电阻率,深部为高电阻率,而松潘地区刚好反过来.这与本次层析结果比较相似,四川盆地呈现浅部低速深部高速的特征,而松潘地区刚好相反,特别是S波速的几条纵剖面均能看到这种现象.
(2) 与双差层析结果(邓文泽等,2014)对比发现,穿过汶川的剖面上龙门山断裂带附近存在局部突出的高速体,与本次层析获取的结果基本一致,可能与基底岩石出露或彭灌杂岩体的分布有关.
(3) 关于下地壳流的数值模拟(Li et al., 2017)表明龙门山地区可能存在下地壳粘弹性物质的快速运动,从层析结果看,松潘地区中下地壳存在大量的低速物质,存在向东扩展的趋势,但由于与上地壳低速体区别不大,还不能完全肯定是下地壳流引起了松潘地区低速物质的分布,需要更多的证据来证明上下地壳运动的解耦作用.
(4) 穿过龙门山断裂带的人工地震反射剖面(Lu et al., 2016)推测20 km深度存在断裂带统一收束的滑脱面,从速度比剖面可以明显看到20 km深度存在上地壳速度比由浅部的低值向深部高值的变化,表明该滑脱带可能的确存在,由于数据范围所限在层析结果上松潘地区滑脱层的分布不清楚,与下地壳低速体的关系需要进一步研究.
(5) 球坐标解析定位求解导数矩阵的过程公式比较复杂,但比起差分方法多次模型计算速度更快,精度更高,是值得尝试的快速导数优化的方式.类似于气象模式计算的伴随模式,计算因变量时同时计算因变量对所有自变量的导数是必然要求,目前已经发展为根据复杂的主程序计算逻辑来自动生成伴随矩阵程序.
(6) 层析成像是解决深部速度结构分布有效的方法之一,然而受到地震和台站分布的限制,其分辨能力一直有待提高,特别是较小范围地区的层析很难获得深部的射线分布.考虑到使用相对残差的诸多缺陷,如何利用穿过深部的远震射线是值得探讨的课题.主要的思路是远震和近震网格剖分采用同一体系,远震网格较大,由于采用解析计算并不影响精度,而近震网格较精细,可以采用数值计算.
综合层析图像结果和以上的讨论,获得了关于龙门山地区盆山结构框架的主要结论包括:
(1) 龙门山地区阻挡了青藏高原向东的扩展,导致松潘地区地壳增厚,速度比降低,从而导致中下地壳存在大量的低速物质活动,并且有一部分物质在汶川一带的中地壳进入盆地.
(2) 龙门山断裂带的两条主干断裂可能收束在20 km深度的滑脱带上,汶茂断裂倾角稍高,控制了彭灌杂岩的位置,彭灌杂岩的根部在收束的滑脱带上尖灭.
(3) 四川盆地深部为扬子板块高速物质控制,高原内部深部则以低速物质为主.浅部受盆地沉积和松潘地体影响,断裂带两侧均呈现低速特征,而断裂带浅部由于夹持了古老的杂岩体而呈现高速特征.
由于本次工作台站平均间距超过15 km,虽然密集的近震分布弥补了龙门山断裂带地区的分辨率,但在山区和平原地区分辨率仍有待提高,还需要更为密集的台站观测和跨过断裂带的人工地震探测剖面予以证实.
致谢感谢GFZ地震仪器联合会和北京港震机电公司提供的探测仪器和技术支持,感谢冯梅研究员关于本台阵的噪音成像处理工作,其获得的面波三维速度分布对层析处理有一定的指导作用,感谢李海兵研究员对项目工作的支持和帮助,同时也感谢评审提出的宝贵意见.
Bratt S R, Bache T C.
1988. Locating events with a sparse network of regional arrays. Bulletin Seismological Society America, 78(2): 780-798.
|
|
Bullen K E, Bot B A. 1985.
An Introduction to the Theory of Seismology. New York: Cambridge University Press.
|
|
Cook K L, Royden L H, Burchfiel B C, et al.
2013. Constraints on Cenozoic tectonics in the southwestern Longmen Shan from low-temperature thermochronology. Lithosphere, 5(4): 393-406.
DOI:10.1130/L263.1 |
|
Crotwell H P, Owens T J, Ritsema J.
1999. The TauP toolkit: flexible seismic travel-time and ray-path utilities. Seismological Research Letters, 70(2): 154-160.
DOI:10.1785/gssrl.70.2.154 |
|
Deng W Z, Che J H, Guo B, et al.
2014. Fine velocity structure of the Longmenshan fault zone by double-difference tomography. Chinese Journal of Geophysics (in Chinese), 57(4): 1101-1110.
DOI:10.6038/cjg20140408 |
|
Guo X, Gao R, Keller G R, et al.
2013. Imaging the crustal structure beneath the eastern Tibetan Plateau and implications for the uplift of the Longmen Shan range. Earth Planet Science Letter, 379(5): 72-80.
|
|
Kissling E, Ellsworth W L, Eberhart-Phillips D, et al.
1994. Initial reference models in local earthquake tomography. Journal of Geophysical Research: Solid Earth, 99(B10): 19635-19646.
DOI:10.1029/93JB03138 |
|
Koulakov I.
2009. LOTOS code for local earthquake tomographic inversion: benchmarks for testing tomographic algorithms. Bulletin of the Seismological Society of America, 99(1): 194-214.
DOI:10.1785/0120080013 |
|
Li H, Wang H, Xu Z, et al.
2012. Characteristics of the fault-related rocks, fault zones and the principal slip zone in the Wenchuan Earthquake Fault Scientific Drilling Project Hole-1 (WFSD-1). Tectonophysics, 584(1): 23-42.
|
|
Li H B, Xu Z Q, Niu Y X, et al.
2014a. Structural and physical property characterization in the Wenchuan earthquake Fault Scientific Drilling project—hole 1 (WFSD-1). Tectonophysics, 619-620: 86-100.
DOI:10.1016/j.tecto.2013.08.022 |
|
Li Y J, Liu S F, Chen L W, et al.
2017. Mechanism of crustal deformation in the Sichuan-Yunnan region, southeastern Tibetan Plateau: Insights from numerical modeling. Journal of Asian Earth Sciences, 146: 142-151.
DOI:10.1016/j.jseaes.2017.05.018 |
|
Li Y Q, Jia D, Wang M M, et al.
2014b. Structural geometry of the source region for the 2013 MW6.6 Lushan earthquake: Implication for earthquake hazard assessment along the Longmen Shan. Earth and Planetary Science Letters, 390: 275-286.
DOI:10.1016/j.epsl.2014.01.018 |
|
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 Geosciences, 7: 361-365.
DOI:10.1038/ngeo2130 |
|
Lu R Q, He D F, Xu X W, et al.
2016. Crustal-scale tectonic wedging in the central Longmen Shan: Constraints on the uplift mechanism in the southeastern margin of the Tibetan Plateau. Journal of Asian Earth Sciences, 117: 73-81.
DOI:10.1016/j.jseaes.2015.11.019 |
|
Qian H, Mechie J, Li H B, et al.
2016. An approach to jointly invert hypocenters and 1D velocity structure and its application to the Lushan earthquake series. Journal of Seismology, 20(1): 213-232.
DOI:10.1007/s10950-015-9521-0 |
|
Rawlinson N, Urvoy M.
2006. Simultaneous inversion of active and passive source datasets for 3-D seismic structure with application to Tasmania. Geophysical Research Letters, 33(24): L24313.
DOI:10.1029/2006GL028105 |
|
Robert A, Zhu J, Vergne J, et al.
2010. Crustal structures in the area of the 2008 Sichuan earthquake from seismologic and gravimetric data. Tectonophysics, 491(1-4): 205-210.
DOI:10.1016/j.tecto.2009.11.010 |
|
Wang C Y, Han WB, Wu J P, et al.
2007. Crustal structure beneath the eastern margin of the Tibetan Plateau and its tectonic implications. Journal of Geophysical Research: Solid Earth, 112(B7): B07307.
DOI:10.1029/2005JB003873 |
|
Wang S J, Wang F Y, Zhang J S, et al.
2015. A study of deep seismogenic environment in Lushan MS7.0 earthquake zone by wide-angle seismic reflection/refraction profile. Chinese Journal of Geophysics, 58(5): 474-485.
DOI:10.1002/cjg2.2015.58.issue-5 |
|
Wang X B, Zhang G, Fang H, et al.
2014. Crust and upper mantle resistivity structure at middle section of Longmenshan, eastern Tibetan plateau. Tectonophysics, 619-620: 143-148.
DOI:10.1016/j.tecto.2013.09.011 |
|
Xu J, Shao Z G, Liu J, et al.
2017. Analysis of interaction between great earthquakes in the eastern Bayan Har block based on changes of coulomb stress. Chinese Journal of Geophysics (in Chinese), 60(10): 4056-4068.
DOI:10.6038/cjg20171031 |
|
Xu X, Keller G R, Gao R, et al.
2015. Uplift of the Longmen Shan area in the eastern Tibetan plateau: an integrated geophysical and geodynamic analysis. International Geology Review, 58(1): 14-31.
|
|
Yang Y H, Fan J, Hua Q, et al.
2017. Inversion for the focal mechanisms of the 2017 Jiuzhaigou M7.0 earthquake sequence using near-field full waveforms. Chinese Journal of Geophysics (in Chinese), 60(10): 4098-4104.
DOI:10.6038/cjg20171034 |
|
Yi G X, Long F, Vallage A, et al.
2016. Focal mechanism and tectonic deformation in the seismogenic area of the 2013 Lushan earthquake sequence, southwestern China. Chinese Journal of Geophysics (in Chinese), 59(10): 3711-3731.
DOI:10.6038/cjg20161017 |
|
Zhang H, Thurber C H.
2003. Double-difference tomography: the method and its application to the Hayward fault, California. Bulletin of the Seismological Society of America, 93(5): 1875-1889.
DOI:10.1785/0120020190 |
|
Zhang Z J, Wang Y H, Chen Y, et al.
2009. Crustal structure across Longmenshan fault belt from passive source seismic profiling. Geophysical Research Letters, 36(17): L17310.
DOI:10.1029/2009GL039580 |
|
Zhang Z J, Yuan X H, Chen Y, et al.
2010. Seismic signature of the collision between the east Tibetan escape flow and the Sichuan Basin. Earth and Planetary Science Letters, 292(3-4): 254-264.
DOI:10.1016/j.epsl.2010.01.046 |
|
Zhao D P, Hasegawa A, Horiuchi S.
1992. Tomographic imaging of P and S wave velocity structure beneath northeastern Japan. Journal of Geophysical Research: Soild Earth, 97(B13): 19909-19928.
DOI:10.1029/92JB00603 |
|
邓文泽, 陈九辉, 郭飚, 等.
2014. 龙门山断裂带精细速度结构的双差层析成像研究. 地球物理学报, 57(4): 1101–1110.
DOI:10.6038/cjg20140408 |
|
徐晶, 邵志刚, 刘静, 等.
2017. 基于库仑应力变化分析巴颜喀拉地块东端的强震相互关系. 地球物理学报, 60(10): 4056–4068.
DOI:10.6038/cjg20171031 |
|
杨宜海, 范军, 花茜, 等.
2017. 近震全波形反演2017年九寨沟M7.0地震序列震源机制解. 地球物理学报, 60(10): 4098–4104.
DOI:10.6038/cjg20171034 |
|
易桂喜, 龙锋, VallageA, 等.
2016. 2013年芦山地震序列震源机制与震源区构造变形特征分析. 地球物理学报, 59(10): 3711–3731.
DOI:10.6038/cjg20161017 |
|