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

何永鋒,李 鍇,劉炳燦,姚國政,趙克常,張獻兵,曾樂貴

(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。