地球物理学报  2015, Vol. 58 Issue (2): 453-462   PDF    
CAP方法反演震源机制的误差分析:以胶东半岛两次显著中等地震为例
郑建常, 林眉, 王鹏, 徐长朋    
山东省地震局, 济南 250014
摘要:利用区域波形数据使用CAP方法反演中强地震的震源机制正逐渐得到广泛应用.本文以胶东半岛近期发生的两次显著中等地震为例,讨论了使用CAP方法反演震源机制时的误差估计,展示了反演结果的不确定性分析过程.2013年11月23日和2014年1月7日在山东莱州和乳山分别发生了M4.6和M4.3级中等地震,两次事件均造成了较大影响.我们基于CAP方法,使用自助抽样(bootstrap)技术多次重复反演过程,得到大样本量的震源机制解数据;基于这些数据,使用粒子群算法和聚类分析技术给出了优化解,估计了震源机制解的误差范围,并利用震源机制解的P、T轴给出了震源球上的概率密度分布.
关键词莱州地震     乳山地震     波形反演     聚类分析     不确定性    
Error analysis for focal mechanisms from CAP method inversion: An example of 2 moderate earthquakes in Jiaodong Peninsula
ZHENG Jian-Chang, LIN Mei, WANG Peng, XU Chang-Peng    
Earthquake Administration of Shandong Province, Ji'nan 250014, China
Abstract: As an effective focal mechanism inversion method for regional earthquakes, CAP (Cut and Paste) is widely used in China in recent years. Its quality and error level need to be evaluated for such solutions are increasingly retrieved. On the other hand, when using the CAP method, those phases which fitted well are usually chosen for inversion, and the other phases which are thought 'bad’ or fitted not so good are ignored. It has been found that different station combinations will lead to varied results with unneglectable discrepancies. Objectively speaking, in a scientific perspective, this artificially selected process will increase uncertainties in final inversion results, especially under the present instrument status in China. Furthermore, because the grid search scheme used in the CAP method which is not evenly distributed on focal sphere, we can not give a convincing proof to illustrate that whether the cause of badly-fitted waveforms comes from data error or from un-sufficient searched solutions.
Two earthquakes of M4.6 and M4.3 occurred in Shandong Peninsula on 23 November 2013 and 7 January 2014, respectively. The former is the largest event in the Shandong area since 1995. For convenience, we label the former event as the Laizhou earthquake, and the latter as the Rushan earthquake, according to their epicenters. Taking these two events for example, this paper discusses error estimation for focal mechanism inversion using the CAP method. The paper also presents an uncertainty analysis process for the inversion results.
Briefly, a bootstrap technique is adopted, waveforms are randomly sampled with equal probability from origin dataset, and then used as data for a repeating inversion procedure. After a large number of inversions, e.g., 1000 times, we finally get bootstrap results consisting of 1000 focal mechanisms.
Based on these focal mechanisms, we conduct the following work. (1) We employ a PSO (particle swarm optimization) algorithm to search a solution, of which the Kagan angle is minimal to all the double couple models, and use its standard deviation as the inversion results uncertainty range. The uncertainty of solutions for the Laizhou event is ±23.7°; for Rushan event, is ±6.4°. These two solutions can be evaluated as ‘A’ level according to the Hardebeck's indicator of mechanisms quality.
(2) A clustering analysis is used for bootstrap results. For the Rushan event, the clustering center is coincident to the PSO optimized solution; for the Laizhou event, several centers are found in clustering, despite isolate solutions. There are two clustering centers, of which the corresponding data proportion is about 98.7 percent.
(3) Projecting P and T axes of bootstrap results onto a focal sphere, calculating its probability density, we get the probability density distribution of focal mechanisms on the focal sphere. Then we can give the confidence interval on different levels for mechanism solutions.
The method demonstrated in this paper is not confined to achieving more accurate focal mechanism and obtaining rational inversion error, while it can also be used to exclude isolate and incorrect solutions effectively, and avoid the effect of data from stations with larger disturbances. Therefore, this method can be used to invert focal mechanisms automatically immediately after moderate earthquakes occur. Sustaining by powerful computational capabilities, we can get more accurate and reliable focal mechanism results without manual work.
Key words: Laizhou earthquake     Rushan earthquake     Waveform inversion     Clustering analysis     Uncertainty    
1 引言

据山东台网测定,2013年11月23日13时44分在山东省莱州市(37.10°N,120.02°E)发生M4.6级地震,这次地震是山东陆地地区自1995年苍山5.2级地震后发生的最大地震,影响范围广,山东东部市地普遍有感;2014年1月7日22时24分在山东乳山(36.80°N,121.70°E)发生M4.3级地震,这次地震也造成胶东地区大面积有感.这两次事件是1970年以来胶东半岛陆地及近海地区发生的最强烈的地震活动,其中莱州地震震中区在1970年以来的小震目录上属于典型的少震、弱震区,活动水平不高,很少有ML≥3.0地震发生,仅在1991年2月以及2012年7月分别发生最大ML3.8级小震序列各一次;乳山地震震中区历史上曾发生公元1046年岠嵎山51 / 2级和1939年乳山下初51 / 2级地震.虽然这两次地震的震中区历史上没有强烈地震活动,但胶东半岛北部近海曾发生多次6、7级强震,如1548年渤海海峡7.0、1948年威海近海6.0以及1969年渤海7.4级等.因此确定这两次显著中等地震的震源机制对于研究区域地质构造的活动特征,以及研判该地区的地震危险性等具有重要的科学价值.

我们使用近年来在国内得到广泛使用的CAP(Cut and Paste)方法反演这两次地震的震源机制.为了得到更准确的解,并且合理地估计反演结果的不确定性,我们使用自助抽样法(bootstrap)对反演过程随机重复,在大样本量的反演结果基础上,使用粒子群算法搜索优化解,利用动态聚类技术对结果进行聚类分析,从而得到了更加稳定可靠的断层面解,给出了可能的误差范围,并进一步给出了震源球上P、T轴的概率密度分布. 2 理论与方法 2.1 CAP方法反演震源机制

震源机制和传播效应决定了观测波形的变化.如果地壳模型已知,可以准确地计算波形传播过程中的效应,因此我们可以通过理论波形s(t)和观测波形u(t)的拟合来估计震源的断层面参数.

双力偶震源产生的理论位移s(t)可以表示为(Zhu and Helmberger,1996):

其中,i=1,2,3 对应三种基本断层响应,即:垂直走滑、垂直倾滑以及倾角为45°的倾滑;Gi为格林函数,Ai是辐射系数,φ是台站方位角,M0为标量地震矩.θ,δ,λ分别为断层的走向、倾角、滑动角.系数Ai由6个矩张量分量和台站方位角表示如下(Jost and Herrmann,1989):

走向θ、倾角δ、滑动角λ,以及标量地震矩M0等可以通过求解以下方程进行估计:

波形反演可以使用全波形数据,也可以单独使用体波或面波震相进行拟合.CAP方法是一种联合使用体波和面波进行反演的方法,近年来在国内得到了广泛的应用(吕坚等,2008; 黄建平等,2009; 郑勇等,2009; 龙锋等,2010; 韩立波等,2012),由 于该方法分别截取波形的Pnl部分和面波部分分别 拟合(Zhao and Helmberger,1994; Zhu and Helmberger,1996),并在反演的过程中允许它们在适当的时间变化范围内相对移动,在一定程度上避免了因为地壳模型不准确而引起的震相到时的误差因素,对速度模型和地壳横向变化的依赖性较小,因此在实际的区域地震震源机制求解中有明显的优势.CAP方法使用频率F-波数K法(Zhu and Rivera,2002)计算格林函数,使用网格搜索方法搜寻最优震源机制参数和震源深度.考虑到波形随震中距的衰减,方法定义误差函数如下:

式中,r为台站震中距,r0为选定的参考震中距,p为指数因子.参考有关研究,对体波p=1,面波p=0.5(韩立波和蒋长胜,2012). 2.2 CAP方法的优化解及其不确定性估计

地球内部的任意震源可以表示为6个独立分量的矩张量,由于CAP方法限制震源为双力偶模型,并且无需发震时刻的对齐,因此只需对震源模型的三个角度,即走向θ、倾角δ、滑动角λ,以及标量地震矩M0进行搜索,理论上而言,仅需要2个台站的三分向波形就可以求解;虽然研究显示,对于大多数3个三分向台的组合,使用波形反演就可以得到相对准确的震源机制,但实际情况也显示,不同的台站组合波形反演得到的解之间仍然存在一定的差异(Godano et al.,2009; 郑建常和陈运泰,2012).目前国内台网密度已经达到相当水平,在东部地区,一个中等地震通常有数十甚至上百个宽频带台能够记录到清晰的波形,以此次莱州地震为例,通过对原始波形进行去均值、去趋势、积分等简单变换后,根据直观的观察,震中距300 km范围内,采样率100 Hz的宽频带三分向波形有近40个台站的资料可用.在使用CAP方法求解震源机制时,一般的做法是选择部分波形拟合较好的台进行反演,删去拟合不好的台或者震相;有些情况下,甚至仅使用面波部分而删除体波震相,需知面波尤其是径向和切向分量,很容易受到台站下方浅层地壳结构的影响.由于CAP方法是采用网格搜索的方法,因而这种人为的选择,必然会为反演结果增加主观的不确定性因素.我们无法令人信服地说明,拟合不好的波形究竟是数据本身确实存在干扰,还是说搜索到的解无法满足该条数据.另外,CAP方法虽然可以在最后的输出结果中给出断层面参数的不确定性,但该估计值只是面向所使用的台站数据的结果,在上述的人为选择下,该不确定性估计能够在多大程度上客观地反映最终解的整体不确定性,是无法说明的.

为了求得更加稳定可靠的解并且合理客观地给出解的误差估计,我们在相对丰富的观测数据基础上,采用自助抽样统计方法进行分析.具体方法是在可用的观测台站中可重复地随机抽取一定数量的台站组成新的台站组合,使用该台站组合的观测数据重复反演过程.在大量的重复计算后(例如,超过1000次),可以有效地排除观测质量不高或存在较大干扰误差的数据的影响,从而得到更加接近真实解的结果,并且可以有效地给出解的不确定性.

另外,由于CAP方法在搜索断层面解时采用的是网格搜索的方法,然后通过插值计算误差函数e的最小值,并且由于固定步长的走向、倾角、滑动角的尝试位置在震源球上的分布是不均匀的(许向彤等,1995),因此在最终解中可能会有空缺(gap)的存在.为了求解优化解,我们进一步使用Kagan(1991)定义的双力偶模型最小空间旋转角,对上面自助抽样得到的大量满足条件的震源机制结果进行分析,定义与所有解的空间偏转角度和为目标函数,使用粒子群非线性优化方法搜索该目标函数最小的结果,视为最优解. 2.3 聚类分析

在震源机制求解中,常见的情况是在震源球上存在几簇相对集中分布的解,对这些可能的解直接取数学平均是不甚合理的,并且在数据存在较大误差或干扰的情况下,满足条件的可能解的分布范围也许会相当大.因此针对这一现象,刁桂苓等(1992)俞春泉等(2009)分别使用系统聚类和动态聚类技术,对所有的可能解进行聚类分析,求取聚类中心作为反演的优化解,数值试验和实际应用都有很好的效果.聚类分析可以很好地排除孤立解和错误解,从而在大量的数据中获取更加接近真实解的结果.本文在使用不同台站组合重复进行波形反演后,同样得到了大量的震源机制解数据,受台站布局和数据误差的影响,这些解或多或少存在差别,因此对这些结果进行聚类分析是很有必要的. 3 数据与资料

本文使用了山东台网提供的波形资料,其中还包括了邻省如辽宁、河北、江苏等省交换资料的部分台站.图 1给出了本项研究使用的台站分布,其中个 别台如JIM、ZSL、HUD等为短周期台,在波形反演 中没有使用.本文研究中,首先由观测记录直接读取初动符号,用于约束波形反演;然后将观测数据扣除仪器响应,经过去均值、零漂等预处理后积分至位移 记录,旋转到Z-R-T坐标系,对观测波形和理论波形同样进行带通滤波,然后用于反演.

图 1 本文研究的两次地震震中及山东台网台站分布图Fig. 1 Map of stations in Sh and ong Network and two earthquakes studied in this paper. Red circles denote epicenter,triangles are stations,and solid black lines are faults.

使用Chang等(2006)给出的朝鲜半岛南部至黄海地区的中上地壳速度结构模型用于本文的震源机制反演.相关地质资料显示,胶东半岛、南黄海以及朝鲜半岛南部在大地构造分区上都属于下扬子地块,地质构造属性相对较为一致(Ree et al.,1996). 4 结果与分析 4.1 乳山M4.3震源机制

选择震中距在250 km以内的15个台站的宽频带波形记录进行反演,Pnl和面波的反演波段分别选择0.05~0.15 Hz和0.033~0.067 Hz频段.图 2给出了不同深度的最佳双力偶解,及拟合误差随不同深度变化的关系,由图可见,震源深度在4 km时观测波形和理论波形的错配值最小,说明事件的震源深度较浅.

图 2 2014年1月7日乳山M4.3级地震不同震源深度的波形拟合误差及最佳震源机制解Fig. 2 Waveform fit errors and best focal mechanisms as function of depth for Jan.7,2014 Rushan M4.3 event

由CAP方法反演得到的最佳震源机制:节面A的参数为:走向202°、倾角75°、滑动角153°;节面B的参数为:走向299.5°、倾角64°、滑动角16.7°;参考乳山序列的双差定位结果(李冬梅和郑建常,2014)分析认为,节面B可能是乳山地震的发震断层;震源机制显示为左旋走滑型,反演得到此次地震的矩震级MW=4.2.

图 3给出了对应最优解的理论波形和观测波形的拟合情况.15个台一共75个震相,其中理论波形与观测波形相关系数大于0.9的有39个,超过50%;相关系数大于0.6(相关性较好)的有67个,约占89.3%;最佳解的方差减少(variance reduction)为70.3%,说明理论波形很好地拟合了观测波形,反演结果是可靠的.个别台(如WEH)平均相关系数较差,可能与台站位于震源机制解的节面线附近,振幅相对较小所致;另外如CHD台的拟合程度不好,可能与该台处于海域、噪声干扰较大有关.

图 3 2014年1月7日乳山M4.3地震最优解的理论波形(红)与观测波形(黑)
波形图下方第一行数字为各段理论地震波形相对实际观测波形的移动时间,正值表示理论波形相对观测波形超前.第二行数字为理论波形与观测波形的相关系数(百分比). 波形图左侧字母为台站,其下数字分别为台站震中距(km)和方位角(°). 图左侧的震源球上红色 区域代表压缩区,白色代表拉张区,震源球采用下半球投影.震源球上标注的“+”和“-”表示反演使用台站的P波初动.
Fig. 3 Comparison between synthetics(red) and observed(black)seismograms of Jan.7,2014 Rushan M4.3 event
The numbers on the lower left side of each seismogram are the time shifts(upper) and cross-correlation coefficient in percent(lower). Positive time shifts mean that the observed data have been delayed. The letters on the left side are stations,the numbers below it are epicentral distance(in km) and azimuth(degree). The red color in beach-ball denotes compression area,while white is extension. The ‘+’ and ‘-’ signs on beach-ball indicate polarities on inversion used stations. Lower hemisphere projection is used.

使用自助抽样的统计方法,对乳山地震震源机制解的不确定性进行估计.选用震中距在300 km以内的采样率为100 Hz的22个三分向宽频带台的观测波形组成原始数据集,为了保证用于反演的数据的样本量,设用于反演的台站数为20个,对原始数据集进行每个台站等概率、可重复地随机抽取,抽取出的台站波形组成新的数据集,然后用于CAP方法的波形反演.对上述的抽取台站反演过程重复1000次,将反演得到的震源机制的断层节面解和P、T轴绘制在一个震源球上,见图 4.可以看出,反演中除去个别反演过程的断层面解出现一定程度的偏离外,其余结果集中分布,均显示为近走滑的机制;图 4中的P、T轴位置和断层节面线集中成丛,大致显示出断层面解的误差范围.

图 4 自助抽样得到的1000次乳山地震震源机制解及粒子群最优解(下半球极射投影)震源球上黑色细线条表示自助抽样结果的断层节面线,红色线条表示粒子群优化解的节面线.Fig. 4 1000 focal mechanisms of Rushan M4.3 event retrieved by a bootstrapping process,all nodal lines(black) and P,T axes(blue and red points,respectively)are plotted on one beach-ball. The red lines on the beach-ball show the optimized solution given by a Particle Swarm Optimization method. Its corresponding P,T axes are also displayed on the beach-ball(yellow and green point,respectively). Lower hemisphere projection is used.

使用粒子群非线性优化方法,以与自助抽样给出的1000个机制解(图 4)的Kagan角之和为目标函数,搜索最优解.结果显示最优解为,节面A:走向 208.4°,倾角89.7°,滑动角154.3°;节面B:走向 298.5°,倾角64.3°,滑动角0.3°;最优解与图 4所示1000个解的平均夹角4.37°,以其与所有解Kagan角的2倍标准差为震源机制解的误差范围,结果显示不确定性为6.44°(图 5).

图 5 乳山地震自助抽样结果与粒子群优化解的Kagan角分布Fig. 5 Kagan angles of bootstrap results to the PSO solution for Rushan event

对自助抽样结果进行动态聚类分析,结果显示最优解为:节面A:走向208.0°,倾角89.3°,滑动角 154.1°,节面B:走向298.4°,倾角64.1°,滑动角1.8°,与粒子群优化解非常一致(见表 1).

表 1 波形反演乳山M4.3地震震源机制解结果 Table 1 Parameters of focal mechanism results from waveform inversion for Jan.7,2014 Rushan M4.3 event

将自助抽样结果中的P、T轴投影到震源球上(图 4),对其进行概率密度统计分析,结果见图 6.

图 6 乳山地震自助抽样结果的震源球概率密度分布俯视图(未进行极射投影)
色标中正值表示T轴的概率密度分布,负值表示P轴的概率密度分布.
Fig. 6 Probability density distribution of solutions on beach-ball(top view of lower hemisphere,without projection)
Positive values on the color scale(corresponding to the red area on beach-ball)indicate probability of T axis,while negative values(corresponding to blue area)mean probability of P axis.
4.2 莱州M4.6震源机制

使用CAP方法对2013年11月23日莱州M4.6 地震进行反演(郑建常等,2015),同样进行CAP反演情况的自助抽样统计分析.选用震中距在270 km以内的采样率为100 Hz的22个三分向宽频带台的观测波形组成原始数据集,采用全样本随机抽取方法,自助抽样反演1000次,图 7给出了反演得到的震源机制的断层节面解和P、T轴在震源球上的分布情况.结果显示,莱州地震的自助抽样结果同样很好地显示出了反演得到震源机制解的误差范围,相对于乳山地震,出现了极个别反演过程的结果偏离较大的情况.

图 7 自助抽样得到的莱州地震震源机制解及粒子群最优解(下半球极射投影)
黑色节面线为自助抽样得到的震源机制解;红色节面线为粒子群最优解.
Fig. 7 1000 Focal mechanisms of Laizhou M4.6 event retrieved by a bootstrapping process,all nodal lines(black) and P,T axes(blue and red points,respectively)are plotted on one beach-ball. The red lines on the beach-ball show the optimized solution given by a PSO method. Its corresponding P,T axes are also displayed on the beach-ball(cyan and yellow point,respectively). Lower hemisphere projection is used.

使用粒子群非线性优化方法,搜索与自助抽样结果的旋转角最小的解.结果显示最优解为,节面A:走向236.9°,倾角76.2°,滑动角-169.3°;节面B:走向144.3°,倾角79.6°,滑动角-14.0°;与所有自助抽样解的平均偏转角17.4°(图 8),以其与所有解Kagan角的2倍标准差(图 8红色虚线所示)为误差范围,结果显示震源机制解的不确定性为23.7°.

图 8 莱州地震自助抽样结果与粒子群优化解的Kagan角分布Fig. 8 Kagan angles of bootstrap results to the PSO solution for Laizhou event

从自助抽样得到的所有机制解在震源球上的分布情况(图 7)可以直观地看出,断层节面线尤其是北西向节面呈现出两组集中.由于我们定义的粒子群优化的目标函数是搜索与所有自助抽样解的空间旋转最小,因此从图 7可以看出,最优解的节面位置处于其中一组的边缘位置,在此情况下,对自助抽样结果进行聚类分析是有意义的.

以两个震源机制解之间的最小空间旋转角(即Kagan角)为距离的定义,对1000次自助抽样结果进行聚类分析,图 9给出了聚类谱系图(由于完整的聚类树过于密集和庞大,因此我们只显示了Kagan>7°的部分),以50°为阈值,可以将结果分为5类.图 10给出了聚类分析的结果,属于Ⅰ类的数据占 32.1%,Ⅱ类66.6%,其余三类数据合计仅有1.3%. 由图 10可以看出,其余三类的断层节面线和P、T轴位置明显偏离集中区且机制解类型与绝大部分结果(走滑型)不一致,是典型的孤立解.孤立解(或错误解)的出现,可能说明我们使用的数据中个别台站(或分向)存在较大干扰.使用俞春泉等(2009)的方 法求取了四类解的聚类中心,其中I类解的聚类中心的断层面参数(设北东向节面为发震断层面)为:走向231.6°,倾角88.7°,滑动角-168.2°;Ⅱ类解的 聚类中心为:走向238.5°,倾角74.1°,滑动角-164.8°.

图 9 莱州地震自助抽样结果的聚类谱系图Fig. 9 Dendrogram plot of the hierarchical binary cluster tree for Laizhou event

图 10 莱州地震自助抽样结果的聚类分析
(a)断层节面线的分类显示;(b)机制解P(+)、T(⊙)轴位置的分类显示.断层节面线和P、T轴颜色表示分类,与图 9分类颜色一致.
Fig. 10 Clustering results of focal mechanisms from a bootstrap process for Laizhou event
The different colors of nodal lines and P(“+” sign in right panel),T(“⊙” in right panel)axes denote different classes,which are corresponding to colors shown in Fig. 9.

将自助抽样结果中的P、T轴投影到震源球上,对其进行概率密度统计分析,结果见图 11.可见,同聚类分析的结果一致,P轴位置的概率密度在震源球上出现了两个极值区,分别对应Ⅰ类和Ⅱ类两个聚类中心.

图 11 莱州地震自助抽样结果的震源球概率密度分布
(a)下半球俯视图(未进行极射投影);(b)东南方向,45°三维侧视图.色标说明同图 6.
Fig. 11 Probability density distribution of solutions on beach-ball
(a)Top view of lower hemisphere,without any projection;(b)Side view of whole beach-ball,from south-east 45° position. Positive values on the color scale(corresponding to the red area on beach-ball)indicate probability of T axis,while negative values(corresponding to blue area)mean probability of P axis.
5 讨论与结论

基于山东省宽频带数字地震波形资料,本文首先使用CAP方法反演了近期胶东半岛地区发生的两次显著中等地震活动的震源机制,讨论了如何合理地估计CAP方法反演震源机制的误差范围以及如何确定优化解的问题.我们首先使用自助抽样方法,对原始数据进行等概率随机抽样,多次重复波形反演过程,排除了人为选择数据的干扰,得到大样本量的震源机制解数据;在此基础上,我们

(1)使用了粒子群优化算法从中搜索与这些机制解空间偏转角最小的解当作优化解,以Kagan角的二倍标准差作为反演结果的不确定性范围,结果显示:乳山地震的粒子群优化解为:走向298.5°,倾角64.3°,滑动角0.3°,不确定性为±6.4°;莱州地震 的优化解为:走向236.9°、倾角76.2°、滑动角-169.3°,不确定性为±23.7°.

(2)对自助抽样结果进行聚类分析,其中:乳山地震结果的聚类中心与粒子群优化解基本一致;莱州地震结果存在多个聚类,排除孤立解后,有两个聚 类中心,其对应两类数据合计占结果的98.7%,说明此次地震的真实解在这两类数据范围内.

(3)将自助抽样结果中的P、T轴投影到震源球,对其进行概率密度统计,给出了机制解在震源球上的概率密度分布图.

本文方法不单可以得到更准确的震源机制优化解、给出科学合理的误差估计,而且可以有效地排除孤立解和错误解,克服存在较大干扰台站的数据的影响,因此在震后应急的震源机制自动化求解中可以发挥作用.在强大计算能力的支持下,无须人工干预即可得到准确可靠的震源机制结果,从而为震害评估、趋势分析等提供重要的科学依据.

自助抽样结果显示乳山地震震源机制解的不确定性要小于莱州地震,笔者推测可能有莱州地震使用的台站中个别台的干扰较大的原因,另外也无法排除莱州地震的震源破裂过程可能更加复杂的可能.CAP方法中用于计算格林函数的F-K方法使用狄拉克-Delta函数作为震源时间函数(Zhu and Rivera,2002),对于小震级的事件该简化方案更为适用,莱州地震(M4.6)与乳山地震(M4.3)震级相差不大,但莱州地震的自助抽样结果出现了两个概率较高的聚类中心,在使用大部分相同台站的情况下,这可能意味着莱州地震的震源破裂随时间的变化可能与狄拉克-Delta函数存在一定的偏离.Rodríguez-Lozoya等(2008)的研究显示,区域中等地震也可能有复杂的震源破裂过程,在该问题上的深入研究需要更进一步的工作.

致谢 感谢两位匿名审稿专家提出的宝贵意见.聚类分析中使用了俞春泉等(2009)提供的部分开放代码,粒子群搜索使用了S Chen给出的Matlab软件包,在此一并表示感谢!
参考文献
[1] Chang S J, Baag C E. 2006. Crustal structure in Southern Korea from joint analysis of regional broadband waveforms and travel times. Bull. Seism. Soc. Amer., 96(3): 856-870.
[2] Diao G L, Yu L M, Li Q Z. 1992. Hierarchical clustering analysis of the focal mechanism solution-Taking the Haicheng Earthquake Sequences for example. Earthquake Research in China (in Chinese), 8(3): 86-92.
[3] 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. Amer., 99(4): 2243-2264.
[4] Han L B, Jiang C S, Bao F. 2010. Source parameter determination of 2010 Taikang MS4.6 earthquake sequences. Chinese J. Geophys. (in Chinese), 55(9): 2973-2981, doi: 10.6038/j.issn.0001-5733.2012.09.016.
[5] Han L B, Jiang C S. 2012. Focal mechanism inversion of 8 Jun 2011 Toksun MS5.3 earthquake. Acta Seismologica Sinica (in Chinese), 34(3): 415-422.
[6] Huang J P, Ni S D, Fu R S, et al. 2009. Source mechanism of the 2006 MW5.1 Wen'an earthquake determined from a joint inversion of local and teleseismic broadband waveform data. Chinese J. Geophys. (in Chinese), 52(1): 120-130.
[7] Jost M L, Hermann R B. 1989. A student's guide to and review of moment tensors. Seism. Res. Lett., 60(2): 37-57.
[8] Kagan Y Y. 1991. 3-D rotation of double-couple earthquake sources. Geophys. J. Int., 106(3): 709-716.
[9] Li D M, Zheng J C. 2014. Relocation analysis of Oct. 1, 2013 Rushan swarms. Qilu Earthquake Sciences (in Chinese), 10(1): 26-28.
[10] Long F, Zhang Y J, Wen X Z, et al. 2010. Focal mechanism solutions of ML≥4.0 events in the MS6.1 Panzhihua-Huili earthquake sequence of Aug 30, 2008. Chinese J. Geophys. (in Chinese), 53(12): 2852-2860, doi: 10.3969/j.issn.0001-5733.2010.12.008.
[11] Lü J, Zheng Y, Ni S D, et al. 2008. Focal mechanisms and seismogenic structures of the MS5.7 and MS4.8 Jiujiang-Ruichang earthquakes of Nov. 26, 2005. Chinese J. Geophys. (in Chinese), 51(1): 158-164, doi: 10.3321/j.issn:0001-5733.2008.01.020.
[12] Ree J-H, Cho M, Kwon S-T, et al. 1996. Possible eastward extension of Chinese collision belt in South Korea: The Imjingang belt. Geology, 24(12): 1071-1074.
[13] Rodríguez-Lozoya H E, Quintanar L, Ortega R, et al. 2008. Rupture process of four medium-sized earthquakes that occurred in the Gulf of California. J. Geophys. Res., 113(B10301), doi: 10.1029/2007JB005323.
[14] Xu X T, Xu Z H, Zhang D N. 1995. A probabilistic grid test method of determining earthquake focal mechanisms using P-wave onset polarity data. Seismological and Geomagnetic Observation and Research (in Chinese), 16(4): 34-42.
[15] 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, doi: 10.3969/j.issn.0001-5733.2009.05.030.
[16] Zhao L S, Helmberger D V. 1994. Source estimation from broadband regional seismograms. Bull. Seism. Soc. Amer., 84(1): 91-104.
[17] Zheng J C, Li D M, Wang P, et al. 2015. Focal mechanisms and seismic tectonic features of the 2013 Laizhou M4.6 earthquake sequence. Seismology and Geology (in Chinese), in press.
[18] Zheng J C, Chen Y T. 2012. Stability of sparse station data inversion for deviatoric moment tensor solution of regional earthquakes. Acta Seismologica Sinica (in Chinese), 34(1): 31-43.
[19] Zheng Y, Ma H S, Lü J, et al. 2009. Source mechanism of strong aftershocks (Ms≥5.6) of the 2008/05/12 Wenchuan earthquake and the implication for seismotectonics. Science in China Series D: Earth Sciences, 52(6): 739-753.
[20] Zhu L P, Helmberger D V. 1996. Advancement in source estimation techniques using broadband regional seismograms. Bull. Seism. Soc. Amer., 86(5): 1634-1641.
[21] Zhu L P, Rivera L A. 2002. A note on the dynamic and static displacements from a point source in multilayered media. Geophys. J. Int., 148(3): 619-627.
[22] 刁桂苓, 于利民, 李钦祖. 1992. 震源机制解的系统聚类分析——以海城地震序列为例. 中国地震, 8(3): 86-92.
[23] 韩立波, 蒋长胜, 包丰. 2012. 2010年河南太康Ms4.6地震序列震源参数的精确确定. 地球物理学报, 55(9): 2973-2981, doi: 10.6038/j.issn.0001-5733.2012.09.016.
[24] 韩立波, 蒋长胜. 2012. 2011年6月8日新疆托克逊Ms5.3地震震源机制解反演. 地震学报, 34(3): 415-422.
[25] 黄建平, 倪四道, 傅容珊等. 2009. 综合近震及远震波形反演2006文安地震(Mw5.1)的震源机制解. 地球物理学报, 52(1): 120-130.
[26] 李冬梅, 郑建常. 2014. 2013年10月1日乳山震群重新定位结果分析. 齐鲁地震科学,10(1):26-28.
[27] 龙锋, 张永久, 闻学泽等. 2010. 2008年8月30日攀枝花—会理6.1级地震序列ML≥4.0事件的震源机制解. 地球物理学报, 53(12): 2852-2860, doi: 10.3969/j.issn.0001-5733.2010.12.008.
[28] 吕坚, 郑勇, 倪四道等. 2008. 2005年11月26日九江-瑞昌Ms5.7, Ms4.8地震的震源机制解与发震构造研究. 地球物理学报, 51(1): 158-164, doi: 10.3321/j.issn:0001-5733.2008.01.020.
[29] 许向彤, 许忠淮, 张东宁. 1995. 求震源机制P波初动解的格点尝试概率法. 地震地磁观测与研究, 16(4): 34-42.
[30] 俞春泉, 陶开, 崔效锋等. 2009. 用格点尝试法求解P波初动震源机制解及解的质量评价. 地球物理学报, 52(5): 1402-1411, doi: 10.3969/j.issn.0001-5733.2009.05.030.
[31] 郑建常, 王鹏, 李冬梅等. 2015. 2013年莱州M4.6地震序列震源机制与发震构造初探. 地震地质, 待刊.
[32] 郑建常,陈运泰. 2012. 稀疏台网反演区域地震偏量矩张量解的稳定性. 地震学报, 34(1): 31-43.
[33] 郑勇, 马宏生, 吕坚等. 2009. 汶川地震强余震(Ms≥5.6)的震源机制解及其与发震构造的关系. 中国科学D辑: 地球科学, 39(4): 413-426.