地震学报  2007, Vol. 23 Issue (3): 229-239
引用本文
许力生, 蒋长胜, 陈运泰, 李春来, 张天中. 2007. 2004年首都圈地区中小地震的矩张量反演. 地震学报, 23 (3): 229-239.  
Xu Lisheng, Jiang Changsheng, Chen Yuntai, Li Chunlai, Zhang Tianzhong. 2007. MOMENT TENSOR INVERSION OF THE 2004 SMALL- MODERATE SIZE EARTHQUAKES IN THE CAPITAL REGION. Acta Seismologica Sinica , 23 (3): 229-239.   
2004年首都圈地区中小地震的矩张量反演
许力生 , 蒋长胜, 陈运泰, 李春来, 张天中    
(中国北京100081中国地震局地球物理研究所)
2006-08-16收到初稿， 2007-02-06决定采用修改稿.
基金项目： 国际973项目（2001CB711005,2004DFA06000）和中国地震局重点攻关项目（2004BA601B）资助.中国地震局地球物理研究所论著07AC1007.
通讯作者,E-mail:xuls@cea-igp.ac.cn
摘要：摘要 通过首都圈地区宽频带数字地震波形反演， 获得了2004年发生在该地区的51次中小地震的矩张量解， 由此确定了这些地震的标量地震矩、 矩震级、 双力偶分量和补偿线性矢量偶极分量的大小以及断层面参数和应力轴参数， 并通过数值试验对反演结果进行了评价.
关键词： 首都圈地区    中小地震    矩张量反演    
MOMENT TENSOR INVERSION OF THE 2004 SMALL- MODERATE SIZE EARTHQUAKES IN THE CAPITAL REGION
Xu Lisheng , Jiang Changsheng, Chen Yuntai, Li Chunlai, Zhang Tianzhong    
(Institute of Geophysics, China Earthquake Administration, Beijing 100081, China)
Abstract: The moment tensor solutions of 51 small-moderate size earthquakes occurred in the Capital Region in the year of 2004 are obtained by inverting the broadband waveform data. Accordingly, other source parameters, such as scalar seismic moments, moment magnitudes, double-couple components and compensated-linear-vector-dipole components, are determined as well as fault parameters and stress-axis parameters. The inverted results are evaluated by groups of numerical tests.
Key words: Capital Region    small-moderate size earthquakes    moment tensor inversion    
引言

地震的震源机制解在地球科学的发展中一直在发挥重要作用，而且必将发挥更大的作用. 早在20世纪60年代，地震的断层面解就被用来研究应力场方向(McKenzie，1969). 数十年来，有大量的研究都与地震的震源机制解有关(Tapponnier，Molnar，1976; Zoback et al，1987; Xu et al，1992; Xu，2001； Kubo，Fukuyama，2003). 早期的震源机制大都是由P波初动符号确定的，由于台站的分布所限，只有少部分地震的震源机制才能得以确定. 1970年代以来，波形资料可以用来确定地震的矩张量解及震源机制(Gilbert，Dziewonski，1975; Patton，Aki，1979)，而且可用的波形资料(如地球的自由振荡、 面波、 体波等各种震相)使能够确定震源机制的地震数目大大增加了.

当然，通过地震矩张量的反演，除地震的断层面解以外，还可以得到地震矩张量、 标量地震矩、 矩震级、 各向同性分量和补偿线性矢量偶极成分等其它震源参数，而这些参数正在被其它研究所用. 例如，地震矩-频度关系的研究(Frohlich，Davis，1993； Sornette et al， 1996，; Kagan，1997)，非双力偶源的全球分布的研究(Kawakatsu，1991)，地震分类的研究(Reasenberg，1999)，地壳势能变化的研究(Tanimoto，Okamoto，2000)，以及震源机制与时间关系的研究(Kagan，2000)等等.

1995—2000年，我国建立了许多区域或地方的数字地震台. 这些台站的建立使我们能够利用它们记录到的波形资料，通过矩张量反演方法确定中小地震的矩张量解、 震源机制和其它震源参数. 本文通过地震矩张量反演方法，确定了首都圈地区发生于2004年的中小地震的矩张量及其相关参数.

1 方法

在频率域里，用矩张量描述的位于坐标原点的地震点源在某一观测点r引起的地震位移可以表示为

式中，ω为角频率，ui(r，ω)为地震的观测位移谱，Mjk(ω)为地震矩张量谱，Gij，k(r，ω)为格林函数谱. 即，地震位移的谱等于格林函数的谱和矩张量的谱的乘积； 反过来，矩张量谱等于地震位移谱和格林函数谱的商.

对于中小地震而言，由于震源持续时间较短，我们总可以找到一个合适的频段，以致在这个频段内实际地震的震源时间函数可以当作一个简单的脉冲，从而可以用狄拉克-δ函数来代表震源时间函数. 在这种情况下，式(1)可以写成

即描述震源的地震矩张量与时间或频率无关.

在实际工作中，我们可以选择适当的带通滤波器，使观测资料和格林函数位于相同的频段，并满足地震的震源时间函数为狄拉克-δ函数的条件. 如图 1所示，在滤波之前，观测波形与格林函数的谱呈现出不同的形状，观测波形的谱的拐角频率不同于格林函数的谱的拐角频率(图 1b). 然而，在滤波之后，观测波形的谱与格林函数的谱的形状非常相似，它们的拐角频率也基本一致(图 1c). 在这种情况下，实际地震的震源时间函数即可近似为狄拉克-δ函数. 与此同时，我们注意到，滤波后的观测波形(图 1a)和格林函数的各分量(图 1d)仍然保留着它们的基本特征，没有丢失用以确定矩张量解的长周期信息.

图 1 滤波前后的观测震相与相应的格林函数及其振幅谱 (a) 滤波前(细线)与滤波后(粗线)的观测震相； (b) 滤波前观测震相的振幅谱(粗线)与滤波前格林函数的 6个分量的振幅谱(细线)； (c) 滤波后观测震相的振幅谱(粗线)与滤波后格林函数6个分量的 振幅谱(细线)； (d) 滤波前(细线)与滤波后(粗线)的格林函数
2 观测资料与格林函数

1995—2000年，中国地震局在北京、 河北和天津地区建立了107个数字地震台. 其中48个台架设了宽频带数字地震仪，59个台站架设了短周期数字地震仪. 在这项工作中，我们处理了2004年度首都圈地区(113°E~121°E，37°N~41.5°N)48个宽频带地震仪记录到的地震资料，挑选了51个地震事件，其中震级最小的地震为ML2.0，它们满足初至的直达波同时被3个或3个以上台站记录清楚的条件. 这些宽频带地震仪的频率范围为0.05~20 Hz，但是，我们对观测资料使用了0.2~5 Hz的带通滤波器，以剔除较低频和较高频噪声而不破坏用于矩张量反演的信息.

格林函数是利用平面分层介质的反射、 折射率方法计算的(Kennett，1983).当地的地壳模型是在中国地区平均地壳速度模型的基础上(中国地震局震害防御司，1992)，通过观测地震图和合成地震图到时的比较建立的. 为了得到一个平均的分层地壳模型，我们选择了不同震中距的实际记录，然后计算相应震中距的合成地震图，并与观测地震图进行比较. 如此反复，我们可以得到一个模型，使基于它的合成地震图的P波到时和S波到时与观测地震图的基本一致. 经过调试，我们最终采用了3层地壳模型： 第一层为200 m的覆盖层，P波速度为2.40 km/s，S波速度为1.0 km/s，密度为1.81 g/cm3，QP=0.03，QS=0.06； 第二层的厚度为15 km，P波速度为5.71 km/s，S波速度为3.20 km/s，密度为2.7 g/cm3，QP=0.001，QS=0.002； 第三层的厚度为25 km，P波速度为6.83 km/s，S波速度为3.70 km/s，密度为3.0 g/cm3，QP=0.001，QS=0.002. 由于实际的地壳分层界面未必是水平的，所以，不可能得到一个唯一的水平分层速度模型可以在所有的震中距上达到观测走时与合成走时的一致. 由于水平分层模型只是对实际复杂模型的简化，所以，不可能考虑次生震相的到时，而只能考虑直达P波与S波的到时. 由于在计算合成地震图时，我们不考虑震源机制对波形的影响，所以我们对观测地震图和合成地震图做了如下处理，以同时突出P波和S波： 以观测地震图或合成地震图的垂直分量为实部、 切向分量为虚部建立复时间序列，并计算此复时间序列的模(模时间序列)，然后进行比较. 图 2a是陡河台(DOH)记录到的震中距为21 km的一次地震的模时间序列与响应的合成图的模时间序列，图 2b是兴隆台(XIL)记录到的震中距为90 km的一次地震的观测图模时间序列与响应的合成图的模时间序列. 从图中可以看出，观测到时与合成到时可以得到相当好地吻合. 不过，也可以看到次生震相的到时很难得到吻合，而且震中距越大，次生震相越复杂，其到时也越难以吻合.

图 2 观测地震图与合成地震图到时的比较 (a) 震中距为21 km的观测地震图(实线)与合成地震图(虚线)的模时间序列； (b) 震中距为90 km的观测地震图(实线)与合成地震图(虚线)的模时间序列
3 地震矩张量反演

地震矩张量的反演是通过软件QUICKMT(许力生，陈运泰，2004)来实现的. 该软件适合于全波的反演，也适合于某一震相的反演. 在本工作中，我们只使用了直达P波. 首先，根据地震的大小和仪器的频率特性给定一个频率范围. 当程序开始运行时，带通滤波器的频率上限会由大到小逐渐变化，生成不同频段的带通滤波器. 这样，地震矩张量反演的过程在不同的频段进行，分别得到不同的结果. 最后，我们选择波形拟合得最好的频段，把这个频段的反演结果作为最终的结果. 例如，北京时间2004年10月11日22点3分30秒发生在39.598°N、 118.180°E的一次ML2.0地震的反演输出结果如下：

发震时刻: 2004-10-11，22:03:30；

事件纬度39.598°N、 经度118.180°E、 深度9.9 km；

台站数目为6、 震相数目为6、 分向为Z；

获得矩张量解的资料的频率上限为0.9375 Hz;

地震矩数量级为1.0×1011 N·m;

标量地震矩为10.08 N·m；

矩震级MW1.9；

矩张量m11，m12，m13，m22，m23，m33分别为-3.11，2.97，3.23，-2.08，-7.88，5.19 N·m；

爆炸成分为0；

补偿线性矢量偶极(CLVD)成分为-0.25 N·m；

最佳双力偶(DC)成分为10.08 N·m；

节面Ⅰ： 走向、 倾角、 滑动角分别为27°，75°，77°；

节面Ⅱ： 走向、 倾角、 滑动角分别为249°，20°，130°;

T轴： 本征值、 方位角、 倾角分别为10.33，280°，58°;

B轴： 本征值、 方位角、 倾角分别为-0.51，31°，12°;

P轴： 本征值、 方位角、 倾角分别为： -9.83，128°，29°.

台站 纬度 经度 震中距/km 震中距/° 方位角/°
DOH 39.752 118.328 21 0.19 36.55
QIX 40.199 118.311 67 0.61 9.49
CHL 39.758 119.095 80 0.72 76.96
JIX 40.076 117.496 78 0.71 312.41
XLD 40.402 118.102 89 0.81 355.76
TLK 40.146 119.046 95 0.86 50.33

以上的反演输出结果不但包含震源参数信息，如地震矩张量解、 标量地震矩、 矩震级以及双力偶解等，而且包括反演所用的资料信息，如资料的数量、 台站参数以及资料的频率上限等.

图 3为北京时间2004年10月11日22点3分30秒发生在39.598°N、 118.180°E的ML2.0地震的震源机制解的“海滩球”表示及观测波形与合成波形的比较. 可以看出，该震源机制解很好地解释了频率低于0.937 5 Hz的直达P波.

图 3 2004年10月11日22点3分30秒ML2.0地震的震源机制解的“海滩球”表示(震源球下半球投影) (a)及观测波形与合成地震波形(b)的比较. 图(b)中左边的数字分别表示观测波形的最大振幅、 相关系数和合成波形的最大振幅， 右边分别表示台站名称、 分向和震相

2004年1月—2005年1月在首都圈地区发生了200多次ML2.0以上的地震，但满足初至波为直达P波、 波形清楚且同时被3个以上台站记录到的地震事件只有51次. 我们对这51次地震的资料使用了同样的反演方法，得到了它们的地震矩张量和其它相关参数. 图 4给出了这些地震震源机制的“海滩球”表示，在表 1中列出了这些地震的震源参数.

图 4 首都圈地区的宽频带地震台分布以及2004年1月—2005年1月发生的ML＞2.0 的51次地震的震中分布及其震源机制解的“海滩球”表示(震源球下半球投影)

表 1 首都圈地区中小地震矩心矩张量解列表
4 台站方位覆盖与台站数目对反演结果的影响

众所周知，地震矩张量的反演结果受观测台站方位覆盖的影响，也受观测台站数目的影响. 为此，我们进行了必要的数值试验，考察了台站的方位覆盖以及台站数目对矩张量反演结果震源机制的影响，同时也考察了在10％随机噪声存在的情况下，台站的方位覆盖或台站的数目对结果的影响.

我们假设断层走向和倾角及滑动角都是45°、 震源深度为10 km的地震点源，计算震中距为80 km，方位角分别为0°，23°，45°，67°，90°，113°，135°，158°，180°，203°，203°，225°，247°，270°，293°，315°和337°的“观测点”的地震图，同时为每个观测点计算格林函数，然后在不同的台站方位覆盖、 不同的台站数目以及为观测地震图添加随机噪声的情况下实施矩张量反演. 表 2给出了该试验结果. 其中，试验1是无噪声情况下的结果，试验2是有噪声情况下的结果.

表 2 台站方位覆盖等因素对矩张量反演结果影响的评价试验结果

试验s010203，s050607，s091011和s131415是在3个台站相对于事件的方位角不同而张角都为45°情况下的结果; 试验s010407，s050811，s091215和s131603是在3个台站相对于事件的方位角不同而张角都为135°情况下的结果; 试验s010511，s050915，s091303和s130107是在3个台站相对于事件的方位角不同而张角都为225°情况下的结果； 试验s010814，s051202，s091607和s130412是在3个台站相对于事件的方位角不同而张角都为293°情况下的结果. 从试验结果可以看出，无论观测点相对于事件的方位如何，无论观测点相对于事件的张角是多大(45°，135°，225°或293°)，3个观测点的资料足以确定事件的震源机制(误差不超过3°)，即使在10％的随机噪声存在的情况下也是如此.

试验s01020304，s03040506，s05060708和s07080910是在4个台站相对于事件的方位角不同而张角都为67°情况下的结果； 试验s0102030405，s0304050607，s0506070809和s0708091011是在5个台站相对于事件的方位角不同而张角都为90°情况下的结果. 从试验结果可以看出，增加观测点可以更好地稳定反演结果； 而且，观测点越多，随机噪声对结果的影响也越小.

由于即使在3个观测点且张角为45°情况下，反演结果依然可靠稳定，所以在4个和5个观测点情况下，我们只测试了张角为67°和90°的情形.

5 讨论和结论

本研究的主要目的是通过反演首都圈地区宽频带数字地震波形资料，确定2004年发生在该地区的中小地震的矩张量解，由此确定这些地震的标量地震矩、 矩震级、 双力偶分量和补偿线性矢量偶极分量的大小以及断层面参数和应力轴参数.

为了避免首波的影响，我们只选择了初至震相为直达P波的记录资料. 同时，为了反演结果的稳定与可靠，我们只挑选了同时被3个以上台站记录清楚的地震事件. 在2004—2005年发生的众多地震事件中，只有51次地震满足上述条件.

通过矩张量反演，我们得到了这51次地震的矩张量解. 通过进一步的分解与计算，得到了这些地震事件最佳双力偶解、 应力轴参数、 标量地震矩、 补偿线性矢量偶极成分、 矩震级.

为了检验反演结果的可信度，针对观测台站相对于地震事件的方位分布、 观测台站相对于地震事件的方位覆盖、 以及随机噪声对反演结果的影响，我们设计了必要的数值试验. 结果表明: ① 张角大于45°的3个或3个以上的台站足以确定震源机制； ② 台站张角越大，结果误差越小； ③ 台站数目越多，结果误差越小； ④ 台站数目越多，随机噪声对结果的影响越小. 因此，我们选择同时被3个或3个以上的台站记录到的地震事件进行矩张量反演是正确的.

尽管同时被3个或3个以上的台站记录的地震事件的矩张量反演结果是稳定可靠的，但考虑到台站张角越大结果误差越小，台站数目越多结果误差越小，以及台站数目越多随机噪声对结果的影响越小的数值试验结果，我们根据台站相对于地震事件的方位角张角θ的大小，对反演结果进行了分类: A类： θ≥180°； B类： 90°≤θ<180°； C类： θ<90°. 在表 1中，我们对每次地震的类别也做了标注，供使用时参考.

本工作所用的资料系由北京数字遥测地震台网提供，特此致谢.

参考文献
[1] 许力生, 陈运泰. 2004. 格林函数库技术与快速地震矩张量反演[G]//陈运泰, 滕吉文, 阚荣举, 等主编. 中国大陆地震学与地球内部物理学研究进展: 庆贺曾融生院士八十寿辰. 北京: 地震出版社: 625-630.(1)
[2] 中国地震局震害防御司. 1992. 地震工作手册[M]. 北京: 地震出版社: 1-633.(1)
[3] 许忠淮. 2001. 东亚地区现今构造应力图的编制[J]. 地震学报, 23(5): 492-501.
[4] Frohlich C, Davis S D.1993. Teleseismic b values； or, much ado about 1.0[J]. J Geophys Res, 98: 631-644.(1)
[5] Gilbert F, Dziewonski A M. 1975. An application of normal mode theory to the retrieval of structureal parameters and source mechanisms from seismic data[J]. Phil Trans R Soc London, Ser A, 278: 187-269.(1)
[6] Kawakatsu H. 1991. Enigma of earthquakes at ridge-transform-fault plate boundaries-distribution of non-double couple parameter of Harvard CMT solutions[J]. Geophys Res Lett, 18: 1 103-1 106.(1)
[7] Kagan Y Y. 1997. Seismic moment-frequency relation for shallow earthquakes: Regional comparison[J]. J Geophys Res, 102: 2 835-2 852. Kagan Y Y. 2000. Temporal correlations of earthquake focal mechanisms[J]. Geophys J Int, 143: 881-897.(1)
[8] Kennett B L N. 1983. Seismic Wave Propagation in Stratified Media[M]. Cambridge: Cambridge University Press: 1-339.(1)
[9] Kubo A, Fukuyama E. 2003. Stress field along the Fyukyu Arc and the Okinawa Trough inferred from moment tensors of shallow earthquakes [J]. Earth Planet Sci Lett, 210: 305-316.(1)
[10] McKenzie D P. 1969. The relationship between fault plane solutions for earhtqkaes and directions of principal stresses[J]. Bull Seism Soc Amer, 59: 591-601.(1)
[11] Patton H, Aki K. 1979. Bias in the estimate of seismic moment tensor by the linear inversion method[J]. Geophys J R astr Soc, 59: 479-495.(1)
[12] Reasenberg P A. 1999. Foreshock occurrence before large earthquakes[J]. J Geophy Res, 104: 4 755-4 768. (1)
[13] Sornette D, Knopoff L, Kagan Y Y, et al. 1996. Rank-ordering statistics of extreme events: Application to the distribution of large earthquakes[J]. J Geophys Res, 101: 13 883-13 893.(1)
[14] Tanimoto T, Okamoto T. 2000. Change of crustal potential energy by earthquakes: An indicator for extensional and compressional tectonics[J]. Geophys Res Lett, 27: 2 313-2 316.(1)
[15] Tapponnier P, Molnar P. 1976. Slip-line field theory and large scale continental tectonics[J]. Nature, 264: 319-324.(1)
[16] Xu Z H, Wang S Y, Huang Y R, et al. 1992. Tectonic stress field of China inferred from a large number of small earthquakes[J]. J Geophy Res, 97(B8): 11 867-11 877.(1)
[17] Zoback M D, Zoback M L, Mount V S, et al. 1987. New evidence on the state of stress of the San Andreas fault system[J]. Science, 238: 1 105-1 111.(1)