李 劍，韓 焱，徐利娜，馬亮星

（1.中北大學(xué)信息探測(cè)與處理技術(shù)研究所，山西太原030051；2.中國(guó)兵器工業(yè)試驗(yàn)測(cè)試研究院，陜西華陰714200；3.山西北方機(jī)械制造有限責(zé)任公司，山西太原030009）

提取初至波到時(shí)的改進(jìn)能量因子算法

李 劍1，韓 焱1，徐利娜2，馬亮星3

（1.中北大學(xué)信息探測(cè)與處理技術(shù)研究所，山西太原030051；2.中國(guó)兵器工業(yè)試驗(yàn)測(cè)試研究院，陜西華陰714200；3.山西北方機(jī)械制造有限責(zé)任公司，山西太原030009）

針對(duì)目前能量因子算法存在初至波到時(shí)提取精度低、誤判率高的問(wèn)題，提出了提取初至波到時(shí)的改進(jìn)能量因子算法。該算法采用三分量能量合成信息作為特征函數(shù)，引入高階統(tǒng)計(jì)量因子提高波形幅值突變特征，利用廣義S變換對(duì)應(yīng)的瞬時(shí)優(yōu)勢(shì)頻率構(gòu)建自適應(yīng)時(shí)窗，增強(qiáng)初至波檢測(cè)因子的時(shí)變特性。試驗(yàn)結(jié)果表明，提出的改進(jìn)算法，提取初至波到時(shí)誤檢率降低，精度提高，在地下震源定位、層析成像、地質(zhì)勘探等地下空間領(lǐng)域具有廣泛的應(yīng)用價(jià)值。

初至波到時(shí)；地下震源定位；自適應(yīng)時(shí)窗；廣義S變換

0 引言

初至波是一種特殊類型的波，它是從震源激發(fā)直接傳到接收點(diǎn)的彈性波，具有起跳時(shí)間早、能量強(qiáng)的特點(diǎn)，是判斷震動(dòng)波到達(dá)的重要依據(jù)，同時(shí)也是實(shí)現(xiàn)地下震源定位、走時(shí)層析成像的前提和保障。

由于震源近場(chǎng)尤其是壓實(shí)區(qū)和破碎區(qū)電磁干擾嚴(yán)重，震動(dòng)波的瞬時(shí)屬性受噪聲的影響較大，造成初至波到時(shí)提取精度低、誤判率高。為了實(shí)現(xiàn)高精度的震源定位，如何在大噪聲背景下提取可靠的初至波到時(shí)是目前地下空間定位領(lǐng)域亟待解決的問(wèn)題。

目前初至波到時(shí)的提取算法可以分為以下幾類：基于地震記錄整體特征的方法，如相關(guān)法［1］、線性最小平方預(yù)測(cè)法等［2］。這類方法對(duì)噪聲有一定的抑制作用，但受地震道相關(guān)性等因素影響，對(duì)于地表情況復(fù)雜地區(qū)的資料效果較差；基于地震初至波多維信息的人工智能方法，如神經(jīng)網(wǎng)絡(luò)拾取技術(shù)［3］、BP網(wǎng)格法等［4］。這類方法實(shí)現(xiàn)復(fù)雜，效率較低，工程實(shí)現(xiàn)難度大；基于地震記錄瞬時(shí)特征的方法，如極值法［5］、能量因子法等［6］。這類方法工程應(yīng)用性較強(qiáng)，但震動(dòng)信號(hào)信噪比較低時(shí)，難以準(zhǔn)確拾取。

目前時(shí)窗長(zhǎng)度采用過(guò)零點(diǎn)法或三次零點(diǎn)法［7］，由于爆炸近場(chǎng)彈性波特性與遠(yuǎn)場(chǎng)相比，波形混疊現(xiàn)象嚴(yán)重，高頻信號(hào)無(wú)法在基線上下實(shí)時(shí)穿插，因此時(shí)窗內(nèi)的樣本不能最優(yōu)地表征當(dāng)前時(shí)刻的模態(tài)信息，不具有時(shí)變特性，造成初至波的到時(shí)拾取精度低，誤判率高。針對(duì)上述問(wèn)題，本文提出了一種提取初至波到時(shí)的改進(jìn)能量因子算法。

1 能量因子算法基本原理

能量因子法是一種經(jīng)典的初至波到時(shí)拾取算法，其原理如圖1所示。設(shè)第i個(gè)時(shí)間點(diǎn)的長(zhǎng)/短時(shí)窗識(shí)別因子p（i）公式如下：

式（1）中，x（j）（j=1，2，…，N）表示震動(dòng)數(shù)據(jù)，M和N分別表示長(zhǎng)、短時(shí)間窗內(nèi)的樣本數(shù)。

圖1 能量因子法滑窗原理圖Fig.1 The schematic diagram of sliding window of the energy factorization method

為了在不等窗長(zhǎng)識(shí)別因子p（i）基礎(chǔ)上精確地確定變化點(diǎn)的位置，李勇等構(gòu)建了等窗長(zhǎng)識(shí)別因子q（i）［8］

式（2）中，x（j）（j=1，2，…，N）表示震動(dòng)數(shù)據(jù)，N表示等時(shí)間窗內(nèi)的樣本數(shù)。

結(jié)合特征函數(shù)及這兩種窗口能量比的優(yōu)點(diǎn)，給出初至波檢測(cè)因子R（i）公式

式（3）中，CF（i）（i=1，2，…，N）表示震動(dòng)數(shù)據(jù)對(duì)應(yīng)的特征函數(shù)，M和N分別表示長(zhǎng)、短時(shí)間窗內(nèi)的樣本數(shù)。

通過(guò)逐個(gè)采樣點(diǎn)滑窗的方式，求取初至波檢測(cè)因子峰值對(duì)應(yīng)的采樣點(diǎn)即為初至波到達(dá)時(shí)間。

2 改進(jìn)的能量因子算法

為了獲取更為精準(zhǔn)的初至波到時(shí)，本文對(duì)公式（3）進(jìn)行了改進(jìn)，利用高階統(tǒng)計(jì)量增強(qiáng)突變點(diǎn)的特征信息，通過(guò)S變換將時(shí)窗長(zhǎng)度自適應(yīng)調(diào)整到時(shí)窗內(nèi)信號(hào)的最小周期，從而增加初至波檢測(cè)因子的時(shí)變特性，提高初至波提取的精度。初至波檢測(cè)因子R（i）'由三個(gè)因子組成

其中P（i）'為不等窗長(zhǎng)下瞬時(shí)能量識(shí)別因子，Q（i）'為等窗長(zhǎng)下瞬時(shí)能量識(shí)別因子，K（i）為對(duì)應(yīng)瞬時(shí)時(shí)窗內(nèi)四階統(tǒng)計(jì)量峰度。

1）不等窗長(zhǎng)下瞬時(shí)能量識(shí)別因子

在式（1）不等窗長(zhǎng)識(shí)別因子的基礎(chǔ)上，將時(shí)窗長(zhǎng)度自適應(yīng)調(diào)整到當(dāng)前時(shí)刻對(duì)應(yīng)的瞬時(shí)優(yōu)勢(shì)頻率的最小周期，將三軸合成的能量信息作為特征函數(shù)輸入式（1），得到了不等窗長(zhǎng)下瞬時(shí)能量識(shí)別因子P（i）'

其中長(zhǎng)窗長(zhǎng)M=5N，短窗長(zhǎng)N=k/Ωi（t），k為權(quán)值，Ωi（t）為三軸震動(dòng)信號(hào)經(jīng)廣義S變換后第i時(shí)刻的瞬時(shí)優(yōu)勢(shì)頻率，它是指經(jīng)廣義S變換得到的時(shí)頻譜中，每一瞬時(shí)時(shí)刻所對(duì)應(yīng)的寬頻譜中，能量最大的頻率點(diǎn)。在最優(yōu)阻抗匹配條件下，在炸點(diǎn)近場(chǎng)，初至波到達(dá)時(shí)三軸的瞬時(shí)優(yōu)勢(shì)頻率一致［11］：

2）等窗長(zhǎng)下瞬時(shí)能量識(shí)別因子

在式（2）等窗長(zhǎng)識(shí)別因子的基礎(chǔ)上，同上將時(shí)窗長(zhǎng)度自適應(yīng)調(diào)整到當(dāng)前時(shí)刻對(duì)應(yīng)的瞬時(shí)優(yōu)勢(shì)頻率的最小周期，將三軸合成的能量信息作為特征函數(shù)輸入式（2）中得到了等窗長(zhǎng)下瞬時(shí)能量識(shí)別因子Q（i）'

其中等窗長(zhǎng)N=k/Ωi（t），Ωi（t）為三軸震動(dòng)信號(hào)經(jīng)廣義S變換后第i時(shí)刻的瞬時(shí)優(yōu)勢(shì)頻率，其中X（i）為三維能量合成后的數(shù)據(jù)。

3）高階統(tǒng)計(jì)量因子

為了加強(qiáng)初至波到時(shí)的波形突變特征，采用四階統(tǒng)計(jì)量峰度［12］構(gòu)建識(shí)別因子K（i）

其中X（i）為時(shí)刻i對(duì)應(yīng)的瞬時(shí)時(shí)窗內(nèi)的三分量信號(hào)能量合成數(shù)據(jù)。

改進(jìn)算法的流程圖如圖2所示。

主要的算法步驟：

1）利用廣義S變化獲取瞬時(shí)優(yōu)勢(shì)頻率，設(shè)定滑窗步長(zhǎng)Δt；

2）將節(jié)點(diǎn)獲取的三分量信號(hào)進(jìn)行能量合成；

3）計(jì)算當(dāng)前時(shí)刻的長(zhǎng)、短時(shí)窗長(zhǎng)度；

4）根據(jù)式（5）計(jì)算不等窗長(zhǎng)識(shí)別因子P（i）'；

5）根據(jù)式（6）計(jì)算等窗長(zhǎng)識(shí)別因子Q（i）'；

6）根據(jù)式（7）計(jì)算峰度Q（i）'；

7）根據(jù)式（4）計(jì)算當(dāng)前時(shí)刻初至波識(shí)別因子R；

8）X（i）以Δt進(jìn)行滑窗；

9）判斷是否滑窗結(jié)束，如結(jié)束執(zhí)行10），否則轉(zhuǎn)向3）；

10）計(jì)算R最大值對(duì)應(yīng)的時(shí)間點(diǎn)，確定達(dá)到時(shí)間，運(yùn)行結(jié)束。

圖2 算法流程圖Fig.2 Algorithm flow chart

3 應(yīng)用及效果分析

為了驗(yàn)證本算法的可行性，進(jìn)行了地下人工爆破試驗(yàn)，將3 kg TNT作為震源埋設(shè)在地下5 m采用中北大學(xué)信息探測(cè)與處理技術(shù)研究所研發(fā)的全向震動(dòng)傳感器獲取三維震動(dòng)信號(hào)。采樣率為20 k Hz，采樣時(shí)間為：10 s。圖3為地下炸點(diǎn)起爆后，獲取的局部三分向震動(dòng)信號(hào)圖。利用廣義S變換，求取X軸的時(shí)頻圖（如圖4（a））及瞬時(shí)優(yōu)勢(shì)頻率譜（如圖4（b））所示。

圖3 三分向時(shí)域信號(hào)圖Fig.3 The time domain graph of signal of three-component

圖4 自適應(yīng)瞬時(shí)窗圖Fig.4 The graph of adaptive instantaneous window

由圖4可知，通過(guò)分析廣義S變換求取的時(shí)頻域，可以明顯地區(qū)分初至波區(qū)、S波區(qū)和面波區(qū)。通過(guò)計(jì)算每一點(diǎn)的瞬時(shí)優(yōu)勢(shì)頻率，得到了瞬時(shí)優(yōu)勢(shì)頻率譜。由圖4（b）可知，瞬時(shí)優(yōu)勢(shì)頻率譜可以有效地反映非平穩(wěn)信號(hào)的時(shí)變特性，利用瞬時(shí)優(yōu)勢(shì)頻率可以構(gòu)建自適應(yīng)時(shí)窗。

設(shè)置滑窗步進(jìn)為Δt=1/F s，特征函數(shù)CF（i）= x（i）2+x（i）·x（i+1），按照式（3），得到了基于三次零點(diǎn)法的初至波識(shí)別因子；同時(shí)按照?qǐng)D2拾取流程，得到了改進(jìn)后的初至波識(shí)別因子，如圖5所示。

圖5 識(shí)別因子圖Fig.5 The graph of identification factors

表1 初至波到時(shí)相關(guān)參數(shù)表Tab.1 The relevant parameter table of first-break wave

結(jié)合表1和圖5可得，改進(jìn)后能量因子法，識(shí)別因子曲線上升速度更快，峰值更高，初至波時(shí)間提取更加準(zhǔn)確。由于增加了廣義S變換和高階統(tǒng)計(jì)量因子，因此程序運(yùn)行速度慢了一倍。但對(duì)于地下震源定位而言，20 s的處理速度可以滿足后續(xù)定位的要求。

4 結(jié)論

本文提出了提取初至波到時(shí)的改進(jìn)能量因子算法，該算法引入高階統(tǒng)計(jì)量因子提高波形幅值突變特征，利用廣義S變換對(duì)應(yīng)的瞬時(shí)優(yōu)勢(shì)頻率構(gòu)建時(shí)窗，使時(shí)窗長(zhǎng)度自適應(yīng)調(diào)整到時(shí)窗內(nèi)信號(hào)的最小周期，增強(qiáng)初至波檢測(cè)因子的時(shí)變特性。試驗(yàn)結(jié)果表明，改進(jìn)后的能量因子算法，提取直達(dá)初至波到時(shí)誤檢率降低，精度提高，在地下震源定位、層析成像等地下空間領(lǐng)域具有一定的工程應(yīng)用價(jià)值。但由于本算法增加了廣義S變換和高階統(tǒng)計(jì)量因子，使程序運(yùn)算速度降低，在保證直達(dá)初至波提取精度的前提下，如何優(yōu)化算法，降低算法的復(fù)雜度是下一步研究的重點(diǎn)。

［1］陳愛(ài)萍，鄒文，何光明，等.基于分維和相關(guān)性的自動(dòng)初至拾取技術(shù)及應(yīng)用［J］.物探化探計(jì)算技術(shù)，2009，31（2）：100-108.

［2］王彩霞，白超英，王馨.地震震相初至自動(dòng)檢測(cè)技術(shù)綜述［J］.地球物理學(xué)進(jìn)展，2013（10）：2363-2366.

［3］徐鈺，段衛(wèi)星.高精度初至自動(dòng)拾取綜合方法研究，物探與化探，2010，34（5）：596-599.

［4］王金峰，羅省賢.BP神經(jīng)網(wǎng)絡(luò)的改進(jìn)及其在初至波拾取中的應(yīng)用［J］.物探化探計(jì)算技術(shù)，2006，28（1）：67-69.

［5］劉希強(qiáng).用于三分向記錄震相識(shí)別的小波變換法［J］.地震研究，2000（3）：125-131.

［6］左國(guó)平，王彥春，隋榮亮.利用能量比法拾取地震初至的一種改進(jìn)方法［J］.石油物探，2004，43（4）.346-349.

［7］張偉，王彥春，李洪臣，等.地震道瞬時(shí)強(qiáng)度比法拾取初至波［J］.地球物理學(xué)進(jìn)展，2009，24（1）：110-114.

［8］張軍華，趙勇.用小波變換與能量比方法聯(lián)合拾取初至波［J］.物探化探計(jì)算技術(shù)，2002，11（5）：310-313.

［9］李洪林，宋玲玲，潘樹(shù)林.時(shí)窗方法與改進(jìn)的能量比方法聯(lián)合拾取初至波比［J］.內(nèi)蒙古石油化工，2007（8）：380-384.

［10］潘樹(shù)林，高磊，陳輝，等.可控震源地震記錄初至拾取方法研究［J］.石油物探，2010（2）：210-214.

［11］詹毅，鐘本善.利用小波變換提高地震波初至拾取的精確度［J］.成都理工大學(xué)學(xué)報(bào)（自然科學(xué)版），2004（6）：704-708.

［12］程仲平，陳鷹鵬，何超群.空變時(shí)窗約束地震波初至的拾取方法［J］.地球物理學(xué)進(jìn)展，2011，26（5）：1626-1627.

Modified Energy Factor Algorithm to Extract First-break Wave Arrival Time

LI Jian1，HAN Yan1，XU lina2，MA liangxing3
（1.Institute of Signal Capture&Processing Technology，North University Of China，Taiyuan 030051 2.Test And Measuring Academy of Norinco Group，Huayin 714200 3.Northern shanxi machinery manufacturing co.，LTD，Taiyuan 030009）

Considering the problem that the extraction accuracy of first-break wave arrival time is low，and the rate of false is high，a modified method based on generalized S transform and higher order statistics was proposed to effectively extract direct first break arrival in this paper.The algorithm introduced higher-order statistics to improve the mutation characteristics of the waveform amplitude，used synthesis information of three-component energy as a characteristic function，and took advantage of the instantaneous dominant frequency corresponding to generalized Stransform to build an adaptive time window and to strengthen the time-varying properties of first break detection factor.Experimental verification showed that，compared with the energy factor algorithm，the proposed adaptive instantaneous energy factor method reduced the false detection rate and improved the accuracy when extracting the primary wave.Besides，it could improve the precision of underground source localization and had wide application in the underground source location，dimensional physical field reconstruction.

first-break wave；distributed source location；adaptive window；S transform

P315.6；TN911.6

A

1008-1194（2015）05-0094-04

2015-03-30

國(guó)家973計(jì)劃前期研究專項(xiàng)基金資助（2011CB311804）；國(guó)家自然科學(xué)基金科學(xué)儀器基礎(chǔ)專項(xiàng)基金資助（61227003）

李劍（1985—），男，山西太原人，博士，研究方向：新型傳感技術(shù)及多維信息處理。E-mail：liiian851208 @126.com