微地震信號(hào)到時(shí)自動(dòng)拾取方法

劉勁松，王 赟，姚振興

1中國(guó)科學(xué)院地質(zhì)與地球物理研究所，北京 100029 2中國(guó)科學(xué)院地球化學(xué)研究所，貴陽(yáng) 550002

1 引 言

近年來(lái)美國(guó)大規(guī)模的頁(yè)巖氣開采引發(fā)了各國(guó)對(duì)頁(yè)巖氣資源的關(guān)注.我國(guó)致密砂巖油氣藏和頁(yè)巖氣藏資源豐富，在松遼、鄂爾多斯、四川及南方諸多盆地均有分布［1］，這些致密砂巖油氣藏和頁(yè)巖氣藏的天然氣資源將成為我國(guó)重要的能源支柱.對(duì)致密砂巖油氣藏和頁(yè)巖氣藏主要采用水平井分段水力壓裂技術(shù)方式開采.水力壓裂使地層形成新的裂縫，產(chǎn)生微震.監(jiān)測(cè)微地震的分布和特征可為水力壓裂設(shè)計(jì)和壓裂效果的監(jiān)測(cè)提供依據(jù)，提高開采效益［2－6］.

微地震記錄的特點(diǎn)是頻率高、信噪比低，因此微地震事件的自動(dòng)識(shí)別和初至到時(shí)拾取對(duì)實(shí)現(xiàn)海量微地震數(shù)據(jù)的自動(dòng)處理有重要意義［7］.對(duì)于天然地震事件，已提出了多種自動(dòng)識(shí)別方法.這些方法主要根據(jù)在地震記錄中地震事件到達(dá)前后質(zhì)點(diǎn)振動(dòng)性質(zhì)差別構(gòu)建特征函數(shù)加以判斷.例如，根據(jù)在時(shí)間域能量和 能 量 變 化 構(gòu) 建 特 征 函 數(shù)，Allen［8－9］、Baer和Kradolfer［10］提出了長(zhǎng)短時(shí)均值比方法（STA／LTA法）；根據(jù)在時(shí)間－頻率域能量變化，Gibbons等［11］提出了頻譜比方法；根據(jù)地震信號(hào)到達(dá)前后地震波形數(shù)據(jù)統(tǒng)計(jì)性質(zhì)的差別，Maeda［12］提出了AIC方法，Saragiotis［13］提出了基于地震波形偏斜度和峰度的PAI－S／K方法；以及常旭和劉伊克［14］討論的分形分維方法等.王繼等［15］討論了用AIC方法對(duì)流動(dòng)地震臺(tái)陣觀測(cè)記錄初至震相的自動(dòng)檢測(cè)問題.與天然地震的記錄相比，微震的震級(jí)更小，通常在1級(jí)以下，而且信噪比更低，有必要深入地討論這些方法對(duì)微地震事件的識(shí)別問題.

本文利用STA／LTA、AIC和PAI－S／K方法，采用兩種時(shí)窗，并利用全局優(yōu)化搜索方法選取拾取的控制參數(shù)，對(duì)在我國(guó)西部某地觀測(cè)到的震中距范圍在0.3至20km、震級(jí)范圍－0.5至1級(jí)的約13359個(gè)微震記錄進(jìn)行初至到時(shí)拾取，比較了上述三種方法的優(yōu)缺點(diǎn)，提出了相應(yīng)的改進(jìn)方法，以利于今后的微震監(jiān)測(cè)工作.

2 地震信號(hào)到時(shí)自動(dòng)拾取的幾種方法

2.1 長(zhǎng)短時(shí)均值比（STA／LTA）方法

Allen［8－9］提出的STA／LTA方法，是目前廣泛使用的一種地震信號(hào)到時(shí)自動(dòng)拾取方法，其主要原理是根據(jù)地震波形特征函數(shù)的長(zhǎng)短時(shí)均值，以及地震波形過(guò)零軸的次數(shù)統(tǒng)計(jì)等特征拾取初至，原理如下：

設(shè)di，i＝0，1，2，…，n－1為離散地震波形數(shù)據(jù)，fi為高通濾波后的數(shù)據(jù)：

根據(jù)質(zhì)點(diǎn)振動(dòng)的能量和能量變化定義特征函數(shù)Ei為

其中f′i為fi對(duì)時(shí)間的導(dǎo)數(shù)，C2為常數(shù).特征函數(shù)的短時(shí)均值函數(shù)αi和長(zhǎng)時(shí)均值函數(shù)βi定義為

振幅絕對(duì)值的平均值eabsi定義為：

判定初至的準(zhǔn)則是同時(shí)滿足下列條件［8］：

（1）αi≥C5βi，eabsi≤C7，并且后續(xù)的波形同時(shí)滿足下列條件；

（2）總過(guò)零次數(shù)（mzc）大于等于I4，或者，連續(xù)小振幅過(guò)零次數(shù)大于等于I3＋mzc／I3；

（3）大振幅過(guò)零次數(shù)大于等于I6并且總過(guò)零次數(shù)大于等于I4，或者，當(dāng)前時(shí)刻距滿足條件1的時(shí)刻大于等于D5.

以上C1－C8以及I3、I4、I5等為常數(shù).SAC軟件系統(tǒng)［16］含有STA／LTA方法的功能模塊（apk命令，pkdet等相關(guān)函數(shù)），SAC軟件建議的缺省值如表1所示.

表1 STA／LTA方法參數(shù)意義及取值Table 1 Parameters of STA／LTA method

2.2 AIC方法

該方法的基本原理是求取地震信號(hào)AIC函數(shù)局部最小值的位置.AIC函數(shù)定義為［12，15］：

其中｛xi｝，i＝1，2，3，…，N為地震信號(hào)的離散波形序列.

在時(shí)窗正好包含了地震信號(hào)前后一段波形的情況下，AIC函數(shù)在初至點(diǎn)峰值尖銳，準(zhǔn)確率較高.但AIC函數(shù)受時(shí)窗的影響較大，在不同的時(shí)窗下，其局部最小值出現(xiàn)的位置往往不同，如圖1所示.因此該方法適用于已知初至的大致位置的情況.

圖1 （a）一個(gè)地震信號(hào)的垂直分量波形；（b）時(shí)窗采用初至前1s至初至后3s的AIC函數(shù)；（c）時(shí)窗采用初至前1s至初至后9s的AIC函數(shù)Fig.1 Vertical component of an earthquake signal（a）and its AIC function of time window from 1second before onset to 3second after onset（b），and AIC function of time window from 1second before onset to 9seconds after onset（c）

2.3 基于高階統(tǒng)計(jì)量偏斜度和峰度的PAI－S／K方法

2002年，Saragiotis提出了基于地震波形峰度和偏斜度的拾取初至方法，稱之為PAI－S／K方法［13］.對(duì)于一個(gè)有限長(zhǎng)度的離散實(shí)數(shù)序列｛xi｝，其偏斜度S（skewness）和峰度K（kurtosis）定義為：

其中為｛xi｝序列的平均值，σx為其標(biāo)準(zhǔn)差.

利用地震波形峰度和偏斜度檢測(cè)地震波到時(shí)的基本原理是判斷峰度曲線和偏斜度曲線的極值點(diǎn)，以極值點(diǎn)前曲線斜率最大的位置作為地震波到時(shí)的位置，該方法具有計(jì)算穩(wěn)定強(qiáng)健、實(shí)現(xiàn)簡(jiǎn)單等優(yōu)點(diǎn).圖2為一個(gè)地震記錄及其峰度和偏斜度曲線.

3 實(shí)際數(shù)據(jù)的拾取結(jié)果

本文應(yīng)用上述方法，對(duì)我國(guó)西部某地的13359條微地震觀測(cè)數(shù)據(jù)記錄進(jìn)行了初至到時(shí)自動(dòng)拾取，同時(shí)人工拾取了P波到時(shí)和S波到時(shí).數(shù)據(jù)涉及1600多個(gè)微地震事件，絕大多數(shù)的地震震級(jí)小于1，震中距范圍從0.3km到20km左右.自動(dòng)拾取采用兩種時(shí)窗，時(shí)窗1為［tp－3，tp＋2tps］，即P波初至前3s至初至后2倍縱橫波走時(shí)差，該時(shí)窗包含了P波初至和S波初至.時(shí)窗2為［tp－3，tp＋tps］，該時(shí)窗起始于P波初至前3s，終止于S波開始時(shí)刻，未包含S波.時(shí)窗1對(duì)應(yīng)于一般情形，實(shí)際數(shù)據(jù)處理中多數(shù)為這種情況.時(shí)窗2對(duì)應(yīng)已知初至大致位置的情況，例如流動(dòng)地震觀測(cè)中的遠(yuǎn)震.拾取前對(duì)數(shù)據(jù)進(jìn)行了5～30Hz的帶通濾波.

對(duì)于STA／LTA方法，首先采用SAC軟件的缺省參數(shù)進(jìn)行拾取，自動(dòng)拾取的初至與人工拾取的初至?xí)r差在0.3s以內(nèi)的記錄約占55%，時(shí)差在0.1s以內(nèi)的約占54%.由于SAC軟件中建議的缺省參數(shù)是針對(duì)天然地震的，其震級(jí)相對(duì)較大，對(duì)于微地震的情況，這些參數(shù)可能不是最佳的.為得到最佳的參數(shù)，本文采用DE全局搜索方法［17］，以人工拾取的初至作為參照，以時(shí)差在0.3s以內(nèi)的記錄所占百分比作為目標(biāo)函數(shù)，自動(dòng)搜索最佳的拾取參數(shù).得到的結(jié)果是，小于0.3s的記錄可增加到77.1%，如表2所示.最佳參數(shù)見表1的第3列.從所求得的最佳參數(shù)表明，至少對(duì)于本文所用的數(shù)據(jù)，在特征函數(shù)中加入導(dǎo)數(shù)項(xiàng)，其效果不如僅采用振幅平方作為特征函數(shù).

AIC方法的拾取結(jié)果如表2所示.結(jié)果表明，在時(shí)窗2的情況下，AIC方法是本文討論的幾種方法中準(zhǔn)確率最高的方法.由AIC函數(shù)的計(jì)算公式可知，某個(gè)點(diǎn)上的AIC函數(shù)值，依賴于搜索時(shí)窗內(nèi)每個(gè)點(diǎn)的值，當(dāng)搜索時(shí)窗變化時(shí)，對(duì)應(yīng)波形某特定點(diǎn)上的AIC函數(shù)值也是不同的，因此AIC方法不適合沿連續(xù)地震記錄滑動(dòng)進(jìn)行.

PAI－S／K方法的拾取結(jié)果見表2的第3～4行.Saragiotis［13］用該方法對(duì)44個(gè)區(qū)域地震范圍的地震信號(hào)記錄進(jìn)行了自動(dòng)拾取，并與STA／LTA方法的結(jié)果進(jìn)行了對(duì)比，其結(jié)論是該方法優(yōu)于STA／LTA方法.但從本文的實(shí)際數(shù)據(jù)拾取結(jié)果看，在時(shí)窗1的情況下該方法拾取的準(zhǔn)確率不如采用優(yōu)化參數(shù)后的STA／LTA方法.PAI－S／K方法通過(guò)判斷搜索時(shí)窗內(nèi)的最大值點(diǎn)確定初至位置，因而拾取結(jié)果也受搜索時(shí)窗范圍的影響.該方法只有一個(gè)控制參數(shù)，即計(jì)算峰度和偏斜度的滑動(dòng)時(shí)窗長(zhǎng)度.同樣采用DE全局搜索方法，可得到對(duì)于本文實(shí)際數(shù)據(jù)最佳的滑動(dòng)時(shí)窗長(zhǎng)度.對(duì)于偏斜度，最佳滑動(dòng)時(shí)窗長(zhǎng)度為2.70s，對(duì)于峰度，最佳滑動(dòng)時(shí)窗長(zhǎng)度為0.72s.滑動(dòng)時(shí)窗長(zhǎng)度對(duì)結(jié)果也有一定影響，但在一定變化范圍內(nèi)，對(duì)拾取結(jié)果的影響不大.例如，對(duì)于峰度法，當(dāng)滑動(dòng)時(shí)窗長(zhǎng)度在［0.5，1.5］之間變化時(shí)，與人工拾取結(jié)果在0.3s以內(nèi)的結(jié)果所占百分比的變化在1個(gè)點(diǎn)左右.

圖2 一個(gè)地震記錄波形（a）及其峰度（b）和偏斜度曲線（c）Fig.2 An earthquake signal（a）and its kurtosis（b）and skewness（c）curves

表2 各種方法初至拾取情況一覽表（表中計(jì)算耗時(shí)未計(jì)入硬盤IO時(shí)間）Table 2 Picking results of each method（the IO operation time is not included in time consumption）

4 對(duì)PAI－S／K方法的改進(jìn)

4.1 移動(dòng)時(shí)窗峰度的快速算法

對(duì)于一個(gè)長(zhǎng)度為N的地震波形序列｛xi｝，設(shè)求峰度的時(shí)窗長(zhǎng)度為M，按照公式（7）沿地震波形滑動(dòng)求取峰度，要調(diào)用（N－M）次公式，計(jì)算量為O（M·（N－M）），對(duì)于海量的多道連續(xù)觀測(cè)記錄，其計(jì)算效率較低.本文將上述公式展開，將公式變換為求取離散實(shí)數(shù)序列的算術(shù)和、平方和、立方和、4次方和的線性組合形式，減少重復(fù)計(jì)算，計(jì)算量減少為O（N），極大地提高了計(jì)算效率.

將（7）式展開并變換得到：

將時(shí)窗內(nèi)離散序列的算術(shù)和、平方和、立方和、4次方和分別記為smj，j＝1，2，3，4，得到：

采用（10）—（12）式計(jì)算的優(yōu)點(diǎn)是，當(dāng)時(shí)窗移動(dòng)時(shí)計(jì)算新時(shí)窗的smj不必重復(fù)計(jì)算時(shí)窗內(nèi)所有點(diǎn)的值，只須在上一個(gè)時(shí)窗各值的基礎(chǔ)上，加上新進(jìn)入時(shí)窗的值，并減去離開時(shí)窗的值即可.經(jīng)實(shí)際數(shù)據(jù)計(jì)算測(cè)試，在時(shí)窗長(zhǎng)度取0.5s、1.0s、1.5s的情況下，采用快速算法的計(jì)算速度分別提高了22倍、32倍、43倍.對(duì)于偏斜度以及AIC函數(shù)的計(jì)算，同樣可采用類似的方法提高速度.

4.2 到時(shí)判別算法

原方法采用峰度最大值點(diǎn)左側(cè)斜率最大的點(diǎn)作為初至位置，但實(shí)際數(shù)據(jù)的結(jié)果顯示，該點(diǎn)往往滯后于真實(shí)的初至位置，而峰度曲線開始起跳的位置更接近初至位置.在某些情況下，初至附近的峰度極值點(diǎn)也不是最大值點(diǎn).針對(duì)這些情況，本文將STA／LTA方法與PAI－S／K法結(jié)合，以地震波形的峰度K作為特征函數(shù)，以特征函數(shù)的長(zhǎng)短時(shí)均值比和峰度值作為判斷初至位置的依據(jù).具體方法如下：

定義特征函數(shù)Ei為

其中Ki為第i個(gè)樣點(diǎn)處的峰度，Ki的值由第i個(gè)樣點(diǎn)前M個(gè)地震波形數(shù)據(jù)計(jì)算.與STA／LTA的方法類似，短時(shí)均值函數(shù)αi和長(zhǎng)時(shí)均值函數(shù)βi定義為：

滿足下列條件的位置判定為初至：

（1）αi≥C5βi并且αi≥C6；

（2）滿足條件1的位置記為n，由n開始向前（時(shí)間減小方向）搜索Ki的值，直到滿足Ki≤Ki－1.

圖3為一個(gè)地震記錄的波形和峰度、以及峰度的長(zhǎng)短時(shí)均值曲線，圖4為初至附近的局部放大圖.

上述公式中的C3—C6為常數(shù).將上述方法應(yīng)用于本文的實(shí)際數(shù)據(jù)，同樣采用DE全局搜索法搜索最佳參數(shù)，保持C3＝0.6，C4＝0.03不變，以人工拾取的初至作為參照，以時(shí)差在0.3s以內(nèi)的記錄所占百分比作為目標(biāo)函數(shù)，自動(dòng)搜索最佳的C5和C6的值以及滑動(dòng)時(shí)窗長(zhǎng)度.對(duì)于本文的數(shù)據(jù)當(dāng)C5＝2.71，C6＝1.43，滑動(dòng)時(shí)窗長(zhǎng)度為0.79s時(shí)，自動(dòng)識(shí)別的效果最佳，得到的結(jié)果如表2中所示.從表2可以看出，在時(shí)窗1的情況下，識(shí)別準(zhǔn)確率比原PAI－S／K方法提高了約10個(gè)百分點(diǎn)，比STA／LTA方法提高了約5個(gè)百分點(diǎn).

本方法主要有兩個(gè)參數(shù)C5和C6，分別是長(zhǎng)短時(shí)均值比閥值和峰度值短時(shí)均值的閥值.在P波初至附近的極值點(diǎn)處，長(zhǎng)短時(shí)均值比和極值都隨地震信號(hào)的信噪比增加而增加，峰度的極值尤其明顯.根據(jù)（7）式和（8）式，偏斜度和峰度均為無(wú)量綱量，因而對(duì)于不同類型的地震記錄（速度、加速度、不同的儀器增益等），該方法的拾取參數(shù)具有一定的普遍意義.

5 結(jié) 論

本文通過(guò)地震信號(hào)初至拾取STA／LTA方法、AIC方法和PAI－S／K方法的討論，以及我國(guó)西部某地13359條微震實(shí)際數(shù)據(jù)的自動(dòng)拾取與人工拾取結(jié)果對(duì)比分析，可以獲得如下認(rèn)識(shí)：

STA／LTA法歷史悠久，算法穩(wěn)定強(qiáng)健，目前仍被廣泛使用.該方法可沿地震記錄波形滑動(dòng)進(jìn)行，因而可用于地震信號(hào)的檢測(cè)，但該方法涉及參數(shù)較多，不易掌控.

圖3 （a）—（d）分別為：垂直分量地震記錄波形，地震波形的峰度曲線，峰度曲線的長(zhǎng)時(shí)均值，峰度曲線的短時(shí)均值Fig.3 From（a）to（d）are an earthquake signal and its kurtosis，long time average of kurtosis，and short time average of kurtosis curves，respectively

圖4 圖3中初至附近的波形（實(shí)線）、峰度曲線長(zhǎng)時(shí)均值（虛線）、短時(shí)均值（點(diǎn)線）.波形曲線按照峰度曲線做了歸一化處理.Fig.4 The enlarged part of Fig.3near onset.Solid line，dashed line and dotted line represent waveform，LTA curve and STA curve of waveform respectively.The waveform curve is rescaled according to kurtosis curve.

AIC方法在時(shí)窗正好包含了地震信號(hào)前后一段波形的情況下，能夠得到非常好的結(jié)果，適用于已有震源先驗(yàn)信息的情形，例如流動(dòng)地震觀測(cè)中記錄到的遠(yuǎn)震信號(hào)和區(qū)域地震信號(hào)，以及在已經(jīng)用其它方法檢測(cè)出地震信號(hào)所在時(shí)段后，再應(yīng)用AIC方法判斷到時(shí).

PAI－S／K方法原理簡(jiǎn)單實(shí)現(xiàn)容易，算法同樣穩(wěn)定強(qiáng)健，由于該方法采用了高階統(tǒng)計(jì)量判斷信號(hào)初至，具有一定的抗干擾能力.

本文提出的改進(jìn)后的PAI－S／K法，可方便地沿連續(xù)地震記錄波形滑動(dòng)進(jìn)行，因而可用于地震信號(hào)的檢測(cè)，在時(shí)窗1的情況下，小于0.3s的記錄占比達(dá)到83.8%，是上述方法中最高的；該方法涉及的參數(shù)較少，同時(shí)由于峰度的無(wú)量綱特性，因而所涉及參數(shù)的選取具有一定的普遍性.

（

）

［1］ 張金川，徐波，聶海寬等.中國(guó)頁(yè)巖氣資源勘探潛力.天然氣工業(yè)，2008，28（6）：136－142.

Zhang J C，Xu B，Nie H K，et al.Exploration potential of shale gas resources in China.NaturalGasIndustry（in Chinese），2008，28（6）：136－142.

［2］ 趙金洲，王松，李勇明.頁(yè)巖氣藏壓裂改造難點(diǎn)與技術(shù)關(guān)鍵.開發(fā)工程，2012，32（4）：1－4.

Zhao J Z，Wang S，Li Y M.Difficulties and key techniques in the fracturing treatment of shale gas reservoirs.NaturalGas Industry（in Chinese），2012，32（4）：1－4.

［3］ 劉振武，撒利明，楊曉等.頁(yè)巖氣勘探開發(fā)對(duì)地球物理技術(shù)的需求.石油地球物理勘探，2011，46（5）：810－820.

Liu Z W，Sa L M，Yang X，et al.Needs of geophysical technologies for shale gas exploration.OilGeophysical Prospecting（in Chinese），2011，46（5）：810－820.

［4］ 張山，劉清林，趙群等.微地震監(jiān)測(cè)技術(shù)在油田開發(fā)中的應(yīng)用.石油物探，2002，41（2）：226－231.

Zhang S，Liu Q L，Zhao Q，et al.Application of microseismic monitoring technology in development of oil field.GeophysicalProspectingforPetroleum（in Chinese），2002，41（2）：226－231.

［5］ 王愛國(guó)，周瑤琪，陳勇等.基于微地震技術(shù)的油田裂縫監(jiān)測(cè)及模擬.中國(guó)海洋大學(xué)學(xué)報(bào)，2008，38（1）：116－120.

Wang A G，Zhou Y Q，Chen Y，et al.Monitoring and simulation of fracture in oil field based on micro－seismic technology.PeriodicalofOceanUniversityofChina（in Chinese），2008，38（1）：116－120.

［6］ Duncan P M，Eisner L.Reservoir characterization using surface microseismic monitoring.Geophysics，2010，75（5）：139－146.

［7］ Warpinski N.Microseismic monitoring：inside and out.JournalofPetroleumTechnology，2009，61（11）：80－85.

［8］ Allen R V.Automatic earthquake recognition and timing from single traces.Bull.Seismol.Soc.Amer.，1978，68（5）：1521－1532.

［9］ Allen R V.Automatic phase pickers：their present use and future prospects.Bull.Seismol.Soc.Amer.，1982，72（6）：225－242.

［10］ Baer M，Kradolfer U.An automatic phase picker for local and teleseismic events.Bull.Seismol.Soc.Amer.，1987，77（4）：1437－1445.

［11］ Gibbons S J，Ringdal F，Kv?rna T.Detection and characterization of seismic phases using continuous spectral estimation on incoherent and partially coherent arrays.Geophys.J.Int.，2008，172（1）：405－421.

［12］ Maeda N.A method for reading and checking phase times in auto－processing system of seismic wave data.Zisin＝Jishin，1985，38（3）：365－379.

［13］ Saragiotis C D.PAI－S／K：A robust automatic seismic P phase arrival identification scheme.IEEETransactionson GeoscienceandRemoteSensing，2002，40（6）：1395－1404.

［14］ 常旭，劉伊克.地震記錄的廣義分維及其應(yīng)用.地球物理學(xué)報(bào)，2002，45（6）：839－846.

Chang X，Liu Y K.The generalized fractal dimension of seismic records and its application.ChineseJ.Geophys.（in Chinese），2002，45（6）：839－846.

［15］ 王繼，陳九輝，劉啟元等.流動(dòng)地震臺(tái)陣觀測(cè)初至震相的自動(dòng)檢測(cè).地震學(xué)報(bào)，2006，28（1）：42－51.

Wang J，Chen J H，Liu Q Y，et al.Automatic onset phase picking for portable seismic array observation.Acta SeismologicaSinica（in Chinese），2006，28（1）：42－51.

［16］ 彼得·鮑曼.新地震觀察實(shí)踐手冊(cè).中國(guó)地震局監(jiān)測(cè)預(yù)報(bào)司譯.金嚴(yán)，陳陪善，許忠淮校.北京：地震出版社，2006.

Peter B.New Manual of Seismological of Observatory Practice（in Chinese）.Beijing：Seismological Press，2006.

［17］ Storn R，Price K.Differential evolution－a simple and efficient adaptive scheme for global optimization over continuous spaces.Technical Report TR－95－012，ICSI，1995.