地球上有多种引起地面震动的事件，例如天然地震、爆破、陨石撞击和矿山塌陷等，震源深度是区别爆破等特殊浅源事件与天然地震的关键信息(Dreger and Woods, 2002).一般来说，大部分天然构造地震的孕震层通常分布在中上地壳的脆性区，深度范围约为5~15 km (Scholz, 1998).但是一些爆破等非天然地震往往位于浅地表.如果震源深度大于10 km时，则该事件极大概率为天然构造地震(Blandford, 1977).在爆破研究中，震源深度是爆炸当量研究关键的参数(Bowers et al., 2001;Murphy et al., 2013)，在不同震源深度假定条件下，估计的当量结果存在较大差异(Murphy et al., 2013;Zhang and Wen, 2013; Zhao et al., 2008, 2017;赵连锋等，2017；谢小碧和赵连锋，2018).核试验震源深度还是评价地面沉降的关键参数，地面沉降范围随着震源深度的增加而减小，在未成岩的沉积土层中，当深度大于215 m时，直径为15 m试验空腔造成的地面沉降几乎可以忽略(Houser, 1969).

基于地震学方法，针对浅源地震震源深度的准确测定研究，前人已开展了大量的工作.如利用地震波走时可以反演地震深度，但是在台网稀疏情况下，仅利用远震P波到时获取的震源深度与发震时刻存在折中，测得的结果误差较大(Fisk, 2002).而利用pP与P的到时差则能提高震源深度的测定精度(Cohen, 1970;Bakun and Johnson, 1973;Douglas et al., 1999)，Gitterman等(2016)根据远震P波的零频率特征和pP深度震相测定了2006—2009年朝鲜三次核试验深度(~2 km).Murphy等(2011)通过远震P波振幅比研究了2006年和2009年朝鲜核试验震源深度，结果分别为200 m和500 m.基于区域台站波形数据，也可以反演核试验震源深度(Minson and Dreger, 2008). Ford等(2009a)基于区域地震资料反演了多个爆炸源以及塌陷和地震的震源参数，通过分析反演结果，识别了天然地震和爆炸源以及塌陷事件，并发现在数据信噪比较高时，可以获取可靠的爆炸源震源机制.基于波形数据，学者们利用gCAP以及ISOLA方法反演了2017年朝鲜M6.3事件震源参数，并对该事件成因进行了分析(Han et al., 2017; Liu et al., 2018; Wang et al., 2018).

在利用地震波形反演震源深度时，大多数方法没有考虑地形的影响.然而Meunier等(2008)研究发现山脊、山谷起伏地形会对S波形造成明显放大效应.Wang等(2016)发现地形起伏可产生与P振幅相当尾波信号，利用P波以及尾波信息可以提高地震定位精度.地形起伏还会造成两次事件震相的时间移动以及面波振幅比异常(Rodgers et al., 2010;Ichinose et al., 2017).在地形起伏较大的地区，如本文研究区域(图 1)，震中高程约为2000 m，而台站分布在40~1800 m之间，台站高程相差约1700 m.浅源地震深度与区域台站的高程差量级相近，地形起伏对震源深度反演所造成的影响需要进行研究.

图 1

Fig. 1

图 1 研究区域地形图以及台站(三角)分布 (a)SEM正演计算区域，红色圆点为2017年9月3日M6.3事件震中.(b)为(a)中方框放大情形. Fig. 1 The topography and stations (triangle) in the study region (a) The simulation region, the red dot denotes epicenter of the M6.3 event on 3 September 2017.(b)Magnified view for the box in (a).

本文进行了一系列正演模拟以及测试，对比不同地形条件下的波形差异，进而探讨地形起伏对浅源地震深度反演造成的偏差.然后基于不同地形条件下的格林函数，利用波形方法反演朝鲜2017年9月3日M6.3事件震源深度.并分析了速度模型以及震源复杂性对理论测试结果的影响，还进一步讨论了速度模型误差对震源深度测定的影响.

地震图数值模拟方法主要有有限差分(FDM)、有限元(FEM)和谱元法(SEM)等多种方法，其中SEM结合了有限元的灵活性与伪谱法的准确性，可处理复杂的地形(Komatitsch and Tromp, 1999;Komatitsch et al., 2004).而针对中小地震震源参数反演，也发展了多种基于波形的方法.如体波振幅以及振幅比(Hardebeck and Shearer, 2003)，近场全波形MTINV(Dreger and Kaverina, 2000;林向东等, 2013)，近场体波与面波Cut and Paste(CAP)(Zhu and Helmberger, 1996).其中，CAP方法在九江—瑞昌、汶川、芦山、九寨沟和南北地震带等中强地震震源参数研究中得到了广泛的应用，并取得了丰富的研究成果(吕坚等，2008；郑勇等，2009；Zheng et al., 2010; Xie et al., 2013;李志伟等，2015；易桂喜等，2017). Zhu和Ben-Zion (2013)在CAP的基础上加入了各向同性(ISO)和补偿线性单力偶(CLVD)分量，从而提出了gCAP方法，并在震源矩张量反演中得到广泛应用(Ross et al., 2015；李圣强等, 2013；张广伟和雷建设, 2015；Wang et al., 2018).本文将使用SEM和gCAP方法来研究地形对震源深度反演的影响.

在利用SEM方法计算合成理论地震波形时，避免吸收边界对波形模拟的影响，将震中(纬度41°，经度129°)放在计算区域的中心，设置经纬度范围(纬度30°—52°，经度117°—139°).为了消除三维复杂结构的影响，在波形方法反演中常用的滤波周期为5 s以上(Dreger and Woods, 2002; Zhu and Helmberger, 1996).在SEM计算中应当包含上述频段，将地震波准确计算的最小周期设置为3 s左右.根据Komatitsch等(2004)计算准则来划分网格，网格步长D应该小于1倍最小波长，若最小速度为3.6 km·s^-1，波长即为10.8 km，为了得到准确的结果，在实际计算中D取值为5 km.台站震中距均小于500 km，计算时长设置为300 s.震源深度以0.5 km为间隔，在0.5~2 km的范围内取值，包含浅源地震事件的分布深度.首先基于简化均匀速度模型，P和S波速度分别为6.4 km·s^-1和3.6 km·s^-1，并选用各向同性爆炸源作为震源机制解来合成理论地震图，在讨论部分我们将测试多层速度模型以及震源复杂性的影响.

基于同样的一维速度模型，利用频率-波数(FK)方法(Zhu and Rivera, 2002)计算不同震源深度的格林函数，以SEM计算的理论地震图作为观测值，然后使用垂向、径向分量的体波和面波数据，其滤波频带范围为0.05~0.16 Hz，并采用gCAP方法进行反演.从测试结果来看，在没有添加地形时，一维反演均能很好恢复正演输入的结果.以震源深度1 km为例，波形拟合系数均为0.99(图 2).其他震源深度情形(0.5 km，1.5 km, 2 km)，反演结果以及波形拟合情况与震源深度1 km时类似.我们还测试了天然地震震源深度5~20 km，间隔5 km的情形，其输入震源深度均能得到准确的恢复(表 1 Test 1).

图 2

Fig. 2

图 2 震源深度1 km时，层状速度模型理论测试结果 (a)波形拟合误差随震源深度分布情况；(b) FK理论(黑色虚线)和SEM正演计算(红色)波形拟合图，波形左侧台站名左下角数字为方位角，右下角数字为震中距，波形下方的两行数字分别表示理论地震图相对观测地震图的移动时间及二者的相关系数. Fig. 2 The test result for the case with 1 km focal depth based on flat layer model (a) Waveform mismatch versus centroid depth; (b) Comparison between synthetic waveforms calculated based on FK (black dashed line) and SEM program (red). The letters in the left of the figure represent name of stations and the numbers below the stations are their azimuths and epicentral distances. Numbers below each trace denote relative time shift between synthetics and data (upper numbers) and their correlation coefficient in percentage (lower numbers).

表 1 (Table 1)

表 1 平层与起伏地形条件下正演测试结果 Table 1 The results of focal depth based on the flat model and model with topography 震源深度输入值 (km) 相对于海平面深度 (km) Test 1 (km) Test 2 (km) 0.5 -1.326 0.5 0.4 1.0 -0.826 1.0 0.8 1.5 -0.326 1.5 1.3 2.0 0.174 2.0 1.8 5.0 3.174 5.0 4.7 10.0 8.174 10.0 9.5 15.0 13.174 15.0 14.4 20.0 18.174 20.0 19.5

表 1 平层与起伏地形条件下正演测试结果 Table 1 The results of focal depth based on the flat model and model with topography

在平坦地形速度模型的基础上，我们加入实际高程数据，来研究地形对于震源深度反演的影响.地形数据来自ETOPO15(http://www.ngdc.noaa.gov/mgg/global/)，为了检验地形添加的正确性，我们将SEM中恢复的地形形态与topo15地形数据进行了对比(图 3)，可以看出两者整体吻合情况较好.在地形网格划分中，使用5 km步长能较好再现研究区域大致地形.接下来，我们对比了基于不同地形条件下合成的理论地震图.以震源埋深1 km为例，图 4中展示GUS(震中距490 km)和JCT(震中距246 km)两个台站波形对比情况，滤波范围为0.05~0.2 Hz，包含实际地震数据处理频带.相对于平坦地表，加入地形之后，主要震相(P、S)形态较为相似，振幅有一定改变.并且地形使得地震波形复杂化，产生一些尾波，这与前人的研究结果较为一致(Wang et al., 2016).

图 3

Fig. 3

图 3 (a) 研究区域地形图，红色圆点和黑色三角分别代表震中位置和示例台站；(b) AA′剖面SEM输出地形(蓝色)和topo15地形(红色)对比 Fig. 3 (a) The topography in the study region, red dot and black triangles denote the epicenter and stations, respectively; (b) Comparison of topography between output data of SEM (blue) and topo15 (red) on AA′ profile

图 4

Fig. 4

图 4 平坦地形(黑色)和起伏地形(红色)条件对地震波形的影响：(a) GUS台站波形图; (b) JCT台站波形图 Fig. 4 Comparison of synthetic waveforms based on flat layer model (black) and the model with topography (red). Waveforms for GUS station (a) and JCT station (b)

加入地形之后，可以看到，震源区地形较高(2000 m左右)，而台站高程均小于2000 m.天然地震的震源深度一般分布在5~20 km，为此我们测试了震源深度为5，10，15和20 km四种情形下，地形对震源深度反演的影响.在反演时FK使用一维速度模型，深度间隔为0.1 km.以输入震源埋深为10 km为例，反演得到震源深度为9.5 km，比真实震源浅0.5 km.震源埋深5, 15和20 km均能得到类似的结果(表 1 Test 2)，这表明在天然地震震源深度反演中，基于平层假设反演得到的震源深度与基于震源区自由地表的深度接近.震源深度反演结果与输入值偏差在0.5 km左右.

针对研究区域的浅源地震来说，震源有可能在台站的上方，为此还测试了震源埋深为0.5，1，1.5和2 km四种情形，反演得到的结果比输入值小0.2 km左右(图 5).在震源较浅时，基于平层速度模型假设也可以得到近似震源埋深的结果，其结果比实际值一般要浅.对于深度在500~2000 m之间的浅源事件来说，利用一维速度模型反演得到的结果基本可靠，能为核试验以及塌陷等浅源事件成因分析提供参考.

图 5

Fig. 5

图 5 地形起伏条件下，震源埋深0.5~2.0 km时深度反演测试结果 (a)震源埋深0.5 km时反演测试结果；(b, c, d)分别为震源埋深1 km、1.5 km和2 km时反演测试结果. Fig. 5 The test results of cases inverted based on topographical model (a) The test result for the case with 0.5 km focal depth; (b, c, d) Same as (a) but for 1 km, 1.5 km and 2 km.

2017年9月以来，朝鲜半岛北部发生多次有感地震.其中，2017年9月3日11时30分(北纬41.3°，东经129.1°)发生6.3级地震，震源深度0 km.由其产生的地震波在中国东北和华北的部分地区产生了强烈震感.朝鲜政府宣称这是氢弹试验，并取得了成功，这是朝鲜进行的第六次核试验，赵连锋等(2017)利用纵横波振幅比值确认这是一次爆破事件.前人对该地区浅震震源机制进行了大量研究，得到了丰富的研究成果，在计算理论地震图时，MDJ2速度模型被广泛应用(Ford et al., 2009b; Cesca et al., 2017; Liu et al., 2018; Wang et al., 2018).

基于MDJ2速度结构模型，利用FK频率-波数方法(Zhu and Rivera, 2002)计算了深度0.1~3 km之间，步长为0.1 km的格林函数.根据最大震中距范围，设置体波和面波最大时窗长度分别为40 s和80 s，为了降低噪声以及地壳各向异性的影响，将P波和面波分别进行0.05~0.12 Hz和0.05~0.1 Hz带通滤波.结果显示，地震震源深度为0.8 km.最佳拟合情况波形如图 6a所示，大部分台站波形拟合较好，其中大部分Rayleigh面波拟合系数在0.95以上，而切向分量波形拟合较差，这是因为切向能量较弱，信噪比较差导致的.这也反映了震源存在较大的各向同性分量.震源机制反演结果显示：ISO以及CLVD反演结果分别为0.93和-0.44，根据Zhu和Ben-Zion(2013)计算公式显示，该地震含有约86%各向同性分量和少量的CLVD以及双力偶分量，与Han等(2017)和Wang等(2018)研究结果一致.

图 6

Fig. 6

图 6 2017年9月3日M6.3事件理论(红色)与观测(黑色)波形拟合情况 (a)平层速度模型; (b)带地形的速度模型，详细说明参照图 2. Fig. 6 Comparison between observed (black) and synthetic (red) waveforms of the M6.3 event on 3 September 2017 (a) The case based on the flat model; (b) The same as (a) but for model with topography, see caption of Fig. 2 for details.

为了考察地形对朝鲜事件的震源深度反演影响，基于MDJ2模型，我们利用SEM计算了实际地形条件下的格林函数.Wang等(2018)通过InSAR数据研究了该次地震震源深度为450 m，Tian等(2018)根据相对定位后地震的精确位置以及地形资料，推断该地震的震源深度约760 m，Liu等(2018)得到震源深度最优深度为1.5 km，赵连锋等(2017)通过精确定位获取震中高程, 利用其与遂道洞口的高程差估计出最小埋深，推断该次地震震源深度约1~2.4 km.综合不同的研究结果，我们计算了震源深度为0.5、1 km和2 km三种情形，每个震源深度的格林函数有18个(Zhu and Zhou, 2016).滤波参数以及台站分布与不带地形情形一致.反演结果显示震源深度1 km处波形拟合误差最小，加入地形之后，大部分台站波形拟合情况有所提高.加入地形之后，基于更加准确的格林函数，可以提高地震震源参数反演的准确性.

为了考察层状速度结构对震源深度反演的影响，我们以MDJ2分层模型为基础进行测试.以埋深1 km为例，在不考虑地形时，震源埋深能很好的恢复.加入地形之后，震源深度反演结果为0.8 km，与输入值相差0.2 km，这与单层速度模型情况下地形对震源深度反演结果影响大小一致.在理论数据测试中，我们使用了各向同性源，而在实际案例中，震源往往包含CLVD以及双力偶分量.为此，基于多层速度模型，我们讨论了震源复杂性的影响.以M6.3地震实际反演的震源矩张量结果作为输入，震源深度反演结果也浅于输入值，与各向同性源情形类似.我们还测试了三类典型断层情况下，地形对震源深度反演的影响，以震源深度1 km为例，机制解分别为断层倾角为45°的倾滑地震(45°和90°)、垂直倾滑地震(90°和90°)和垂直走滑地震(90°和0°).不同震源机制解正演测试显示，在有地形起伏情形下，震源深度反演结果均比输入值略浅.为了考察地形起伏大小的影响，我们还以爆炸源为例，测试了1和2倍真实地形情况下的情形，发现2倍地形条件下反演得到的震源深度比输入值偏浅约500 m左右.起伏地形会对S波形造成明显放大效应以及产生与P波振幅相当尾波信号(Meunier et al., 2008；Wang et al., 2016).地形起伏加剧时，波形改变程度也随之增大.在反演时，振幅放大的S波以及尾波信号包含在反演时窗之内，可能会对震源反演结果造成影响.

对浅源地震来说，浅层速度结构对其深度反演结果影响可能较大.MDJ2速度模型来自于面波信息，该模型在研究区域有较好适应性.为了考察该速度模型的准确性，我们提取了M6.3地震在不同台站的1~10 s群速度频散曲线，并计算了MDJ2速度模型浅层(< 5 km)10%速度扰动的理论频散曲线.结果显示，理论计算值能很好的预测观测结果，这表明MDJ2模型浅层5 km处的速度不确定性在10%之内(图 7).为此，我们将MDJ2模型浅层5 km处的V_P和V_S速度增减10%，形成两个模型.基于这两个模型重新反演了该地震的震源深度，台站分布和数据处理流程与前文一致.震源参数反演结果如表 2所示，震源深度不确定性在1 km左右.

图 7

Fig. 7

图 7 MDJ2速度模型不确定分析 (a) MDJ2(黑色)及其浅层5 km增加(蓝色)以及减少(绿色)10%速度结构，(b)观测(灰色)以及不同模型下理论面波群速度曲线 Fig. 7 Uncertainty analysis of the MDJ2 velocity model (a) MDJ2 (black), the model with both VP and VS perturbations of +10% (blue) and -10% (green); (b) Observed (gray) and synthetic group velocities calculated based on the three models

表 2 (Table 2)

表 2 2017年9月3日M6.3事件地震参数反演误差分析 Table 2 Uncertainty analysis of the focal depth of the M6.3 event on 3 September 2017 Velocity Model ζ ISO(%) Centroid depth (km) +10% VP and VS 0.90 81% 1.2 -10% VP and VS 0.94 88% 1.9

表 2 2017年9月3日M6.3事件地震参数反演误差分析 Table 2 Uncertainty analysis of the focal depth of the M6.3 event on 3 September 2017

除了震源深度，震源矩张量信息也可以帮助理解震源的类型(Ford et al., 2009a).为了更好地展示反演得到的震源机制解，我们将结果投影在Hudson plot图上(Hudson et al., 1989).从图 8中可以看出，震源位于爆炸源区域，基于不同地形条件得到的震源矩张量比较接近，同时与Liu等(2018)反演的结果一致.相对于Wang等(2018)反演结果中负的CLVD，而本文反演结果为正的CVLD分量，这可能是使用了不同台站波形数据以及地壳模型造成的.Cesca等(2017)发现矩张量中水平分量(M_xx和M_yy)与垂直分量(M_zz)之间存在折中，会造成CLVD含量发生变化.我们还假定纯爆炸源震源机制对实际波形进行拟合，发现Love面波拟合情况较差.爆炸源中的Love波成分可能是由约14%的双力偶源分量、地壳各向异性以及爆炸源的非均匀性产生(Ford et al., 2009b; Vavryčuk and Kim, 2014).

图 8

Fig. 8

图 8 CLVD和ISO分量Hudson示意图 黑色与红色圆点分别代表基于平层和带地形速度模型反演的结果.绿色和黄色圆点分别代表Liu等(2018)和Wang等(2018)反演的结果. Fig. 8 The Hudson plots for the CLVD and ISO components The black and red dots indicate the preferred moment tensor inverted using flat and topographical velocity model. The green and yellow dots display the moment tensors from independent inversion result of Liu et al. (2018) and Wang et al. (2018).

近年来，现代观测技术的发展丰富了核试验等浅源地震震源深度的研究.例如，Wei(2017)基于InSAR观测的地表形变，反演得到2016年1月6日朝鲜核试验的深度范围为420~700 m；Zhang和Wen(2013)根据高精度的卫星影像估算了2009年和2013年的核试验震源深度分别为610 m和430 m.此外，利用区域震相的尾波波形特征，通过单一地震台站拟合区域波形全波列的包络，可以同时获得核试验的当量和埋深(Pasyanos et al., 2012;林鑫和姚振兴，2016).将上述方法与区域波形反演方法相结合，快速反演得到震源深度以及震源性质，可以为核试验以及塌陷等浅源地震事件研究提供更准确的约束.

本文研究了地形对浅源震源深度反演的影响.首先进行了正演计算与理论测试，在没有地形时，基于层状均匀模型可以得到正确的结果.当存在地形起伏时，基于层状均匀模型也可以反演得到与震源埋深较为一致的结果，但会系统偏浅0.2 km左右.基于MDJ2速度结构模型，我们反演了2017年9月3日朝鲜M6.3事件震源参数，其深度约为0.8 km.加入地形之后，波形拟合度提高，震源深度为1 km左右，与均匀层状模型假设反演结果相差不大.案例研究表明，基于均匀层状速度模型，利用近震资料反演的震源深度基本可靠，能为浅源地震成因分析提供关键约束.

但是本文仅使用了一维速度模型，没有考虑地壳各向异性的影响，而且本文测试台站高程均小于震源区海拔.而针对于地形起伏更大的地区，如龙门山断裂带，周围台站高程变化较大，介质横向变化剧烈，Zhu和Zhou(2016)等发展了gCAP3D方法，并基于三维速度结构模型反演了该地区2013年芦山地震序列的震源参数.在龙门山断裂带东北缘还存在一些震源深度为3 km左右的极浅源地震(Luo et al., 2010; Zheng et al., 2010)，较大地形起伏以及三维速度结构对这些地震震源深度反演的影响需要进一步研究.