地球物理学报  2014, Vol. 57 Issue (8): 2555-2572   PDF    
一种地方与区域地震震源机制反演技术:广义极性振幅技术(一)——原理与数值实验
严川, 许力生    
中国地震局地球物理研究所, 北京 100081
摘要:基于数十年来已有的研究进展,提出了一种地方和区域地震震源机制反演技术——广义极性振幅技术(GPAT),并通过一系列数值实验检验了这种技术的可行性和抗干扰能力.首先,从地震波场概念出发,利用P波初动极性与广义震相振幅构建矢量,建立反演系统,并给出求解技术;然后,考虑影响反演结果的各种因素,包括台站布局、台站数目、随机噪声、震中位置误差、震源深度误差和速度模型误差,分别进行了单一因素影响测试;最后,同时考虑各种因素进行了综合测试.实验结果表明,GPAT是可行的,具有良好的抗干扰能力.需要强调的是,在众多影响因素中,速度模型误差对反演结果的影响最大.
关键词地方与区域地震     震源机制反演     广义极性振幅技术     数值实验    
An inversion technique for the mechanisms of local and regional earthquakes:generalized polarity and amplitude technique (Ⅰ)——Principle and numerical tests
YAN Chuan, XU Li-Sheng    
Institute of Geophysics, China Earthquake Administration, Beijing 100081, China
Abstract: An inversion technique for focal mechanisms of local and regional earthquakes, named Generalized Polarity and Amplitude Technique(GPAT), is proposed based on the progress already made over decades, and its feasibility and capability of resisting disturbance are tested by a series of numerical experiments. At first, an inversion equation system is constructed with the vectors consisting of the polarities of the P first motion and the maximal amplitudes of generalized phases based on the concept of wave field, and a technique for solving the equation system is specified. Then, a series of experiments are conducted for the effect of various factors on inversion results, including distribution of observation stations, number of the stations, random noises, errors in epicenter location, errors in focal depth and errors in velocity model. Finally, a special experiment is designed for comprehensive effect of all the factors. The experiments suggest that the GPAT is feasible, and has rather good capability of resisting disturbance. It must be stressed that the error in velocity model has the strongest effect on inversion result compared with other factors.
Key words: Local and regional earthquakes     Focal mechanism inversion     Generalized polarity and amplitude technique     Numerical test    
1 引言

地震的发生是应力作用的结果,震源机制类型取决于应力作用方式.因此,借助于震源机制可以认识震区的应力状态.所以,确定震源机制的工作显得十分重要.

地震矩张量的引入使地震产生的远场位移与描 述震源机制的地震矩张量之间成为线性关系(Gilbert,1973; Aki and Richards, 1980; Lay and Wallace, 1995).如果地震产生的位移与描述地震波传播路径效应的格林函数已知,则很容易通过波形反演得到地震矩张量,即震源机制(Aki and Richards, 1980; Lay and Wallace, 1995).然而,波形反演技术只适用于远场记录或较大的地震(Dziewonski et al., 1981),最多应用到区域地震记录或中等大小的地震(Helmberger and Engen, 1980; Zhu and Helmberger, 1996).

尽管波形反演方法可以推广到区域地震记录或中等大小的地震,但应用的条件是比较苛刻的. Helmberger和Engen(1980)最早发现,在区域距离上(震中距为2°~12°的范围)存在一种广义体波(P 波之后面波之前的波列,称之为Pnl波)对地壳结构 相对稳定,可用于区域地震震源机制的反演.事实表 明,的确如此(Wallace et al., 1981; Liu and Helmberger, 1985). 不过,Dreger和Helmberger(1991)后来发现,当所用台站记录之间的Pn波和Sn波到时变化分别在2%和3%以内且波形相似时,地壳模型可以被认为相对平直,此时利用Pnl波确定震源机制才是可行的.据此,Dreger和Helmberger(1993)正式提出一种区域距离上的广义体波反演方法(下文称Pnl波方法),并称可以利用稀疏台网和单台反演震源机制.作为对Dreger和Helmberger(1993)方法的发展,Zhao和Helmberger(1994)提出了称之为CAP(Cut and Paste)的方法,与前者相比,后者增加了面波的使用,即将整个地震记录分为Pnl波和面波两部分并进行分离(即Cut),然后根据相位校正再拼接(即Paste).现阶段使用较为广泛的是经Zhu和Helmberger(1996)改进后CAP的新版,通常被用来确定区域和地方范围内3~5级地震的震源机制.由此可见,虽然波形反演方法可以推广到区域地 震记录或中等大小的地震,但前提是当地的地壳结构要相对平直.另外,CAP方法中的拼接环节容易引入人 为误差.同时可以看出,无论是Pnl波方法还是CAP方法,他们对于小震或微震的震源机制仍无能为力.

自地震产生的初动极性被观测证实呈现四象限分布以来,P波初动符号就成为确定震源机制的重 要信息(Kasahara,1963; Lindh et al., 1978; Warren,1979; 许忠淮等,1983;Reasenberg and Oppenheimer, 1985; Saikia and Herrmann, 1985; Natale et al., 1991; Schwartz,1995; Hardebeck and Shearer, 2002; 俞春泉等,2009).但是,利用P波初动确定震源机制需要方位覆盖好且稠密的台站分布,因此,仅利用初动方向确定震源机制的方法(初动法)的应用受到很大限制.为了克服初动法的缺点,有人引入了 振幅信息(Kisslinger,1980; Kisslinger et al., 1982; 梁尚鸿等,1984; 吴大铭等,1989; Snoke et al., 1984; Snoke,1989; Ebel and Bonjer, 1990; Nakamura et al., 1999; Hardebeck and Shearer, 2003; 刘杰等,2004; Cesca et al., 2006).不妨称为初动-振幅方法.有些人把振幅比作为观测资料(Kisslinger,1980; Kisslinger et al., 1982; Julian and Foulger, 1996),有些人把绝对振幅作为观测资料(Ebel and Bonjer, 1990; Godano et al., 20092011).与初动方法相比,初动-振幅法使用的信息量较大,因此,大大降低了对台站密度和方位覆盖的要求.不过,有研究表明,只有振幅信息没有受到“污染”的情况下,它才对改善震源机制具有积极意义(Hardebeck and Shearer, 2003).但不管怎样,Snoke 等(1984)Snoke(1989)作了一件非常重要的工作,即考虑P、SV、SH的极性和振幅比SV/P,SH/P,SV/SH等信息的任意组合,编写了一个程序FOCMEC,并一直在不断改进.这是至今对初动-振幅法的最全面的总结.统略初动-振幅方法,我们发现,由于问题本身的复杂性,至今没有一种公认的方法能像远场波形反演方法那样被广泛接受和使用.

综合考虑现有方法的优点和缺点,我们尝试一种新的非线性反演方法,不妨称之为广义极性振幅法(GPAT).之所以称之为广义极性振幅法,是对现有初动-振幅类方法的一般化,因为我们不但考虑P波初动的极性,同时也考虑最大振幅的极性;更为重要的是,我们所说的震相不再局限于P波、S波或是面波,只要在记录上表现出明显的震相特征,无论它叫什么,其最大振幅即可被采用.这一点对于利用震相极其复杂的地方和区域地震记录反演震源机制非常关键.同时,与现有方法不同,我们把观测波场与合成波场的相似性作为求解的目标,这一点对于充分利用最大振幅的极性及其大小的作用至关重要. 2 原理 2.1 目标函数的表达

一般地,地震产生的地表位移可以表示为

这里,j=1,2,3,分别代表东西、北南和垂直(上下)三个分向,而uj可以是P波、S波、面波和其他震相.不失一般性,我们可以将uj表示为

上式中,上角标P、S、F和O分别表示P波、S波、面波和其他震相.尽管GPAT并不局限于某个或某几个震相,但为了讨论问题的方便,我们仍只考虑P波、S波和面波的情形.

地震波场在时间和空间上都是连续的,但GPAT只关心有限的观测点记录的P波、S波和面波带有极性的最大振幅以及P波初动极性,并据此分别建立两个行矢量

式(3)和式(4)中,i=1,2,…,N,表示不同的观测点,j=1,2,3,分别代表东西、北南和垂直(上下)三个分向,k=1,2,3分别代表直达P波、S波和面波; a为矢量 v 1的元素,表示绝对值最大但带有极性的振幅,b为矢量 v 2的元素,代表P波初动极性(b=+1代表初动向东、北和上,b=-1代表初动向西、南和下,而b=0则代表初动方向不明确).

仿照式(3)和(4),对于给定的震源产生的波场可以写出类似的矢量

为了建立目标函数的方便,令

上式中,w为振幅信息与P波初动信息的相对权重. αm为观测矢量元素,而βm为合成矢量元素,T表示行矢量向列矢量的转置.

可以想见,如果给定震源的标量地震矩、震源位置和震源机制与实际震源的完全相同,则有

而如果仅有标量地震矩不同但其他参数相同,则二者相似,那么 η ′和η 线性相关.令ρ为二者的相关系数,则

上式中,和分别表示αm和βm的平均值.

相对权重w依赖于aijk和bijk,但如何为w取值成为一个问题.为了阐明这个问题,我们定义2个列向量 X和Y,

并评价其相关系数,

不难想见,κ的值依赖于w.如果w=0,则κ=1,表明只有振幅发挥作用; 如果w→+∞,则κ→-1,表明振幅几乎不发挥作用;那么,只要知道振幅和极性发挥相同作用时w的取值,便可以根据需要为w赋值.从(13)可以看出,当κ=0时,振幅和极性发挥相同作用,即w0=w(κ=0),这个值的大小取决于具体的数据矢量 X和Y .例如,图 1展示了某种情况下权重w与相关系数κ的关系曲线.在这种情况下,w0=4.045,若要极性发挥较大作用,则w取值要大于4.045;若要振幅发挥较大作用,则w取值要小于4.045.

根据上文,βm依赖于待定事件的震源机制和震源位置.若震中位置确定,则依赖于震源机制和震源 深度,因此,ρ必然是震源机制和震源深度的函数,即

上式中,φ,θ,λ,d分别表示断层走向、倾角、滑动角和震源深度.因此,我们的反演问题就是ρ→1或者Δρ=1-ρ→0时的非线性求解问题. 2.2 求解技术

由(10)和(14)式可以看出,目标函数ρ是待求地震的震源机制和震源深度的非线性函数.求解这个非线性问题的技术多种多样,如流行的遗传算法、模拟退火法等,但经试验后我们认为网格搜索法仍是解决这个问题的最佳选择.这也许是为什么很多人在解决这类问题时都采用这种技术的原因(Langston 1982;; 许忠淮等,1983; Reasenberg and Oppenheimer, 1985; Zhao and Helmberger, 1994; Zhu and Helmberger, 1996; Hardebeck and Shearer, 2002; 俞春泉等,2009).遗传算法、模拟退火法等似乎能提高求解效率,但这个高效率是以恰当的控制参数为前提的.然而,这些控制参数如何设置通常是使用者棘手的问题.这些控制参数不当,常常出现的局面是效率高但解是伪解,或者解是真解但效率低.所以,网格搜索技术被认为仍是最稳妥的技术.当然,使用网格搜索技术也需权衡效率与精度之间的问题.比如,就震源机制的求解问题而言,若搜索步长为10°,则需迭代12960次;若步长为1°,则需迭代11793600次;后者是前者的1000倍左右.大量实际 经验表明,5°的搜索步长既可以保证精度也可以保 证 效率(许忠淮等,1983; Reasenberg and Oppenheimer, 1985; Zhao and Helmberger, 1994; Zhu and Helmberger 1996; 俞春泉等,2009). 如果5°的精度仍嫌不足,那么,简单的变步长搜索即可解决这个问题,例如,GPAT便采用两步搜索,第一步步长设为5°,搜索范围为全空间;第二步步长为1°,搜索范围为第一步解为中心的10°空间.至于震源深度的步长,根据具体情况,设置为不小于1 km.即使设置为1 km,一次求解过程所需的机时通常大约为60 s(基于本文实验). 因此,变步长网格搜索技术是我们求解这个非线性问题的最终选择.

图 1 用以阐述如何为w取值的示意图 横轴为κ,纵轴为w.曲线为w随κ的变化曲线,该线在纵轴的截距为振幅与初动极性起相同作用的w值. Fig. 1 Illustration of how the w is given a value Horizontal axis is for κ while vertical axis is for w. The curve shows the variation of w with change of κ,whose intercept on vertical axis is the w value indicating that amplitudes and polarities of first motion make the same contribution.
2.3 解的不确定性

在实际工作中,我们通常面临的是非线性的近似问题.本文讨论的震源机制的求解问题就是这样的问题.由于震源位置、记录噪声、台站分布等原因,反演结果通常具有不确定性.线性问题的不确定性描述似乎存在一些被认可的方法(Flinn,1965; Anderson,1982),但是,非线性问题解的不确定性描述却非然.对于求解震源机制这个问题,已有的工 作提供了重要的启示,即用解的分布范围描述解的不 确定性(许忠淮等,1983; Reasenberg and Oppenheimer, 1985; Hardebeck and Shearer, 2002). 然而,众所周知,这种解的不确定性边界强烈依赖于目标函数的类型及其值的大小,不同的目标函数或者取值 都可能导致不同的不确定边界.对于GPAT采用的目标函数,我们尝试利用数值实验的方法确定解的不确定边界对应的目标函数边界,反过来,用目标函数的边界确定解的不确定边界.

假设在观测资料完备的情况下,解的不确定范围不超过我们感兴趣的精度1°,我们针对各种可能的震源机制类型,通过大量数值实验测量获得,当震源机制各参数改变1°时目标函数的变化δρ,那么这个δρ对应的震源参数的变化范围就是解的不确定范围.如图 2所示,假设完美的台网布局如图 3a 所示,观测资料不带有任何误差,随机产生震源机制,重复20000次实验.结果发现,δρ的均值渐近线逼近0.5,即0.5%.这意味着,解的1°不确定范围对应着目标函数值1与99.5%之间的目标函数值的范围.反过来,目标函数值1与99.5%的范围对应的震源参数的不确定性范围为1°.

图 2 随机震源机制参数改变1°引起的 目标函数值的改变与实验次数的关系 横轴为随机实验次数,纵轴为目标函数值的改变. Fig. 2 The relation between change of the objective function value caused by change of 1° on parameters of focal mechanisms Horizontal axis represents number of experiments, and vertical axis represents the change of the objective function value.

图 3 震中与台站分布.图中六角星为震中,三角形为台站 (a)台站相对于震中等间隔360°覆盖,共18个台,分两组,一组的震中距为0.6°,另一组为1°;(b)台站相对于震中等间隔180°覆盖,其他同(a);(c)台站相对于震中等间隔90°覆盖,其他同(a);(d)台站相对于震中0°覆盖,共18个台,等间距线型分布于0.6°与1°之间.图中纵横 坐标是以度为单位的长度单位. Fig. 3 Distributions of stations relatively to epicenter. David-star are epicenters and triangles are stations (a)Relatively to the epicenter,18 stations cover 360° azimuth at the same interval, and are divided into 2 groups having 0.6° and 1° epicentral distances,respectively;(b)The stations cover 180° azimuth, and others are the same as(a);(c)The stations cover 90° azimuth, and others are the same as(a);(d)The stations st and in the same line as the epicenter and between the epicentral distances 0.6° and 1° at the same interval.

在实际中,我们假设震源参数与目标函数值在最大值附近的较小范围内成线性关系,并将目标函数值范围与解的不确定性范围的对应关系加以推广,用目标函数的最大值与它的99.5%之间的值对应的解的范围描述解的不确定性,即

其中φ0、δ0、λ0和d0为目标函数最大值对应的断层走向、倾角、滑动角和震源深度,而φi、δi、λi和di为目标函数最大值与其99.5%之间的值对应的解. 2.4 标量地震矩与矩震级的确定

由于小震或微震的震源时间函数通常可近似为迪拉克函数(Ebel and Bonjer, 1990),因此,当待求地震的震源机制和震源深度确定后,待求地震的标量地震矩可以直接由观测振幅和待求地震的合成振幅确定(Godano et al., 2009),即

式中,N为采用的振幅数,an和a′ n分别表示观测及其对应的合成振幅.相应地,矩震级由下式确定(Kanamori,1977):

其中,M0是标量地震矩,其单位为N·m. 3 数值实验

为了检验GPAT的可行性与抗干扰能力,我们进行了必要的数值实验.考虑到观测台站相对于震中的方位覆盖总是有限的现实情况,设计了台站分布对反演结果的影响测试;考虑到观测点的数目时多时少的现实情况,设计了台站数目对反演结果的影响测试;考虑到观测资料总是或多或少受到噪声干扰的现实情况,设计了随机噪声对反演结果的影响测试;考虑到震源机制反演是在地震定位之后进行且定位总是存在误差的现实情况,设计了定位误差对反演结果的影响测试;考虑到计算格林函数的速度模型总是与实际模型存在差异的现实情况,设计了速度模型误差对反演结果的影响测试;最后,综合上述情况,设计了综合因素对反演结果的影响测试.

数值实验中使用的速度模型为表 1所示的水平分层模型,波形数据均利用反射-折射率方法计算(Kennett,1983; Kennett and Engdahl, 1991),波形的采样频率为10 Hz.数值实验中采用的测试地震的震源深度为5 km,标量地震矩为1015 N·m,相当于Mw3.9地震.考虑到震源机制相对于台站分布的变化可能对反演结果的影响,测试的震源机制由100组随机数确定.

表 1 数值实验使用的速度模型 Table 1 Velocity model used in numerical experiments
3.1 台站分布的影响

观测台站相对于震中的方位分布,最理想的应该是360°全方位覆盖,其次是180°和90°,最糟糕莫过于0°张角分布.为此,我们构建了如图 3所示的4种台网布局,分别呈360°、180°、90°及0°张角分布.在本实验中,为了排除台站数目过少可能带来的影响,台站的数目均设置为18.

图 4展示了测试机制(输入机制)与反演机制(反演结果)的比较.图 5展示了反演机制相对于测试机制的偏差的统计分析结果.我们分别计算走向、倾角和滑动角之间的偏差,然后挑选偏差最大者进行统计,并考虑小于1°、1°~5°、5°~10°、10°~20°、20°~40°以及大于40°的6个区间出现偏差的反演结果占总数的比例.由图可见,所有偏差均小于1°.

图 4 在台站布局的影响测试中,反演机制与测试机制的比较 红色填充区域为反演机制,而黑色迹线为测试机制.(a)、(b)、(c)和(d)依次对应于图 3中的四种情形. Fig. 4 Comparison of the inverted mechanisms with the test ones in the test for effect of the station distribution Red shading represents inverted mechanisms, and black nodal line represents test ones. (a),(b),(c) and (d)correspond to those in Fig. 1 in the same order.

图 5 在台站布局的影响测试中,反演机制参数相对于测试机制参数最大偏差统计分析 横轴为偏差,纵轴为占比.(a)、(b)、(c)和(d)分别对应图 4所示的四种情形. Fig. 5 Statistical analysis of the maximal bias of the parameters of the inverted mechanism relatively to those of the test mechanisms in the test for effect of the station distribution Horizontal axis represents the maximal bias, and vertical axis represents the ratio over the total. (a),(b),(c) and (d)correspond to the 4 cases shown in Fig. 4,respectively.

这个实验表明,只要有足够的观测点,观测点相对于震中的方位覆盖并不影响反演结果,即便台站和震中都分布在一条直线上. 3.2 台站数目的影响

就震源机制反演问题而言,为了结果的稳定与 可靠,观测点通常越多越好.但是,实际上观测点往 往是非常有限的.究竟用多少观测点才能获得稳定 而可靠的结果与使用的具体方法有关.对于本文的方法,我们通过本数值实验予以认识.

上一个数值实验表明,18个观测点,无论如何分布,足以获得稳定而可靠的反演结果.考虑到理想情况下能够确定震源机制的台站数至少为2,同时考虑到台站呈线型分布的情况过于极端,所以,在本实验中我们采用如图 6所示的张角呈90°的台站分 布,并对台站数目分别为10、6、4、2四种情况进行测试.

图 6 震中与台站分布(纵横坐标是以度为单位的长度单位) 台站相对于震中等间隔90°覆盖,分两组,一组震中距为0.6°,另一组为1°.(a)10个台;(b)6个台;(c)4个台;(d)2个台. Fig. 6 Distributions of stations relatively to epicenter The epicenter is surrounded by stations covering 90° azimuth at the same interval, and are divided into 2 groups with 0.6° and 1° epicentral distances,respectively.(a)10 stations;(b)6 stations;(c)4 stations;(d)2 stations only.

图 7展示了反演机制与测试机制的比较,图 8展示了反演机制相对于测试机制的偏差的统计结果.可以看出,所有的反演机制与测试机制之间的偏差均在1°以内.

图 7 震中与台站分布(纵横坐标是以度为单位的长度单位) (a)、(b)、(c)和(d)依次对应于图 6中的四种情形. Fig. 7 Comparison of the inverted mechanisms with the test ones in the test for effect of the station number (a),(b),(c) and (d)correspond to those in Fig. 5 in the same order.

这个实验表明,呈90°张角的2个观测点也能获得稳定可靠的反演结果. 3.3 随机噪声的影响

实际的观测资料中总是包含各种各样的干扰,随机噪声就是其中之一.通过上面两个实验,关于台站分布和台站数目对反演结果的影响有了清楚的认识,本实验旨在测试不同水平的随机噪声对反演结果的影响.

尽管以上实验表明,具有90°张角的2个观测点就能够获得可靠的结果,但是,考虑到2个台站的情形过于极端,所以,我们选择常见的6个台站、90°张角的情形进行测试,即图 6b所示的情形.

我们分别为合成记录添加最大振幅的5%、10%、20%及30%的随机噪声形成观测记录,然后利用这些观测记录进行实验.注意,我们只是添加了幅度为最大振幅5%、10%、20%及30%的噪声信号,对噪声没有其他任何限制;同时,我们注意到,当噪声水平达到30%时,记录已经达到了我们使用资料能够容忍的极限——初动较难辨认.

图 9展示了反演机制与测试机制的比较.从图中可以看出,噪声添加后反演机制已经开始偏离测试机制,并且随着噪声水平的增加,这种偏离呈增大趋势.图 10展示了反演机制相对于测试机制偏离的统计结果.可以看出,当噪声水平为5%时,偏差小 于1°的比例约为85%,在1°~5°范围内约为14%,在5°~10°范围内约为1%,大于10°的比例为零.当噪声水平为10%时,偏差小于1°的比例约为66%,在1°~5°范围内约为29%,在5°~10°范围内约为1%,在10°~20°范围内约为3%,在20°~40°范围内约为1%,在大于40°的范围内为零.当噪声水平增加至20%时,偏差小于1°的比例约为39%,在1°~5°范围内约为46%,在5°~10°范围内约为10%,在10°~20°范围内约为4%,在20°~40°范围内约为1%,在大于40°的范围内为零.当噪声水平增大至30%时,偏差小于1°的比例约为31%,在1°~5°范 围内约为44%,在5°~10°范围内约有19%,在10°~20° 范围内约为5%,在20°~40°范围内约为1%,在大于40°的范围内为零.

图 8 在台站数目的影响测试中,反演机制参数相对于测试机制参数最大偏差统计分析 (a)、(b)、(c)和(d)分别对应图 6所示的四种情形. Fig. 8 Statistical analysis of the maximal bias of the parameters of the inverted mechanism relatively to those of the test mechanisms in the test for effect of the station number (a),(b),(c) and (d)correspond to the 4 cases shown in Fig. 6,respectively.

图 9 在随机噪声的影响测试中,反演机制与测试机制的比较 (a)、(b)、(c)和(d)依次来自于噪声水平为5%、10%、20%和30%的数值实验. Fig. 9 Comparison of the inverted mechanisms with the test ones in the test for effect of the r and om noise (a),(b),(c) and (d)come from the experiments of 5%,10%,20% and 30% noise level in order.

图 10 在随机噪声误差的影响测试中,反演机制参数相对于测试机制参数最大偏差统计分析 (a)、(b)、(c)和(d)分别对应图 9中的四种情形. Fig. 10 Statistical analysis of the maximal bias of the parameters of the inverted mechanism relatively to those of the test mechanisms in the test for effect of the r and om noise (a),(b),(c) and (d)correspond to the 4 cases shown in Fig. 9,respectively.

这个实验表明,随机噪声的确对反演结果有影响,噪声水平越高,出现差异的几率越高,且出现的差异也越大.例如,当噪声水平为5%时,偏差大于1°的几率为15%,且最大偏差不超过10°; 当噪声水平达到30%时,偏差大于1°的几率达到69%,且最大偏差已达40°.因此,30%似乎是可以容忍的随机噪声上限. 3.4 定位误差的影响

当利用本文的方法反演震源机制时,我们认为震源位置是已知的.然而,事实上,震源位置总有一定的不确定性,这种不确定性必然会影响反演结果.定位误差对反演结果的影响主要是通过射线的方位角和离源角(Hardebeck and Shearer, 2002)及路径的改变导致的格林函数的极性和振幅变化引起的.这种影响的大小依赖于位置不确定性的大小.为此,我们设计了两组实验,分别考察震中位置误差与震源深度误差对反演结果的影响.第一组,震源深度不变,震中位置发生改变;第二组,震中位置不变,震源深度发生改变.

图 11展示了实验使用的观测点与地震事件的位置.在实验中,我们分别给震中0.5 km、1 km、2 km及3 km的改变.实验结果展示于图 12图 13.可以看出,当震中误差为0.5 km时,偏差小于1°的比例约为24%,在1°~5°范围内为58%,在5°~10°范围内约为17%,在10°~20°范围内约为1%,20°以上为零.当震中误差为1 km时,偏差小于1°的比例约为23%,在1°~5°范围内约为55%,在5°~10°范围内约为20%,在10°~20°范围内约为2%,20°以上为零.当震中误差为2 km时,偏差小于1°的比例约为13%,在1°~5°范围内约为58%,在 5°~10°范围内约为21%,在10°~20°范围内约为8%,20°以上为零.当震中误差为3 km时,偏差小于1°的比例约为18%,在1°~5°范围内约为75%,在5°~10°范围内约为6%,在10°~20°范围内约为1%,20°以上为零.

图 11 震中与台站分布(纵横坐标是以度为单位的长度单位) 黑色六角星为“真实”震中,白色六角星为“观测”震中.台站相对于震中等间隔90°覆盖,共6个台站,分两组, 一组震中距为0.6°,另一组为1°.在(a)、(b)、(c)和(d)中,震中分别具有0.5 km、1 km、2 km和3 km的误差. Fig. 11 Distributions of stations relatively to epicenter Solid david-stars are “real” epicenters, and empty ones are “observed” epicenters. The epicenter is surrounded by 6 stations covering 90°azimuth at the same interval, and are divided into 2 groups with 0.6° and 1° epicentral distances,respectively. In subplots(a),(b),(c) and (d),the epicenters have errors of 0.5 km,1 km,2 km and 3 km,respectively.

这个实验表明,震中误差对反演结果确有影响,误差越大,出现偏差的几率越高,但出现的偏差似乎都不超过20°.例如,当误差为0.5 km时,偏差大于1°的几率为76%,但最大偏差不大于20°; 当误差达到3 km时,偏差大于1°的几率达到82%,且最大偏差也不大于20°.由此可见,震中误差对反演结果存在影响,但小于3 km的误差导致大于20°的偏差出现的几率几乎为零.

图 12 在震中位置误差的影响测试中,反演机制与测试机制的比较 (a)、(b)、(c)和(d)依次对应于图 11中的四种情形. Fig. 12 Comparison of the inverted mechanisms with the test ones in the test for effect of the location errors (a),(b),(c) and (d)correspond to those in Figure 11 in the same order.

第二组实验仅考虑震源深度误差对反演结果的影响,所以震中位置固定不变.台站布局与第一组实验情况相同.我们将准确的深度设在5 km,深度误差从1 km到5 km以1 km间距发生变化.前面的实验表明,本实验采用的台站布局和台站数目足以确定任何一种类型的震源机制,所以,考虑到实验效率,在本实验中我们并没有对每个深度进行100次随机实验,而是具有代表性地选取走向正北的纯走滑断层、正断层、逆断层以及走向、倾角和滑动角都是45°的走滑兼逆冲机制四种类型.

图 13 在震中位置误差的影响测试中,反演机制参数相对于测试机制参数最大偏差统计分析 (a)、(b)、(c)和(d)分别对应图 12所示的四种情形. Fig. 13 Statistical analysis of the maximal bias of the parameters of the inverted mechanism relatively to those of the test mechanisms in the test for effect of the epicentral location errors (a),(b),(c) and (d)correspond to the 4 cases shown in Fig. 12,respectively.

图 14和15展示了实验结果.可以看出,当震源深度为5 km时,目标函数值均为1,观测资料与合成资料达到最佳吻合,反演机制与测试机制之间的差异为0°;当震源深度偏离5 km时,反演机制与测试机制出现差别,且随着深度差别的增大,偏差表现出趋势性增大,其中走滑机制对于震源深度最敏感,倾滑机制相对迟钝.例如,对于走滑机制,1 km的深度误差很可能引起30°的偏差,而对于倾滑机制,3 km的深度误差引起的偏差仅5°左右.另外,总体而言,震源深度向深处偏离时,这个误差对反演结果影响较小.

图 14 震源深度偏差对震源机制的影响 横轴为震源深度,纵轴为目标函数.红色震源球为真实深度对应的反演结果.(a)、(b)、(c)和(d)分别为 走向正北的走滑断层、正断层、逆断层及走向、倾角和滑动角都为45°的走滑兼逆冲断层四种情形. Fig. 14 Effects of the depth errors on focal mechanisms Horizontal axis is for focal depth and vertical axis is for objective function value. The red beach balls represent the inverted results at correct depth.(a),(b),(c) and (d)are exhibitions for cases of pure strike-slip,normal-slip,thrust-slip with northward strike and thrust-strike-slip with 45° of strike,dip and rake,respectively.

图 15 在震源深度误差的影响测试中,反演机制参数相对于测试机制参数最大偏差统计分析 横轴为震源深度,纵轴为最大偏差.(a)、(b)、(c)和(d)分别对应图 14所示的四种情形. Fig. 15 Statistical analysis of the maximal bias of the parameters of the inverted mechanism relatively to those of the test mechanisms in the test for effect of the depth errors Horizontal axis is for focal depths, and vertical axis is for the maximal biases.(a),(b),(c) and (d)correspond to the 4 cases shown in Fig. 14,respectively.
3.5 速度模型误差的影响

真实速度结构与理论模型之间总会有差异,这种差异势必影响格林函数,进而影响反演结果.为了考察速度模型误差对反演结果的影响,我们假定“真实”的速度模型为表 1给出的模型,然后对该模型中每层的P波及S波速度添加1%、3%、5%和10%的 随机扰动后形成新的模型,最后利用这些模型进行实验.实验的台网分布和台站数目与噪声试验相同.

图 16展示了反演机制与测试机制的比较.从图中可以看出,与噪声和定位误差造成的影响相比,速度模型误差对反演结果的影响更加明显,且随着误差水平的增加,反演机制与测试机制间的差别也呈增大趋势.图 17展示了不同误差水平时的统计结果.当模型误差为1%时,偏差小于1°的比例约为18%,在1°~5°之间约为68%,在5°~10°之间约为10%,在10°~20°之间约为2%,在20°~40°之间约为1%,大于40°的约为1%.当模型误差为3%时,偏差小于1°的比例为零,在1°~5°之间约为34%,在5°~10°之间约为36%,在10°~20°之间约为22%,在20°~40°之间约为6%,大于40°的约为2%.当模型误差为5%时,偏差小于1°的比例约为2%,在1°~5°之间约为32%,在5°~10°之间约为37%,在10°~20°之间约为12%,在20°~40°之间约为11%,大于40°的约为6%.当模型误差为10%时,偏差小于1°的比例为零,在1°~5°之间约为14%,在5°~10°之间约为10%,在10°~20°之间约为34%,在20°~40°之间约为21%,大于40°的约为21%.

图 16 在模型误差影响测试中,反演机制与测试机制的比较 (a)、(b)、(c)和(d)依次来自于模型误差水平为1%、3%、5%和10%的实验. Fig. 16 Comparison of the inverted mechanisms with the test ones in the test for effect of the model errors (a),(b),(c) and (d)come from the experiments with velocity-model-error level of 1%、3%、5% and 10%,respectively.

这个实验表明,模型误差对反演结果影响非常明显,误差越大,出现偏差的几率越高,且出现的偏差也越大.例如,当误差为1%时,偏差大于20°的几率为2%,且大于40°的几乎为零; 当误差达到10%时,偏差大于20°的几率达到42%,且大于40°的已超过20%.由此可见,5%似乎是模型误差上限.

以上实验测试了台站分布、台站数目、随机噪声、震中误差、震源深度误差以及模型误差对反演结果的影响.总体而言,台站分布、台站数目、随机噪声和震中误差对反演结果影响较小,而震源深度误差与速度模型误差对反演结果的影响较大.如果将偏差小于10°的结果作为一等结果,将偏差小于20°但大于10°的结果作为二等结果,将偏差小于40°但大于20°的结果作为三等结果,那么,只有模型误差或 者震源深度误差可能导致三等结果或更糟糕的结果. 3.6 综合因素的影响

真实的观测数据不可能受单一因素的干扰,因此,我们对以上因素的综合影响进行测试至关重要. 在这个实验中,我们考虑30%的随机噪声、2 km的震中误差、3 km的深度误差以及5%的速度模型误差.注意,30%为能够容忍的噪声上限,2 km的震中误差和3 km的深度误差对地方震而言也是较大的误差,而5%的速度模型误差也是可以容忍的上限. 实验用的台站分布和台站数目与噪声试验完全相同.

图 17 在模型误差影响测试中,反演机制参数相对于测试机制参数最大偏差统计分析 (a)、(b)、(c)和(d)分别对应图 16所示的四种情形. Fig. 17 Statistical analysis of the maximal bias of the parameters of the inverted mechanism relatively to those of the test mechanisms in the test for effect of the model errors (a),(b),(c) and (d)correspond to the 4 cases shown in Fig. 16,respectively.

图 18展示了反演机制与测试机制的比较以及统计分析结果.从图 18a中可以看出,多重因素对反演结果的影响非常明显.从图 18b中可以看出,偏差小于1°的比例为零,在1°~5°之间约为9%,在5°~10°之间约为27%,在10°~20°之间约为29%,在20°~40°之间约为20%,大于40°的约为15%.

图 18 在综合因素影响的测试中,反演机制与测试机制的比较(a) 以及反演机制参数相对于测试机制参数最大偏差统计分析(b) Fig. 18 Comparison of the inverted mechanisms with the test ones(a) and statistical analysis of the maximal bias of the parameters of the inverted mechanism relatively to those of the test mechanisms(b)in the test for comprehensive effect of various factors

这个实验表明,若考虑30%的随机噪声、2 km的震中误差、3 km的深度误差、5%的速度模型误差,并在台站相对于震中呈90°张角分布和台站数目为6的情况下,反演结果总是包含或多或少的偏 差.其中,65%的结果属于一类或二类结果,35%的结果属于三类或更糟糕的结果.

各种因素的组合多种多样,无法一一考虑.我们只希望通过这个实验能够对各种因素对反演结果的综合影响有所认识,或者说,仅希望这个测试结果能 够对从实际资料中获得的结果的合理解释提供参考. 4 讨论与结论

GPAT是对已有初动-振幅类反演方法的一般化,也是对包括远震震源机制反演在内的所有地震的震源机制反演方法的一般化.之所以说是对已有初动-振幅类方法的一般化,是因为GPAT所使用的极性并不局限于初动极性,它包括任何一个或多个振幅的极性,振幅也不局限于P波、S波和面波的振幅,它可以是任何一个或多个震相的振幅.之所以说是对震源机制反演方法的一般化,是因为GPAT不但适合于地方地震和区域地震,也适合于远场地震.也就是说,根据GPAT的技术原理,三分向记录任何极性信息和振幅信息都可以作为资料,而且,地震不分大小,都可以用GPAT求解震源机制.

GPAT方法以合成地震波场与观测地震波场的相似性为准则,以观测矢量和合成矢量的相关系数为目标函数,而目标函数是地震震源机制参数和震源深度的非线性函数.我们采用变步长网格搜索技术求解这个非线性方程,既可以保证解的精度,也可以保证解的稳定.在目前的计算机技术条件下,通常完成一个求解过程的时间也不过60s左右.需要说明的是,为了提高求解效率,实际上我们也尝试过遗传算法.但是,遗传算法需要许多控制参数,不合理的参数设置可能导致伪解,也可能导致效率过低,这些控制参数的调整无形中增加了工作量,减低了 工作效率,所以,最终放弃了遗传算法.而变步长网格搜索技术因其稳定可靠的优点成为我们的最终选择.

关于描述非线性方程解的不确定性问题,至今还没有一种被广泛认可的方法,GPAT方法借鉴了已有的用解的分布空间描述其不确定性的方法(许忠淮等,1983; Reasenberg and Oppenheimer, 1985; Hardebeck and Shearer, 2002).但不同的是,基于GPAT方法的特点,我们通过大量数值实验定量给出了确定解的不确定性范围的技术途径.我们相信这种方法不是最好的方法,但却是合理的方法.根据这个方法,当相关系数达到最大值ρmax 时,必然存在另一个值0.995 ρmax .在这二者之间δρ的范围内,存在一个对应的解空间R,这个解空间的边界就是解的不确定范围.可以想见,ρmax δρ与观测资料的质量有关,所以,观测资料的质量直接影响解空间R的大小.资料质量越好,R越小,即不确定性越小;反之,R越大,即不确定性越大.

GPAT方法把确定标量地震矩放在确定震源机制和震源深度之后进行,主要原因是它采用观测波矢量和合成波矢量的相似性作为目标函数.当然,这样做也减少了解空间的自由度.需要说明的是,这样计算标量地震矩的前提是事件的震源时间函数为迪拉克-δ函数.对于小震来说这个前提是成立的,但是,对于大震便不再恰当.因此,随着地震的增大,标量地震矩可能被高估.不过,如果根据具体情况为格林函数添加合适的震源时间函数,GPAT同样可以给出合适的标量地震矩.

为了检验GPAT的可行性以及抗干扰能力,我们考虑了常见的干扰因素,设计了不同实验进行检验.几乎所有数值实验表明,GPAT是完全可行的.关于GPAT的抗干扰能力,可以归纳如下:(1)台站相对于震中的方位覆盖对反演结果几乎无影响;(2)相对于震中张角呈90°的、大于或等于2个以上的台站数目对反演结果几乎无影响;(3)相对于震中张角呈90°的6个台站情况,30%的随机噪声可能会导致约1%~2%的三类结果;(4)相对于震中张角呈90°的6个台站情况,3 km的震中偏差不会导致三类结果;但是,反演结果对震源深度比较敏感,尤其是走滑机制对震源深度最为敏感,1 km的偏差可能会导致三类或更糟糕的结果;(5)相对于震中张角呈90°的6个台站情况,反演结果对模型误差极为敏感,例如,当模型误差达到10%时,已经有42%的结果为三类或更糟糕的情况.考虑到 震源深度误差直接与速度模型误差有关,所以,归根到底,影响反演结果的主要因素是速度模型误差.因此,除速度模型误差外,GPAT具有良好的抗干扰能力.

为了考察GPAT的实用性,我们只设计了一个综合因素对反演结果影响的测试.这个实验结果只从某个侧面或视角反映了GPAT的实用性.但毕竟只考虑了综合因素组合的一种情况.更多的实用性的实证应该来自对实际资料的应用.

根据本文的数值实验,我们认为GPAT是可行的,具有较强的抗干扰能力.

致谢 非常感谢两位评审专家的悉心评点.
参考文献
[1] Aki K, Richards P G. 1980. Quantitative Seismology Theory and Methods, vol. I, W. H. San Francisco, California: Freeman and Company.
[2] Anderson K R. 1982. Robust earthquake location using M-estimates. Phys. Earth Planet. Interiors., 30(2-3): 119-130.
[3] Cesca S, Bufornand E, Dahm T. 2006. Amplitude spectra moment tensor inversion of shallow earthquakes in Spain. Geophys. J. Int., 166(2): 839-854.
[4] de Natale G, Ferraro A, Virieux J. 1991. A probability method for local earthquake focal mechanisms. Geophys. Res. Lett., 18(4): 613-616.
[5] Dreger D, Helmberger D. 1991. Source parameters of the sierra-madre earthquake from regional and local body waves. Geophys. Res. Lett., 18(11): 2015-2018.
[6] Dreger D S, Helmberger D V. 1993. Determination of source parameters at regional distances with three-component sparse network data. J. Geophys. Res., 98(5): 8107-8125.
[7] Dziewonski A M, Chou T A, Woodhouse J H. 1981. Determination of earthquake source parameters from waveform data for studies of global and regional seismicity. J. Geophys. Res., 86(B4): 2825-2852.
[8] Ebel J E, Bonjer K P. 1990. Moment tensor inversion of small earthquakes in southwestern Germany for the fault plane solution. .Geophys. J. Int, 101(1): 133-146.
[9] Flinn E A. 1965. Confidence regions and error determinations for seismic event location. Rev. Geophys., 3(1): 157-185.
[10] Francis T W, Wang P D, Chen Y T. 1989. Determination of focal mechanism using SH to P amplitude ratio. Acta Seismologica Sinica (in Chinese), 11(3): 275-281.
[11] Gilbert F. 1973. Derivation of source parameters from low-frequency spectra. Philos. Trans. R. Soc. Lond., 274(1239): 369-371.
[12] Godano M, Regnier M, Deschamps A, et al. 2009. Focal mechanisms from sparse observations by nonlinear inversion of amplitudes: method and tests on synthetic and real data. Bull. Seism. Soc. Am., 99(4): 2243-2264.
[13] Godano M, Bardainne T, Regnier M, et al. 2011. Moment-tensor determination by nonlinear inversion of amplitudes. Bull. Seism. Soc. Am., 101(1): 366-378.
[14] Hardebeck J L, Shearer P M. 2002. A new method for determining first-motion focal mechanisms. Bull. Seism. Soc. Am., 92(6): 2264-2276.
[15] Hardebeck J L, Shearer P M. 2003. Using S/P amplitude ratios to constrain the focal mechanisms of small earthquakes. Bull. Seism. Soc. Am., 93(6): 2434-2444.
[16] Helmberger D V, Engen G R. 1980. Modeling the long-period body waves from shallow earthquakes at regional distances. Bull. Seism. Soc. Am., 70(5): 1699-1714.
[17] Julian B R, Foulger G R. 1996. Earthquake mechanisms from linear-programming inversion of seismic-wave amplitude ratios. Bull. Seism. Soc. Am., 86(4): 972-980.
[18] Kanamori H. 1977. The energy release in great earthquakes. J. Geophys. Res., 82(20): 2981-2987.
[19] Kasahara K. 1963. Computer program for a fault-plane solution. Bull. Seism. Soc. Am., 53(1): 1-13.
[20] Kennett B L N. 1983. Seismic Wave Propagation in Stratified Media. Cambridge: Cambridge University Press.
[21] Kennett B L N, Engdahl E R. 1991. Traveltimes for global earthquake location and phase identification. Geophys. J. Int., 105(2): 429-465.
[22] Kisslinger C. 1980. Evaluation of s to p amplitude ratios for determining focal mechanisms from regional network observations. Bull. Seism. Soc. Am., 70(4): 999-1014.
[23] Kisslinger C, Bowman J R, Koch K. 1982. Determination of focal mechanism from SV/P amplitude ratios at small distances. Phys. Earth Planet. Inter., 30(2-3): 172-176.
[24] Langston C A. 1982. Single-station fault plane solutions. Bull. Seism. Soc. Am., 72(3): 729-744.
[25] Lay T, Wallace T C. 1995. Modern Global Seismology. New York: Academic Press.
[26] Liang S H, Li Y M, Shu P Y, et al. 1984. On the determining of source parameters of small earthquakes by using amplituderatios of P and S from regional network observations. Chinese J. Geophys. (in Chinese), 27(3): 247-257.
[27] Lindh A, Fuis G, Mantis C. 1978. Seismic amplitude measurements suggest foreshocks have different focal mechanisms than aftershocks. Science, 201(4350): 56-59.
[28] Liu H L, Helmberger D V. 1985. The 23: 19 aftershock of the 15 October 1979 Imperial Valley earthquake: more evidence for asperity. Bull. Seism. Soc. Am., 75(3): 689-708.
[29] Liu J, Zheng S H, Kang Y, et al. 2004. The focal mechanism determinations of moderate 2 small earthquakes using the first motion and amplitude ratio of P and S wave. Earthquake (in Chinese), 24(1): 19-26.
[30] Nakamura A, Horiuchi S, Hasegawa A. 1999. Joint focal mechanism determination with source-region station corrections using short-period body-wave amplitude data. Bull. Seism. Soc. Am., 89(2): 373-383.
[31] Reasenberg P, Oppenheimer D H. 1985. FPFIT, FPPLOT, and FPPAGE, FORTRAN computer programs for calculating and displaying earthquake fault-plane solutions. U. S. Geological Survey Open-File Report 85-739, 109p.
[32] Saikia C K, Herrmann R B. 1985. Application of waveform modeling to determine focal mechanisms of four 1982 Miramichi aftershocks. Bull. Seism. Soc. Am., 75(4): 1021-1040.
[33] Schwartz S Y. 1995. Source parameters of aftershocks of the 1991 Costa Rica and 1992 Cape Mendocino, California, earthquakes from inversion of local amplitude ratios and broadband waveforms. Bull. Seism. Soc. Am., 85(6): 1560-1575.
[34] Snoke J A, Munsey J W, Teague A C, et al. 1984. A program for focal mechanism determination by combined use of polarity and SV-P amplitude ratio data. Earthquake Notes, 55(3): 15.
[35] Snoke J A. 1989. Earthquake Mechanisms.//James D E ed. Encyclopedia of Geophysics. New York: Van Nostrand Reinhold Company, 239-245.
[36] Wallace T C, Helmberger D V, Mellman G R. 1981. A technique for the inversion of regional data in source parameter studies. J. Geophys. Res., 86(B3): 1679-1685.
[37] Warren D H. 1979. Fault-plane solutions for microearthquakes preceding the Thanksgiving Day, 1974, Earthquake at Hollister, California. Geophys. Res. Lett., 6(8): 633-636.
[38] Xu Z H, Yan M, Zhao Z H. 1983. Evaluation of the direction of tectonic stress in north China from recorded data of a large number of small earthquake. Acta Seismologica Sinica (in Chinese), 5(3): 268-279.
[39] Yu C Q, Tao K, Cui X F, et al. 2009. P- wave first-motion focal mechanism solutions and their quality evaluation. Chinese J. Geophys. (in Chinese), 52(5): 1402-1411.
[40] Zhao L S, Helmberger D V. 1994. Source estimation from broadband regional seismograms. Bull. Seism. Soc. Am., 84(1): 91-104.
[41] Zhu L P, Helmberger D V. 1996. Advancement in source estimation techniques using broadband regional seismograms. Bull. Seism. Soc. Am., 86(5): 1634-1641.
[42] 梁尚鸿, 李幼铭, 束沛镒等. 1984. 利用区域地震台网P、S 振幅比资料测定小震震源参数. 地球物理学报, 27(3): 249-256.
[43] 刘杰, 郑斯华, 康英等. 2004. 利用P波和S波的初动和振幅比计算中小地震的震源机制解. 地震, 24(1): 19-26.
[44] 吴大铭, 王培德, 陈运泰. 1989. 用SH 波和P 波振幅比确定震源机制解. 地震学报, 11(3): 275-281.
[45] 许忠淮, 阎明, 赵仲和. 1983. 由多个小地震推断的华北地区构造应力场的方向. 地震学报, 5(3): 268-279.
[46] 俞春泉, 陶开, 崔效锋等. 2009. 用格点尝试法求解P波初动震源机制解及解的质量评价. 地球物理学报, 52(5): 1402-1411.