用地震波波形拟合方法研究中小地震的震源机制

1、第 24 卷第 3 期地 震 地 磁 观 测 与 研 究SEISMOLOGICAL AND GEOMAGNETICVol 24 No 32003 年6 月OBSERVATION AND RESEARCHJun 2003用地震波波形拟合方法研究 中小地震的震源机制明跃红吴建平楼海（中国北京 100081 中国地震局地球物理研究所）摘要 采用网格搜索方法，利用区域地震三分量全波形资料，对中小地震震源机制的反演方法进 行了研究。反演采用初动约束、分段求相关系数，以及幅值加权的混合方法计算误差函数，通过 格点搜索法确定震源深度。理论模型试验表明，反演结果具有较好的稳定性和可靠性。云南地区 区域数字地震

2、台网的观测资料反演表明，该方法可以合理地拟合观测波形资料。将 2001 年 7 月 9日（°N，°E）发生的 MS 5.3 地震的震源机制解与哈佛大学 CMT 进行了对比，结果表 明，Z 分量波形拟合略比 CMT 拟合的好，R 分量和 T 分量的波形拟合明显好于 CMT 的对应分量。 我们对 2001 年 4 月 5 日在（°N，°E）发生的 mb 地震进行了震源机制反演。结果 表明，该地震主要表现为走滑断层，震源机制解为（180°，80°，11°）。根据这两个地震震源机 制推断的主应力轴，与该地区现代构造应力场特征基本一致

3、。关键词：震源机制；反演；中小地震中图分类号：+3文献标识码：A文章编号：1003-3246(2003)03-0020-07引言震源机制解是研究现代构造应力场的基础资料之一，它在全球应力场、活应力图等方 面得到了广泛的应用。中小地震震源机制解的一致性程度也是区分地震序列类型和预测后 续强震的重要指标。研究表明，大震前孕震区内小震震源的主压应力轴存在重新取向，震 区小震综合断层面解呈现出矛盾符号比下降，不同地震的波形趋于一致的现象。由于受观测台站分布的限制，对于较小的地震，远台记录的振幅较弱，基于长周期远 震波形拟合的震源机制反演，很难得到比较准确的中小地震的震源机制。20 世纪 90 年代以

4、前，由于我国区域地震台网主要为模拟纪录，在震源机制的研究中主要采用模拟记录资料 中的 P 波初动及振幅比信息（梁尚鸿，1986）。对于较小的地震，由于震级小、能量弱，通 常只有少量台站观测到清晰的记录，传统的仅依靠 P 波初动资料和振幅比的方法很难得到作者简介：眀跃红，女，1983 年毕业于北京大学地球物理系，1987 年获中国地震局地球物理研究所硕士学位。现主要从事 地球内部结构、火山地震学等方面的研究本课题为地震科学联合基金资助项目 中国地震局地球物理研究所论著编号 03AC2015 本文收到日期：2003-02-10第 3 期明跃红等：用地震波波形拟合方法研究中小地震的震源机制21可靠的

5、震源机制解，不能适应对大震震前应力场动态变化监测的需要。近年来，随着区域 数字地震台网的建立、非线性反演理论的发展以及在时间域中对震源时间函数的合理表示（Sileny，1992），使得利用少量宽频带近震资料反演震源机制的方法得到了快速发展（Zhao， 1994；Zhu，1996），为中小地震震源机制的研究提供了新的途径。本文采用宽频带数字地 震资料，对云南地区中小地震的震源机制研究进行了探讨。1理论与方法Aki 和 Richards（1980）以及 Jost 和 Herrmann（1989）对地震点源的矩张量进行了较为 详细的分析。矩张量可以在一级近似下完全描述一般点源的等效力。等效力能够与物

6、理震 源模型，如断层面上突然的相对位移、亚稳态相变的快速传播、相变引起的突然体积塌陷 或爆炸引起的体积膨胀。代表断层面上突然位移的等效力构成了人们熟知的双力偶。矩张 量是一个广义的概念，它能描述各种震源模型，双力偶源只是其中的一种。本文采用 Wang 和 Herrmann（1980）以及 Herrmann 和 Wang（1985）给出的频率-波数 积分法（FK 方法），计算层状介质中任意点位错源在地表产生的位移。位移场可表示为uz (r, z = 0,s ) = M xx ZSS (cos 2a ) / 2 - ZDD / 2 + ZEP / 3+ M yy -ZSS ( cos 2a )/

7、2 - ZSS / 2 + ZEP / 3+ M zz ZEP / 3 + M xy ZSS sin 2a+ M xz ZDS cosa + M yz ZDS sinaur (r, z = 0,s) = M xxRSS (cos 2a ) / 2 - ZDD / 2 + REP / 3+ M yy -RSS(cos 2a ) / 2 - ZDD / 2 + REP / 3+ M zz REP / 3 + M xy RSS sin 2a+ M xz RDS cosa+ M yz RDS sin aut (r, z = 0,s) = M xxTSS (sin 2a )/ 2 - M yy TSS

8、(sin 2a ) / 2- M xy TSS sin 2a + M xzTDS sina - M yzTDS sin a（1）式中 r 为震中距， a 为台站相对震中位置的方位角，函数 ZDD、ZDS、ZSS、RDD、RDS、RSS、 TDS、TSS，加上 ZEP 和 REP，是 10 个用于计算层状介质中任意点位错源或爆炸源产生的位移 场所需的格林函数。其中 ZSS 是走滑断层产生的 Z 分量位移，ZDS 是倾滑断层产生的 Z 分量 位移，ZDD 为 45°倾向滑动断层产生的 Z 分量位移，ZEP 表示由 3 个相互垂直的单位偶极子 共同产生的 Z 分量位移。R 和 T 分别表示

9、径向和切向分量的位移。对于位错源引起的地震， ZEP 和 REP 为零，矩张量分量与断层走向、倾向、滑动角之间的关系为M xx = -M 0 (sin 6 cosZ sin 2j + sin 26 sin Z sin 2 j) M yy = M 0 (sin 6 cos Z sin 2j - sin 26 sin Z cos2 j) M zz = M 0 (sin 26 sin Z )M xy = M 0 (sin 6 cosZ cos 2j + sin 26 sin Z sin 2j)M xz = -M 0 (cos6 cos Z cosj + cos 26 sin Z sin j) M y

10、z = -M 0 (cos6 cosZ sinj - cos 26 sin Z cosj) M yx = M xyM zx = M xz（2）式中 $ 、 6 、 Z 分别是断层的走向、倾角和滑动角。22地 震 地 磁 观 测 与 研 究24 卷我们利用地方震和区域地震资料，采用时间域内的波形拟合确定中小地震的震源机制。 对于震源时间函数相对简单的小地震，给定梯形震源时间函数，采用网格搜索法确定断层 的走向、倾角和滑动角和震源的深度。震源时间函数可以是梯形震源时间函数或任意形状 的震源时间函数。目标函数的选择对反演结果具有重要的影响。对它的选取可根据不同的问题和要求采 用多种不同的形式。在矩张

11、量反演中，最初人们通常采用的目标函数为ne(x) = åò2Toi (t) - gi (t, x) dt（3）i=1 0式中 gi (t, x) 为理论地震图，x 为与矩张量有关的参数。oi (t) 为第 i 个台站的观测记录。后来 人们逐渐认识到，采用这种形式的目标函数，由于震源在不同方向的地震波辐射强度不同，强辐射方向由于地震波振幅比较大，因而对误差的影响较大，在振幅较小的方向地震波波 形的影响变得不太重要。另外，由于场地条件等都对地震波振幅产生影响，使得采用这种 目标函数反演得到的地震矩张量，具有相对较大的误差。Helmberger 等（1980）认为，利 用理论地震

12、图与观测波形的相关函数作为误差函数，强调两者之间波形的拟合，对振幅的 绝对大小不太敏感，在矩张量反演中能取得较好的结果。e1(x)ìïM= åí1 -TüTò0 oi (t)gi (t, x)dtïT1 ý（4）i=1 ïîò0oi (t)oi (t)dtò0gi (t, x)gi (t, x)dt 2 ïþ姚振兴、郑天愉（1994）认为，由于 P 波在节线附近 P 波初动振幅往往很小，而它又 是确定矩张量的重要信息，提出除了利用上述误差函数外，还应引入初

13、动方向作为重要的 约束条件。我们认为，当资料的方位覆盖较好时，采用上述目标函数的构建是一种非常有效手段。 但对于某些观测资料的方位覆盖有限的地震，如果在某个角度范围内存在较多的观测资料， 我们有必要尽量利用振幅随方位变化的信息，考虑振幅的影响。为了减少几何扩散和震源 辐射的方向性对误差函数计算的不良影响，我们让不同方向和震中距的地震记录在误差计 算中具有相同的权重，即TNo (t)g (t, x)w 2 dte2 (x) =å ò0i- i i( o (t)imax)2i=1（5）wi = omax / g max式中 oi (t) 和 gi (t) 分别为第 i 个观测记

14、录和理论地震图， oi (t)max（6）为第 i 个观测记录在 0T时间段内绝对振幅的最大值。 omax 和 gmax 分别是所有观测记录和所有理论地震图中绝对振幅的最大值。该目标函数的特点是兼顾了远近不同距离的地震记录，使它们具有相同的权 重；同时以观测资料为依据考虑了振幅的影响。即当理论地震图与观测记录具有相同的波 形，且各理论地震图最大振幅之间的比例关系与观测记录间振幅的相互比例关系相同时，第 3 期明跃红等：用地震波波形拟合方法研究中小地震的震源机制23能最好地拟合观测值。在反演过程中，引入了初动方向约束，同时将全波形资料分为 P 波 段和 S 波与面波段两个部分，在相关系数的计算中

15、具有相同的权重。2 资料与结果云南数字地震台网是“九五”期间我国较早完成数字化改造的区域数字地震台网，台 网内共 23 个台站，台站分布见图 1。1999 年正式运行以来，该区域数字地震台网积累了大 量的中小地震的地震资料。图 1 云南地区数字地震台站分布及构造云南地区是多地震地区，该地区地壳厚度变化强烈，速度结构存在较大的横向变化。 为了减小速度结构变化对震源机制反演的影响，我们根据波形反演方法确定了该地区的平 均速度结构（表 1），以此为模型计算格林函数。为了验证本方法的可靠性，我们选取 2001 年 7 月 9 日 23 时 51 分 秒在（°N，°E）发生的 MS

16、地震，将所得到的震源机制解与哈佛大学 CMT 进行了对比。表 2 中给出了哈佛大学测定的震源参数和震源机制解以及中国地震台网给出的震源参数和我们 反演得到的震源机制解。反演采用了中国地震台网的震源参数。由于震源深度的可靠性较差，反演过程中通过对不同深度进行搜索，将波形拟合误差 最小的深度作为地震的震源深度，相应的震源机制解作为最佳拟合结果，并对 z 分量的初动 加以约束。该地震计算结果见表 3。表 1 平均速度结构模型层厚度/kmP 波速度/km· s-1S 波速度/km· s-1QQ aþ 24地震 地 磁 观 测 与研 究24 卷表 2 震源参数及震源机制日期

17、时间jN /°ZE /°深度MS断层面 1断层面 22001-07-0923：51：233，50，2142，89，1402001-07-0923：51：315，76，17047，80，14表 3 不同深度的震源机制解深度/km断层面 1断层面 2相关系数130，90，18040，90，00.289 523 440，90，10310，80，1800.442 938 6310，80，17040，80，100.522 183 2310，80，17040，80，100.540 450 7310，80，17040，80，100.555 948 2310，80，17040，80，100

18、.547 112 240，80，-10130，80，-1700.441 508 86从表 3 中可以看出，理论波形与观测波形之间相关系数最大的震源深度为 11 km，我们 认为该深度为地震震源位置的最佳深度，相应的震源机制解为：走向 310°，倾角 80°，错 动角 170°。表 3 中给出的结果是以角度 10°为间隔搜索得到的。为了提高震源机制解的精 度，在此基础上，将角度搜索间隔由原来的 10°提高到 1°，得到的震源机制最终解为走向 315°，倾角 76°，错动角 170°，相关系数为 ，地震表现为

19、以走滑为主。图 2 给出了 部分台站（8 个台站）的观测波形和理论波形之间的拟合结果，图中观测值为实线，理论值 为虚线，从左到右依次为 Z 方向、R 方向、T 方向。从图 2 中可以看出各台的初动及波形都 拟合得较好，由该地震的震源机制解计算的应力场为：最大压应力轴的方位 180°，仰角 3°， 最小压应力轴的方位 271°，仰角 17°，中等压应力轴的方位 81°，仰角 73°，近似垂直。图 2 2001 年 7 月 9 日的地震。（a）Z 方向；（b）R 方向；（c）T 方向为了进一步验证方法的可靠性，我们将该震源机制解与哈佛大学

20、 CMT 进行了对比。图3 给出了根据哈佛大学震源机制（233，50，2）计算的理论波形和观测波形之间的对比结果，相应的相关系数为 ，比我们在反演过程中得到的相关系数 低得多。图 3 给出了第 3 期明跃红等：用地震波波形拟合方法研究中小地震的震源机制25与上述 8 个台站相对应的观测波形和理论波形之间的拟合结果。对比图 2 与图 3 可知：Z 分量：图 2 总体好于图 3，图 2 的黑龙潭、云县、通海、云龙台的波形拟合明显好于图 3 所对 应台的波形拟合；R 分量：也是图 2 好于图 3，图 2 的黑龙潭、通海、云龙、昭通台的波形 拟合明显好于图 3 所对应台的波形拟合；T 分量：图 2 的

21、云龙台波形拟合的好，图 3 的畹町台波形拟合的好。从三分量总体看：图 2 的拟合结果好于图 3 的拟合结果。图 3 2001 年 7 月 9 日的地震 CMT。（a）Z 方向；（b）R 方向；（c）T 方向图 4 2001 年 4 月 5 日的地震。（a）Z 方向；（b）R 方向；（c）T 方向（图中实线为观测值，虚线为理论值，横坐标为时间，单位为 s）文中，利用该方法对 2001 年 4 月 5 日 18 时 12 分 秒在（°N，°E）发 生的 mb 地震，进行了反演，得到了这个地震的震源机制解，为走向 180°，倾角 80°， 错动角 11

22、6;，相关系数为 ，该地震主要表现为以走滑为主的断层。另一个断层面是（88， 79，170）。图 4 给出了该地震观测值和理论值的比较，可以看出各台各分量的波形拟合的 都很好，各台的初动也都复合的很好。用该地震的震源机制解计算该地区的应力场，可得26地 震 地 磁 观 测 与 研 究24 卷最大压应力轴的方位 314°，仰角 1°，最小压应力轴的方位 44°，仰角 15°，中等压应力轴 的方位 222°，仰角 75°。断层面近直立时，断层一般沿水平方向滑动，这与我们得到的两个地震的震源机制主 要表现为高角度走滑断层相一致。根据两个地震

23、的震源机制得到，中等主应力轴的方位角 和仰角分别为（81°，73°）和（222°，75°），方向基本垂直。全球地应力测量结果表明， 对于走滑断层来说，垂直向应力通常是中等主应力，这与我们所得的应力方向一致。我们得出的这两个地震的应力场方向都表现为最大应力轴和最小应力轴呈水平方向、 中等应力轴为垂直方向。这个结果也符合西南地区现代构造应力场的特征。参考文献梁尚鸿等. 1986. 攀西地区地震分布和构造应力场特征 J. 地球物理学报，29(6)：557566姚振兴等. 1994. 用 P 波波形资料测定中强地震震源过程的地震距张量反演方法 J. 地球物理学报

24、，37(1)：3644 Aki K and P G Richards. 1980. Quantitative seismology: theory and methods M. W H Freeman and Company, San Francisco, CalifHelmberger D V and G R Engen. 1980. Modeling long-period body waves from shallow earthquakes at regional distances J. Bull Seism Soc Am, 70, 1 6991 714Herrmann R B an

25、d C Y Wang. 1985. A comparison of synthetic seismograms J. Bull Seism Soc Am 75, 4156Jost M L and Herrmann R B. 1989. A student s guide and review of moment tensors J. Seism Res Lett, 60, 3557Sileny J, F Panza and P Campus. 1992. Waveform inversion for pinot source moment tensor retrieval with varia

26、ble hypocentral depth and structural model J. Geophys J Int, 109, 259274Wang C Y and R B Herrmann. 1980. A numerical study of P-, SV-, and SH-wave generation in a plane layered medium J. Bull Seism Soc Am, 70, 1 0151 036Zhao and Helmberger. 1994 Source estimation from broadband regional seismograms

27、J. Bull Seism Soc Am, 84, 91104Zhu and Helmberger, 1996. Advancement in source estimation techniques using broadband regional seismograms J. Bull Seism Soc Am, 86, 1 6341 641The source mechanism inversion for middle and small earthquakes by using seismic waveform fitting techniqueMing Yuehong，Wu Jia

28、nping and Lou Hai（Institute of Geophysics, China Seismological Bureau, Bejing 100081, China）AbstractA source mechanism inversion technique for middle and small earthquakes was studied by using complete regional waveforms of all three components. The grid search algorithm is used to find the best solution of mechanism and focal depth. In inversion process the misfit function is calculated by using first motion constrain and coherent coefficients of waveforms for diffe