一、模型研制的目的、必要性和需求分析
根据中国地震科学实验场2016年制定的研究规划，系统开展川滇地区公共速度模型构建的基础性研究工作，包括结构模型构建、震源参数反演及结构模型验证三个方面。中国地震科学实验场覆盖川甘交界到云南南部区域，位于欧亚板块与印度板块互相碰撞挤压强烈变形地区，地质构造复杂多样，涵盖了川滇菱形块体、滇南、滇西、巴颜喀拉块体东段，包括龙门山、鲜水河、安宁河、则木河、小江、红河、小金河等重要断裂带。实验场区域地震活动频繁、地震类型丰富，有长期积累的观测资料，是开展地震科学实验的理想场所。
川滇地区由于其独特的构造背景，中小地震频发，有着丰富的中小地震的观测资料。而地震的震源参数是研究地震和认识地震的重要组成部分，通过研究地震的震源参数，可以深入的了解地震性质、发震机理、构造应力场等，能够为区域构造、岩石流变特性、应力状态等地震学研究提供关键信息。此外，研究中小地震的震源参数可以认识区域应力场小尺度的时空变化，可以在更广的范围更加全面的认识区域应力状态以及介质的流变性质（罗艳，2010a）；强震余震的中小地震震源参数可以为主震发震机理研究提供重要的信息（Wu and Takeo,2004;郑勇等，2009）。
因此，我们对川滇实验场区域的4.0级以上的地震震源机制解进行研究，给出研究区域的地震震源性质和深度，以及起断层特征，可以为建立统一的川滇实验场模型提供关键性的支撑。
二、研究内容
收集2009－2017年间川滇实验场区的地震的波形数据，对级的地震，采用波形反演方法及序列发展方法计算其震源机制解、震源深度和矩震级；对于速度模型比较好的情况，采用矩张量反演，确定双偶源和非双偶源的成分。对一些显著地震，结合精细的速度结构，利用CAP、gCAP等方法和深度震相波形对比方法对震源深度进行精确确定，将这些地震作为参考震例，进一步提高地震定位的精度。
三、模型预期指标和技术路线
模型预期指标：以ASCII文件的形式和统一的格式，给出可查看和方便读取的川滇地区2009年以来4.0级以上地震震源机制解和矩震级。
技术路线：①收集实验场区2009年以来的地震的波形数据，进行整理和除错。②对原始地震波形数据进行去仪器响应、去均值、去倾斜等预处理，转化为速度纪录。③对不同的震级采用不同的滤波频段，对ML≥4.0级的地震，采用CAP和gCAP波形反演方法及序列发展方法计算其震源机制解、震源深度和矩震级；对ML≥5.5的较强地震，在远震记录上也有较高信噪比的记录，拟采用CAP的发展方法CAPjoint方法对地震的震源参数进行反演，得到更可靠的震源参数。④对反演结进行筛选，确保可靠性。
为了保证结果的可靠性，本项目组将分成三个独立的研究组同时开展研究，采用同样的数据，分别采用不同的方法或者技术，对震源机制、矩心深度等信息进行研究，并确定其不确定度。对比分析结果的可靠性。
四、模型数据来源和质量分析
1. 分析数据和台站分布
地震目录来自中国地震台网中心，收集了川滇地区2009-2017年地震事件共2468个，其中可用地震事件共634个，地震事件共1663个。地震分布如图1所示，其中红色点表示地震事件。台站采用国家固定台网，分布如下图2，实际计算每个事件时，选择地震事件350km以内的台站。
为了更好的保证数据质量，我们除了采用地震目录之外，利用AIC-autopick方法自动拾取初至，然后利用人工检测初至拾取的好坏。最终选取信噪比较高，初至清晰的地震事件的波形数据进行反演，共计选出了400多个地震事件进行了震源机制反演。
图1 地震事件分布图
图2 台站分布图
2. 速度结构模型
由于川滇地区地下结构复杂，因此采用了三种速度模型计算格林函数：云南-贵州地区速度模型来自云南地震局；四川盆地速度模型以及青藏高原东南缘速度模型来自郑勇等（2013）的噪声成像结果，具体如下图3所示。
图3反演采用的三种速度模型
3. 数据和结果的质量分析：
⑴对原始数据：①检视地震波形，去除无地震和低质量数据；②查看震源深度和发震时刻，去除错误信息。
⑵对震级机制解反演结果：①检视理论波形和观测波形的拟合程度，确保两者方向一致，振幅接近，相关系数足够大；②参与反演的台站分布尽可能的均匀，数量至少为6个。
按以上原则对原始数据和反演结果进行筛选，确保震级机制解的可靠性。
五、模型研制关键技术方法
对于ML≥4级的中强地震，分别采用gCAP、CAP同步源全矩张量反演方法对其震源机制解和矩震级、震源深度进行反演。
1. CAP方法。
将近震宽频带数字波形记录分为体波（Pnl）和面波部分，在双力偶源的震源假设下，分别计算它们的理论地震图和实际观测波形之差的目标函数，在给定参数空间范围内进行网格搜索，同时反演震源机制解、矩震级和震源深度。由于CAP方法充分利用了全波形的信息，反演时一定程度上允许理论和观测的Pnl和面波到时分别相对平移，降低了对速度模型和地壳横向变化的依赖性。许多学者的研究结果都充分证明了CAP方法在震源机制解研究方面的有效性与可靠性。对ML≥5.5的较强地震，在远震记录上也有较高信噪比的记录，采用CAP的发展方法CAPjoint方法对地震的震源参数进行反演，得到更可靠的震源参数。CAPjoint方法基于CAP的思想，采用近震资料和远震资料，联合反演，在空间上更大范围覆盖，以更准确地了解地震的震源性质。
2. gCAP方法。
由于川滇区域的地壳速度结构比较复杂，震源机制解可能不仅限于双力偶源，CLVD分量也可能存在。因此，我们也采用gCAP方法对震源机制进行反演。gCAP（Zhu & Helmberger 1996;Zhu &Ben-Zion 2013）方法基于波形拟合，该方法将地震数据分为体波部分（Pnl）和面波部分，分别计算合成波形与真实记录波形的误差函数，网格搜索得到最佳的震源参数、震源深度、地震矩等相关参数，通过分解地震矩，能够得到各向同性、双力偶源以及CLVD成分，是很普遍的用来求解震源机制解的方法。
3. 同步源全矩张量反演算法。
CAP和gCAP方法属于比较常规的反演方法，对于中小地震比较合适。对于地震震级偏大的中强地震，在一些情况下存在一定的偏差。为了进一步分析较大地震的震源性质，并确定本课题的震源机制解的可靠性，我们同时发展了同步源全矩张量反演方法，对实验场的地震震源性质进行研究。
3.1同步源全矩张量反演算法的技术思路
根据表示定理和点源近似，在台站记录到的位移可以写成震源矩张量时间序列和格林函数时间序列的卷积。
（1）
理论上可以直接求解公式（1）获得矩张量时间序列，但这个反演过程极不稳定。在传统的震源机制反演过程中，一般假设各个矩张量独立分量的时间序列相同，即假设：
（2）
且在反演前一般事先给定震源时间函数。这样就只需要反演地震矩张量或者断层的走向、倾角和滑动角。三角函数或梯形函数是比较常用的近似震源时间函数。这样就可以写成线性方程组的形式，直接求解地震矩张量。然而，的设定人为性较强，给的不合适就会影响反演结果。
同步源全矩张量反演方法将同时反演全地震矩张量（包含有ISO，CLVD和DC）和震源时间函数。反演步骤如下：
（1） 给定第k步震源时间函数，然后利用以下公式（3）反演出
（3）
（2） ）利用第k步的，使用以下公式（4）反演出第k+1步的
（4）
（3）重复第1步和第2步，直到波形残差足够小或超过预设的迭代次数。最终反演得到和。
（4）利用矩张量和断层面解（走向、倾角、滑动角）的关系求得断层面解。
3.2 数值算例
为了验证算法的可靠性和稳定性，设定了如图4所示的震源台站分布模型，并使用如表所示的速度结构模型合成理论地震记录然后用于测试。假设震源的两组断层面解分别为（332,57，-105）和（178.1961 35.8949 -68.2704）。我们进行了两组测试，算例一的震源时间函数为一等腰三角形；算例二的震源时间函数为两个部分重叠的三角形。测试反演用到的速度模型列在表1中。
图4. 同步源全矩张量反演测试模型震源台站分布图
算例一测试比较如图5所示，算例二的测试比较如图6所示。测试结果列在表2和表3中。测试结果表明，同步源全矩张量反演算法能够很好的反演得到真实的地震矩张量和震源时间函数。其中断层面解的误差在1°以内。每个台站的反演波形和观测波形相似度高达0.97以上。
表1 测试用速度模型
深度（km）
密度（g/cm3）
Vs(km/s)
Vp(km/s)
Qs
Qp
2.1
1.2
2.5
650
300
2.4
2.1
4.0
650
300
2.8
3.6
6.2
650
300
19
2.9
3.7
6.6
650
300
35
3.1
4.0
7.3
650
300
46
3.4
4.7
8.2
650
300
(a)
(b)
(c)
图5 算例一反演结果，（a）观测（蓝线）和反演（红色虚线）得到的波形记录对比；（b）拟合残差随迭代次数的变化；（c）反演得到的震源时间函数（蓝线）和真实震源时间函数（红色虚线）及滤波后的真实震源时间函数（黑线）的对比。
表2 算例一反演得到的震源机制和真实震源机制对比
真实震源机制
断层面解一
断层面解二
走向
倾角
滑动角
走向
倾角
滑动角
0.15
332
57
-105
178.1960
35.8949
-68.2704
反演得到的震源机制
0.1427
331.224
57.346
-105.917
179.0821
35.9364
-66.8311
(a)
(b)
(c)
图6 算例二反演结果，（a）观测（蓝线）和反演（红色虚线）得到的波形记录对比；（b）拟合残差随迭代次数的变化；（c）反演得到的震源时间函数（蓝线）和真实震源时间函数（红色虚线）及滤波后的真实震源时间函数（黑线）的对比。
表3 算例二反演得到的震源机制和真实震源机制对比
真实震源机制
断层面解一
断层面解二
走向
倾角
滑动角
走向
倾角
滑动角
0.27
332
57
-105
178.1960
35.8949
-68.2704
反演得到的震源机制
0.2434
331.424
57.114
-105.719
178.8243
36.0655
-67.2659
六、模型完成情况与成果
根据项目要求，我们分三个组分别对实验场的震源机制进行了反演。已完成实验场区2009~2017年来ML≥4级的中强地震的震源机制解目录。
6.1 郑勇课题组研究结果
郑勇课题组利用gCAP方法总共得到了634个地震事件的震源机制解，但考虑台站的方位角分布以及台站数量不少于8个，最终得到592个结果比较好的震源机制解结果。采用马文涛等（2004,2008）震源机制解分类方法，结果展示如图7所示，其中走滑断层地震223个，正断层地震112个，逆断层地震257个。
图7 2009-2017震源机制解gCAP反演结果。黑色为走滑地震，蓝色为拉张地震，红色为逆冲地震。
为了更好的了解川滇地区近几十年的震源机制解的时空变化，我们收集了1970-2008年Harvard CMT 的震源机制解；谢祖军等（2013），吕坚等（2012b）的结果，得到了1970-2017年川滇地区震源机制解分布图，如图8。
图8 1970-2017震源机制解结果。2009年以前的结果摘自Harvard GCMT结果。
为了展示拟合的情况和结果的可靠性，图9、图10具体展示了201606161512267、200901021100104两个地震事件的计算结果，图中左边是波形拟合结果，右边两个图是台站分布以及由残差随深度变化确定的震源深度。可以看到，波形拟合度很高，绝大多数波形的拟合相关性系数在0.9以上。同时，深度的拟合情况也很好，满足单一极小，可以将矩心深度误差控制在1km左右。
图9 201606161512267地震事件的波形拟合情况。左边为波形拟合情况，黑色为观测波形，红色为拟合波形。右上图为台阵分布图和地震位置，右下为拟合误差随深度变化情况，可以看到地震的最佳拟合深度为3km左右。
图10 200901021100104地震事件的波形拟合情况。左边为波形拟合情况，黑色为观测波形，红色为拟合波形。右上图为台阵分布图和地震位置，右下为拟合误差随深度变化情况，可以看到地震的最佳拟合深度为18km左右。
6.2 姚华建课题组研究结果
根据数据的信噪比和初至的清晰程度，姚华建学科组对446个地震事件进行了震源机制反演。图11给出了事件200901021100104的反演结果。图12给出了所有446个地震事件的震源沙滩球。
(a)
(b)
(c)
图11 事件200901021100104的反演结果，（a）台站分布和震源沙滩球；（b）为波形对比；（c）各分量的震源速率函数。
图12 姚华建学科组446个地震事件反演得到的震源机制沙滩球。
6.3 储日升学科组结果
储日升学科组利用CAP方法，收集2009-2017年Ml>=4.0地震数据进行反演，共获得536个地震的震源机制解、矩震级和深度。
六、模型验证（测试）与精度评价
为了验证震源机制反演的结果，我们分成了三个独立的小组，使用不同方法进行震源机制反演。我们首先将Mw≥5级的地震的震源机制解与Global CMT的结果进行了对比，两者的震级、走向、倾角和滑动角都非常接近。Global CMT给出了更大的震级和震源深度。这是由于震级主要由振幅决定，越深的事件需要越大的震级才能在地表产生相同的观测振幅。随机挑选了一些地震，选用不同速度模型对其进行可测试，得到的结果仅有微小的变化，表明震源机制解释可靠的。
针对小于5级地震，我们以T轴夹角作为评判两组结果的差异的标准，假定T轴夹角<40两组震源机制基本接近，表4给出了三个学科组两两对比结果接近的事件数目。
表4 三个小组的震源机制反演结果对比
储日升结果vs.郑勇结果
郑勇结果vs.文健结果
储日升结果vs.文健结果
118
67
90
通过分析，造成三个小组的震源机制反演结果差异的原因可能有：
1） 速度模型不一致，储日升和文健均采用相同的水平层状模型，而郑勇采用的是不同事件不同水平层状模型；
2）震源深度不一致，郑勇和储日升反演时对深度进行了搜索，而文健未做深度搜索；
3）反演方法不同，储日升和郑勇均基于双力偶源假设，直接搜索断层走向、倾角和滑动角，而文健未要求地震事件是双力偶源，采用的是全矩张量反演获得地震矩张量然后获得走向、倾角和滑动角。
总体而言， 三家的反演结果在较大地震（Ml>4.5级）上一致性很好，4.0~4.5级地震总体上比较一致。表明本课题的震源机制解及其矩心深度比较可靠。
七、模型使用说明
震源机制解为双力偶解，其深度为震源的矩心深度，与日常基于震相到时定位得到的震源深度的含义不同，不宜混用。
八、应用案例
基于震源机制解产出的深度结果，将川滇地区震源矩心深度的分布模式与南加州地区的震源矩心深度的分布模式进行了对比：川滇地区具有更多的浅源地震。深度分布表明川滇地区脆性-韧性过渡带较南加州地区浅，这可能与新生断层的存在有关。
参考文献
Zhu, L., & Rivera, L. A. (2002). A note on the dynamic and static displacements from a point source in multilayered media. Geophysical Journal International, 148(3), 619-627.
Zhu, L., & Ben-Zion, Y. (2013). Parametrization of general seismic potency and moment tensors for source inversion of seismic waveform data. Geophysical Journal International, 194(2), 839-843.
Wu, C., & Takeo, M. (2004). An intermediate deep earthquake rupturing on a dip‐bending fault: Waveform analysis of the 2003 Miyagi-ken Oki earthquake. Geophysical research letters, 31(24).
罗艳, 倪四道, 曾祥方, 郑勇, 陈棋福, & 陈颙. 2010a. 汶川地震余震区东北端一个余震序列的地震学研究. 中国科学: 地球科学(6): 677-687
郑勇, 马宏生, 吕坚, 倪四道, 李迎春, & 韦生吉. 2009. 汶川地震强余震 (Ms≥ 5.6) 的震源机制解及其与发震构造的关系. 中国科学: D 辑(4): 413-426.
郑勇,葛粲,谢祖军,YingJie YANG,熊熊,许厚泽.芦山与汶川地震震区地壳上地幔结构及深部孕震环境[J].中国科学:地球科学,2013,43(06):1027-1037.
谢祖军，金笔凯，郑勇，葛粲，熊熊，熊诚，许厚泽.近远震波形反演2013 年芦山地震震源参数.中国科学：地球科学，2013,43:1010-1019
吕坚, 王晓山, 苗春兰, 苏金蓉, & 靳玉科. 2012b. 汶川余震序列北川以北段的震源分布特性与破裂复杂性研究. 大地测量与地球动力学, v.32(05): 17-21.
马文涛,徐锡伟,曹忠权,于贵华,李海鸥,徐长朋.震源机制解分类与川滇及邻近地区最新变形特征[J].地震地质,2008,30(04):926-934.
马文涛,徐锡伟,徐平,于贵华.地震三角形分类图解法与华北地区地震成因分析[J].地球物理学进展,2004(02):379-385.