CAP方法反演震源機(jī)制的誤差分析：以膠東半島兩次顯著中等地震為例

鄭建常，林眉，王鵬，徐長朋

山東省地震局，濟(jì)南 250014

1 引言

據(jù)山東臺網(wǎng)測定，2013年11月23日13時(shí)44分在山東省萊州市（37.10°N，120.02°E）發(fā)生M4.6級地震，這次地震是山東陸地地區(qū)自1995年蒼山5.2級地震后發(fā)生的最大地震，影響范圍廣，山東東部市地普遍有感；2014年1月7日22時(shí)24分在山東乳山（36.80°N，121.70°E）發(fā)生M4.3級地震，這次地震也造成膠東地區(qū)大面積有感.這兩次事件是1970年以來膠東半島陸地及近海地區(qū)發(fā)生的最強(qiáng)烈的地震活動(dòng)，其中萊州地震震中區(qū)在1970年以來的小震目錄上屬于典型的少震、弱震區(qū)，活動(dòng)水平不高，很少有ML≥3.0地震發(fā)生，僅在1991年2月以及2012年7月分別發(fā)生最大ML3.8級小震序列各一次；乳山地震震中區(qū)歷史上曾發(fā)生公元1046年岠嵎山51／2級和1939年乳山下初51／2級地震.雖然這兩次地震的震中區(qū)歷史上沒有強(qiáng)烈地震活動(dòng)，但膠東半島北部近海曾發(fā)生多次6、7級強(qiáng)震，如1548年渤海海峽7.0、1948年威海近海6.0以及1969年渤海7.4級等.因此確定這兩次顯著中等地震的震源機(jī)制對于研究區(qū)域地質(zhì)構(gòu)造的活動(dòng)特征，以及研判該地區(qū)的地震危險(xiǎn)性等具有重要的科學(xué)價(jià)值.

我們使用近年來在國內(nèi)得到廣泛使用的CAP（Cut and Paste）方法反演這兩次地震的震源機(jī)制.為了得到更準(zhǔn)確的解，并且合理地估計(jì)反演結(jié)果的不確定性，我們使用自助抽樣法（bootstrap）對反演過程隨機(jī)重復(fù)，在大樣本量的反演結(jié)果基礎(chǔ)上，使用粒子群算法搜索優(yōu)化解，利用動(dòng)態(tài)聚類技術(shù)對結(jié)果進(jìn)行聚類分析，從而得到了更加穩(wěn)定可靠的斷層面解，給出了可能的誤差范圍，并進(jìn)一步給出了震源球上P、T軸的概率密度分布.

2 理論與方法

2.1 CAP方法反演震源機(jī)制

震源機(jī)制和傳播效應(yīng)決定了觀測波形的變化.如果地殼模型已知，可以準(zhǔn)確地計(jì)算波形傳播過程中的效應(yīng)，因此我們可以通過理論波形s（t）和觀測波形u（t）的擬合來估計(jì)震源的斷層面參數(shù).

雙力偶震源產(chǎn)生的理論位移s（t）可以表示為（Zhu and Helmberger，1996）：

其中，i＝1，2，3對應(yīng)三種基本斷層響應(yīng)，即：垂直走滑、垂直傾滑以及傾角為45°的傾滑；Gi為格林函數(shù)，Ai是輻射系數(shù)，φ是臺站方位角，M0為標(biāo)量地震矩.θ，δ，λ分別為斷層的走向、傾角、滑動(dòng)角.系數(shù)Ai由6個(gè)矩張量分量和臺站方位角表示如下（Jost and Herrmann，1989）：

走向θ、傾角δ、滑動(dòng)角λ，以及標(biāo)量地震矩M0等可以通過求解以下方程進(jìn)行估計(jì)：

波形反演可以使用全波形數(shù)據(jù)，也可以單獨(dú)使用體波或面波震相進(jìn)行擬合.CAP方法是一種聯(lián)合使用體波和面波進(jìn)行反演的方法，近年來在國內(nèi)得到了廣泛的應(yīng)用（呂堅(jiān)等，2008；黃建平等，2009；鄭勇等，2009；龍鋒等，2010；韓立波等，2012），由于該方法分別截取波形的Pnl部分和面波部分分別擬合（Zhao and Helmberger，1994；Zhu and Helmberger，1996），并在反演的過程中允許它們在適當(dāng)?shù)臅r(shí)間變化范圍內(nèi)相對移動(dòng)，在一定程度上避免了因?yàn)榈貧つＰ筒粶?zhǔn)確而引起的震相到時(shí)的誤差因素，對速度模型和地殼橫向變化的依賴性較小，因此在實(shí)際的區(qū)域地震震源機(jī)制求解中有明顯的優(yōu)勢.CAP方法使用頻率F－波數(shù)K法（Zhu and Rivera，2002）計(jì)算格林函數(shù)，使用網(wǎng)格搜索方法搜尋最優(yōu)震源機(jī)制參數(shù)和震源深度.考慮到波形隨震中距的衰減，方法定義誤差函數(shù)如下：

式中，r為臺站震中距，r0為選定的參考震中距，p為指數(shù)因子.參考有關(guān)研究，對體波p＝1，面波p＝0.5（韓立波和蔣長勝，2012）.

2.2 CAP方法的優(yōu)化解及其不確定性估計(jì)

地球內(nèi)部的任意震源可以表示為6個(gè)獨(dú)立分量的矩張量，由于CAP方法限制震源為雙力偶模型，并且無需發(fā)震時(shí)刻的對齊，因此只需對震源模型的三個(gè)角度，即走向θ、傾角δ、滑動(dòng)角λ，以及標(biāo)量地震矩M0進(jìn)行搜索，理論上而言，僅需要2個(gè)臺站的三分向波形就可以求解；雖然研究顯示，對于大多數(shù)3個(gè)三分向臺的組合，使用波形反演就可以得到相對準(zhǔn)確的震源機(jī)制，但實(shí)際情況也顯示，不同的臺站組合波形反演得到的解之間仍然存在一定的差異（Godano et al.，2009；鄭建常和陳運(yùn)泰，2012）.目前國內(nèi)臺網(wǎng)密度已經(jīng)達(dá)到相當(dāng)水平，在東部地區(qū)，一個(gè)中等地震通常有數(shù)十甚至上百個(gè)寬頻帶臺能夠記錄到清晰的波形，以此次萊州地震為例，通過對原始波形進(jìn)行去均值、去趨勢、積分等簡單變換后，根據(jù)直觀的觀察，震中距300km范圍內(nèi)，采樣率100Hz的寬頻帶三分向波形有近40個(gè)臺站的資料可用.在使用CAP方法求解震源機(jī)制時(shí)，一般的做法是選擇部分波形擬合較好的臺進(jìn)行反演，刪去擬合不好的臺或者震相；有些情況下，甚至僅使用面波部分而刪除體波震相，需知面波尤其是徑向和切向分量，很容易受到臺站下方淺層地殼結(jié)構(gòu)的影響.由于CAP方法是采用網(wǎng)格搜索的方法，因而這種人為的選擇，必然會(huì)為反演結(jié)果增加主觀的不確定性因素.我們無法令人信服地說明，擬合不好的波形究竟是數(shù)據(jù)本身確實(shí)存在干擾，還是說搜索到的解無法滿足該條數(shù)據(jù).另外，CAP方法雖然可以在最后的輸出結(jié)果中給出斷層面參數(shù)的不確定性，但該估計(jì)值只是面向所使用的臺站數(shù)據(jù)的結(jié)果，在上述的人為選擇下，該不確定性估計(jì)能夠在多大程度上客觀地反映最終解的整體不確定性，是無法說明的.

為了求得更加穩(wěn)定可靠的解并且合理客觀地給出解的誤差估計(jì)，我們在相對豐富的觀測數(shù)據(jù)基礎(chǔ)上，采用自助抽樣統(tǒng)計(jì)方法進(jìn)行分析.具體方法是在可用的觀測臺站中可重復(fù)地隨機(jī)抽取一定數(shù)量的臺站組成新的臺站組合，使用該臺站組合的觀測數(shù)據(jù)重復(fù)反演過程.在大量的重復(fù)計(jì)算后（例如，超過1000次），可以有效地排除觀測質(zhì)量不高或存在較大干擾誤差的數(shù)據(jù)的影響，從而得到更加接近真實(shí)解的結(jié)果，并且可以有效地給出解的不確定性.

另外，由于CAP方法在搜索斷層面解時(shí)采用的是網(wǎng)格搜索的方法，然后通過插值計(jì)算誤差函數(shù)e的最小值，并且由于固定步長的走向、傾角、滑動(dòng)角的嘗試位置在震源球上的分布是不均勻的（許向彤等，1995），因此在最終解中可能會(huì)有空缺（gap）的存在.為了求解優(yōu)化解，我們進(jìn)一步使用Kagan（1991）定義的雙力偶模型最小空間旋轉(zhuǎn)角，對上面自助抽樣得到的大量滿足條件的震源機(jī)制結(jié)果進(jìn)行分析，定義與所有解的空間偏轉(zhuǎn)角度和為目標(biāo)函數(shù)，使用粒子群非線性優(yōu)化方法搜索該目標(biāo)函數(shù)最小的結(jié)果，視為最優(yōu)解.

2.3 聚類分析

在震源機(jī)制求解中，常見的情況是在震源球上存在幾簇相對集中分布的解，對這些可能的解直接取數(shù)學(xué)平均是不甚合理的，并且在數(shù)據(jù)存在較大誤差或干擾的情況下，滿足條件的可能解的分布范圍也許會(huì)相當(dāng)大.因此針對這一現(xiàn)象，刁桂苓等（1992）、俞春泉等（2009）分別使用系統(tǒng)聚類和動(dòng)態(tài)聚類技術(shù)，對所有的可能解進(jìn)行聚類分析，求取聚類中心作為反演的優(yōu)化解，數(shù)值試驗(yàn)和實(shí)際應(yīng)用都有很好的效果.聚類分析可以很好地排除孤立解和錯(cuò)誤解，從而在大量的數(shù)據(jù)中獲取更加接近真實(shí)解的結(jié)果.本文在使用不同臺站組合重復(fù)進(jìn)行波形反演后，同樣得到了大量的震源機(jī)制解數(shù)據(jù)，受臺站布局和數(shù)據(jù)誤差的影響，這些解或多或少存在差別，因此對這些結(jié)果進(jìn)行聚類分析是很有必要的.

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

本文使用了山東臺網(wǎng)提供的波形資料，其中還包括了鄰省如遼寧、河北、江蘇等省交換資料的部分臺站.圖1給出了本項(xiàng)研究使用的臺站分布，其中個(gè)別臺如JIM、ZSL、HUD等為短周期臺，在波形反演中沒有使用.本文研究中，首先由觀測記錄直接讀取初動(dòng)符號，用于約束波形反演；然后將觀測數(shù)據(jù)扣除儀器響應(yīng)，經(jīng)過去均值、零漂等預(yù)處理后積分至位移記錄，旋轉(zhuǎn)到Z－R－T坐標(biāo)系，對觀測波形和理論波形同樣進(jìn)行帶通濾波，然后用于反演.

圖1 本文研究的兩次地震震中及山東臺網(wǎng)臺站分布圖Fig.1 Map of stations in Shandong Network and two earthquakes studied in this paper.Red circles denote epicenter，triangles are stations，and solid black lines are faults.

使用Chang等（2006）給出的朝鮮半島南部至黃海地區(qū)的中上地殼速度結(jié)構(gòu)模型用于本文的震源機(jī)制反演.相關(guān)地質(zhì)資料顯示，膠東半島、南黃海以及朝鮮半島南部在大地構(gòu)造分區(qū)上都屬于下?lián)P子地塊，地質(zhì)構(gòu)造屬性相對較為一致（Ree et al.，1996）.

4 結(jié)果與分析

4.1 乳山M4.3震源機(jī)制

選擇震中距在250km以內(nèi)的15個(gè)臺站的寬頻帶波形記錄進(jìn)行反演，Pnl和面波的反演波段分別選擇0.05～0.15Hz和0.033～0.067Hz頻段.圖2給出了不同深度的最佳雙力偶解，及擬合誤差隨不同深度變化的關(guān)系，由圖可見，震源深度在4km時(shí)觀測波形和理論波形的錯(cuò)配值最小，說明事件的震源深度較淺.

由CAP方法反演得到的最佳震源機(jī)制：節(jié)面A的參數(shù)為：走向202°、傾角75°、滑動(dòng)角153°；節(jié)面B的參數(shù)為：走向299.5°、傾角64°、滑動(dòng)角16.7°；參考乳山序列的雙差定位結(jié)果（李冬梅和鄭建常，2014）分析認(rèn)為，節(jié)面B可能是乳山地震的發(fā)震斷層；震源機(jī)制顯示為左旋走滑型，反演得到此次地震的矩震級MW＝4.2.

圖2 2014年1月7日乳山M4.3級地震不同震源深度的波形擬合誤差及最佳震源機(jī)制解Fig.2 Waveform fit errors and best focal mechanisms as function of depth for Jan.7，2014Rushan M4.3event

圖3給出了對應(yīng)最優(yōu)解的理論波形和觀測波形的擬合情況.15個(gè)臺一共75個(gè)震相，其中理論波形與觀測波形相關(guān)系數(shù)大于0.9的有39個(gè)，超過50%；相關(guān)系數(shù)大于0.6（相關(guān)性較好）的有67個(gè)，約占89.3%；最佳解的方差減少（variance reduction）為70.3%，說明理論波形很好地?cái)M合了觀測波形，反演結(jié)果是可靠的.個(gè)別臺（如 WEH）平均相關(guān)系數(shù)較差，可能與臺站位于震源機(jī)制解的節(jié)面線附近，振幅相對較小所致；另外如CHD臺的擬合程度不好，可能與該臺處于海域、噪聲干擾較大有關(guān).

使用自助抽樣的統(tǒng)計(jì)方法，對乳山地震震源機(jī)制解的不確定性進(jìn)行估計(jì).選用震中距在300km以內(nèi)的采樣率為100Hz的22個(gè)三分向?qū)掝l帶臺的觀測波形組成原始數(shù)據(jù)集，為了保證用于反演的數(shù)據(jù)的樣本量，設(shè)用于反演的臺站數(shù)為20個(gè)，對原始數(shù)據(jù)集進(jìn)行每個(gè)臺站等概率、可重復(fù)地隨機(jī)抽取，抽取出的臺站波形組成新的數(shù)據(jù)集，然后用于CAP方法的波形反演.對上述的抽取臺站反演過程重復(fù)1000次，將反演得到的震源機(jī)制的斷層節(jié)面解和P、T軸繪制在一個(gè)震源球上，見圖4.可以看出，反演中除去個(gè)別反演過程的斷層面解出現(xiàn)一定程度的偏離外，其余結(jié)果集中分布，均顯示為近走滑的機(jī)制；圖4中的P、T軸位置和斷層節(jié)面線集中成叢，大致顯示出斷層面解的誤差范圍.

使用粒子群非線性優(yōu)化方法，以與自助抽樣給出的1000個(gè)機(jī)制解（圖4）的Kagan角之和為目標(biāo)函數(shù)，搜索最優(yōu)解.結(jié)果顯示最優(yōu)解為，節(jié)面A：走向208.4°，傾 角89.7°，滑 動(dòng) 角154.3°；節(jié) 面B：走 向298.5°，傾角64.3°，滑動(dòng)角0.3°；最優(yōu)解與圖4所示1000個(gè)解的平均夾角4.37°，以其與所有解Kagan角的2倍標(biāo)準(zhǔn)差為震源機(jī)制解的誤差范圍，結(jié)果顯示不確定性為6.44°（圖5）.

表1 波形反演乳山M4.3地震震源機(jī)制解結(jié)果Table 1 Parameters of focal mechanism results from waveform inversion for Jan.7，2014Rushan M4.3event

對自助抽樣結(jié)果進(jìn)行動(dòng)態(tài)聚類分析，結(jié)果顯示最優(yōu)解為：節(jié)面 A：走向208.0°，傾角89.3°，滑動(dòng)角154.1°，節(jié)面B：走向298.4°，傾角64.1°，滑動(dòng)角1.8°，與粒子群優(yōu)化解非常一致（見表1）.

將自助抽樣結(jié)果中的P、T軸投影到震源球上（圖4），對其進(jìn)行概率密度統(tǒng)計(jì)分析，結(jié)果見圖6.

4.2 萊州M4.6震源機(jī)制

使用CAP方法對2013年11月23日萊州M4.6地震進(jìn)行反演（鄭建常等，2015），同樣進(jìn)行CAP反演情況的自助抽樣統(tǒng)計(jì)分析.選用震中距在270km以內(nèi)的采樣率為100Hz的22個(gè)三分向?qū)掝l帶臺的觀測波形組成原始數(shù)據(jù)集，采用全樣本隨機(jī)抽取方法，自助抽樣反演1000次，圖7給出了反演得到的震源機(jī)制的斷層節(jié)面解和P、T軸在震源球上的分布情況.結(jié)果顯示，萊州地震的自助抽樣結(jié)果同樣很好地顯示出了反演得到震源機(jī)制解的誤差范圍，相對于乳山地震，出現(xiàn)了極個(gè)別反演過程的結(jié)果偏離較大的情況.

使用粒子群非線性優(yōu)化方法，搜索與自助抽樣結(jié)果的旋轉(zhuǎn)角最小的解.結(jié)果顯示最優(yōu)解為，節(jié)面A：走向236.9°，傾角76.2°，滑動(dòng)角－169.3°；節(jié)面B：走向144.3°，傾角79.6°，滑動(dòng)角－14.0°；與所有自助抽樣解的平均偏轉(zhuǎn)角17.4°（圖8），以其與所有解Kagan角的2倍標(biāo)準(zhǔn)差（圖8紅色虛線所示）為誤差范圍，結(jié)果顯示震源機(jī)制解的不確定性為23.7°.

從自助抽樣得到的所有機(jī)制解在震源球上的分布情況（圖7）可以直觀地看出，斷層節(jié)面線尤其是北西向節(jié)面呈現(xiàn)出兩組集中.由于我們定義的粒子群優(yōu)化的目標(biāo)函數(shù)是搜索與所有自助抽樣解的空間旋轉(zhuǎn)最小，因此從圖7可以看出，最優(yōu)解的節(jié)面位置處于其中一組的邊緣位置，在此情況下，對自助抽樣結(jié)果進(jìn)行聚類分析是有意義的.

圖3 2014年1月7日乳山M.3地震最優(yōu)解的理論波形（紅）與觀測黑波形圖下方第一行數(shù)字為各段理論地震波形相對實(shí)際觀測波形的移動(dòng)時(shí)間，正值表示理論波形相對觀測波形超前.第二行數(shù)字為理論波形與觀測波形的相關(guān)系數(shù)（百分比）.波形圖左側(cè)字母為臺站，其下數(shù)字分別為臺站震中距（km）和方位角（°）.圖左側(cè)的震源球上紅色區(qū)域代表壓縮區(qū)，白色代表拉張區(qū)，震源球采用下半球投影.震源球上標(biāo)注的“＋”和“－”表示反演使用臺站的P波初動(dòng).Fig.3 Comparison between synthetics（red）and observed（black）seismograms of Jan.7，2014Rushan M4.3event The numbers on the lower left side of each seismogram are the time shifts（upper）and cross－correlation coefficient in percent（lower）.Positive time shifts mean that the observed data have been delayed.The letters on the left side are stations，the numbers below it are epicentral distance（in km）and azimuth（degree）.The red color in beach－ball denotes compression area，while white is extension.The‘＋’and‘－’signs on beach－ball indicate polarities on inversion used stations.Lower hemisphere projection is used.

圖4 自助抽樣得到的1000次乳山地震震源機(jī)制解及粒子群最優(yōu)解（下半球極射投影）震源球上黑色細(xì)線條表示自助抽樣結(jié)果的斷層節(jié)面線，紅色線條表示粒子群優(yōu)化解的節(jié)面線.Fig.4 1000focal mechanisms of Rushan M4.3event retrieved by a bootstrapping process，all nodal lines（black）and P，Taxes（blue and red points，respectively）are plotted on one beach－ball.The red lines on the beach－ball show the optimized solution given by a Particle Swarm Optimization method.Its corresponding P，Taxes are also displayed on the beach－ball （yellow and green point，respectively）.Lower hemisphere projection is used.

圖5 乳山地震自助抽樣結(jié)果與粒子群優(yōu)化解的Kagan角分布Fig.5 Kagan angles of bootstrap results to the PSO solution for Rushan event

以兩個(gè)震源機(jī)制解之間的最小空間旋轉(zhuǎn)角（即Kagan角）為距離的定義，對1000次自助抽樣結(jié)果進(jìn)行聚類分析，圖9給出了聚類譜系圖（由于完整的聚類樹過于密集和龐大，因此我們只顯示了Kagan＞7°的部分），以50°為閾值，可以將結(jié)果分為5類.圖10給出了聚類分析的結(jié)果，屬于Ⅰ類的數(shù)據(jù)占32.1%，Ⅱ類66.6%，其余三類數(shù)據(jù)合計(jì)僅有1.3%.由圖10可以看出，其余三類的斷層節(jié)面線和P、T軸位置明顯偏離集中區(qū)且機(jī)制解類型與絕大部分結(jié)果（走滑型）不一致，是典型的孤立解.孤立解（或錯(cuò)誤解）的出現(xiàn)，可能說明我們使用的數(shù)據(jù)中個(gè)別臺站（或分向）存在較大干擾.使用俞春泉等（2009）的方法求取了四類解的聚類中心，其中I類解的聚類中心的斷層面參數(shù)（設(shè)北東向節(jié)面為發(fā)震斷層面）為：走向231.6°，傾角88.7°，滑動(dòng)角－168.2°；Ⅱ類解的聚類中心為：走向238.5°，傾角74.1°，滑動(dòng)角－164.8°.

圖6 乳山地震自助抽樣結(jié)果的震源球概率密度分布俯視圖（未進(jìn)行極射投影）色標(biāo)中正值表示T軸的概率密度分布，負(fù)值表示P軸的概率密度分布.Fig.6 Probability density distribution of solutions on beach－ball（top view of lower hemisphere，without projection）Positive values on the color scale（corresponding to the red area on beach－ball）indicate probability of T axis，while negative values（corresponding to blue area）mean probability of Paxis.

圖7 自助抽樣得到的萊州地震震源機(jī)制解及粒子群最優(yōu)解（下半球極射投影）黑色節(jié)面線為自助抽樣得到的震源機(jī)制解；紅色節(jié)面線為粒子群最優(yōu)解.Fig.7 1000Focal mechanisms of Laizhou M4.6event retrieved by a bootstrapping process，all nodal lines（black）and P，Taxes（blue and red points，respectively）are plotted on one beach－ball.The red lines on the beachball show the optimized solution given by a PSO method.Its corresponding P，Taxes are also displayed on the beach－ball（cyan and yellow point，respectively）.Lower hemisphere projection is used.

圖8 萊州地震自助抽樣結(jié)果與粒子群優(yōu)化解的Kagan角分布Fig.8 Kagan angles of bootstrap results to the PSO solution for Laizhou event

圖9 萊州地震自助抽樣結(jié)果的聚類譜系圖Fig.9 Dendrogram plot of the hierarchical binary cluster tree for Laizhou event

圖10 萊州地震自助抽樣結(jié)果的聚類分析（a）斷層節(jié)面線的分類顯示；（b）機(jī)制解P（＋）、T（⊙）軸位置的分類顯示.斷層節(jié)面線和P、T軸顏色表示分類，與圖9分類顏色一致.Fig.10 Clustering results of focal mechanisms from a bootstrap process for Laizhou event The different colors of nodal lines and P （“＋”sign in right panel），T （“⊙”in right panel）axes denote different classes，which are corresponding to colors shown in Fig.9

將自助抽樣結(jié)果中的P、T軸投影到震源球上，對其進(jìn)行概率密度統(tǒng)計(jì)分析，結(jié)果見圖11.可見，同聚類分析的結(jié)果一致，P軸位置的概率密度在震源球上出現(xiàn)了兩個(gè)極值區(qū)，分別對應(yīng)Ⅰ類和Ⅱ類兩個(gè)聚類中心.

5 討論與結(jié)論

基于山東省寬頻帶數(shù)字地震波形資料，本文首先使用CAP方法反演了近期膠東半島地區(qū)發(fā)生的兩次顯著中等地震活動(dòng)的震源機(jī)制，討論了如何合理地估計(jì)CAP方法反演震源機(jī)制的誤差范圍以及如何確定優(yōu)化解的問題.我們首先使用自助抽樣方法，對原始數(shù)據(jù)進(jìn)行等概率隨機(jī)抽樣，多次重復(fù)波形反演過程，排除了人為選擇數(shù)據(jù)的干擾，得到大樣本量的震源機(jī)制解數(shù)據(jù)；在此基礎(chǔ)上，我們

（1）使用了粒子群優(yōu)化算法從中搜索與這些機(jī)制解空間偏轉(zhuǎn)角最小的解當(dāng)作優(yōu)化解，以Kagan角的二倍標(biāo)準(zhǔn)差作為反演結(jié)果的不確定性范圍，結(jié)果顯示：乳山地震的粒子群優(yōu)化解為：走向298.5°，傾角64.3°，滑動(dòng)角0.3°，不確定性為±6.4°；萊州地震的優(yōu)化解為：走向236.9°、傾角76.2°、滑動(dòng)角－169.3°，不確定性為±23.7°.

（2）對自助抽樣結(jié)果進(jìn)行聚類分析，其中：乳山地震結(jié)果的聚類中心與粒子群優(yōu)化解基本一致；萊州地震結(jié)果存在多個(gè)聚類，排除孤立解后，有兩個(gè)聚類中心，其對應(yīng)兩類數(shù)據(jù)合計(jì)占結(jié)果的98.7%，說明此次地震的真實(shí)解在這兩類數(shù)據(jù)范圍內(nèi).

圖11 萊州地震自助抽樣結(jié)果的震源球概率密度分布（a）下半球俯視圖（未進(jìn)行極射投影）；（b）東南方向，45°三維側(cè)視圖.色標(biāo)說明同圖6.Fig.11 Probability density distribution of solutions on beach－ball（a）Top view of lower hemisphere，without any projection；（b）Side view of whole beach－ball，from south－east 45°position.Positive values on the color scale（corresponding to the red area on beach－ball）indicate probability of Taxis，while negative values（corresponding to blue area）mean probability of Paxis.

（3）將自助抽樣結(jié)果中的P、T軸投影到震源球，對其進(jìn)行概率密度統(tǒng)計(jì)，給出了機(jī)制解在震源球上的概率密度分布圖.

本文方法不單可以得到更準(zhǔn)確的震源機(jī)制優(yōu)化解、給出科學(xué)合理的誤差估計(jì)，而且可以有效地排除孤立解和錯(cuò)誤解，克服存在較大干擾臺站的數(shù)據(jù)的影響，因此在震后應(yīng)急的震源機(jī)制自動(dòng)化求解中可以發(fā)揮作用.在強(qiáng)大計(jì)算能力的支持下，無須人工干預(yù)即可得到準(zhǔn)確可靠的震源機(jī)制結(jié)果，從而為震害評估、趨勢分析等提供重要的科學(xué)依據(jù).

自助抽樣結(jié)果顯示乳山地震震源機(jī)制解的不確定性要小于萊州地震，筆者推測可能有萊州地震使用的臺站中個(gè)別臺的干擾較大的原因，另外也無法排除萊州地震的震源破裂過程可能更加復(fù)雜的可能.CAP方法中用于計(jì)算格林函數(shù)的F－K方法使用狄拉克－Delta函數(shù)作為震源時(shí)間函數(shù)（Zhu and Rivera，2002），對于小震級的事件該簡化方案更為適用，萊州地震（M4.6）與乳山地震（M4.3）震級相差不大，但萊州地震的自助抽樣結(jié)果出現(xiàn)了兩個(gè)概率較高的聚類中心，在使用大部分相同臺站的情況下，這可能意味著萊州地震的震源破裂隨時(shí)間的變化可能與狄拉克－Delta函數(shù)存在一定的偏離.Rodríguez－Lozoya等（2008）的研究顯示，區(qū)域中等地震也可能有復(fù)雜的震源破裂過程，在該問題上的深入研究需要更進(jìn)一步的工作.

致謝感謝兩位匿名審稿專家提出的寶貴意見.聚類分析中使用了俞春泉等（2009）提供的部分開放代碼，粒子群搜索使用了S Chen給出的Matlab軟件包，在此一并表示感謝！

Chang S J，Baag C E.2006.Crustal structure in Southern Korea from joint analysis of regional broadband waveforms and travel times.Bull.Seism.Soc.Amer.，96（3）：856－870.

Diao G L，Yu L M，Li Q Z.1992.Hierarchical clustering analysis of the focal mechanism solution－Taking the Haicheng Earthquake Sequences for example.EarthquakeResearchinChina（in Chinese），8（3）：86－92.

Godano M，Regnier M，Deschamps A，et al.2009.Focal mechanisms from sparse observations by nonlinear inversion of amplitudes：method and tests on synthetic and real data.Bull.Seism.Soc.Amer.，99（4）：2243－2264.

Han L B，Jiang C S，Bao F.2010.Source parameter determination of 2010TaikangMS4.6earthquake sequences.ChineseJ.Geophys.（in Chinese），55（9）：2973－2981，doi：10.6038／j.issn.0001－5733.2012.09.016.

Han L B，Jiang C S.2012.Focal mechanism inversion of 8Jun 2011 ToksunMS5.3earthquake.ActaSeismologicaSinica（in Chinese），34（3）：415－422.

Huang J P，Ni S D，F(xiàn)u R S，et al.2009.Source mechanism of the 2006MW5.1Wen′an earthquake determined from a joint inversion of local and teleseismic broadband waveform data.ChineseJ.Geophys.（in Chinese），52（1）：120－130.

Jost M L，Hermann R B.1989.A student′s guide to and review of moment tensors.Seism.Res.Lett.，60（2）：37－57.

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

Li D M，Zheng J C.2014.Relocation analysis of Oct.1，2013 Rushan swarms.QiluEarthquakeSciences（in Chinese），10（1）：26－28.

Long F，Zhang Y J，Wen X Z，et al.2010.Focal mechanism solutions ofML≥4.0events in theMS6.1Panzhihua－Huili earthquake sequence of Aug 30，2008.ChineseJ.Geophys.（in Chinese），53（12）：2852－2860，doi：10.3969／j.issn.0001－5733.2010.12.008.

LüJ，Zheng Y，Ni S D，et al.2008.Focal mechanisms and seismogenic structures of theMS5.7andMS4.8Jiujiang－Ruichang earthquakes of Nov.26，2005.ChineseJ.Geophys.（in Chinese），51（1）：158－164，doi：10.3321／j.issn：0001－5733.2008.01.020.

Ree J－H，Cho M，Kwon S－T，et al.1996.Possible eastward extension of Chinese collision belt in South Korea：The Imjingang belt.Geology，24（12）：1071－1074.

Rodríguez－Lozoya H E，Quintanar L，Ortega R，et al.2008.Rupture process of four medium－sized earthquakes that occurred in the Gulf of California.J.Geophys.Res.，113（B10301），doi：10.1029／2007JB005323.

Xu X T，Xu Z H，Zhang D N.1995.A probabilistic grid test method of determining earthquake focal mechanisms using P－wave onset polarity data.SeismologicalandGeomagnetic ObservationandResearch（in Chinese），16（4）：34－42.

Yu C Q，Tao K，Cui X F，et al.2009.P－wave first－motion focal mechanism solutions and their quality evaluation.ChineseJ.Geophys.（in Chinese），52（5）：1402－1411，doi：10.3969／j.issn.0001－5733.2009.05.030.

Zhao L S，Helmberger D V.1994.Source estimation from broadband regional seismograms.Bull.Seism.Soc.Amer.，84（1）：91－104.

Zheng J C，Li D M，Wang P，et al.2015.Focal mechanisms and seismic tectonic features of the 2013LaizhouM4.6earthquake sequence.SeismologyandGeology（in Chinese），in press.

Zheng J C，Chen Y T.2012.Stability of sparse station data inversion for deviatoric moment tensor solution of regional earthquakes.ActaSeismologicaSinica（in Chinese），34（1）：31－43.

Zheng Y，Ma H S，LüJ，et al.2009.Source mechanism of strong aftershocks（Ms≥5.6）of the 2008／05／12Wenchuan earthquake and the implication for seismotectonics.ScienceinChinaSeriesD：EarthSciences，52（6）：739－753.

Zhu L P，Helmberger D V.1996.Advancement in source estimation techniques using broadband regional seismograms.Bull.Seism.Soc.Amer.，86（5）：1634－1641.

Zhu L P，Rivera L A.2002.A note on the dynamic and static displacements from a point source in multilayered media.Geophys.J.Int.，148（3）：619－627.

附中文參考文獻(xiàn)

刁桂苓，于利民，李欽祖.1992.震源機(jī)制解的系統(tǒng)聚類分析——以海城地震序列為例.中國地震，8（3）：86－92.

韓立波，蔣長勝，包豐.2012.2010年河南太康MS4.6地震序列震源參數(shù)的精確確定.地球物理學(xué)報(bào)，55（9）：2973－2981，doi：10.6038／j.issn.0001－5733.2012.09.016.

韓立波，蔣長勝.2012.2011年6月8日新疆托克遜MS5.3地震震源機(jī)制解反演.地震學(xué)報(bào)，34（3）：415－422.

黃建平，倪四道，傅容珊等.2009.綜合近震及遠(yuǎn)震波形反演2006文安地震（MW5.1）的震源機(jī)制解.地球物理學(xué)報(bào)，52（1）：120－130.

李冬梅，鄭建常.2014.2013年10月1日乳山震群重新定位結(jié)果分析.齊魯?shù)卣鹂茖W(xué)，10（1）：26－28.

龍鋒，張永久，聞學(xué)澤等.2010.2008年8月30日攀枝花—會(huì)理6.1級地震序列ML≥4.0事件的震源機(jī)制解.地球物理學(xué)報(bào)，53（12）：2852－2860，doi：10.3969／j.issn.0001－5733.2010.12.008.

呂堅(jiān)，鄭勇，倪四道等.2008.2005年11月26日九江－瑞昌MS5.7，MS4.8地震的震源機(jī)制解與發(fā)震構(gòu)造研究.地球物理學(xué)報(bào)，51（1）：158－164，doi：10.3321／j.issn：0001－5733.2008.01.020.

許向彤，許忠淮，張東寧.1995.求震源機(jī)制P波初動(dòng)解的格點(diǎn)嘗試概率法.地震地磁觀測與研究，16（4）：34－42.

俞春泉，陶開，崔效鋒等.2009.用格點(diǎn)嘗試法求解P波初動(dòng)震源機(jī)制解及解的質(zhì)量評價(jià).地球物理學(xué)報(bào)，52（5）：1402－1411，doi：10.3969／j.issn.0001－5733.2009.05.030.

鄭建常，王鵬，李冬梅等.2015.2013年萊州M4.6地震序列震源機(jī)制與發(fā)震構(gòu)造初探.地震地質(zhì)，待刊.

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

鄭勇，馬宏生，呂堅(jiān)等.2009.汶川地震強(qiáng)余震（Ms≥5.6）的震源機(jī)制解及其與發(fā)震構(gòu)造的關(guān)系.中國科學(xué)D輯：地球科學(xué)，39（4）：413－426.