基于面波数据的地下核爆炸的全元素矩张量反演方法*

何永锋,李 锴,刘炳灿,姚国政,赵克常,张献兵,曾乐贵

(1.陆军装甲兵学院，北京 100072; 2.北京大学地球与空间科学学院，北京 100871)

基于面波数据的地下核爆炸的全元素矩张量反演方法*

何永锋1,李 锴1,刘炳灿1,姚国政1,赵克常2,张献兵2,曾乐贵1

(1.陆军装甲兵学院，北京 100072; 2.北京大学地球与空间科学学院，北京 100871)

区域分层介质模型下, 可以将地震波场描述为矩张量各分量作为权重的基本格林函数的线性组合，利用该理论地震波场可以反演实际天然地震或地下核爆炸的震源机制，反演结果中不同震源机制成分的比重，可以用来识别地下核爆炸，该系统方法越来越受到关注。给出了基于广义反射-透射系数方法的水平分层介质模型的地震波场正演公式，并对基于该公式的单台反演结果的准确性、稳定性、可靠性进行了理论分析，为利用该公式对实际地下核爆炸进行反演提供了理论基础，该方法对利用区域少量甚至是单站记录数据检测、识别地下核爆炸具有重要参考意义。

地下核爆炸;面波;格林函数;全元素矩张量反演方法;补偿线性偶极矢量源

通常采用忽略源区非线性效应的等效力模型来描述震源机制[1-2]，震源机制的定量描述主要包括震源强度和震源的断层面解(倾角、滑动角、走向)。F.Gilbert首先引进了矩张量的概念[3], 定义为作用在一点上的等效体力的一阶矩,这种等效力不仅可以描述典型天然地震的位错模型，还可以描述由爆炸导致的源区体积的快速膨胀及由于介质相变导致的源区体积的快速坍塌[4]。

介质相变的力学性质表现为剪切模量的突然改变，这个物理机制可以用补偿线性偶极矢量源(compensate linear vector dipole, CLVD)表示[5-7]。用矩张量表示震源，能将记录波形数据、震源机制和传播路径三者之间的关系归结为线性关系，如果已知震源位置和相应的介质模型下的格林函数，那么由记录波形数据就可以线性地反演出震源机制矩张量。

利用初动波形数据反演震源机制矩张量需要分布较好的近场台站数据，当这个条件无法得到满足时，从区域少量台站甚至是单台记录数据反演稳定的矩张量解就显得十分有意义[8-9]，由于区域的长周期地震波对速度结构的横向变化相对不敏感，在震源定位较准确、波形数据信噪比较高的情况下，利用面波数据也可以反演得到稳定的矩张量解。

地下核爆炸通常是在人烟稀少且严格保密的条件下进行的，获取实施地下核爆炸的本土近场台站数据很困难，因此利用区域少量甚至是单站记录数据来检测、识别地下核实验具有重要现实意义。

1 方法原理

近年来，理论地震图计算能力及效率得到了进一步的提高，基于复杂的分层介质模型的反演方法越来越受到人们的重视。

天然地震的反演研究以M.L.Jost等[10]的工作为代表，C.Y.Wang等[8]在前人工作基础上，给出了任意一个位错点源和爆炸源产生的地震波场表达式及所需的10个基本Green函数，其中包含爆炸源Green函数，理论上可以得到含有爆炸源成分的全元素矩张量解(full moment tensor)。

S.E.Minson等[11]对M.L.Jost等[10]给出的含有爆炸源格林函数的地震波场公式进行了修正，并采用C.K.Saikia[12]给出的离散波数积分方法来计算理论格林函数。本文中，采用X.F.Chen[13]、Z.X.Yao等[14]的基于广义反、透射系数方法的水平层状地球模型中理论地震图的计算方法，本质上两种方法是一致的，矩张量权重表现形式略有不同，经过对两种算法对比分析，得到本文中采用的公式。在圆柱坐标系下，理论位移u在垂向、径向和切向分量分别为：

(1)

式中：θ是台站到震源的方位角，GSS是纯走滑断层(倾角90°，滑动角0°)的格林函数，GDS是纯倾滑断层(倾角90°，滑动角90°)的格林函数，GDD是45°的斜滑断层(倾角45°，滑动角90°)的格林函数[15]，GEP是纯爆炸源的格林函数，z、r、t分别表示垂向、径向和切向。与S.E.Minson等[11]和D.Dreger等[9]给出的公式不同之处，主要体现在GDD,z、GDD,r权重系数上。由于考虑了纯爆炸源的格林函数，由式(1)可以反演出含有爆炸成分的震源机制矩张量，对所反演矩张量没有任何约束条件，不仅能反演力偶(double couple, DC)成分，还能准确地反演对角线(isotropic, ISO)成分。

2 数值模拟

2.1理论震源机制矩张量反演

B.Romanowicz等[16]的研究结果表明，在速度结构比较准确的情况下，利用单台三分量的区域震相进行反演能得到可靠的反演结果[17]。为验证基于式(1)反演方法的准确性，利用该公式对不同震源机制进行理论数值反演, 采用的模型为K.L.Mclaughlin等[18]给出的适于东哈萨克斯坦地区的地壳速度模型(见表1)。

QS、QP分别为横波、纵波品质因数，震源深度统一取d=10 km，震源时间函数为δ函数，理论格林函数的计算结果如图1所示，利用Butterworth带通滤波器进行滤波，滤波周期为20～50 s。

表1 理论地壳模型Table 1 Theoretical crustal model

2.1.1纯爆炸源的反演

先计算基于球对称源的理论地震图，并进行周期范围为20～50 s滤波处理， 将各格林函数代入式(1)，利用最小二乘法在时间域进行反演[19]。得到EXP源准确矩张量解：

2.1.2CLVD源的反演

先计算基于CLVD源的理论地震图，并进行滤波，将各格林函数代入式(1)进行反演。得到CLVD源的准确矩张量解：

2.1.3任意位错源的反演

媳妇一跟我吵架，就哭着跑出去逛街购物，以发泄心中的不满。今天媳妇哭着对我说：“这日子没法过了，你已经一个星期没跟我吵架了。”

采用S.E.Minson等[11]给出的DC源为任意位错源(倾角67°，滑动角45°，走向23°)，先计算该DC源的理论地震图并进行相应的滤波处理，利用式(1)进行反演。得到该DC源的准确矩张量解：

2.1.4混合源的反演

混合源为含有多种成分的震源X(EXP+CLVD+DC)，先计算混合源的理论地震图，并进行相应滤波处理，利用式(1)进行反演。得到该混合源的准确矩张量解：

在介质速度结构和震源深度等信息准确的条件下，基于式(1)进行反演，可以准确得到震源中EXP、CLVD和DC成分，根据各种成分的比重，理论上可以达到区分天然地震和地下核爆炸的目的。

2.2噪声对反演结果的影响

为验证反演方法在噪声干扰下的有效性，采用理论地壳模型(见表1)计算X源(EXP+CLVD+DC)的理论地震图，并分别叠加10%、30%、50%的噪声干扰。相应的反演结果分别为：

2.3震相到时误差对反演结果的影响

实际波形的到时可能与一维地壳模型的理论到时存在一定的偏差[20]。为评估震相到时误差对反演结果的影响，本文中对台站设置±3 s的随机震相到时误差，对X源分别计算在10%、30%和50%的噪声干扰和±3 s的随机震相到时误差的情况下的理论地震图，并进行矩张量反演。数值计算结果表明，存在震相到时误差的情况下，随着噪声干扰水平的增大，反演结果开始出现波动，但是整体结果均在较准确的范围，说明在震相到时误差的情况下，反演方法还是比较很稳定的，与郑建常等[20]得到的结论是一致的。

2.4速度模型对反演结果的影响

为了考察介质速度模型对反演结果的影响，将表1的速度模型进行修改，对其相邻层进行合并计算，对密度、波速和Q分别求平均，合并后的模型称为平均模型。采用平均速度模型计算各格林函数，然后反演基于表1计算的爆炸源、CLVD源、DC源和X源的理论地震图，得到4种震源机制矩张量：

计算结果表明，采用平均模型的反演结果，与基于表1的速度模型的反演结果非常一致。这说明，可以用平均模型来代替较复杂的模型，与郑建常等[20]的结论一致。

2.5震源深度误差对反演的影响

由于不同深度的震源对理论格林函数计算结果影响较大，因此会对最后的矩张量反演结果产生一定范围的误差。本文中利用不同深度震源的理论格林函数反演得到源矩张量，计算其理论地震波形，并与实际观测数据进行比较，计算方差缩减RV，取方差缩减最大结果为最佳解[21-22]。采用此方法，分别求得不同深度h(8 km≤h≤12 km )下EXP源、CLVD源、DC源和X源的RV，其中模拟观测位移的爆炸源、CLVD源、DC源和X源的深度均为10 km，略去数值计算过程，给出最终结果，见表2。

计算结果表明，当理论格林函数的震源深度越接近真实深度时，RV越大，反演结果越接近真实情况。另外，当理论震源深度大于实际震源深度时，震源深度误差对反演的影响较小，如在理论震源深度为10.5 km的效果比深度为9.5 km的反演效果好，与许力生等[23]的结论一致。

表2 不同深度下的方差缩减Table 2 Variance reduction at different source depths

3 结 论

不同的震源机制具有不同的矩张量形式，利用反演得到的矩张量的特征值及特征向量可以对震源机制进行分析，震源矩张量可以分解为对角线部分及偏量部分，地下核爆炸的震源与天然地震震源的矩张量中的对角部分和偏量部分具有不同的表现形式。本文中利用数值方法，验证了基于单台数据的震源矩张量反演方法，分析了存在噪声干扰、震相到时误差、速度模型误差和震源深度定位误差的情况下，反演方法的稳定性。各种误差的分析结果表明，震源深度误差对反演结果的影响较大。地震波的实际传播路径和地壳的三维结构对记录波形影响较大，同时台站的记录数据可能会受到其他的干扰因素，数据品质会受到影响，因此由一维地壳速度模型反演得到准确的结果是困难的，即使利用对分层结构敏感的面波。本文的研究结果对于从震源的角度来了解地下核爆炸的物理机制具有较好的参考意义，同时也为利用单台记录数据反演震源矩张量、并进一步进行识别提供了理论支持。

[1] Stump B W, Johnson L R. The determination of source properties by the linear inversion of seismograms[J]. Bulletin of the Seismological Society of America, 1977,67(6):1489-1502.

[2] Aki K, Richards P G. Quantitative seismology: Theory and methods[M]. San Francisco: Freeman W H and Company, 1980.

[3] Gilbert F. Excitation of the normal modes of the earth by earthquake sources[J]. Geophysical Journal International, 1971,22(2):223-226.

[4] Kennett B L N. Seismic wave propagation in stratified media[M]. Cambridge: Cambridge University Press, 1983.

[5] Knopoff L, Randall M J. The compensated linear-vector dipole: A possible mechanism for deep earthquakes[J]. Journal of Geophysical Research, 1970,75(26):4957-4963.

[6] 何永锋,陈晓非,张海明.地下核爆炸Lg波的激发机制[J].地球物理学报,2005,48(2):367-372.

He Yongfeng, Chen Xiaofei, Zhang Haiming. The excitation of Lg wave by underground nuclear explosion[J]. Chinese Journal of Geophysics, 2005,48(2):367-372.

[7] 何永锋,赵克常,张献兵,等.地下核爆炸地震波二次源特征[J].地球物理学,2012,55(5):1742-1748.

He Yongfeng, Zhao Kechang, Zhang Xianbing, et al. The characteristic of the waveform from the second source induced by underground explosion[J]. Chinese Journal of Geophysics, 2012,55(5):1742-1748.

[8] Wang C Y, Herrmann R B. A numerical study of P-, SV-, and Sh-wave generation in a plane layered medium[J]. Bulletin of the Seismological Society of America, 1980,70(4):1015-1036.

[9] Dreger D, Helmberger D. Determination of source parameters at regional distances with three-components sparse network data[J]. Journal of Geophysical Research, 1993,98(B5):8107-8125.

[10] Jost M L, Herrmann R B. A student's guide to and review of moment tensor[J]. Seismological Research Letters, 1989,60(2):37-57.

[11] Minson S E, Dreger D S. Stable inversions for complete moment tensors[J]. Geophysical Journal International, 2010,174(2):585-592.

[12] Saikia C K. Modified frequency-wavenumber algorithm for regional seismograms using Filon’s quadrature: Modeling of Lg waves in eastern North America[J]. Geophysical Journal International, 1994,118(1):142-158.

[13] Chen X F. A systematic and efficient method of computing normal modes for multilayered half-space[J]. Geophysical Journal International, 1993,115(2):391-409.

[14] Yao Z X, Harkrider D G. A generalized reflection-transmission coefficient matrix and discrete wavenumber method for synthetic seismograms[J]. Bulletin of the Seismological Society of America, 1983,73(6):1685-1699.

[15] Langston C A. Source inversion of seismic waveforms: The Koyna, India, earthquakes of 13 September 1967[J]. Bulletin of the Seismological Society of America, 1981,71(1): 1-24.

[16] Romanowicz B, Dreger D, Pasyanos M, et al. Moniting of strain release in central and northern California using broadband data[J]. Geophysical Research Letters, 1993,20(15):1643-1646.

[17] 赵翠萍.1997-2003年新疆伽师震源区特征的地震学方法研究[D].北京:中国地震局地球物理研究所,2006.

[18] McLaughlin K L, Barker T G, Day S M, et al. Effects of depth of burial on explosion and earthquake regional seismograms: Regional discrimination and yield estimation[R]. Jolla, California, 1988.

[19] Fukuyama E, Dreger D S. Performance test of an automated moment tensor determination system for the future “Tokai” earthquake[J]. Earth, Planets and Space, 2000,52(6):383-392.

[20] 郑建常,陈运泰.基于Langston分解和Hilbert变换约束的区域偏量矩张量反演方法及应用[J].地震学报,2012,34(2):171-190.

Zheng Jianchang, Chen Yuntai. Regional deviatoric moment tensor inversion based on Langston's decomposition and Hilbert transform constraints and its application[J]. Acta Seismologica Sinica, 2012,34(2):171-190.

[21] 林向东，葛洪魁，徐平, 等.近场全波形反演：芦山7.0级地震及余震矩张量解[J].地球物理学报,2013,56(12):4037-4047.

[22] Templeton D C, Dreger D S. Non-double-couple earthquakes in the long valley volcanic region[J]. Bulletin of the Seismological Society of America, 2006,96(1):69-79.

[23] 许力生，陈运泰.震源深度误差对矩张量反演的影响[J].地震学报,1997,19(5):462-470.

Abstract: Powerful techniques have been developed for calculating the plane wave response of horizontally layered models. This method is quite general and is widely used in synthetic wave algorithms. Using this method, we can describe the displacements in terms of a linear combination of the moment tensor elements, and the moment tensor for an arbitrarily oriented dislocation can be given by this method. The moment tensor can be used to distinguish natural earthquakes and underground nuclear experiments according to its different elements. In this paper we rewrite the formula and estimate the reliability of the non-double-couple solutions on the basis of error analysis that includes the variance of modeling and of the noise in the data. Our analysis of synthetic data shows that this method is robust and can be used in the real data analyses. The result is significant for monitoring nuclear explosions by using data from just a few monitoring stations or even from a single station.

Keywords: underground nuclear experiment; surface wave; Green functions; full moment tensor inversion; compensate linear vector dipole source

(责任编辑 丁 峰)

Fullmomenttensorinversionmethodofundergroundnuclearexplosionsbasedonsurfacewavesdata

He Yongfeng1, Li Kai1, Liu Bingcan1, Yao Guozheng1, Zhao Kechang2, Zhang Xianbing2, Zeng Legui1

(1.ArmyArmoredForcesAcademy,Beijing100072,China; 2.SchoolofEarthandSpaceSciences,PekingUniversity,Beijing100871,China)

O382.1国标学科代码1303520

A

10.11883/1001-1455(2017)05-0945-06

2016-01-27；

2016-08-29

国家自然科学基金项目(41374068)

何永锋(1966— )，男，博士研究生，教授,heyfeng@sina.com。