基于面波數(shù)據(jù)的地下核爆炸的全元素矩張量反演方法*

何永鋒,李 鍇,劉炳燦,姚國(guó)政,趙克常,張獻(xiàn)兵,曾樂貴

(1.陸軍裝甲兵學(xué)院，北京 100072; 2.北京大學(xué)地球與空間科學(xué)學(xué)院，北京 100871)

基于面波數(shù)據(jù)的地下核爆炸的全元素矩張量反演方法*

何永鋒1,李 鍇1,劉炳燦1,姚國(guó)政1,趙克常2,張獻(xiàn)兵2,曾樂貴1

(1.陸軍裝甲兵學(xué)院，北京 100072; 2.北京大學(xué)地球與空間科學(xué)學(xué)院，北京 100871)

區(qū)域分層介質(zhì)模型下, 可以將地震波場(chǎng)描述為矩張量各分量作為權(quán)重的基本格林函數(shù)的線性組合，利用該理論地震波場(chǎng)可以反演實(shí)際天然地震或地下核爆炸的震源機(jī)制，反演結(jié)果中不同震源機(jī)制成分的比重，可以用來(lái)識(shí)別地下核爆炸，該系統(tǒng)方法越來(lái)越受到關(guān)注。給出了基于廣義反射-透射系數(shù)方法的水平分層介質(zhì)模型的地震波場(chǎng)正演公式，并對(duì)基于該公式的單臺(tái)反演結(jié)果的準(zhǔn)確性、穩(wěn)定性、可靠性進(jìn)行了理論分析，為利用該公式對(duì)實(shí)際地下核爆炸進(jìn)行反演提供了理論基礎(chǔ)，該方法對(duì)利用區(qū)域少量甚至是單站記錄數(shù)據(jù)檢測(cè)、識(shí)別地下核爆炸具有重要參考意義。

地下核爆炸;面波;格林函數(shù);全元素矩張量反演方法;補(bǔ)償線性偶極矢量源

通常采用忽略源區(qū)非線性效應(yīng)的等效力模型來(lái)描述震源機(jī)制[1-2]，震源機(jī)制的定量描述主要包括震源強(qiáng)度和震源的斷層面解(傾角、滑動(dòng)角、走向)。F.Gilbert首先引進(jìn)了矩張量的概念[3], 定義為作用在一點(diǎn)上的等效體力的一階矩,這種等效力不僅可以描述典型天然地震的位錯(cuò)模型，還可以描述由爆炸導(dǎo)致的源區(qū)體積的快速膨脹及由于介質(zhì)相變導(dǎo)致的源區(qū)體積的快速坍塌[4]。

介質(zhì)相變的力學(xué)性質(zhì)表現(xiàn)為剪切模量的突然改變，這個(gè)物理機(jī)制可以用補(bǔ)償線性偶極矢量源(compensate linear vector dipole, CLVD)表示[5-7]。用矩張量表示震源，能將記錄波形數(shù)據(jù)、震源機(jī)制和傳播路徑三者之間的關(guān)系歸結(jié)為線性關(guān)系，如果已知震源位置和相應(yīng)的介質(zhì)模型下的格林函數(shù)，那么由記錄波形數(shù)據(jù)就可以線性地反演出震源機(jī)制矩張量。

利用初動(dòng)波形數(shù)據(jù)反演震源機(jī)制矩張量需要分布較好的近場(chǎng)臺(tái)站數(shù)據(jù)，當(dāng)這個(gè)條件無(wú)法得到滿足時(shí)，從區(qū)域少量臺(tái)站甚至是單臺(tái)記錄數(shù)據(jù)反演穩(wěn)定的矩張量解就顯得十分有意義[8-9]，由于區(qū)域的長(zhǎng)周期地震波對(duì)速度結(jié)構(gòu)的橫向變化相對(duì)不敏感，在震源定位較準(zhǔn)確、波形數(shù)據(jù)信噪比較高的情況下，利用面波數(shù)據(jù)也可以反演得到穩(wěn)定的矩張量解。

地下核爆炸通常是在人煙稀少且嚴(yán)格保密的條件下進(jìn)行的，獲取實(shí)施地下核爆炸的本土近場(chǎng)臺(tái)站數(shù)據(jù)很困難，因此利用區(qū)域少量甚至是單站記錄數(shù)據(jù)來(lái)檢測(cè)、識(shí)別地下核實(shí)驗(yàn)具有重要現(xiàn)實(shí)意義。

1 方法原理

近年來(lái)，理論地震圖計(jì)算能力及效率得到了進(jìn)一步的提高，基于復(fù)雜的分層介質(zhì)模型的反演方法越來(lái)越受到人們的重視。

天然地震的反演研究以M.L.Jost等[10]的工作為代表，C.Y.Wang等[8]在前人工作基礎(chǔ)上，給出了任意一個(gè)位錯(cuò)點(diǎn)源和爆炸源產(chǎn)生的地震波場(chǎng)表達(dá)式及所需的10個(gè)基本Green函數(shù)，其中包含爆炸源Green函數(shù)，理論上可以得到含有爆炸源成分的全元素矩張量解(full moment tensor)。

S.E.Minson等[11]對(duì)M.L.Jost等[10]給出的含有爆炸源格林函數(shù)的地震波場(chǎng)公式進(jìn)行了修正，并采用C.K.Saikia[12]給出的離散波數(shù)積分方法來(lái)計(jì)算理論格林函數(shù)。本文中，采用X.F.Chen[13]、Z.X.Yao等[14]的基于廣義反、透射系數(shù)方法的水平層狀地球模型中理論地震圖的計(jì)算方法，本質(zhì)上兩種方法是一致的，矩張量權(quán)重表現(xiàn)形式略有不同，經(jīng)過對(duì)兩種算法對(duì)比分析，得到本文中采用的公式。在圓柱坐標(biāo)系下，理論位移u在垂向、徑向和切向分量分別為：

(1)

式中：θ是臺(tái)站到震源的方位角，GSS是純走滑斷層(傾角90°，滑動(dòng)角0°)的格林函數(shù)，GDS是純傾滑斷層(傾角90°，滑動(dòng)角90°)的格林函數(shù)，GDD是45°的斜滑斷層(傾角45°，滑動(dòng)角90°)的格林函數(shù)[15]，GEP是純爆炸源的格林函數(shù)，z、r、t分別表示垂向、徑向和切向。與S.E.Minson等[11]和D.Dreger等[9]給出的公式不同之處，主要體現(xiàn)在GDD,z、GDD,r權(quán)重系數(shù)上。由于考慮了純爆炸源的格林函數(shù)，由式(1)可以反演出含有爆炸成分的震源機(jī)制矩張量，對(duì)所反演矩張量沒有任何約束條件，不僅能反演力偶(double couple, DC)成分，還能準(zhǔn)確地反演對(duì)角線(isotropic, ISO)成分。

2 數(shù)值模擬

2.1理論震源機(jī)制矩張量反演

B.Romanowicz等[16]的研究結(jié)果表明，在速度結(jié)構(gòu)比較準(zhǔn)確的情況下，利用單臺(tái)三分量的區(qū)域震相進(jìn)行反演能得到可靠的反演結(jié)果[17]。為驗(yàn)證基于式(1)反演方法的準(zhǔn)確性，利用該公式對(duì)不同震源機(jī)制進(jìn)行理論數(shù)值反演, 采用的模型為K.L.Mclaughlin等[18]給出的適于東哈薩克斯坦地區(qū)的地殼速度模型(見表1)。

QS、QP分別為橫波、縱波品質(zhì)因數(shù)，震源深度統(tǒng)一取d=10 km，震源時(shí)間函數(shù)為δ函數(shù)，理論格林函數(shù)的計(jì)算結(jié)果如圖1所示，利用Butterworth帶通濾波器進(jìn)行濾波，濾波周期為20～50 s。

表1 理論地殼模型Table 1 Theoretical crustal model

2.1.1純爆炸源的反演

先計(jì)算基于球?qū)ΨQ源的理論地震圖，并進(jìn)行周期范圍為20～50 s濾波處理， 將各格林函數(shù)代入式(1)，利用最小二乘法在時(shí)間域進(jìn)行反演[19]。得到EXP源準(zhǔn)確矩張量解：

2.1.2CLVD源的反演

先計(jì)算基于CLVD源的理論地震圖，并進(jìn)行濾波，將各格林函數(shù)代入式(1)進(jìn)行反演。得到CLVD源的準(zhǔn)確矩張量解：

2.1.3任意位錯(cuò)源的反演

媳婦一跟我吵架，就哭著跑出去逛街購(gòu)物，以發(fā)泄心中的不滿。今天媳婦哭著對(duì)我說：“這日子沒法過了，你已經(jīng)一個(gè)星期沒跟我吵架了?！?/p>

采用S.E.Minson等[11]給出的DC源為任意位錯(cuò)源(傾角67°，滑動(dòng)角45°，走向23°)，先計(jì)算該DC源的理論地震圖并進(jìn)行相應(yīng)的濾波處理，利用式(1)進(jìn)行反演。得到該DC源的準(zhǔn)確矩張量解：

2.1.4混合源的反演

混合源為含有多種成分的震源X(EXP+CLVD+DC)，先計(jì)算混合源的理論地震圖，并進(jìn)行相應(yīng)濾波處理，利用式(1)進(jìn)行反演。得到該混合源的準(zhǔn)確矩張量解：

在介質(zhì)速度結(jié)構(gòu)和震源深度等信息準(zhǔn)確的條件下，基于式(1)進(jìn)行反演，可以準(zhǔn)確得到震源中EXP、CLVD和DC成分，根據(jù)各種成分的比重，理論上可以達(dá)到區(qū)分天然地震和地下核爆炸的目的。

2.2噪聲對(duì)反演結(jié)果的影響

為驗(yàn)證反演方法在噪聲干擾下的有效性，采用理論地殼模型(見表1)計(jì)算X源(EXP+CLVD+DC)的理論地震圖，并分別疊加10%、30%、50%的噪聲干擾。相應(yīng)的反演結(jié)果分別為：

2.3震相到時(shí)誤差對(duì)反演結(jié)果的影響

實(shí)際波形的到時(shí)可能與一維地殼模型的理論到時(shí)存在一定的偏差[20]。為評(píng)估震相到時(shí)誤差對(duì)反演結(jié)果的影響，本文中對(duì)臺(tái)站設(shè)置±3 s的隨機(jī)震相到時(shí)誤差，對(duì)X源分別計(jì)算在10%、30%和50%的噪聲干擾和±3 s的隨機(jī)震相到時(shí)誤差的情況下的理論地震圖，并進(jìn)行矩張量反演。數(shù)值計(jì)算結(jié)果表明，存在震相到時(shí)誤差的情況下，隨著噪聲干擾水平的增大，反演結(jié)果開始出現(xiàn)波動(dòng)，但是整體結(jié)果均在較準(zhǔn)確的范圍，說明在震相到時(shí)誤差的情況下，反演方法還是比較很穩(wěn)定的，與鄭建常等[20]得到的結(jié)論是一致的。

2.4速度模型對(duì)反演結(jié)果的影響

為了考察介質(zhì)速度模型對(duì)反演結(jié)果的影響，將表1的速度模型進(jìn)行修改，對(duì)其相鄰層進(jìn)行合并計(jì)算，對(duì)密度、波速和Q分別求平均，合并后的模型稱為平均模型。采用平均速度模型計(jì)算各格林函數(shù)，然后反演基于表1計(jì)算的爆炸源、CLVD源、DC源和X源的理論地震圖，得到4種震源機(jī)制矩張量：

計(jì)算結(jié)果表明，采用平均模型的反演結(jié)果，與基于表1的速度模型的反演結(jié)果非常一致。這說明，可以用平均模型來(lái)代替較復(fù)雜的模型，與鄭建常等[20]的結(jié)論一致。

2.5震源深度誤差對(duì)反演的影響

由于不同深度的震源對(duì)理論格林函數(shù)計(jì)算結(jié)果影響較大，因此會(huì)對(duì)最后的矩張量反演結(jié)果產(chǎn)生一定范圍的誤差。本文中利用不同深度震源的理論格林函數(shù)反演得到源矩張量，計(jì)算其理論地震波形，并與實(shí)際觀測(cè)數(shù)據(jù)進(jìn)行比較，計(jì)算方差縮減RV，取方差縮減最大結(jié)果為最佳解[21-22]。采用此方法，分別求得不同深度h(8 km≤h≤12 km )下EXP源、CLVD源、DC源和X源的RV，其中模擬觀測(cè)位移的爆炸源、CLVD源、DC源和X源的深度均為10 km，略去數(shù)值計(jì)算過程，給出最終結(jié)果，見表2。

計(jì)算結(jié)果表明，當(dāng)理論格林函數(shù)的震源深度越接近真實(shí)深度時(shí)，RV越大，反演結(jié)果越接近真實(shí)情況。另外，當(dāng)理論震源深度大于實(shí)際震源深度時(shí)，震源深度誤差對(duì)反演的影響較小，如在理論震源深度為10.5 km的效果比深度為9.5 km的反演效果好，與許力生等[23]的結(jié)論一致。

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

3 結(jié) 論

不同的震源機(jī)制具有不同的矩張量形式，利用反演得到的矩張量的特征值及特征向量可以對(duì)震源機(jī)制進(jìn)行分析，震源矩張量可以分解為對(duì)角線部分及偏量部分，地下核爆炸的震源與天然地震震源的矩張量中的對(duì)角部分和偏量部分具有不同的表現(xiàn)形式。本文中利用數(shù)值方法，驗(yàn)證了基于單臺(tái)數(shù)據(jù)的震源矩張量反演方法，分析了存在噪聲干擾、震相到時(shí)誤差、速度模型誤差和震源深度定位誤差的情況下，反演方法的穩(wěn)定性。各種誤差的分析結(jié)果表明，震源深度誤差對(duì)反演結(jié)果的影響較大。地震波的實(shí)際傳播路徑和地殼的三維結(jié)構(gòu)對(duì)記錄波形影響較大，同時(shí)臺(tái)站的記錄數(shù)據(jù)可能會(huì)受到其他的干擾因素，數(shù)據(jù)品質(zhì)會(huì)受到影響，因此由一維地殼速度模型反演得到準(zhǔn)確的結(jié)果是困難的，即使利用對(duì)分層結(jié)構(gòu)敏感的面波。本文的研究結(jié)果對(duì)于從震源的角度來(lái)了解地下核爆炸的物理機(jī)制具有較好的參考意義，同時(shí)也為利用單臺(tái)記錄數(shù)據(jù)反演震源矩張量、并進(jìn)一步進(jìn)行識(shí)別提供了理論支持。

[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波的激發(fā)機(jī)制[J].地球物理學(xué)報(bào),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] 何永鋒,趙克常,張獻(xiàn)兵,等.地下核爆炸地震波二次源特征[J].地球物理學(xué),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年新疆伽師震源區(qū)特征的地震學(xué)方法研究[D].北京:中國(guó)地震局地球物理研究所,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] 鄭建常,陳運(yùn)泰.基于Langston分解和Hilbert變換約束的區(qū)域偏量矩張量反演方法及應(yīng)用[J].地震學(xué)報(bào),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] 林向東，葛洪魁，徐平, 等.近場(chǎng)全波形反演：蘆山7.0級(jí)地震及余震矩張量解[J].地球物理學(xué)報(bào),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] 許力生，陳運(yùn)泰.震源深度誤差對(duì)矩張量反演的影響[J].地震學(xué)報(bào),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

(責(zé)任編輯 丁 峰)

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國(guó)標(biāo)學(xué)科代碼1303520

A

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

2016-01-27；

2016-08-29

國(guó)家自然科學(xué)基金項(xiàng)目(41374068)

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