基于最大順序統計量的P波初動極性自動讀取方法的小震震源機制反演及在云南小江地區的應用研究

裴瑋來，周仕勇

北京大學地球與空間科學學院，北京 100871

0 引言

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

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

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

(1)

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

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

(4)以機器學習為代表的利用人工智能搭建神經網絡進行地震波震相拾取的方法.通過模擬人類的學習行為來訓練神經網絡，從而達到完成特定任務的目的(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).基于人工智能的方法雖然能給出較高的準確率，但需要提供大量標記樣本數據來訓練模型，因此在不同地區和數據間的遷移和應用存在一定程度的困難且較為耗時.

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

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

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

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

1 方法原理

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

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

圖1 互信息到時判定利用互信息構建的目標函數Z(t,ε)，取目標函數值最大的時刻tinit作為到時.定義式中不同的事件用不同顏色區分出來.此處閾值約為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)

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

(3)

其中m和n取0或1,Pm 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指代使互信息目標函數Z(t,ε)最大的t的取值，tinit即為閾值已知時認定的真實地震波到時：

(5)

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

(6)

上式中I為邏輯函數，其括號內條件成立則取值為1，否則取值為0.

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

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

我們提出最大順序統計量準則：閾值ε應當與其對應到時tinit前置的背景噪音Asamp，在正態分布假設下的最大順序統計量相近，即

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

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

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

(8)

這里引入正態分布假設是因為真實的噪音往往是服從正態分布的.因此真正的閾值所對應的前置噪音應當也是大致服從正態分布.在正態分布假設下，前置的背景噪聲可以認為是服從如下總體的隨機采樣：

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

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

(10)

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

(11)

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

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

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

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

因此有

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

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

(15)

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

1.2 基于P波初動的震源機制解與區域應力場反演流程

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

圖2 本研究采用的自動化反演流程圖圖中的每一項計算過程都實現了自動化.其中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.

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

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

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

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

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

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

2 研究區域與研究資料

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

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

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

圖3 本研究選取的研究區域，使用的地震臺與選取的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 結果分析

3.1 初動極性判定結果

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

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

圖4 POI算法初動極性判定結果與人工初動極性判定結果的比較圖中的每一行代表一類人工拾取初動結果，而每一行中三種極性的占比該類別數據被POI自動算法拾取為向上、未知與向下3種情況的概率.從圖中可以看出Impulsive質量的數據判定結果一致性高，對人工判定極性為上、下的數據分別有高達96.5%和92.5%的拾取正確率.對于極性質量一般的數據判定的正確率平均超過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 初動極性判定案例人工標記該條記錄極性向上，質量等級為Impulsive.POI算法計算其初動極性100%的概率向上.圖中左上部分展示的是閾值的概率密度分布，右下圖展示的是與之對應的初至到時的概率密度分布.主圖中的紅線和綠線展示兩者之間是如何對應的.本例中到時幾乎百分之百集中在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 初動極性判定案例人工標記該條記錄極性向上，質量等級為Emergent.POI算法計算其初動極性100%的概率向下.圖中左上部分展示的是閾值的概率密度分布，右下圖展示的是與之對應的初至到時的概率密度分布.主圖中的紅線和綠線展示兩者之間是如何對應的.本例中到時幾乎百分之百集中在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 初動極性判定案例人工標記該條記錄極性未知.POI算法計算其初動極性100%的概率未知.圖中左上部分展示的是閾值的概率密度分布，右下圖展示的是與之對應的初至到時的概率密度分布.主圖中的紅線和綠線展示兩者之間是如何對應的.本例中到時幾乎百分之百大于6 s，意味著POI算法判定在本段波形數據中不存在初動.因此極性為未知.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 初動極性判定案例人工標記該條記錄極性未知.POI算法計算其初動極性100%的概率向上.圖中左上部分展示的是閾值的概率密度分布，右下圖展示的是與之對應的初至到時的概率密度分布.主圖中的紅線和綠線展示兩者之間是如何對應的.本例中到時幾乎百分之百集中在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 震源機制解反演結果

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

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

表1 震源機制解質量分級表Table 1 Quality grading of focal mechanisms

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

圖9 發生于2016年9月7日07時17分51.47秒的地震事件的P波初動極性拾取情況與震源機制解該事件的走向、滑動角及傾角分別為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個質量較高解的地震事件的震源機制結果.其中走滑類型、正斷層類型、逆沖斷層類型和混合型地震事件的數量分別為87，84，31，25，它們占總事件數量的百分比分別為38.5%，37%，13.5%.11%.走滑類型是占比最多的地震事件類型，其次是正斷層型，這與前人計算出的震源機制解類型大致一致(嚴川,2015)，也與區域的構造背景和GPS觀測數據相一致.

圖10 4種震源機制解的空間分布情況(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.

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

3.3 應力場反演結果

我們進一步統計分析了這227個地震事件的P軸和T軸的空間方位分布.結果見圖11.P軸的方位分布整體集中在西北-南東方向，T軸的方位角分布整體更為集中，其優勢走向仍為北東-南西.P軸與T軸的優勢展布方向與前人的結果一致(吳建平等，2004)，與區域構造背景和GPS觀測數據相符合(張培震等，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.

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

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

圖12 應力主軸方向分布圖紅、綠、藍依次代表最大至最小應力主軸.(a)南區應力主軸取向;(b)北區應力主軸取向.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中可以看出，在南北兩個區域中，應力主軸的方向幾乎沒有改變.最大壓應力軸與豎直方向有一定的偏離，次大壓應力軸與最大壓應力軸的走向相同，均為北西-南東方向，與上節中P軸的反演結果一致，而最小壓應力軸則是與T軸的展布方向基本一致.最大壓應力軸與次大壓應力軸的傾角不確定度較大，但水平最大主應力方向的不確定非常小.與之對應的，最小主應力軸的方位角和傾角都被約束的非常好，且與前面得到的T軸方向相近，為北東-南西的近似水平方向.

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

橫向對比兩個區域的R值，北區的R值為0.1，南區的R值為0.13，二者非常接近.從R值中可以推出，最大壓應力軸所對應的應力大小與中間壓應力軸相近.這說明P軸方向與豎直方向的應力較為接近，而他們都遠遠大于最小壓應力軸T軸的主應力值.

4 討論與結論

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

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

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

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