基于最大順序統(tǒng)計(jì)量的P波初動(dòng)極性自動(dòng)讀取方法的小震震源機(jī)制反演及在云南小江地區(qū)的應(yīng)用研究

裴瑋來，周仕勇

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

0 引言

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

(1)基于長短窗信號(hào)比的震相到時(shí)檢測方法.通過計(jì)算長短兩個(gè)時(shí)窗的信號(hào)的能量比值來判定信號(hào)的初至到時(shí)(Stevenson,1976).短時(shí)能量(STA)顯著高于長時(shí)能量(LTA)時(shí)，設(shè)定閾值判定地震信號(hào)到達(dá)的時(shí)間.長短窗算法是最經(jīng)典的地震波信號(hào)拾取算法之一，常被拿來與其他方法進(jìn)行對(duì)比.為了增加算法的靈敏度和精確度，許多學(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ù))，并尋找給定時(shí)窗中AIC函數(shù)的最小值的位置來決定時(shí)窗內(nèi)地震信號(hào)的到時(shí)(Maeda,1985).后者通過計(jì)算隨機(jī)變量之間的互信息值，或以其他信息論概念為基礎(chǔ)構(gòu)建的目標(biāo)函數(shù)來判定地震震相到時(shí).互信息的概念如下(Shannon,1948)：

(1)

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

(3)基于高階統(tǒng)計(jì)量如偏度、峰度等參數(shù)開發(fā)的震相拾取方法，通過計(jì)算地震波形的峰度和偏斜度的極值點(diǎn)，尋找地震波震相到時(shí)(Nakamula et al.,2007;Ross et al.,2016;Saragiotis et al.,2002).其中偏斜度描述隨機(jī)變量相對(duì)于均值的不對(duì)稱程度，而峰度描述隨機(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)用存在一定程度的困難且較為耗時(shí).

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

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

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

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

1 方法原理

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

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

圖1 互信息到時(shí)判定利用互信息構(gòu)建的目標(biāo)函數(shù)Z(t,ε)，取目標(biāo)函數(shù)值最大的時(shí)刻tinit作為到時(shí).定義式中不同的事件用不同顏色區(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)性最高的到時(shí)隨機(jī)變量X，標(biāo)識(shí)了震相到時(shí)t.在閾值ε確定時(shí)，到時(shí)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即為閾值已知時(shí)認(rèn)定的真實(shí)地震波到時(shí)：

(5)

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

(6)

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

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

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

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

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

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

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

(8)

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

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

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

(10)

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

(11)

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

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

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

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

因此有

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

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

(15)

將Y=1或Y=-1代入公式(15)中，所求得的概率即為初動(dòng)極性為正或?yàn)樨?fù)的概率.

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

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

圖2 本研究采用的自動(dòng)化反演流程圖圖中的每一項(xiàng)計(jì)算過程都實(shí)現(xiàn)了自動(dòng)化.其中POI極性拾取是本文采用的新的初動(dòng)極性拾取方法.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.

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

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

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

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

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

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

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

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

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

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

圖3 本研究選取的研究區(qū)域，使用的地震臺(tái)與選取的674個(gè)地震空間分布圖倒三角代表固定臺(tái)站(昆明臺(tái)、東川臺(tái))，正三角代表寬頻帶流動(dòng)地震臺(tái)站，圓點(diǎn)代表地震位置.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 初動(dòng)極性判定結(jié)果

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

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

圖4 POI算法初動(dòng)極性判定結(jié)果與人工初動(dòng)極性判定結(jié)果的比較圖中的每一行代表一類人工拾取初動(dòng)結(jié)果，而每一行中三種極性的占比該類別數(shù)據(jù)被POI自動(dòng)算法拾取為向上、未知與向下3種情況的概率.從圖中可以看出Impulsive質(zhì)量的數(shù)據(jù)判定結(jié)果一致性高，對(duì)人工判定極性為上、下的數(shù)據(jù)分別有高達(dá)96.5%和92.5%的拾取正確率.對(duì)于極性質(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 初動(dòng)極性判定案例人工標(biāo)記該條記錄極性向上，質(zhì)量等級(jí)為Impulsive.POI算法計(jì)算其初動(dòng)極性100%的概率向上.圖中左上部分展示的是閾值的概率密度分布，右下圖展示的是與之對(duì)應(yīng)的初至到時(shí)的概率密度分布.主圖中的紅線和綠線展示兩者之間是如何對(duì)應(yīng)的.本例中到時(shí)幾乎百分之百集中在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 初動(dòng)極性判定案例人工標(biāo)記該條記錄極性向上，質(zhì)量等級(jí)為Emergent.POI算法計(jì)算其初動(dòng)極性100%的概率向下.圖中左上部分展示的是閾值的概率密度分布，右下圖展示的是與之對(duì)應(yīng)的初至到時(shí)的概率密度分布.主圖中的紅線和綠線展示兩者之間是如何對(duì)應(yīng)的.本例中到時(shí)幾乎百分之百集中在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 初動(dòng)極性判定案例人工標(biāo)記該條記錄極性未知.POI算法計(jì)算其初動(dòng)極性100%的概率未知.圖中左上部分展示的是閾值的概率密度分布，右下圖展示的是與之對(duì)應(yīng)的初至到時(shí)的概率密度分布.主圖中的紅線和綠線展示兩者之間是如何對(duì)應(yīng)的.本例中到時(shí)幾乎百分之百大于6 s，意味著POI算法判定在本段波形數(shù)據(jù)中不存在初動(dòng).因此極性為未知.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 初動(dòng)極性判定案例人工標(biāo)記該條記錄極性未知.POI算法計(jì)算其初動(dòng)極性100%的概率向上.圖中左上部分展示的是閾值的概率密度分布，右下圖展示的是與之對(duì)應(yīng)的初至到時(shí)的概率密度分布.主圖中的紅線和綠線展示兩者之間是如何對(duì)應(yīng)的.本例中到時(shí)幾乎百分之百集中在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波初動(dòng)極性與SH/P波的振幅比資料，反演了上述674個(gè)地震事件的震源機(jī)制解.

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

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

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

圖9 發(fā)生于2016年9月7日07時(shí)17分51.47秒的地震事件的P波初動(dòng)極性拾取情況與震源機(jī)制解該事件的走向、滑動(dòng)角及傾角分別為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個(gè)質(zhì)量較高解的地震事件的震源機(jī)制結(jié)果.其中走滑類型、正斷層類型、逆沖斷層類型和混合型地震事件的數(shù)量分別為87，84，31，25，它們占總事件數(shù)量的百分比分別為38.5%，37%，13.5%.11%.走滑類型是占比最多的地震事件類型，其次是正斷層型，這與前人計(jì)算出的震源機(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)橹?，展示了?duì)小江主斷裂帶的南北走向特征的刻畫.雖然偶爾有滑動(dòng)方向不嚴(yán)格遵從左旋走滑的地震事件發(fā)生，大多數(shù)沿小江斷裂帶主斷層分布的事件的震源機(jī)制解仍然是左旋走滑，也與整個(gè)區(qū)域的構(gòu)造背景符合.

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

我們進(jìn)一步統(tǒng)計(jì)分析了這227個(gè)地震事件的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個(gè)小震的震源機(jī)制解對(duì)小江斷裂帶的區(qū)域應(yīng)力場進(jìn)行了反演.

由于研究區(qū)域較小，地震活動(dòng)分布不均勻.考慮到程序反演出穩(wěn)定結(jié)果對(duì)資料的數(shù)量有限制，僅將研究區(qū)域劃分成南北兩塊，以25.7°為分界線，分別反演南北兩個(gè)區(qū)域的平均應(yīng)力張量.反演過程中，為了客觀不加入任何先驗(yàn)條件，我們將每個(gè)地震事件的兩個(gè)節(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中可以看出，在南北兩個(gè)區(qū)域中，應(yīng)力主軸的方向幾乎沒有改變.最大壓應(yīng)力軸與豎直方向有一定的偏離，次大壓應(yīng)力軸與最大壓應(yīng)力軸的走向相同，均為北西-南東方向，與上節(jié)中P軸的反演結(jié)果一致，而最小壓應(yīng)力軸則是與T軸的展布方向基本一致.最大壓應(yīng)力軸與次大壓應(yīng)力軸的傾角不確定度較大，但水平最大主應(yīng)力方向的不確定非常小.與之對(duì)應(yīng)的，最小主應(yīng)力軸的方位角和傾角都被約束的非常好，且與前面得到的T軸方向相近，為北東-南西的近似水平方向.

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

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

4 討論與結(jié)論

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

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

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

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