基于粒子群優(yōu)化的雙力偶模型振幅譜反演方法及應(yīng)用

鄭建常 陳運(yùn)泰

1）中國北京100081中國地震局地球物理研究所

2）中國濟(jì)南250014山東省地震局

基于粒子群優(yōu)化的雙力偶模型振幅譜反演方法及應(yīng)用

鄭建常1，2），陳運(yùn)泰1）

1）中國北京100081中國地震局地球物理研究所

2）中國濟(jì)南250014山東省地震局

發(fā)展了一種基于全波形振幅譜的頻率域雙力偶震源機(jī)制反演方法.通過理論振幅譜與觀測振幅譜的擬合搜尋斷層面參數(shù)，基于粒子群優(yōu)化算法可以在較短的時間內(nèi)得到穩(wěn)定可靠的解.數(shù)值試驗(yàn)表明，在定位誤差較大，以及臺站布局較差的情況下，振幅譜反演仍可較為準(zhǔn)確地得到震源機(jī)制，并且由此計(jì)算得到的最優(yōu)震源深度仍比較接近真實(shí)的震源位置.使用該方法用2010年5月17日渤海ML4.0地震的震源機(jī)制進(jìn)行了檢驗(yàn)，結(jié)果與加權(quán)P波初動解非常一致.應(yīng)用該方法對山東半島及近海地區(qū)2003—2010年14次ML≥4.0地震震源機(jī)制進(jìn)行了估計(jì).

矩張量 頻率域 震源機(jī)制 粒子群優(yōu)化 振幅譜 全波

引言

地震震源機(jī)制的測定是地震學(xué)的一項(xiàng)基本工作，它有助于了解震源區(qū)和區(qū)域性構(gòu)造應(yīng)力狀態(tài)以及斷層的構(gòu)造特性，是理解地震孕育過程的重要途徑.傳統(tǒng)的求解震源機(jī)制的方法包括初動符號法、振幅比法等.對于中、小地震，在數(shù)據(jù)量較少以及臺站布局等不理想的情況下，這些方法的應(yīng)用都受到一定的限制（吳大銘等，1989）.區(qū)域性地震的波形記錄包含了有關(guān)震源和區(qū)域地殼的豐富的信息，因此，在對地殼結(jié)構(gòu)已知的情況下，可以利用矩張量反演的方法來求解這些地震的震源機(jī)制（Dreger，Helmberger，1993；?íleny｀，Pseneik，1995）.

矩張量的波形反演既可以在時間域內(nèi)進(jìn)行，也可以在頻率域內(nèi)執(zhí)行.對較大的遠(yuǎn)震，矩張量的頻率域反演方法已經(jīng)有了長足的發(fā)展（Mendiguren，1977；Kanamori，Given，1981，1982；Romanowicz，1982；Beck，Patton，1991；Cotton，Campillo，1995；Dufumier，Trampert，1997；Sasatani，1997；周榮茂，陳運(yùn)泰，1999；Hernandezetal，1999；Nakanoetal，2008；Cescaetal，2010）.但是對于區(qū)域性的中小地震，由于其優(yōu)勢頻率較高，易受到噪聲和地殼結(jié)構(gòu)誤差及橫向不均勻性的影響（Fordetal，2010），因此區(qū)域地震的矩張量反演比較復(fù)雜（鄭建常，陳運(yùn)泰，2012）.由于時間域反演需要對每個臺站進(jìn)行震相校正，人工對齊波形的到時，相對而言，波形的振幅譜對震相的平移不敏感，因此在頻率域內(nèi)利用振幅譜的反演方法在研究區(qū)域性地震時更具優(yōu)勢（卞愛飛等，2010）.

Patton和Zandt（1991）首先對用于遠(yuǎn)震的面波反演算法進(jìn)行修改使之適用于區(qū)域性地震，并對內(nèi)華達(dá)州1988年的一次核爆事件進(jìn)行了研究；吳忠良和陳運(yùn)泰（1994）使用直達(dá)P波、直達(dá)S波和SP轉(zhuǎn)換波的位移地震圖在頻率域內(nèi)對1985年4月18日云南祿勸MS6.1地震的15次余震在簡單介質(zhì)模型下進(jìn)行了近震源記錄的矩張量反演，結(jié)果顯示對ML4.0—5.0的地震可以較好地給出震源機(jī)制解.Thio和Kanamori（1995）在利用區(qū)域相速度模型對觀測數(shù)據(jù)進(jìn)行震相校正，消除橫向不均勻結(jié)構(gòu)影響后，使用短周期面波（10—50s）譜反演了TERRA－scope臺陣記錄的超過180次的MW3.2—6.5南加州地方震事件的矩張量，結(jié)果與初動符號和其它波形反演方法得到的結(jié)果一致；Pasyanos等（1996）比較了區(qū)域面波反演與時間域矩張量反演方法的結(jié)果，并進(jìn)行了中等地震矩張量解的準(zhǔn)實(shí)時處理的嘗試.Zahradník等（2001）發(fā)展了利用振幅譜與初動符號聯(lián)合求解小震震源機(jī)制的方法，并應(yīng)用該方法對希臘科林斯灣（Corinth Gulf）地區(qū)5次MW＜4的中小地震震源機(jī)制進(jìn)行了研究；Cesca等（2006）提出了使用直達(dá)體波振幅譜在頻率域反演偏量矩張量解的方法，使用合成數(shù)據(jù)對該方法進(jìn)行了數(shù)值試驗(yàn)，在區(qū)域震中距范圍內(nèi)對西班牙地區(qū)中低震級的淺震震源機(jī)制的應(yīng)用研究表明這一方法是穩(wěn)定可靠的；許力生等（2007）應(yīng)用頻率域矩張量反演方法獲得了首都圈地區(qū)2004年51次中小地震的矩張量解，并通過數(shù)值試驗(yàn)對反演結(jié)果進(jìn)行了評價；Lashak等（2010）使用頻率域矩張量反演方法，研究了伊朗地區(qū)9次中等強(qiáng)度的區(qū)域地震的震源機(jī)制.

近年來，隨著數(shù)字化臺網(wǎng)建設(shè)的加快，尤其是我國東部地區(qū)的地震臺網(wǎng)密度有了很大的提高，可以求得更加可靠的小震震源機(jī)制解.但是由于受場地條件的限制，對于發(fā)生于海域的區(qū)域性地震，多數(shù)臺站距離較遠(yuǎn)，無法得到清晰的初動；而一些近海和島嶼上的臺站，雖然可以記錄到波形，但臺站的方位覆蓋差，并且由于受到海浪和潮汐的影響，噪聲干擾大，因此在使用振幅比或波形反演求解震源機(jī)制的方法時都受到了很大限制.鑒于這種情況，本研究致力于發(fā)展在頻率域內(nèi)使用理論振幅譜來擬合觀測振幅譜的矩張量反演方法，以求解區(qū)域性中、小地震的震源機(jī)制.

1 理論與方法

1.1 基本理論

在點(diǎn)源和遠(yuǎn)場近似的情況下，地震在觀測點(diǎn)r引起的位移可以表示為（許力生，陳運(yùn)泰，2002）

式中，ω為角頻率，ui（r，ω）為地震的觀測位移譜，Mjk（ω）為地震矩張量譜，Gij，k（r，ω）為格林函數(shù)（微商）譜.對于中小地震，由于震源破裂時間很短，可以用狄拉克δ－函數(shù)來表示震源時間函數(shù)（許力生等，2007）.這種情況下，式（1）可以寫成

對于近區(qū)域中小地震，只考慮矩心周圍的格林函數(shù)一階項(xiàng)，位移譜ui（ω）可以表示為對應(yīng)6種基本震源類型的基本地震圖的譜Eij（ω）的加權(quán)和

系數(shù)aj與矩張量Mjk有關(guān)，在地理坐標(biāo)系（x朝北為正，y朝東為正，z朝下為正）中，

由矩張量Mij的特征矢量可得斷層機(jī)制的走向、傾角、滑動角；由特征值可得標(biāo)量地震矩M0；在純雙力偶震源模型的約束條件下，矩張量的跡和行列式為零.

1.2 反演方法

在純剪切位錯的情況下，由于震源機(jī)制中沒有體積分量，因此無需考慮各向同性分量的基本地震圖，即公式（4）中a6＝0.其余5個權(quán)重系數(shù)ai可以表示為斷層的走向φ、傾角δ和滑動角λ的函數(shù).這5個系數(shù)與5種基本理論振幅褶積求和，由于格林函數(shù)譜是復(fù)數(shù)譜，對結(jié)果取模得到理論振幅譜.為了便于比較，反演時對理論和觀測振幅譜進(jìn)行歸一化

式中，上下標(biāo)s和k分別表示臺站和震相.

本文使用基于一維速度模型的離散波數(shù)方法計(jì)算Green函數(shù)（Bouchon，1981）.使用L1范數(shù)來計(jì)算觀測振幅譜與合成振幅譜之間的錯配程度，錯配函數(shù)定義如下：

式中，obs表示觀測振幅譜，syn表示合成振幅譜，n＝臺站數(shù)×分向數(shù)×使用頻率數(shù).

對于約束為純雙力偶模型的矩張量，其震源機(jī)制的求解過程是典型的非線性最優(yōu)化問題.傳統(tǒng)的格點(diǎn)搜索會耗費(fèi)大量的計(jì)算時間，且得到的解不是十分精確.因此在實(shí)際應(yīng)用中，多采用不同的非線性優(yōu)化算法（Sambridge，Drijkoningen，1992；Kobayashi，Nakanishi，1994；Zhouetal，1995；楊文采，2002；?íleny｀，1998；Wuetal，2008；Godanoetal，2009）.我們使用了最近發(fā)展的粒子群算法來進(jìn)行震源機(jī)制解的非線性優(yōu)化搜索.

1.3 粒子群算法

粒子群算法（particle swarm optimization algorithm，簡寫為PSO）是近年來發(fā)展起來的一種新的全局優(yōu)化算法（Kennedy，Eberhart，1995），屬于進(jìn)化算法（evolutionary algorithm）的一種.與遺傳算法相似，粒子群算法也是從隨機(jī)解出發(fā)，通過群體迭代尋找最優(yōu)解，也是通過適應(yīng)度來評價解的品質(zhì).但該算法比遺傳算法規(guī)則簡單，沒有遺傳算法的“交叉”和“變異”操作，而是通過粒子在解空間追隨當(dāng)前的最優(yōu)粒子來尋找全局最優(yōu).該算法基于復(fù)雜適應(yīng)系統(tǒng)（complex adaptive system，簡寫為CAS）理論，通過模擬生物群體遷徙、覓食、躲避獵食者的行為，離散的個體之間通過共享信息和集體協(xié)作使群體達(dá)到最優(yōu)目的，是一種基于群體智能的優(yōu)化方法.

粒子群算法相對于遺傳算法和神經(jīng)網(wǎng)絡(luò)算法，具有一些優(yōu)點(diǎn).Hassan等（2005）以及Yang等（2008）通過一系列基準(zhǔn)問題的測試，證明了與遺傳算法相比，粒子群優(yōu)化算法達(dá)到同樣精度的解所需的計(jì)算用時更少，并且對于連續(xù)變量的非線性問題，粒子群優(yōu)化算法的計(jì)算效率更高.

為了避免反演過程收斂到局部極小值，我們對粒子群算法進(jìn)行了修改，采用斷層的走向、傾角、滑動角的10°×10°×10°步長，對震源機(jī)制參數(shù)的解空間進(jìn)行網(wǎng)格搜索得到的結(jié)果作為初始解，這樣就保證了搜索過程可以穩(wěn)定地收斂到全局最優(yōu)解.數(shù)值試驗(yàn)表明，使用該算法進(jìn)行搜索，反演過程可以較快地收斂到真實(shí)解（圖1）.

2 數(shù)值試驗(yàn)

2.1 用于數(shù)值試驗(yàn)的數(shù)據(jù)與資料

為了檢驗(yàn)方法的可行性，采用合成數(shù)據(jù)進(jìn)行數(shù)值試驗(yàn).用于數(shù)值試驗(yàn)的臺站分布見圖2.為了模擬實(shí)際地震的情況，假設(shè)臺網(wǎng)給出的定位位置存在偏差，假設(shè)定位的震中相對真實(shí)的震中位置在經(jīng)度和緯度上各偏移了0.05°，相當(dāng)于偏離了約8km，這在近區(qū)域或者地方震的定位中屬于偏離較大、精度較差的結(jié)果.用于檢驗(yàn)的“真實(shí)”震源深度設(shè)為9km，標(biāo)量地震矩M0＝1.0×1015N·m，矩震級MW4.0；震源機(jī)制設(shè)為走滑型，其斷層面參數(shù)設(shè)為（該震源機(jī)制來自2010年10月24日周口太康MS4.6地震）：走向224°、傾角89°、滑動角－172°（圖2）.同時，參考實(shí)際需要，為了模擬近海海域地震觀測的實(shí)際情況，給出的臺站分布的方位角張角范圍不到90°，臺站的震中距范圍在100—300km之間.臺站相對真實(shí)震源的震中距和方位角見表1.

圖1 粒子群算法搜索過程（a）迭代10次時解的分布；（b）迭代50次時解的分布；（c）迭代100次時解的分布；（d）迭代結(jié)束時解的分布.（a）—（d）圖中左半部分為迭代中解的錯配值變化，其中藍(lán)色點(diǎn)為迭代中所有解的錯配平均值，黑色點(diǎn)為這些解中的相對最優(yōu)解；右半部分為解在震源球上的分布示意圖，藍(lán)色空心圓表示初始解的P軸在震源球上的投影；綠色實(shí)心圓為每一步迭代過程中解的P軸的投影Fig.1 Convergence of PSO algorithm’s searching procedure on focal sphere（a）Solution after 10iterations；（b）After 50iterations；（c）After 100iterations；（d）Iteration end.Left part of each panel shows misfit values for solutions in each iteration step.Blue points are mean misfit value for all solutions in each step，while black points are misfit value for relative optimized solution in each step.Right part of each panel denotes solution distribution on the focal sphere.Blue circles on the sphere denote Paxes of initial solutions，while green points are that of iterative solutions in each step

圖2 用于數(shù)值試驗(yàn)的臺站分布（a）及震源機(jī)制（b）（紅色實(shí)心圓表示定位結(jié)果，藍(lán)色空心圓表示實(shí)際震中）Fig.2 （a）shows earthquakes（circle）and stations（triangles），and（b）plots the focal mechanism used in numerical experiments.In（a）figure solid circle represents location result，blue circle represents true epicenter of test event

表1 用于數(shù)值試驗(yàn)的臺站的震中距及方位角Table 1 Epicentral distance and azimuth of stations for numerical experiments

使用Bouchon（1981）的離散波數(shù)法，在0.01—1Hz頻段內(nèi)計(jì)算合成波形，作為用于數(shù)值檢驗(yàn)的真實(shí)波形.圖3給出了ST3臺的添加了10%噪聲的合成波形.

2.2 計(jì)算結(jié)果

在與真實(shí)震中位置存在偏離的定位位置上進(jìn)行反演，使用如圖3添加了10%水平均值為零的高斯白噪聲信號的合成波形作為真實(shí)波形，像真實(shí)數(shù)據(jù)一樣進(jìn)行濾波.在深度方向進(jìn)行搜索，起始深度為3km，搜索步長1km，最大深度20km，共18個深度位置.選用了0.10—0.20Hz的頻段進(jìn)行震源參數(shù)的反演，最終計(jì)算得到各個深度的最優(yōu)解以及錯配函數(shù)值見圖4.反演得到的最佳深度h＝10km，與真實(shí)的深度h＝9km相差1km.在該深度位置上得到的標(biāo)量地震矩為1.16×1015N·m，矩震級MW4.0，與真實(shí)值相差不大.反演得到的最佳震源機(jī)制與真實(shí)震源機(jī)制之間的差別不大，解的詳細(xì)對比情況見表2.

表2 真實(shí)震源機(jī)制與反演結(jié)果的斷層面參數(shù)的比較Table 2 Comparison between true solution and inversion result

由表2可見，反演得到的結(jié)果與真實(shí)解相差不大，得到的兩個節(jié)面的走向的差異均在5°范圍內(nèi)，其中P軸方位的差異不到1°，T軸方位差在10°以內(nèi).

由真實(shí)解計(jì)算得到的觀測振幅譜與對應(yīng)搜索得到的最佳解的理論譜的對應(yīng)關(guān)系見圖5.值得說明的是，我們反演使用的頻段為0.10—0.20Hz，分別在0.08—0.10Hz和0.20—0.208Hz頻段使用了余弦邊瓣.由圖5可以看出，除了個別頻率點(diǎn)外，大部分頻率都符合得相當(dāng)好.

對10km深度得到的最佳解計(jì)算合成波形，與真實(shí)解的合成波形進(jìn)行比較，并計(jì)算了相關(guān)系數(shù)，結(jié)果見表3及圖6.可以看出，大部分臺站分向的波形擬合得相當(dāng)好.

表3 最佳解合成波形與真實(shí)解合成波形的相關(guān)系數(shù)比較Table 3 Correlations between observed and synthetic waveforms

3 方法檢驗(yàn)

3.1 數(shù)據(jù)與資料

據(jù)山東臺網(wǎng)測定，北京時間2010年5月17日20時2分13.0秒，在渤海海峽（38.38°N，120.53°E）發(fā)生了ML4.0地震.

這次地震發(fā)生在廟島群島的西北側(cè)，距離遼東半島陸地超過100km，距離山東半島最近不到70km.我們使用了來自于山東臺網(wǎng)的波形數(shù)據(jù)，在震中距150km范圍內(nèi)主要有山東臺網(wǎng)的6個臺站（BHC，CHD，LOK，ZHY，YTA，LZH），主要分布在地震的南側(cè)，以及遼寧臺網(wǎng)用于數(shù)據(jù)交換的DL2和HSH兩個臺，位于地震的北側(cè).地震及臺站分布見圖7.

距離最近的BHC臺存在問題，數(shù)據(jù)無法使用.根據(jù)臺站波形的信噪比情況，選擇了CHD，LOK，DL2，ZHY四個臺站的數(shù)據(jù)用于反演.首先對臺站的波形數(shù)據(jù)進(jìn)行傅里葉變換，得到了觀測譜.圖8給出了CHD臺（震中距55.22km）和DL2臺（震中距111.86km）的垂直向波形的速度振幅譜和噪聲譜.由圖8可見，兩個臺均有較好的信噪比.

圖9 本文使用的地殼速度結(jié)構(gòu)模型Fig.9 Crustal velocity structure used in this paper

反演使用的速度模型是Park等（2006）給出的黃海地區(qū)中上地殼速度結(jié)構(gòu)模型（圖9）.該模型也是韓國地礦研究所（KIGAM，Korea Institute of Geology，Mining and Materials）使用的黃海至朝鮮半島地區(qū)的地震定位速度模型.

由于矩張量反演對震中位置的較小誤差不敏感，因此我們使用中國地震臺網(wǎng)中心給出的全國小震目錄中的定位結(jié)果，在深度方向進(jìn)行搜尋.起始深度為3km，搜尋步長1km，最大深度20 km，共18個深度位置.

3.2 0.10—0.20Hz頻段計(jì)算結(jié)果與分析

使用0.10—0.20Hz頻段（周期5—10s）的振幅譜進(jìn)行反演，圖10給出了不同深度的最優(yōu)解的錯配值隨深度的變化關(guān)系，以及最終得到的全局最優(yōu)解.

反演得到此次地震的最佳深度在9km，該結(jié)果比臺網(wǎng)中心給出的定位結(jié)果（12km）的震源深度要淺.由圖10可以看出，我們使用的周期5—10s的波形對震源深度的分辨能力不是很好，嘗試震源深度在8—13km之間時，解的錯配值的差別不是很大.最佳震源機(jī)制解的錯配函數(shù)misfit＝0.405 7.反演得到此次地震的標(biāo)量地震矩約3.05×1014N·m，相當(dāng)于矩震級MW3.7.

圖10 0.10—0.20Hz頻段反演得到的渤海ML4.0地震不同深度的最優(yōu)解和錯配值（a）以及全局最優(yōu)解（b）Fig.10 Inversion results for Bohai Sea ML＝4.0event based on 0.10—0.20Hz data.（a）Optimized solutions at different focal depths；（b）Global optimized focal mechanism

圖11給出了最佳雙力偶解計(jì)算得到的合成譜與觀測振幅譜.使用反演得到的最佳震源機(jī)制解計(jì)算得到合成波形，圖12給出了在對齊到時的情況下CHD臺的波形擬合情況，3個分向的平均相關(guān)系數(shù)達(dá)到了0.757.由圖可見，合成波形能夠再現(xiàn)觀測波形的基本震幅和主要特征，由理論譜得到的合成波形與實(shí)際的觀測波形擬合的相當(dāng)一致，這表明了方法的可靠性以及反演過程的準(zhǔn)確性.

3.3 0.10—0.25Hz頻段計(jì)算結(jié)果與分析

使用0.10—0.25Hz頻段進(jìn)行反演，得到了各個深度的最佳震源機(jī)制解以及相應(yīng)的觀測振幅譜與合成譜的錯配值（圖13）.由圖13可見，最佳震源深度仍然為9km，與0.10—0.20Hz頻段反演結(jié)果一致.在此深度上，最佳震源機(jī)制解對應(yīng)的錯配函數(shù)最小值為0.397，得到的標(biāo)量地震矩為3.339×1014N·m，相當(dāng)于矩震級MW3.7.與5—10s周期反演結(jié)果相比，最佳震源機(jī)制非常一致，標(biāo)量地震矩也很接近.

由反演得到的最佳雙力偶模型震源機(jī)制計(jì)算的合成波形與觀測波形的平均相關(guān)系數(shù)約0.739（圖略），稍低于0.10—0.20Hz頻段的結(jié)果.由此可見，在反演過程中增加相對高頻的成分，可以提高深度的分辨能力，但受一維速度結(jié)構(gòu)模型的限制，有可能會減弱合成波形的擬合程度.

3.4 與初動符號解的比較

根據(jù)山東臺網(wǎng)的波形資料，讀取了17個臺的初動符號，使用P波初動格點(diǎn)嘗試法求解震源機(jī)制，得到的結(jié)果見表4.

與初動符號的震源機(jī)制相比，二者之間的最小空間旋轉(zhuǎn)角（Kagan，1991）為11.24°，節(jié)面走向的差不到10°，P和T軸的方位位置均比較接近.進(jìn)一步說明了方法的可靠性.

4 山東半島及近海地區(qū)地震的反演結(jié)果

2003年以來，山東半島兩側(cè)近海海域陸續(xù)發(fā)生了多次3—4級中小地震，受臺站布局以及觀測資料的限制，對這些地震一直沒有開展過研究.我們采用上述的基于粒子群優(yōu)化的雙力偶約束的譜相關(guān)反演方法，對該區(qū)的地震震源機(jī)制進(jìn)行了初步分析.為了保證反演質(zhì)量，使用波形的相對低頻部分，在合適的信噪比的情況下，我們共選擇了該區(qū)2003—2010年ML≥4.0的14次地震用于計(jì)算.圖14給出了這14次地震的震源機(jī)制解下半球投影的“海灘球”表示，表5中列出了這些地震的震源參數(shù).

圖13 0.10—0.25Hz頻段反演得到的渤海ML4.0地震不同深度的最優(yōu)解和錯配值（a）以及全局最優(yōu)解（b）Fig.13 Inversion results for Bohai Sea ML＝4.0event based on 0.10—0.25Hz data.（a）Optimized solutions at different focal depths；（b）Global optimized focal mechanism

表4 由P波初動符號得到的斷層面解與振幅譜相關(guān)方法得到解的比較Table 4 Comparison of the fault plane solution obtained by polarity analysis with that of this paper

圖14 研究區(qū)域的臺站分布及本文得到的震源機(jī)制解Fig.14 Stations in studied area and focal mechanisms of earthquakes estimated by this paper

表5 本文得到的14次ML≥4.0地震的震源機(jī)制解參數(shù)Table 5 Focal mechanisms of 14earthquakes studied in this paper

5 討論與結(jié)論

發(fā)展了一種基于全波形振幅譜的頻率域雙力偶震源反演方法.通過理論譜與觀測譜振幅的擬合搜尋斷層面參數(shù)，基于粒子群優(yōu)化算法，可以在較短的時間內(nèi)得到相對穩(wěn)定可靠的解.

數(shù)值試驗(yàn)表明，在定位誤差較大以及臺站布局較差的情況下，通過振幅譜反演可以較為準(zhǔn)確地得到震源機(jī)制，并且計(jì)算得到的最優(yōu)震源深度接近真實(shí)的震源位置.

對于時間域波形反演，振幅譜方法速度更快，扣除儀器響應(yīng)更容易，并且該方法既不需要確切的發(fā)震時刻，又無需對觀測波形與合成波形進(jìn)行人工對齊，減少了中間誤差.此外，由于數(shù)據(jù)量少，觀測譜與合成譜之間的擬合程度的量化也比波形要容易.因此更適合于區(qū)域中等地震的矩張量反演的自動化處理.

使用該方法對2010年5月17日渤海ML4.0地震的震源機(jī)制進(jìn)行了檢驗(yàn)計(jì)算.計(jì)算結(jié)果顯示，此次地震的最佳深度為9km，標(biāo)量地震矩M0＝3.049×1014N·m，相當(dāng)于矩震級MW3.7；不同頻段的反演結(jié)果相對穩(wěn)定，震源機(jī)制結(jié)果與P波初動解的結(jié)果較為一致.

應(yīng)用該方法對2003年以來山東半島及其兩側(cè)近海海域發(fā)生的14次ML≥4.0地震的震源機(jī)制進(jìn)行了研究，給出了最佳的雙力偶震源機(jī)制解.反演得到的最佳深度與臺網(wǎng)給出的地震目錄的深度不盡一致.根據(jù)有關(guān)研究（許力生，陳運(yùn)泰，1997），深度的誤差對震源機(jī)制的結(jié)果影響不大.此外需要說明的是，受臺站布局和早期數(shù)字化觀測的數(shù)據(jù)質(zhì)量較差的影響，對于距離陸地較遠(yuǎn)的海域地震，其結(jié)果的可靠程度相對較低.

感謝審稿專家為本文提出的修改意見和建議.本文使用的粒子群優(yōu)化算法程序?yàn)镾 Chen（2010）給出的基于 Matlab的源代碼（http：∥www.mathworks.com/matlabcentral/fileexchange/25986－another－particle－swarm－toolbox），并進(jìn)行了部分修改，在此謹(jǐn)致謝忱.

卞愛飛，於文輝，周華偉.2010.頻率域全波形反演方法研究進(jìn)展［J］.地球物理學(xué)進(jìn)展，25（3）：982－－993.

吳大銘，王培德，陳運(yùn)泰.1989.用SH波和P波振幅比確定震源機(jī)制解［J］.地震學(xué)報(bào)，11（3）：275－－281.

吳忠良，陳運(yùn)泰.1994.近震源寬頻帶記錄的地震矩張量反演［J］.地震學(xué)報(bào)，16（2）：141－－152.

許力生，陳運(yùn)泰.1997.震源深度誤差對矩張量反演的影響［J］.地震學(xué)報(bào)，19（5）：462－－470.

許力生，陳運(yùn)泰.2002.震源時間函數(shù)與震源破裂過程［J］.地震地磁觀測與研究，23（6）：1－－8.

許力生，蔣長勝，陳運(yùn)泰，李春來，張?zhí)熘?2007.2004年首都圈地區(qū)中小地震的矩張量反演［J］.地震學(xué)報(bào)，29（3）：229－－239.

楊文采.2002.非線性地球物理反演方法：回顧與展望［J］.地球物理學(xué)進(jìn)展，27（2）：255－－261.

鄭建常，陳運(yùn)泰.2012.稀疏臺網(wǎng)反演區(qū)域地震偏量矩張量解的穩(wěn)定性［J］.地震學(xué)報(bào)，34（1）：31－－43.

周榮茂，陳運(yùn)泰.1999.由矩張量反演得到的北部灣地震的震源機(jī)制［J］.地震學(xué)報(bào)，21（6）：561－－569.

Beck L S，Patton H J.1991.Inversion of regional surface wave spectra for source parameters of aftershock from the Loma－Prieta earthquake［J］.BullSeismSocAmer，81（5）：1726－－1736.

Bouchon M.1981.A simple method to calculate Green’s functions for elastic layered media［J］.BullSeismSocAmer，71（4）：959－－971.

Cesca S，Buforn E，Dahm T.2006.Amplitude spectra moment tensor inversion of shallow earthquakes in Spain［J］.GeophysJInt，166（2）：839－－854.

Cesca S，Heimann S，Dahm T.2010.Rapid directivity detection by azimuthal amplitude spectra inversion［J］.JSeismol，15（1）：147－－164.

Cotton F，Campillo M.1995.Frequency domain inversion of strong motions：Application to the 1992Landers earthquake［J］.JGeophysRes，100（B3）：3961－－3975.

Dreger D S，Helmberger D V.1993.Determination of source parameters at regional distances with three－component sparse network data［J］.JGeophysRes，98（B5）：8107－－8125.

Dufumier H，Trampert J.1997.Contribution of seismic tomography in moment－tensor inversions using teleseismic surface－wave spectra［J］.BullSeismSocAmer，87（1）：114－－122.

Ford S R，Dreger D S，Walter W R.2010.Network sensitivity solutions for regional moment tensor inversions［J］.Bull SeismSocAmer，100（5A）：1962－－1970.

Godano M，Regnier M，Deschamps A，Bardainne T，Gaucher E.2009.Focal mechanisms from sparse observations by nonlinear inversion of amplitude：Method and tests on synthetic and real data［J］.BullSeismSocAmer，99（4）：2243－－2264.

Hassan R，Cohanim B，de Weck O.2005.A comparison of particle swarm optimization and the generic algorithm［C］∥Proceedingsofthe46thAIAA/ASME/ASCE/AHS/ASCStructures，StructuralDynamicsandMaterialsConference.Austin，Texas：AIAA：1－－13.

Hernandez B，Cotton F，Campillo M.1999.Contribution of radar interferometry to a two－step inversion of the kinematic process of the 1992Landers earthquake［J］.JGeophysRes，104（B6）：13083－－13099.

Kagan Y Y.1991.3－D rotation of double－couple earthquake sources［J］.GeophysJInt，106（3）：709－－716.

Kanamori H，Given J W.1981.Use of long－period surface waves for rapid determination of earthquake source parameters［J］.PhysEarthPlanetInteri，27（1）：8－－31.

Kanamori H，Given J W.1982.Use of long－period surface waves for rapid determination of earthquake source parameters 2：Preliminary determination of source mechanisms of large earthquakes（MS≥6.5）in 1980［J］.PhysEarthPlanet Interi，30（2－－3）：260－－268.

Kennedy J，Eberhart R.1995.Particle swarm optimization［C］∥ProceedingsoftheIEEEInternationalConferenceon NeuralNetworks.Perth，Australia：IEEE：1942－－1945.

Kobayashi R，Nakanishi I.1994.Application of genetic algorithms to focal mechanism determination［J］.GeophysRes Lett，21（8）：729－－732.

Lashak A B，Zare M，Mortezanejad G，Bayranvand S P.2010.Moment tensor inversion of nine events in Iran using INSN data［J］.JSeismol，14（4）：751－－760.

Mendiguren J.1977.Inversion of surface wave data in source mechanism studies［J］.GeophysResLett，82（5）：889－－894.Nakano M，Kumagai H，Inoue H.2008.Waveform inversion in the frequency domain for the simultaneous determination of earthquake source mechanism and moment function［J］.GeophysJInt，173（3）：1000－－1011.

Pasyanos M E，Dreger D S，Romanowicz B.1996.Toward real－time estimation of regional moment tensors［J］.Bull SeismSocAmer，86（5）：1255－－1269.

Patton H J，Zandt G.1991.Seismic moment tensors of western U S earthquakes and implication for the tectonic stress field［J］.JGeophysRes，96（B11）：18245－－18259.

Romanowicz B A.1982.Moment tensor inversion of long period Rayleigh waves：A new approach［J］.JGeophysRes，87（B7）：5395－－5407.

Sambridge M，Drijkoningen G.1992.Genertic algorithms in seismic waveform inversion［J］.GeophysJInt，109（2）：323－－342.

Sasatani T.1997.Source characteristics of the 1994Hokkaido Toho－Oki earthquake deduced from wide band strong motion records［J］.JourFacSciHokkaidoUniv，SerVII（Geophysics），10（2）：269－－293.

?íleny｀J，Pseneik I.1995.Mechanisms of local earthquakes in 3－D inhomogeneous media determined by waveform inversion［J］.GeophysJInt，121（2）：459－－474.

?íleny｀J.1998.Earthquake source parameters and their confidence regions by agenetic algorithm with a‘memory’［J］.GeophysJInt，134（1）：228－－242.

Thio H K，Kanamori H.1995.Moment tensor inversions for local earthquakes using surface waves recorded at TERRA－scope［J］.BullSeismSocAmer，85（4）：1021－－1038.

Wu Y M，Zhao L，Chang C H，Hsu Y J.2008.Focal mechanism determination in Taiwan by genetic algorithm［J］.Bull SeismSocAmer，98（2）：651－－661.

Yang F，Zhang C，Sun T.2008.Comparison of particle swarm optimization and genetic algorithm for HMM training［C］∥19thInternationalConferenceonPatternRecognition.New York：IEEE，2008：1－－4.

Zahradník J，Jansky J，Paptsimpa K.2001.Focal mechanism of weak earthquakes from amplitude spectra and polarities［J］.PureApplGeophys，158（4）：647－－665.

Zhou R，Tajima F，Stoffa P L.1995.Earthquake source parameter determination using genetic algorithms［J］.Geophys ResLett，22（4）：517－－520.

Amplitude spectrum inversion for double－couple source model with particle swarm optimization algorithm

Zheng Jianchang1，2），Chen Yun－tai1）

1）InstituteofGeophysics，ChinaEarthquakeAdministration，Beijing100081，China
2）EarthquakeAdministrationofShandongProvince，Jinan250014，China

An inversion method is proposed to best fit the observation in frequency domain with synthetic amplitude spectra of full waveform in seismic source mechanism studies.Grid search is performed for finding fault parameters by correlating theoretical spectra with observations.Employing particle swarm optimization algorithm，stable and reliable solutions can be obtained in shorter time.Numerical test shows that，in the case of large location error and worse station coverage，the inversion can result in the solutions close to the given true focal mechanism，and the focal depth near the real focus.The method has been applied to processing local and regional recordings of the 17May 2010ML4.0 earthquake in Bohai Sea region.The result is in agreement with that obtained from P wave polarity analysis.We also estimated focal mechanism solutions of 14ML≥4.0earthquakes in Shandong peninsula and its adjacent sea area from 2003to 2010.

moment tensor；frequency domain；focal mechanism；particle swarm optimization；amplitude spectrum；full waveform

10.3969/j.issn.0253－3782.2012.03.003

P315.3＋3

A

鄭建常，陳運(yùn)泰.2012.基于粒子群優(yōu)化的雙力偶模型振幅譜反演方法及應(yīng)用.地震學(xué)報(bào)，34（3）：308－－322.

Zheng Jianchang，Chen Yun－tai.2012.Amplitude spectrum inversion for double－couple source model with particle swarm optimization algorithm.ActaSeismologicaSinica，34（3）：308－－322.

地震科技星火計(jì)劃項(xiàng)目（XH12028Y）、山東省科技發(fā)展計(jì)劃項(xiàng)目（2008GG10008004）與中國地震局監(jiān)測預(yù)報(bào)司震情跟蹤定向工作任務(wù)項(xiàng)目（2010020102）及山東省地震局博士基金聯(lián)合資助.

2011－05－02收到初稿，2011－08－16決定采用修改稿.

e－mail：zjcmail＠yeah.net