劉勁松，王 赟，姚振興

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

1 引 言

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

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

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

2 地震信號到時自動拾取的幾種方法

2.1 長短時均值比（STA／LTA）方法

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

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

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

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

振幅絕對值的平均值eabsi定義為：

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

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

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

（3）大振幅過零次數(shù)大于等于I6并且總過零次數(shù)大于等于I4，或者，當(dāng)前時刻距滿足條件1的時刻大于等于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方法

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

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

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

圖1 （a）一個地震信號的垂直分量波形；（b）時窗采用初至前1s至初至后3s的AIC函數(shù)；（c）時窗采用初至前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)計量偏斜度和峰度的PAI－S／K方法

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

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

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

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

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

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

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

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

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

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

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

4.1 移動時窗峰度的快速算法

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

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

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

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

4.2 到時判別算法

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

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

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

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

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

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

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

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

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

5 結(jié) 論

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

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

圖3 （a）—（d）分別為：垂直分量地震記錄波形，地震波形的峰度曲線，峰度曲線的長時均值，峰度曲線的短時均值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中初至附近的波形（實線）、峰度曲線長時均值（虛線）、短時均值（點(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方法在時窗正好包含了地震信號前后一段波形的情況下，能夠得到非常好的結(jié)果，適用于已有震源先驗信息的情形，例如流動地震觀測中記錄到的遠(yuǎn)震信號和區(qū)域地震信號，以及在已經(jīng)用其它方法檢測出地震信號所在時段后，再應(yīng)用AIC方法判斷到時.

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

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

（

）

［1］ 張金川，徐波，聶海寬等.中國頁巖氣資源勘探潛力.天然氣工業(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］ 趙金洲，王松，李勇明.頁巖氣藏壓裂改造難點(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］ 劉振武，撒利明，楊曉等.頁巖氣勘探開發(fā)對地球物理技術(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)測技術(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］ 王愛國，周瑤琪，陳勇等.基于微地震技術(shù)的油田裂縫監(jiān)測及模擬.中國海洋大學(xué)學(xué)報，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é)報，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］ 王繼，陳九輝，劉啟元等.流動地震臺陣觀測初至震相的自動檢測.地震學(xué)報，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］ 彼得·鮑曼.新地震觀察實踐手冊.中國地震局監(jiān)測預(yù)報司譯.金嚴(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.