地球物理学进展  2017, Vol. 32 Issue (5): 1930-1935   PDF    
引用本文
王丽霞, 张辉, 李赫. 2017. 基于粒子群算法对岷县两次中等地震震源机制解的误差分析. 地球物理学进展, 2017, 32(5): 1930-1935, doi: 10.6038/pg20170509  
WANG Li-xia, ZHANG Hui, LI He. 2017. Error analysis based on PSO algorithm for focal mechanisms of the two Minxian moderate earthquakes. Progress in Geophysics (in Chinese), 32(5): 1930-1935, doi: 10.6038/pg20170509   
基于粒子群算法对岷县两次中等地震震源机制解的误差分析
王丽霞1, 张辉1, 李赫2     
1. 中国地震局兰州地震研究所, 兰州 730000
2. 天津市地震局, 天津 300201
收稿日期：2017-02-13; 修回日期：2017-08-09.
基金项目：中国地震局地震科技星火计划攻关项目（XH15042，XH14049）和国家自然科学基金项目（41304048）共同资助.
通讯作者：张辉, 副研究员, 主要从事地震学及地震综合预测研究.(E-mail:zhh_gx@163.com)
摘要：本文采用CAP方法反演了2011年两次岷县中等地震的震源机制解，首先讨论了方位角对反演结果的影响，通过对不同方位角组合的比较，得到的结果比较稳定，进一步得出两次地震的震源机制解.然后，采用自助抽样（bootstrap）方法，基于CAP方法进行多次重复反演，得到岷县两次中等地震事件的大样本量震源机制解数据.基于这些大样本数据，以Kagan角为目标函数，利用粒子群非线性优化方法搜索震源机制解的优化解，并估计了震源机制解的误差范围，MS 4.0、MS 4.5的误差分别为6.0°和18.3°.结果分析表明，CAP结果与得到的优化解基本接近，其误差在可接受的误差范围内.
关键词：CAP方法    震源机制解    自助抽样    粒子群算法    
Error analysis based on PSO algorithm for focal mechanisms of the two Minxian moderate earthquakes
WANG Li-xia1 , ZHANG Hui1 , LI He2     
1. Lanzhou Institute of Seismology, CEA, Lanzhou 730000, China
2. Earthquake Administration of Tianjin, Tianjin 300201, China
Abstract: We accept the Cut and Past(CAP) method to inverse the focal mechanisms of the two Minxian moderate earthquakes and discuss error estimation for the inversion in this paper. First, we inverse the focal mechanism with combinations of different azimuth stations used the CAP method and the results show that the azimuth has little effect on the inversion. And we select a dozen stations to calculate the focal mechanisms of the two Minxian moderate earthquakes. Then, We sample by equal probability randomly the waveforms of the two earthquakes recorded by the seismic station and employ these data repeating inversion procedure. We obtain 1000 focal mechanisms by 1000 times inversions. Based on these solutions, we employ the PSO (Particle Swarm Optimization) algorithm to confirm the optimal solution when the Kagan angle is minimal to all the double couple models and evaluate the uncertainty range of the inversion result using its standard deviation. The error of solutions are separately ±6.0° and ±18.3° for the Minxian earthquakes with MS 4.0 and MS 4.5.Meanwhile, We project P and T axes of bootstrap results onto a focal sphere and achieve the probability density distribution of focal mechanisms. We conclude that the method demonstrated in this paper is not confined to achieve more accurate focal mechanisms and obtain inversion error, while it can be used to exclude isolate and incorrect solutions effectively and avoid the effect of data from station azimuth. Therefore, the CAP method can be applied to get reliable solution automatically immediately after a moderate earthquake occur.
Key words: CAP method     focal mechanisms     bootstrap technique     PSO algorithm    
0 引言

甘肃岷县地区2000年以来先后发生了2003年11月13日5.2级、2004年9月7日5.0级、2013年7月22日6.6级地震，中强地震活动频繁，地震集中在临潭—宕昌断裂、礼县—罗家堡断裂和西秦岭北缘断裂的交汇三角区域活动.2011年2月和11月在该区域发生了4.0、4.5级两次中等地震，之后在岷县漳县6.6级地震前中等地震出现了长达21个月的平静，因此确定这两次中等地震的震源机制对于研究该区域地质构造的活动特征，分析区域的地震危险性等具有重要的科学意义.

本文使用近年来国内广泛应用的CAP(Cut and Past)方法分别反演这两次地震的震源机制解.为了得到更为准确的结果，我们使用自助抽样方法(bootstrap)随机可重复地抽取反演所参与的台站样本，得到大量反演结果(抽样反演1000次).基于大样本量的反演结果，使用粒子群算法搜索优化解，得到更加稳定可靠的断层面解，给出可能的误差范围(郑建常等，2015)，并进一步给出了震源球上P、T轴的概率密度分布.

1 理论方法 1.1 CAP方法反演震源机制

反演震源机制解的CAP方法(Zhao and Helmberger, 1994；Zhu and Helmberger, 1996)，其主要思想是将区域近震宽频带波形资料分割为Pnl波和面波部分，分别赋予不同权重，计算理论地震图和实际观测波形之间的目标误差函数，在给定的参数空间中进行网格搜索，得到二者全局误差最小的震源机制解和震源矩心深度.对于绝大多数的构造地震，CAP方法能反演得到最佳的双力偶节面解，但在一些例如火山和地热活动较频繁的地区，矿洞坍塌或核爆等，仍可发现一些非双力偶源的地震，Zhu和Ben-Zion(2013)在CAP方法基础上添加了非双力偶成分进行反演，发展为全矩张量反演方法gCAP(a generalized CAP method).本文研究的岷县两次中等地震属于天然地震事件，因此采用CAP方法便能得到可靠的震源机制双力偶节面解.

1.2 震源机制解误差分析及优化解

国内大量的研究(吕坚等，2008；郑勇等，2009；韩立波等，2012；张辉等，2014；杨萍等，2017)表明CAP方法在震源机制解与地震矩心深度研究方面的有效性与可靠性，且反演结果具有对速度模型和地壳横向变化依赖性相对较小的优点.随着数字化测震台站的建设，对一次地震事件可选用的台站越来越多，因此不同的台站组合波形反演得到的解之间存在一定的差异(Godano et al., 2009；郑建常和陈运泰，2012).在使用CAP方法求解震源机制时，通常是挑选部分波形拟合较好的台，删去拟合不好的台或者震相，有些情况下，可能仅保留个别台站的面波部分进行反演，我们知道面波尤其是径向和切向分量，很容易受到台站下方浅层地壳结构的影响.由于CAP方法采用网格搜索的方法，这种人为的选择也会为反演结果增加主观的不确定性因素；另外，CAP方法在最后的输出结果中给出断层面参数的不确定性，但该估计值只是面向所使用的台站数据的结果，由于加入人为的因素，也增加了解的整体不确定性.

为了得到更加稳定可靠的震源机制解，采取自助抽样方法，重复随机抽取一定数量的台站，用不同的台站组合基于CAP方法进行大量的反演过程，得到大样本量的震源机制解.对自助抽样的大量震源机制结果进行分析，定义两个震源机制解的P、T、B轴坐标架间的最小空间旋转角为Kagan角(1991)，把所有解与其平均值之间的Kagan角的和作为目标函数，利用粒子群非线性优化方法(Partical Swarm Optimization-PSO)搜索该目标函数的最小值为优化解，同时用每个可选解和所有解的平均值之间Kagan角的标准差来表述自助抽样反演结果的离散程度，以2倍Kagan角的标准差来描述震源机制解的误差范围.

2 数据处理与结果分析 2.1 数据处理和速度模型

本文使用了甘肃地震台网及周边邻省(青海、宁夏、四川)部分台站的宽频带波形资料，按照台站方位角覆盖及信噪比要求进行波形数据处理，筛选出的台站分布如图 1中所示.反演时对选取台站的观测数据去除仪器响应，并旋转至大圆路径得到径向、切向和垂向记录，进而分解为Pnl波和面波两部分，赋予不同的权重，通过0.05~0.2 Hz和0.02~0.1 Hz的4阶Butterworth带通滤波器来压制噪声.采用频率-波数法(F-K)(Zhu and Rivera, 2002)计算格林函数，得到理论地震图，并采用与观测波形相同的分解、滤波过程.

图 1 本文研究的两次地震震中及参与反演的台站分布图五角星代表震中，三角形为参与了Ms4.0地震反演的台站，圆点为参与Ms4.5地震反演的台站，黑色线条代表断层. Figure 1 Map of station involving the inversion and two earthquakes studied in this paper Stars denote epicenter; Triangles are stations involved in the inversion of MS 4.0, circles are stations involved in the inversion of MS 4.5; and solid black lines are faults.

在crust 2.0速度和密度分层结构模型的基础上，结合李少华等(2012)人的研究结果，得到甘东南区域一维速度结构模型，速度模型分为5层(详见表 1).

表 1 甘东南区域一维速度结构模型 Table 1 1-D velocity model in southeast area of Gansu province
2.2 CAP方法反演结果

图 1所示两个地震事件周围可以选用的地震台站较多，选用表 1的速度模型，按所在方位角0°~90°、0°~180°、0°~270°和0°~360°组成不同的台站组，采用CAP方法反演得到岷县两次地震事件的震源机制解(表 2).由表 2可以看出，两次地震事件不同台站组合得到节面解参数存在一定的差异，但结果总体上比较稳定.

表 2 不同方位角台站组合反演结果 Table 2 The result of combinations from different azimuth

我们挑选方位角覆盖均匀的12个台站，在深度从1~15 km，每隔1 km的震源深度上进行波形拟合，并得到各深度上拟合误差.图 2为两次地震震源机制解随深度的变化图，图 3为所对应最佳震源深度时的波形拟合情况.

图 2 震源机制解反演误差随深度变化图 (a)2月23日MS 4.0地震；(b)11月2日MS 4.5地震. Figure 2 Error plots as a function of depth in focal mechanism solution (a)MS 4.0;(b)MS 4.5.

对于岷县4.0级地震的震源机制随深度变化不大(图 2a)，深度为7 km时反演方差最小.图 3a为深度7 km时的波形拟合图，除去拟合不好未参与反演的震相外，共有55个震相，相关系数达到0.85以上占81.82%，说明理论波形很好地拟合了观测波形，反演结果是可靠的.最佳双力偶解：节面Ⅰ走向131°、倾角41°、滑动角57°，节面Ⅱ走向352°、倾角57°、滑动角115°，矩震级MW=4.36，震源机制类型为走滑兼逆冲型.

图 3 岷县两次地震最佳解的观测波形(黑色)与理论波形(红色)拟合图 红线表示理论地震图，黑线表示观测地震图，绿线表示不参与反演震相；波形下方的两行数据分别表示理论地震图相对观测地震图的移动时间及二者的相关系数(用百分比表示). Figure 3 Comparison between synthetics(red) and observed (black) seismograms of Minxian The red waveforms denote synthetics seismograms, the black waveforms are observed seismograms; the green waveforms are not involved in inversion. 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.

如图 2b显示岷县4.5级地震CAP反演结果随深度变化比较稳定，震源矩心深度为6 km，其对应的最佳双力偶解：节面Ⅰ走向317°、倾角53°、滑动角52°，节面Ⅱ走向189°、倾角51°、滑动角129°，矩震级为MW=4.59，其观测波形和理论波形的拟合情况如图 3b所示，除去拟合效果不好的震相外，相关系数高于0.80的占78.95%，拟合效果较好.

3 震源机制解最优解 3.1 2月23日MS 4.0震源机制解

基于自助抽样技术，利用粒子群优化方法从大量机制解中得到优化解，并估计震源机制解的不确定性.选择震中距在400 km以内22个台站资料，为保证用于反演台站的样本量，设用于反演的台站数为12个，对原始数据进行每个台站等概率、可重复地随机抽取1000次组成新的数据集，基于CAP方法反演得到1000个震源机制解.图 4a给出了全部震源机制断层节面解及P、T轴在一个震源机制球上分布情况，可以看出，除个别偏离的断层节面解外，断层节面解结果整体呈集中分布，P、T轴的分布集中于成簇，其分布情况大致可以反映出震源机制结果的误差分布情况.

图 4 自助抽样得到的1000次震源机制解及粒子群优化解 (a)2月23日MS 4.0地震；(b)11月2日MS 4.5地震.黑色为节面为1000次的抽样结果，绿色节面为粒子群非线性搜索得到的优化解；蓝色为T轴位置，红色为P轴位置，白色圆圈为粒子群方法得到的P、T轴优化解. Figure 4 1000 focal mechanisms by bootstrapping process earthquake (a)MS 4.0;(b)MS 4.5. All nodal lines(black) and P, T axes(blue and red points, respectively) are plotted on one beach-ball. The green lines on the beach-ball show the optimized solution given by a PSO method. Its corresponding P, Taxes are also displayed on the beach-ball (the hollow points).

对自助抽样后反演得到的震源机制结果，利用粒子群非线性优化方法搜索优化解，其结果为：节面Ⅰ走向129.1°、倾角41.9°、滑动角54.1°，节面Ⅱ走向353.3°、倾角57.2°、滑动角117.8°.其最优解与图 4a 1000次自助抽样解的平均夹角为3.8°，以其与所有解Kagan角的2倍标准差为震源机制解的误差范围，结果显示解的不确定性值为6.05°(图 5a虚线).将自助抽样计算得到的P、T轴结果进行概率密度统计分析(图 6a)，得到了机制解在震源球上的密度分布特征.

图 5 地震自助抽样结果与粒子群优化解的Kagan角分布 (a)2月23日MS 4.0地震；(b)11月2日MS 4.5地震.虚线为Kagan角的2倍标准差：不确定性值. Figure 5 Kagan angles of bootstrap results tothe PSO solution for two events (a)MS 4.0;(b)MS 4.5. The dotted lines denote the uncertainty results.

图 6 P、T轴概率密度分布图 (a)MS 4.0地震；(b)MS 4.5地震.负值代表T轴，正值代表P轴. Figure 6 Probability density distribution of solutions on beach-ball (a)MS 4.0;(b)MS 4.5.The positive values and negative values respectively denote P, T axes.

表 3 反演得到的震源机制解结果 Table 3 Parameters of focal mechanism results fron waveform inversion for two earthquakes
3.2 11月2日MS 4.5震源机制解

对岷县4.5级地震事件进行自助抽样统计并重复反演过程，从反演中的25个台站中随机可重复抽取13个台站组成1000个反演样本组合，得到1000次震源机制解，结果分布见图 4b，节面及T轴分布相对分散，P轴分布相对集中，形成两个集中点.使用粒子群非线性优化方法，搜索与自助抽样结果的旋转角最小的解，结果显示最优解为：节面Ⅰ走向320.9°、倾角48.2°、滑动角59.3°，节面Ⅱ走向182.6°、倾角50.2°、滑动角119.7°，与图 4b 1000次自助抽样解的平均偏转角为7.4°.将Kagan角的2倍标准差作为震源机制解的误差范围，图 5b中虚线所示其不确定性值为18.30°；对P、T轴进行概率密度统计分析得到图 6b结果，正负两个极值对应于P、T轴位置，与得到反演结果P、T轴分布一致.

4 结论与讨论 4.1

本文基于区域数字地震台网资料，首先采用CAP方法反演了2011年岷县两次中等地震的震源机制解，使用自助抽样方法，对原始数据进行等概率随机抽样，利用粒子群非线性优化方法得到了两次地震的优化解，结果表明：

(1) 采用CAP方法反演地震的震源机制解时，台站方位角单一对结果具有一定的影响，但是在分析得到的误差范围之内，还是比较稳定的.因此，在测震台站足够多的情况下应充分考虑台站与地震分布情况，如果台站分布不是很理想，利用该方法也可以得到比较稳定的震源机制解.

(2) 选用方位角覆盖均匀的12个台站，反演得到的两次岷县地震事件的震源参数：MS 4.0地震震源深度7 km，节面Ⅰ走向131°、倾角41°、滑动角57°；MS 4.5地震震源深度6 km，节面Ⅰ走向317°、倾角53°、滑动角52°.

(3) 通过自助抽样技术，基于CAP方法反演了两次中等地震的1000次震源机制解，利用粒子群非线性算法搜索得到最优解，同时由Kagan角的二倍标准差给出反演结果的不确定性范围，得到：MS 4.0地震震源机制优化解节面Ⅰ走向129.1°、倾角41.9°、滑动角54.1°，误差范围为6.0°；MS 4.5地震震源机制优化解节面Ⅰ走向320.9°、倾角48.2°、滑动角59.3°，误差范围为18.3°.

4.2

粒子群非线性优化方法不仅可以得到更准确的震源机制优化解、给出科学合理的误差估计，而且可以有效地排除孤立解和错误解.因此在区域地震震源机制的反演中，采用粒子群非线性优化方法，基于CAP方法反演震源机制解具有很好的应用前景.

致谢 CAP反演程序由美国圣路易斯大学的Zhu L P教授提供，文中部分图件使用GMT绘图，在此一并致谢！
参考文献
[] 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[J]. Bull. Seism. Soc. Amer., 99(4): 2243–2264. DOI:10.1785/0120080210
[] Han L B, Jiang C B, Bao F. 2012. Source parameter determination of 2010 Taikang MS 4..6 earthquake sequences[J]. Chinese J. Geophys. (in Chinese), 55(9): 2973–2981. DOI:10.6038/j.issn.0001-5733.2012.09.016
[] Kagan Y Y. 1991. 3-D rotation of double-couple earthquake sources[J]. Geophys. J. Int., 106(3): 709–716. DOI:10.1111/gji.1991.106.issue-3
[] Li S H, Wang Y B, Liang Z B, et al. 2012. Crustal structure in southeastern Gansu from regional seismic waveform inversion[J]. Chinese J. Geophys. (in Chinese), 55(4): 1186–1197. DOI:10.6038/j.issn.0001-5733.2012.04.015
[] Lü J, Zheng Y, Ni S D, et al. 2008. Focal mechanisms and seismogenic structures of the MS 5..7 and MS 4.8 Jiujiang-Ruichang Earthquakes of Nov. 26, 2005[J]. Chinese J. Geophys. (in Chinese), 51(1): 158–164. DOI:10.3321/j.issn:0001-5733.2008.01.020
[] Yang P, Zhang H, Feng J G. 2017. Preliminary study on focal mechanism solution and seismogenic structure of Qilian.. Qinghai MS5.2 earthquake on Nov. 23, 2015[J]. China Earthq. Eng. J. (in Chinese), 39(1): 150–153.
[] Zhang H, Zhang L P, Feng J G. 2014. Focal mechanism solutions and characteristics of the July 22, 2013 Minxian-Zhangxian MS 6.6 Earthquake sequence[J]. Earthquake (in Chinese), 34(4): 110–117.
[] Zhao L S, Helmberger D V. 1994. Source estimation from broadband regional seismograms[J]. Bull. Seismol. Soc. Amer., 84(1): 91–104.
[] Zheng J C, Chen Y T. 2012. Stability of Sparse Station Data Inversion for Deviatoric moment tensor solution of regional earthquakes[J]. Acta Seismol. Sin (in Chinese), 34(1): 31–43.
[] Zheng J C, Lin M, Wang P, et al. 2015. Error analysis for focal mechanisms from CAP method inversion:An example of 2 moderate earthquakes in Jiaodong Peninsula[J]. Chinese J. Geophys. (in Chinese), 58(2): 453–462. DOI:10.6038/cjg20150209
[] 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[J]. Sci. China Ser. D Earth Sci., 52(6): 739–753. DOI:10.1007/s11430-009-0074-3
[] Zhu L P, Helmberger D V. 1996. Advancement in source estimation techniques using broadband regional seismograms[J]. Bull. Seismol. Soc. Amer., 86(5): 1634–1641.
[] Zhu L P, Ben-Zion Y. 2013. Parametrization of general seismic potency and moment tensors for source inversion of seismic waveform data[J]. Geophys. J. Int., 194(2): 839–843. DOI:10.1093/gji/ggt137
[] Zhu L P, Rivera L A. 2002. A note on the dynamic and static displacements from a point source in multilayered media[J]. Geophys. J. Int., 148(3): 619–627. DOI:10.1046/j.1365-246X.2002.01610.x
[] 韩立波, 蒋长胜, 包丰. 2012. 2010年河南太康MS 4.6地震序列震源参数的精确确定[J]. 地球物理学报, 55(9): 2973–2981. DOI:10.6038/j.issn.0001-5733.2012.09.016
[] 李少华, 王彦宾, 梁子斌, 等. 2012. 甘肃东南部地壳速度结构的区域地震波形反演[J]. 地球物理学报, 55(4): 1186–1197. DOI:10.6038/j.issn.0001-5733.2012.04.015
[] 吕坚, 郑勇, 倪四道, 等. 2008. 2005年11月26日九江-瑞昌MS 5.7、4.8地震的震源机制解与发震构造研究[J]. 地球物理学报, 51(1): 158–164. DOI:10.3321/j.issn:0001-5733.2008.01.020
[] 杨萍, 张辉, 冯建刚. 2017. 2015年11月23日青海祁连MS 5.2地震震源机制解及发震构造初步探讨[J]. 地震工程学报, 39(1): 150–153.
[] 张辉, 张浪平, 冯建刚. 2014. 2013年7月22日岷县漳县6.6级地震序列震源机制解及其特征分析[J]. 地震, 34(4): 110–117.
[] 郑建常, 陈运泰. 2012. 稀疏台网反演区域地震偏量矩张量解的稳定性[J]. 地震学报, 34(1): 31–43.
[] 郑建常, 林眉, 王鹏, 等. 2015. CAP方法反演震源机制的误差分析:以胶东半岛两次显著中等地震为例[J]. 地球物理学报, 58(2): 453–462. DOI:10.6038/cjg20150209
[] 郑勇, 马宏生, 吕坚, 等. 2009. 汶川地震强余震(MS ≥ 5.6) 的震源机制解及其与发震构造的关系[J]. 中国科学D辑:地球科学, 39(4): 413–426.