裴瑋來，周仕勇

北京大學(xué)地球與空間科學(xué)學(xué)院，北京 100871

0 引言

對地震波形記錄中的地震事件進(jìn)行準(zhǔn)確檢測并拾取地震震相的到時與初動極性，是地震資料分析的基本內(nèi)容之一.目前在地震信號采集之中，最常用且被廣泛認(rèn)同的單臺地震事件檢測與震相到時讀取方法主要有以下幾類：

(1)基于長短窗信號比的震相到時檢測方法.通過計算長短兩個時窗的信號的能量比值來判定信號的初至到時(Stevenson,1976).短時能量(STA)顯著高于長時能量(LTA)時，設(shè)定閾值判定地震信號到達(dá)的時間.長短窗算法是最經(jīng)典的地震波信號拾取算法之一，常被拿來與其他方法進(jìn)行對比.為了增加算法的靈敏度和精確度，許多學(xué)者先后采用不同的特征函數(shù)替代地震序列值作為輸入(Allen,1978,1982;Baer and Kradolfer,1987)，來提高震相檢測的精確性.

(2)基于熵的概念進(jìn)行震相檢測的方法.較有代表性的方法包括以赤池信息準(zhǔn)則構(gòu)建目標(biāo)函數(shù)的震相拾取算法以及互信息方法兩類.前者通過熵的概念權(quán)衡模型復(fù)雜度和數(shù)據(jù)擬合優(yōu)良性，構(gòu)建目標(biāo)函數(shù)(AIC函數(shù))，并尋找給定時窗中AIC函數(shù)的最小值的位置來決定時窗內(nèi)地震信號的到時(Maeda,1985).后者通過計算隨機(jī)變量之間的互信息值，或以其他信息論概念為基礎(chǔ)構(gòu)建的目標(biāo)函數(shù)來判定地震震相到時.互信息的概念如下(Shannon,1948)：

(1)

其中p(x,y)表征隨機(jī)變量X與Y的取值的聯(lián)合概率密度分布；p(x)與p(y)表征隨機(jī)變量X和Y的邊緣概率密度分布；Pmi(x;y)表征點互信息，衡量兩個隨機(jī)事件之間的相關(guān)性.一般而言，基于熵的概念進(jìn)行震相檢測的算法在選取的時窗較為合適時能夠給出非常精確的結(jié)果.

(3)基于高階統(tǒng)計量如偏度、峰度等參數(shù)開發(fā)的震相拾取方法，通過計算地震波形的峰度和偏斜度的極值點，尋找地震波震相到時(Nakamula et al.,2007;Ross et al.,2016;Saragiotis et al.,2002).其中偏斜度描述隨機(jī)變量相對于均值的不對稱程度，而峰度描述隨機(jī)變量分布曲線的尖峭程度.

(4)以機(jī)器學(xué)習(xí)為代表的利用人工智能搭建神經(jīng)網(wǎng)絡(luò)進(jìn)行地震波震相拾取的方法.通過模擬人類的學(xué)習(xí)行為來訓(xùn)練神經(jīng)網(wǎng)絡(luò)，從而達(dá)到完成特定任務(wù)的目的(Wang and Teng,1997;Zhao and Takano,1999;Hara et al.,2019;Ross et al.,2018;Zhou et al.,2019;Zhu and Beroza,2019;Wong et al.,2021).基于人工智能的方法雖然能給出較高的準(zhǔn)確率，但需要提供大量標(biāo)記樣本數(shù)據(jù)來訓(xùn)練模型，因此在不同地區(qū)和數(shù)據(jù)間的遷移和應(yīng)用存在一定程度的困難且較為耗時.

(5)各類雜交方法與其他方法.Zhang等(2003)對單分量波形取一系列滑動時窗進(jìn)行離散小波變換,并在每個時窗內(nèi)應(yīng)用AIC準(zhǔn)則自動拾取判斷是否有P波震相；Bai和Kennett(2001)結(jié)合了包括長短窗、平均-短期-長期比值、互相關(guān)、PCA在內(nèi)的多種不同方法進(jìn)行震相拾取；蔣策等(2018)結(jié)合AIC方法與長短窗方法的結(jié)果自動選擇合適的值作為震相到時；趙大鵬等(2012)利用直達(dá)P波信號的峰度和Kurtosis-AIC兩個指標(biāo)，采用峰度觸發(fā)時間窗、AIC挑震相的思想來達(dá)到區(qū)域地震事件實時檢測和直達(dá)P波到時精細(xì)拾取的目的.

以上的幾種方法主要是針對單臺、單分量記錄通過振幅或者能量的變化來檢測P波震相.雖然它們基本實現(xiàn)了快速、自動化的P波震相到時檢測，但是以上的方法大多難以估計拾取到時的不確定度，給定量估計初至到時的誤差帶來了困難.

作為小震震源機(jī)制解反演的重要資料，P波初動的極性是地震資料分析中的又一項基礎(chǔ)性工作.P波初動極性的判定在一定程度上基于P波初至到時的準(zhǔn)確檢測.因此現(xiàn)有的P波初動極性的自動判定方法大多數(shù)使用了上文提到的一種或多種P波震相到時讀取算子.比如Pugh等(2016)依賴上文中提到的長短窗函數(shù)法.Ross等(2018)依賴人工智能拾取P波震相到時的方法.Chen和Holland(2016)依賴基于峰度、偏度與赤池信息準(zhǔn)則的聯(lián)合到時拾取算法.當(dāng)P波到時被準(zhǔn)確挑取時，P波初動的判定本身沒有過多的爭議和可開發(fā)的空間，多數(shù)情況下我們都能得到準(zhǔn)確的判定結(jié)果，但是當(dāng)其所依賴的P波初至到時拾取算子不夠準(zhǔn)確或者出現(xiàn)錯誤時，微小的到時拾取誤差就有可能導(dǎo)致不一樣的結(jié)果.雖然大多數(shù)前人提出的方法都經(jīng)過了檢驗并取得了不錯的結(jié)果，但在將其應(yīng)用到具體問題中時，還是有一定的犯錯的幾率.在目前主流的P波初動極性判定方法中，僅有Pugh等(2016)具備對初動極性誤差進(jìn)行定量分析的能力.然而由于該方法采用的初至拾取算子的誤差定量估計不夠嚴(yán)謹(jǐn)，導(dǎo)致其對P波初動極性的錯誤率分析不夠嚴(yán)謹(jǐn)可靠.

現(xiàn)有方法在P波初至與P波初動極性拾取方面存在的問題，歸根結(jié)底是因為初至到時的不確定度難以估計.因此，一種數(shù)學(xué)上更加嚴(yán)密、能夠考慮結(jié)果的不確定度，且到時拾取更加精確的自動化到時拾取算法是有必要且意義重大的.本文采用的新的自動化到時拾取算法，POI(Probability of arrival time and polarity based on Order statistics and Information theory)算法，是基于最大順序統(tǒng)計量和互信息理論構(gòu)建的.該方法能夠定量計算給定時窗中震相到時的概率密度分布，使P波到時與初動極性的不確定度得到約束.該方法并在自動化反演震源參數(shù)等方面可以得到很好的應(yīng)用.

1 方法原理

1.1 基于最大順序統(tǒng)計量的P波到時與極性讀取方法

POI算法主要由兩步構(gòu)成.第一步是對給定時窗中的地震波數(shù)據(jù)，使用互信息理論，構(gòu)造合適的特征函數(shù)來拾取地震波的初至?xí)r刻.地震波到達(dá)往往導(dǎo)致波形記錄的振幅突然變大.因此通過構(gòu)造一個體現(xiàn)振幅變化的隨機(jī)變量，觀察其與表征地震波到時的隨機(jī)變量的相關(guān)性，就可以準(zhǔn)確地拾取出地震波初至到時.假設(shè)地震波形記錄Ai，定義一個振幅值ε作為閾值(圖1).其中閾值ε是運用互信息算法對波形的振幅值A(chǔ)i進(jìn)行區(qū)分從而作為確認(rèn)初動到時的參數(shù).其值表征噪音和信號的振幅的絕對值大小的分界值.超過該閾值的數(shù)據(jù)點或數(shù)據(jù)段被認(rèn)為其來自地震信號，而小于該閾值的數(shù)據(jù)點或數(shù)據(jù)段來自背景噪音，即有

圖1 互信息到時判定利用互信息構(gòu)建的目標(biāo)函數(shù)Z(t,ε)，取目標(biāo)函數(shù)值最大的時刻tinit作為到時.定義式中不同的事件用不同顏色區(qū)分出來.此處閾值約為300.Fig.1 Arrival time picking using pointwise mutual informationWhen amplitude threshold is given,arrival time is determined by maximizing objective function Z(t,ε).Here threshold is about 300.

(2)

由互信息理論，與表征振幅值變化的隨機(jī)變量Y相關(guān)性最高的到時隨機(jī)變量X，標(biāo)識了震相到時t.在閾值ε確定時，到時t的取值完全由下列目標(biāo)函數(shù)決定：

(3)

其中m和n取0或1,Pm n代表聯(lián)合概率密度分布，其涵義為

P(X=x0,Y=y0)=P00(t,ε)=Pr{i≤t∧|Ai|≤ε},

P(X=x0,Y=y1)=P01(t,ε)=Pr{i≤t∧|Ai|>ε},

P(X=x1,Y=y0)=P10(t,ε)=Pr{i>t∧|Ai|≤ε},

P(X=x1,Y=y1)=P11(t,ε)=Pr{i>t∧|Ai|>ε}.

(4)

用tinit指代使互信息目標(biāo)函數(shù)Z(t,ε)最大的t的取值，tinit即為閾值已知時認(rèn)定的真實地震波到時：

(5)

若已知ε的概率密度分布，則可直接推出到時tinit的概率密度分布：

(6)

上式中I為邏輯函數(shù)，其括號內(nèi)條件成立則取值為1，否則取值為0.

方法的第二步是通過引入最大順序統(tǒng)計量的概念來計算合理的參數(shù)閾值ε的概率密度分布，從而使得到時tinit的概率密度分布可以直接通過等式(6)計算出來.

理論上閾值的取值范圍可以是從0到正無窮.給出不同的閾值，或者閾值的概率密度分布，所推導(dǎo)出的到時，或到時的概率密度分布往往是不同的.我們希望給出的閾值的概率密度分布能夠使得最終推導(dǎo)出的到時的概率密度分布接近真實的到時.為此需要推導(dǎo)能夠計算出真實到時tinit的閾值ε所滿足的條件.

我們提出最大順序統(tǒng)計量準(zhǔn)則：閾值ε應(yīng)當(dāng)與其對應(yīng)到時tinit前置的背景噪音Asamp，在正態(tài)分布假設(shè)下的最大順序統(tǒng)計量相近，即

max(|Asamp|)≈ε,(7)

上式左端表示在正態(tài)分布假設(shè)下，前置噪音樣本Asamp的最大順序統(tǒng)計量.其中前置噪音樣本定義為

Asamp=ξ(ε)={Ai:|Ai|≤ε,i≤η(ε)}.

(8)

這里引入正態(tài)分布假設(shè)是因為真實的噪音往往是服從正態(tài)分布的.因此真正的閾值所對應(yīng)的前置噪音應(yīng)當(dāng)也是大致服從正態(tài)分布.在正態(tài)分布假設(shè)下，前置的背景噪聲可以認(rèn)為是服從如下總體的隨機(jī)采樣：

Asamp～N(0,σ2),(9)

其中σ代表總體的噪聲水平，可以用前置噪音樣本估計：

(10)

此時前置噪音的最大順序統(tǒng)計量的概率密度分布可以由下式給出：

(11)

上式中erf代表誤差函數(shù)，card函數(shù)為返回元素個數(shù)的函數(shù).考慮所有的閾值取值的情況，若閾值取值的概率密度分布已知，則前置噪音樣本的最大順序統(tǒng)計量的概率密度即為邊緣分布：

=x|ε)×p(ε)dε,(12)

式(12)代表參數(shù)ε邊緣化的過程，等式左邊為將閾值邊緣化后最大順序統(tǒng)計量的概率密度分布.根據(jù)最大順序統(tǒng)計量準(zhǔn)則，合理的閾值的概率密度分布應(yīng)該使其推出的最大順序統(tǒng)計量的總體的概率密度分布接近自身，即

p(ε=x)≈p(max(|Asamp|)=x),(13)

因此有

式(14)中，待求的概率密度分布p(ε)為方程中的唯一未知量.因此可以通過解方程的形式將p(ε)解出，再由式(6)即可求得初至到時的概率密度分布.

初至波峰A定義為初至到時后第一個極值點的振幅值，其概率密度分布可以由初至到時的概率密度直接求得.當(dāng)求解的閾值ε高于所有振幅絕對值|Ai|時，在所選時窗內(nèi)無法定義初至到時與初至波峰A，初動的極性被判定為未知.除此之外，初至到時與初至波峰A的概率密度分布存在且可求.Pugh等(2016)中提出了完備的的貝葉斯極性判定公式，公式支持正負(fù)兩種極性的計算.我們沿用其中的公式，計算初動的極性向上(Y=1)或向下(Y=-1)的概率：

(15)

將Y=1或Y=-1代入公式(15)中，所求得的概率即為初動極性為正或為負(fù)的概率.

1.2 基于P波初動的震源機(jī)制解與區(qū)域應(yīng)力場反演流程

圖2顯示了本文工作設(shè)計的震源機(jī)制解與區(qū)域應(yīng)力場反演流程.這套反演流程采用本文提出的POI P波初動極性判定方法，并綜合了已有的HASH算法和SATSI算法.流程包括初動極性拾取、SH/P波振幅比計算、震源機(jī)制反演與應(yīng)力場反演4個主要步驟.由于所有采用的方法本身都是能夠?qū)崿F(xiàn)自動化的，因此這套反演流程不需要人工的介入.

圖2 本研究采用的自動化反演流程圖圖中的每一項計算過程都實現(xiàn)了自動化.其中POI極性拾取是本文采用的新的初動極性拾取方法.Fig.2 Automatic inversion process used in this studyEach step of inversion is independent of human involvement.POI method is the new polarity determination method adopted in the inversion process.

反演流程的第一個步驟是通過原始地震目錄截取波形數(shù)據(jù).對于給定的地震目錄和連續(xù)波形記錄，選取合適震中距的臺站記錄(120 km以內(nèi))，并通過區(qū)域一維速度結(jié)構(gòu)估算P波震相與S波震相的理論到時，而后截取P波到時前10 s至后30 s的時間窗.在經(jīng)過1～15 Hz的濾波后，選取P波初至到時前后各3 s的速度記錄數(shù)據(jù)作為P波震相拾取時窗，為后續(xù)POI極性拾取做準(zhǔn)備.分別選用P波到時前5 s到后2 s的時間窗，和S波理論到時前后共7 s的時間窗，分別計算P波與S波的位移記錄的最大振幅值，為計算振幅比做準(zhǔn)備.

反演流程的第二個步驟是通過波形數(shù)據(jù)分別同步計算P波初動極性和振幅比.P波初動極性的計算采用1.1節(jié)中提出的POI算法，設(shè)置質(zhì)量控制參數(shù)(概率閾值，如：0.99)，保留所有概率閾值之上讀取的極性數(shù)據(jù).反演過程中使用的振幅比資料則需要在初步振幅比值的基礎(chǔ)上進(jìn)行了自由界面校正和臺站校正.

由于放置在地表的地震儀器記錄到的位移記錄為入射波與反射波的疊加，而反演震源機(jī)制解的過程中使用的振幅比資料事實上只有入射波的部分，因此需要對初步振幅比資料進(jìn)行自由界面校正.在眾多種類的振幅比資料中，SH/P振幅比資料穩(wěn)定性較好，在自由界面校正的過程中簡單易算，不需要考慮P與SV的互相轉(zhuǎn)化的問題.而SV/P以及S/P在計算過程中由于SV在自由界面表面反射的過程中存在大于臨界角的情況，需要對入射角度進(jìn)行嚴(yán)格分類計算.因此本研究的振幅比資料主要采用SH/P振幅比資料.盡管吳大銘等(1989)通過研究不同入射角范圍的反透射系數(shù)，認(rèn)為總體而言針對P波的自由界面校正是不必要的，考慮到本研究中入射角的展布范圍，實際反演中還是對振幅比值進(jìn)行了自由界面校正.自由界面校正的公式參考Aki和Richards(2002).由于定位結(jié)果中包含一定的誤差，而這些誤差可能會對振幅比的自由界面校正量帶來影響，我們因此根據(jù)hypoinverse程序給出的深度誤差計算了校正量的上下界.對于深度定位誤差過大導(dǎo)致的振幅比的自由界面校正系數(shù)難以確定的情況，舍棄該振幅比值，僅保留在校正量誤差在以2倍以內(nèi)的數(shù)據(jù)，進(jìn)行了自由界面校正.

臺站校正是指通過計算理論輻射花樣，與實際臺站記錄到的振幅比資料做對比，在一定程度上修正淺層地殼，特別是臺站下方的場地效應(yīng)以及衰減效應(yīng).臺站校正整體上是一個經(jīng)驗性的校正.本研究中計算了所有臺站的臺站校正，采用的處理流程遵循Hardebeck和Shearer(2003)的處理流程.將經(jīng)過自由界面校正與臺站校正后得到的振幅比資料用于震源機(jī)制解反演.

反演流程的第三個步驟是通過HASH算法，使用第二個步驟計算的初動極性和振幅比數(shù)據(jù)對所有地震的震源機(jī)制進(jìn)行反演和篩選.HASH程序的源代碼是計算S波與P波的振幅比.算法考慮了地震深度、一維速度結(jié)構(gòu)、P波初動極性等數(shù)據(jù)資料的不確定度，對同一地震事件采用不同的地震深度與不同的速度結(jié)構(gòu)組合，反演地震震源機(jī)制，并在所有結(jié)果中進(jìn)行聚類，然后返回可能的震源機(jī)制解.我們將HASH源碼中計算S/P的目標(biāo)函數(shù)改為計算SH/P的目標(biāo)函數(shù)，其他流程不變.

反演流程的最后一步是通過SATSI算法，使用第三個步驟計算的震源機(jī)制解作為數(shù)據(jù)資料，反演區(qū)域應(yīng)力場的應(yīng)力主軸方向以及三個主應(yīng)力的大小關(guān)系.SATSI算法是基于法向應(yīng)力反演，默認(rèn)斷層滑動沿斷層面上的剪切應(yīng)力方向發(fā)生的應(yīng)力張量求解技術(shù).SATSI算法可以在保持連續(xù)性的前提下同時求解多個子區(qū)域的應(yīng)力場，從而觀察研究區(qū)域的應(yīng)力主方向的空間分布與時間變化.

2 研究區(qū)域與研究資料

小江斷裂帶位于青藏高原東南緣的川滇菱形塊體東側(cè)，毗鄰華南地塊，走向近南北，全長約400 km，分為北、中、南三段，是川滇地區(qū)的主要活動斷裂之一.多個研究結(jié)果顯示，由于川滇塊體向東南方向遠(yuǎn)離青藏高原，并圍繞東喜馬拉雅構(gòu)造順時針旋轉(zhuǎn)，小江斷裂帶的主要運動形式以左旋走滑運動為主.小江斷裂帶主斷層周邊分布多條次級斷裂，斷層的幾何形態(tài)與運動狀態(tài)均比較復(fù)雜.為了研究小江斷裂帶中北段以及南段的地震活動性及斷層結(jié)構(gòu)，2015年至2019年間，北京大學(xué)周仕勇課題組聯(lián)合中國地震局地球物理所許力生研究員沿小江斷裂帶走向布設(shè)了一個以寬頻帶儀器為主的流動地震臺陣(圖3).臺陣覆蓋小江斷裂帶的主斷層結(jié)構(gòu)，臺站間距平均小于20 km.流動臺陣的全部連續(xù)三分量記錄以及昆明、東川兩個固定臺的連續(xù)波形記錄構(gòu)成了本研究的主要資料來源.

本研究中待反演震源機(jī)制的地震目錄由自動化監(jiān)測算法結(jié)合人工處理得到.利用自動化震相拾取聚類算法對近震的P、S初至震相進(jìn)行初步的拾取和聚類，得到了包含1953個地震事件的初始地震目錄.為了進(jìn)一步提高震相目錄的準(zhǔn)確性和震相數(shù)量，增加求解資料，通過ObsPy(Beyreuther M et al，2010)和crazyseismic軟件(Yu et al.，2017)對初步地震目錄進(jìn)行了人工篩選以及震相增補(bǔ)，去掉誤識別并且補(bǔ)充初步拾取過程中遺漏的P波初至震相.經(jīng)過人工篩選和擴(kuò)充后的震相文件包含約15000個P波震相到時.它們來自1866個地震事件的所有震中距小于120 km的臺站的可識別記錄.平均每一個事件的P波震相記錄數(shù)量大于8.

對于1866個通過篩選的地震事件，采用hypoinverse程序?qū)λ惺录M(jìn)行重定位(Klein，2002).重定位所采用的速度模型參考吳建平等(2006).經(jīng)過重定位后，選取P波記錄數(shù)大于5，到時殘差小于0.5 s，水平定位誤差小于3 km，垂直定位誤差小于5 km的地震事件，共674個符合條件的事件被挑選出來(圖3).空間上看這些地震事件大多數(shù)分布在小江斷裂帶主斷層的附近，少數(shù)事件有叢集特征的現(xiàn)象.

圖3 本研究選取的研究區(qū)域，使用的地震臺與選取的674個地震空間分布圖倒三角代表固定臺站(昆明臺、東川臺)，正三角代表寬頻帶流動地震臺站，圓點代表地震位置.Fig.3 Demonstration of stations and earthquake catalog used in this studyInverted triangles represent permanent stations,i.e.,Dongchuan station and Kunming station.Triangles represent mobile stations.The catalog contains 674 events and they distribute evenly along the main fault.

3 結(jié)果分析

3.1 初動極性判定結(jié)果

為了檢測方法的有效性，以上節(jié)確定的研究區(qū)674個定位精度較高的地震事件為基礎(chǔ)資料，按照1.2節(jié)中的反演流程首先進(jìn)行P波初動的極性判定.研究區(qū)臺陣共獲得了來自這674個地震的6688條地震記錄，將這6688條地震記錄截取的P波到時前3 s至后3 s時窗作為輸入數(shù)據(jù)進(jìn)入POI，開展P波初至到時與極性的概率計算.對于6688個輸入數(shù)據(jù)，計算了每一條的P波初動極性概率以及時窗中的震相到時.圖5至圖8顯示了其中4組計算案例.我們選取不確定度小于0.01的數(shù)據(jù)(初動極性向上的概率大于0.99或初動極性向下的概率大于0.99)作為有效數(shù)據(jù)(impulsive)，用作下一步震源機(jī)制解反演的輸入資料，而將其他數(shù)據(jù)作為可靠性較低(emergent)或者極性未知(unknown)的初動數(shù)據(jù)舍棄.在6688組數(shù)據(jù)中，共有4084條記錄被選中.

另一方面，為了驗證新方法拾取初動極性的準(zhǔn)確性，我們對6688組Z分量數(shù)據(jù)的初動極性進(jìn)行了人工拾取，以便與自動拾取的初動極性做比較.人工拾取時采用傳統(tǒng)的分類方法，將極性分為向上、向下、未知三類，并且對于向上和向下兩類判定出極性的數(shù)據(jù)，按照質(zhì)量又將每一類極性的波形分成了impulsive和emergent兩類.6688組數(shù)據(jù)中，人工辨別極性清晰可讀，質(zhì)量較好的數(shù)據(jù)(impulsive)有2490條，極性可讀但質(zhì)量一般的數(shù)據(jù)(emergent)有849條，而極性未知的數(shù)據(jù)(unknown)有3349條.圖4給出人工讀取與POI自動讀取結(jié)果的比較圖.如果能假設(shè)人工能清晰讀取的2490條極性結(jié)果是正確的，我們可以看到，POI自動讀取的正確率在95%左右.此外POI還能將大量人工難以判斷的極性高概率地識別出來.

圖4 POI算法初動極性判定結(jié)果與人工初動極性判定結(jié)果的比較圖中的每一行代表一類人工拾取初動結(jié)果，而每一行中三種極性的占比該類別數(shù)據(jù)被POI自動算法拾取為向上、未知與向下3種情況的概率.從圖中可以看出Impulsive質(zhì)量的數(shù)據(jù)判定結(jié)果一致性高，對人工判定極性為上、下的數(shù)據(jù)分別有高達(dá)96.5%和92.5%的拾取正確率.對于極性質(zhì)量一般的數(shù)據(jù)判定的正確率平均超過75%.Fig.4 Comparison between automatic determined polarities and manual picking polaritiesEach line represents a category in manual picking.For each category,the probability of up,unknown and down polarities decided by POI algorithm is demonstrated.It′s obvious that our algorithm works well for impulsive polarities,with precision higher than 96.5% and 92.5% for up and down polarities,respectively.Our algorithm also has an accuracy above 75% for emergent polarities.

圖5 初動極性判定案例人工標(biāo)記該條記錄極性向上，質(zhì)量等級為Impulsive.POI算法計算其初動極性100%的概率向上.圖中左上部分展示的是閾值的概率密度分布，右下圖展示的是與之對應(yīng)的初至到時的概率密度分布.主圖中的紅線和綠線展示兩者之間是如何對應(yīng)的.本例中到時幾乎百分之百集中在3.048 s左右.Fig.5 An example of polarity determinationManual polarity picking classification of the seismic trace is up-impulsive.Polarity probability calculated by POI is 100%upward.The top left picture shows the pdf of threshold value.The down right picture shows the pdf of arrival time.Their correspondence is demonstrated in main graph.In this case pdf of arrival time is almost concentrated at 3.048 s.

圖6 初動極性判定案例人工標(biāo)記該條記錄極性向上，質(zhì)量等級為Emergent.POI算法計算其初動極性100%的概率向下.圖中左上部分展示的是閾值的概率密度分布，右下圖展示的是與之對應(yīng)的初至到時的概率密度分布.主圖中的紅線和綠線展示兩者之間是如何對應(yīng)的.本例中到時幾乎百分之百集中在3.271 s左右.Fig.6 An example of polarity determinationManual polarity picking classification of the seismic trace is down-emergent.Polarity probability calculated by POI is 100% downward.The top left picture shows the pdf of threshold value.The down right picture shows the pdf of arrival time.Their correspondence is demonstrated in main graph.In this case pdf of arrival time is almost concentrated at 3.271 s.

圖7 初動極性判定案例人工標(biāo)記該條記錄極性未知.POI算法計算其初動極性100%的概率未知.圖中左上部分展示的是閾值的概率密度分布，右下圖展示的是與之對應(yīng)的初至到時的概率密度分布.主圖中的紅線和綠線展示兩者之間是如何對應(yīng)的.本例中到時幾乎百分之百大于6 s，意味著POI算法判定在本段波形數(shù)據(jù)中不存在初動.因此極性為未知.Fig.7 An example of polarity determinationManual polarity picking classification of the seismic trace is unknown.Polarity probability calculated by POI is 100% unknown.The top left picture shows the pdf of threshold value.The down right picture shows the pdf of arrival time.Their correspondence is demonstrated in main graph.In this case pdf of arrival time is almost 100% above 6 s,which means POI method is almost 100% sure that there is no seismic arrival in this waveform section.Thus the polarity is unknown.

圖8 初動極性判定案例人工標(biāo)記該條記錄極性未知.POI算法計算其初動極性100%的概率向上.圖中左上部分展示的是閾值的概率密度分布，右下圖展示的是與之對應(yīng)的初至到時的概率密度分布.主圖中的紅線和綠線展示兩者之間是如何對應(yīng)的.本例中到時幾乎百分之百集中在2.956 s左右.Fig.8 An example of polarity determinationManual polarity picking classification of the seismic trace is unknown.Polarity probability calculated by POI is 100% upward.The top left picture shows the pdf of threshold value.The down right picture shows the pdf of arrival time.Their correspondence is demonstrated in main graph.In this case pdf of arrival time is almost concentrated at 2.956 s.

3.2 震源機(jī)制解反演結(jié)果

采用HASH程序(Hardebeck and Shearer,2002,2003)，利用POI讀取的P波初動極性與SH/P波的振幅比資料，反演了上述674個地震事件的震源機(jī)制解.

考慮到HASH算法的多解性，需要對所得結(jié)果進(jìn)行質(zhì)量評估.本研究利用多種參數(shù)綜合指標(biāo)判斷結(jié)果的質(zhì)量，并對震源機(jī)制解進(jìn)行分級與篩選.篩選標(biāo)準(zhǔn)包含斷層面角度的不確定度數(shù)(Varest)、該組解占所有可能解的比例(Prob)、符合條件的解的總數(shù)量(Q-Number)、錯誤初動極性百分比(Bad Pol)、臺站程度覆蓋百分比(STA cover)、振幅比誤差(SH/P error)等多個方面(表1).結(jié)合HASH程序中的基礎(chǔ)的質(zhì)量評定標(biāo)準(zhǔn)，只有在兩次質(zhì)量標(biāo)準(zhǔn)評級均在C或C以上，或者第一次評為A，第二次評級為D或D以上的解才會被保留.

表1 震源機(jī)制解質(zhì)量分級表Table 1 Quality grading of focal mechanisms

圖9展示了發(fā)震時刻于2016年9月7日的地震事件經(jīng)過HASH反演算法所得到的震源機(jī)制解.圖中標(biāo)注了該事件在反演過程中使用到的初動波形數(shù)據(jù)以及POI算法計算出的極性向上的概率.從圖中可以看出POI算法對初動極性的拾取是基本正確的，采用自動反演流程得到的震源機(jī)制解也與初動極性本身大致吻合.

圖9 發(fā)生于2016年9月7日07時17分51.47秒的地震事件的P波初動極性拾取情況與震源機(jī)制解該事件的走向、滑動角及傾角分別為330°，75°，-17°.Fig.9 The polarities and focal mechanism of event occurred at 2016-09-07T07:17:51.47The strike,dip and rake angle are 330°，75°，-17°,respectively.

圖10展示了篩選后得到227個質(zhì)量較高解的地震事件的震源機(jī)制結(jié)果.其中走滑類型、正斷層類型、逆沖斷層類型和混合型地震事件的數(shù)量分別為87，84，31，25，它們占總事件數(shù)量的百分比分別為38.5%，37%，13.5%.11%.走滑類型是占比最多的地震事件類型，其次是正斷層型，這與前人計算出的震源機(jī)制解類型大致一致(嚴(yán)川,2015)，也與區(qū)域的構(gòu)造背景和GPS觀測數(shù)據(jù)相一致.

圖10 4種震源機(jī)制解的空間分布情況(a)正斷層事件分布;(b)逆沖事件分布;(c)走滑事件分布;(d)混合類型事件分布.Fig.10 Spatial distribution of four type of focal mechanisms(a)Normal type events;(b)Thrust type events;(c)Strike-slip type events;(d)Hybrid type events.

分類型看，走滑型震源機(jī)制解的斷層面展布仍然以南北和東西方向為主，展示了對小江主斷裂帶的南北走向特征的刻畫.雖然偶爾有滑動方向不嚴(yán)格遵從左旋走滑的地震事件發(fā)生，大多數(shù)沿小江斷裂帶主斷層分布的事件的震源機(jī)制解仍然是左旋走滑，也與整個區(qū)域的構(gòu)造背景符合.

3.3 應(yīng)力場反演結(jié)果

我們進(jìn)一步統(tǒng)計分析了這227個地震事件的P軸和T軸的空間方位分布.結(jié)果見圖11.P軸的方位分布整體集中在西北-南東方向，T軸的方位角分布整體更為集中，其優(yōu)勢走向仍為北東-南西.P軸與T軸的優(yōu)勢展布方向與前人的結(jié)果一致(吳建平等，2004)，與區(qū)域構(gòu)造背景和GPS觀測數(shù)據(jù)相符合(張培震等，2002).

圖11 PT軸展布(a)P軸方位角展布;(b)T軸方位角展布.Fig.11 Strike distribution of P and T axes(a)Strike distribution of P axes;(b)Strike distribution of T axes.

更進(jìn)一步，采用Hardebeck和Michael(2006)提出的基于法向應(yīng)力的SATSI反演算法，使用Martínez-Garzón等(2014)開發(fā)的matlab程序包，利用227個小震的震源機(jī)制解對小江斷裂帶的區(qū)域應(yīng)力場進(jìn)行了反演.

由于研究區(qū)域較小，地震活動分布不均勻.考慮到程序反演出穩(wěn)定結(jié)果對資料的數(shù)量有限制，僅將研究區(qū)域劃分成南北兩塊，以25.7°為分界線，分別反演南北兩個區(qū)域的平均應(yīng)力張量.反演過程中，為了客觀不加入任何先驗條件，我們將每個地震事件的兩個節(jié)面都作為可能的斷層面輸入程序，并將正確斷層面數(shù)量的比例保持在默認(rèn)值0.5，阻尼值采用程序默認(rèn)值2000.反演結(jié)果由圖12顯示.

圖12 應(yīng)力主軸方向分布圖紅、綠、藍(lán)依次代表最大至最小應(yīng)力主軸.(a)南區(qū)應(yīng)力主軸取向;(b)北區(qū)應(yīng)力主軸取向.Fig.12 Principal axes of stress fieldOrientation of maximum,medium and minimum principal axes are represented by red,green and blue dots,respectively.(a)Orientation of axes of south region;(b)Orientation of axes of north region.

從圖12中可以看出，在南北兩個區(qū)域中，應(yīng)力主軸的方向幾乎沒有改變.最大壓應(yīng)力軸與豎直方向有一定的偏離，次大壓應(yīng)力軸與最大壓應(yīng)力軸的走向相同，均為北西-南東方向，與上節(jié)中P軸的反演結(jié)果一致，而最小壓應(yīng)力軸則是與T軸的展布方向基本一致.最大壓應(yīng)力軸與次大壓應(yīng)力軸的傾角不確定度較大，但水平最大主應(yīng)力方向的不確定非常小.與之對應(yīng)的，最小主應(yīng)力軸的方位角和傾角都被約束的非常好，且與前面得到的T軸方向相近，為北東-南西的近似水平方向.

R值表征了三個方向主應(yīng)力的相對大小，其定義為：

橫向?qū)Ρ葍蓚€區(qū)域的R值，北區(qū)的R值為0.1，南區(qū)的R值為0.13，二者非常接近.從R值中可以推出，最大壓應(yīng)力軸所對應(yīng)的應(yīng)力大小與中間壓應(yīng)力軸相近.這說明P軸方向與豎直方向的應(yīng)力較為接近，而他們都遠(yuǎn)遠(yuǎn)大于最小壓應(yīng)力軸T軸的主應(yīng)力值.

4 討論與結(jié)論

傳統(tǒng)的人工震相到時讀取與P波初動極性判斷不僅費時，且結(jié)果嚴(yán)重依賴于人工的經(jīng)驗性專業(yè)知識，效率低、主觀性強(qiáng)，已不適合處理基于密集臺陣的豐富的小震記錄.本文采用的POI算法是從互信息理論出發(fā)的.互信息方法本身是成熟的信息學(xué)理論，在數(shù)學(xué)、物理等眾多學(xué)科中有廣泛的應(yīng)用.已有科研工作者在地震信號的拾取中使用互信息理論，以往的互信息方法往往使用互信息值或者其他基于互信息的函數(shù)值作為目標(biāo)函數(shù)，然后設(shè)立判據(jù)尋找初至到時(李輝等，2007).而本文的方法提供了一個完全不同的解決問題的思路.該方法首次在使用互信息方法，將到時拾取問題轉(zhuǎn)化為模型參數(shù)抉擇問題后，引入最大順序統(tǒng)計量準(zhǔn)則的概念，來解算參數(shù)的概率密度分布，并最終推導(dǎo)出初至到時的概率密度.POI算法沒有局限于互信息本身，原理簡單，使用的參數(shù)少，數(shù)學(xué)推理嚴(yán)密，同時在計算的過程中自然地引入了正態(tài)分布假設(shè)，減少了結(jié)果的隨機(jī)性和主觀設(shè)定先驗信息可能造成的誤識別.同時POI算法還可以給出不確定度的估計.對于信噪比較低的數(shù)據(jù)，其他方法往往只能通過設(shè)定閾值，對數(shù)據(jù)進(jìn)行取舍，而本研究采用的POI算法能夠判斷該段數(shù)據(jù)包括P波初動的概率，從而給出概率化的解，這在目前的初動極性檢測算法中是沒有的.

目前流行的地震目錄完備震級檢測和基于海量樣本訓(xùn)練的人工智能方法(Zhou et al.,2018;Zhou et al.,2019;Zhu and Beroza,2019)雖然有效地解決了地震事件的自動檢測與震相到時的自動讀取問題，但由于地震記錄的區(qū)域特性，通常需要針對研究區(qū)尋找包含研究區(qū)區(qū)域特征的新地震記錄樣本集進(jìn)行重新訓(xùn)練，給基于海量樣本訓(xùn)練的人工智能方法的移植性造成困難.相對人工智能類的方法，本文提出的方法無需搜集海量訓(xùn)練樣本，因此有更高的工作效率與可移植性.

本文的第二部分主要內(nèi)容是依據(jù)POI算法，聯(lián)合已有的成熟的反演程序HASH和MSATSI，開發(fā)了一套自動化反演震源機(jī)制解及區(qū)域應(yīng)力場的工作流程.該工作流程所需要的輸入數(shù)據(jù)包括經(jīng)過定位的地震目錄，區(qū)域速度結(jié)構(gòu)與連續(xù)波形記錄.本文采用北京大學(xué)與中國地震局地球物理研究所在小江區(qū)域布置的2015—2019年臺陣數(shù)據(jù)，以及674個定位誤差較小的可靠事件目錄作為資料，運用POI算法，對P波初動到時與極性進(jìn)行了自動拾取，計算了振幅比資料，并反演了小震的震源機(jī)制解以及區(qū)域平均應(yīng)力場.反演過程中得到了4084條P波初動極性記錄以及227個震源機(jī)制解.經(jīng)過與不同結(jié)果的比對，P波初動極性的拾取結(jié)果、震源機(jī)制的反演結(jié)果以及應(yīng)力場推斷結(jié)果均與前人研究較為吻合，符合基本的區(qū)域構(gòu)造特征和其他數(shù)據(jù)觀測資料.因此POI自動拾取算法，以及其聯(lián)合HASH與SATSI構(gòu)建的自動化震源機(jī)制與應(yīng)力場反演流程具有較高的可靠性和應(yīng)用性和可推廣性.我們希望更多地震學(xué)家能夠試用本文提出的方法與反演流程，將其應(yīng)用到不同區(qū)域和不同觀測條件下，不斷完善該方法，尋找最優(yōu)反演流程，并拓展其應(yīng)用范圍.

致謝中國地震局地球物理研究所許力生課題組慷慨地給我們提供他們的地震臺陣記錄資料.3位匿名評審專家給本文提出了許多很好的建設(shè)性修改意見，在此衷心感謝他們.