高 磊,羅關(guān)鳳,劉 蕩,閔 帆,2*
(1.西南石油大學(xué)計(jì)算機(jī)科學(xué)學(xué)院,成都 610500;2.西南石油大學(xué)人工智能研究所,成都 610500)
在地震數(shù)據(jù)處理過程中,確定初至波的到達(dá)時(shí)間是非常重要的,它將直接影響后續(xù)步驟的精度,例如,動(dòng)校正[1]、靜校正[2]和速度分析[3]等。除此之外,錯(cuò)誤的初至波拾取結(jié)果也會(huì)對(duì)后續(xù)子波的分析造成干擾,使專家不能正確地進(jìn)行事件識(shí)別、震源反演,以及震源判斷等。然而,在強(qiáng)烈的背景噪聲和復(fù)雜近地表?xiàng)l件地干擾下,現(xiàn)有算法的準(zhǔn)確率會(huì)降低,如何準(zhǔn)確地拾取初至波是本領(lǐng)域目前亟待解決的問題。
早期的初至波拾取由專家根據(jù)經(jīng)驗(yàn)手動(dòng)完成,這種方式過分依賴于專家經(jīng)驗(yàn),耗時(shí)且主觀[4]。后來發(fā)展為半自動(dòng)拾取,由專家先拾取一部分初至波,然后通過軟件確定剩下的初至波,這種方式雖然在一定程度上提升了效率,但仍不能滿足實(shí)際需求。隨著自動(dòng)化拾取算法開始興起,出現(xiàn)了能量比法、相關(guān)法、聚類法等各種初至波自動(dòng)拾取算法,其中聚類法是目前較為流行的初至波拾取算法。本文通過分析各種聚類算法的優(yōu)缺點(diǎn),選取了k均值(k-means)聚類和基于密度的噪聲應(yīng)用空間聚類(Density-Based Spatial Clustering of Applications with Noise,DBSCAN)這兩種算法來拾取初至波。
由于DBSCAN 對(duì)噪聲不敏感,聚類效果好,但對(duì)于龐大的地震數(shù)據(jù),它的時(shí)間復(fù)雜度高;而k-means 對(duì)噪聲比較敏感,但其時(shí)間復(fù)雜度低。因此,本文提出將這兩種聚類相結(jié)合,先用k-means 技術(shù)找到初至波簇;再用DBSCAN 技術(shù)在初至波簇中拾取初始初至波。這種方式充分利用了兩種聚類算法的優(yōu)點(diǎn),DBSCAN 不直接在整炮數(shù)據(jù)上運(yùn)行,而是在初至波簇中進(jìn)行拾取,這種方式可以提高運(yùn)行速度,準(zhǔn)確率也得到了很大的提升。
本文將DBSCAN 用于初至波拾取,取得了較好的效果,主要工作包括以下幾個(gè)方面:
1)基于k-means 技術(shù)選取初至波簇;
2)基于DBSCAN 技術(shù)在初至波簇中拾取初始初至波;
3)通過局部線性回歸來補(bǔ)齊初始初至波中的缺失值;
4)通過能量比值對(duì)初始初至波中的錯(cuò)誤值進(jìn)行調(diào)整。
初至波拾取在地震勘探領(lǐng)域中發(fā)揮著非常重要的作用,是地震數(shù)據(jù)分析的預(yù)處理階段,被廣泛地應(yīng)用到各領(lǐng)域中,具有重要的研究意義。接下來本文將對(duì)初至波拾取的發(fā)展歷程進(jìn)行一個(gè)簡單的介紹。
目前,國內(nèi)外學(xué)者進(jìn)行了很多研究,提出了多種初至波拾取方法,本文主要從三個(gè)方面進(jìn)行介紹:
1)以單道為基礎(chǔ)的能量比法,主要以Coppens 算法[5]及其變體為代表。Coppens 使用兩個(gè)窗口內(nèi)的能量比來處理數(shù)據(jù);Sabbione 和Velis[6]改進(jìn)了Coppens 算法,以熵和分形維數(shù)算法來拾取初至波;在此基礎(chǔ)上,Allen[7]利用短期平均值與長期平均值之比(Short-Term Average/Long-Term Average,STA/LTA)來提高準(zhǔn)確率;Lee 等[8]用多個(gè)窗口代替了傳統(tǒng)的兩個(gè)窗口以獲得更可靠的結(jié)果。
2)以多道為基礎(chǔ)的互相關(guān)法。Peraldi 等[9]考慮了相鄰道之間的能量,提出了互相關(guān)算法;Senkaya 等[10]進(jìn)一步考慮了每一道數(shù)據(jù)和源小波之間的互相關(guān);Molyneux 等[11]提出使用最大互相關(guān)值作為評(píng)價(jià)準(zhǔn)則的直接相關(guān)算法來拾取初至波。
3)以整個(gè)數(shù)據(jù)集為基礎(chǔ)的聚類算法。Zhu 等[12]提出了一種基于聚類的初至波拾取算法;周竹生等[13]提出了基于加權(quán)k均值聚類的多屬性初至拾取算法;Gao 等[14]將滑動(dòng)窗口與模糊c均值聚類結(jié)合,采用滑動(dòng)窗口獲得范圍帶,模糊c均值只需在范圍帶里面檢測初至波,可以獲得更加精確的效果。
這些算法在高質(zhì)量數(shù)據(jù)集中能夠取得較好的效果,但在信噪比(Signal-to-Noise Ratio,SNR)較低的數(shù)據(jù)集中,都有一定的缺陷。例如以單道為基礎(chǔ)的算法未考慮初至波之間的關(guān)系,拾取精度不穩(wěn)定[15];以多道為基礎(chǔ)的互相關(guān)法在相關(guān)性較弱的地震數(shù)據(jù)中,準(zhǔn)確率有所下降;以整個(gè)數(shù)據(jù)集為基礎(chǔ)的大多數(shù)聚類算法易受到噪聲干擾,聚類效果差,但這種算法考慮了整個(gè)數(shù)據(jù)的分布,并且可以根據(jù)初至波的特性將其最大限度地劃分為一簇。因此,本文將繼續(xù)沿用基于聚類的初至波拾取方法,提出了基于聚類和局部線性回歸的初至波自動(dòng)拾取算法(First-arrival automatic Picking algorithm based on Clustering and Local linear regression,F(xiàn)PCL),取得了較好的效果。
在本章中,將建立一個(gè)數(shù)據(jù)模型,并介紹初至波拾取問題及其相關(guān)定義。
在地震勘探中,野外數(shù)據(jù)的采集主要是通過炮點(diǎn)激發(fā),并由不同的檢波器接收。其中,每個(gè)檢波器接收到的波形記錄對(duì)應(yīng)一道地震數(shù)據(jù),所有地震道的數(shù)據(jù)集合組成了共炮點(diǎn)道集。本文采用共炮點(diǎn)道集作為數(shù)據(jù)集來進(jìn)行初至波拾取。
圖1 是輸入的原始數(shù)據(jù)和拾取結(jié)果。它表示一個(gè)包含了210 道數(shù)據(jù),每道數(shù)據(jù)含有700 個(gè)樣本點(diǎn)的共炮點(diǎn)道集,水平坐標(biāo)表示地震道的數(shù)量,垂直坐標(biāo)表示每道的樣本數(shù)量。對(duì)于每一道數(shù)據(jù),都能找到一個(gè)位于波峰的初至波,其中,白點(diǎn)是每道的初至波位置。
圖1 共炮點(diǎn)道集和初至波Fig.1 Common shot gather and first-arrivals
在地震勘探中,由炮點(diǎn)激發(fā)并最先到達(dá)檢波器的地震波稱為初至波[16]。初至波拾取是地震數(shù)據(jù)分析的重要預(yù)處理階段,它可以定義如下:
問題1 初至波拾取問題。
輸入 共炮點(diǎn)道集A=[aij]是一個(gè)m×n的矩陣,n是地震道數(shù),m是每道的樣本數(shù),aij是第j道的第i個(gè)樣本能量值。
輸出 初至波位置數(shù)組F=[f1,f2,…,fn],fj為第j道的初至波位置。
定義1初至波簇是一個(gè)l×3 的矩陣D=[di′j′],其中,l=|Sk*|代表初至波簇的樣本點(diǎn)個(gè)數(shù),di′1=aij,(i,j)∈Sk*代表能量值,di′2存儲(chǔ)行下標(biāo)i,di′3存儲(chǔ)列下標(biāo)j。
定義2分割線Zk=,其中Zk代表第k簇的分割線,
FPCL 算法由預(yù)拾取和微調(diào)兩個(gè)階段來實(shí)現(xiàn)。預(yù)拾取階段通過k-means 和DBSCAN 技術(shù)來找到初始初至波。由于DBSCAN 具有傳遞性且抗噪能力強(qiáng),適合用于拾取初至波,但它直接在整炮數(shù)據(jù)上運(yùn)行時(shí)間復(fù)雜度較高,因此,本文利用k-means 時(shí)間復(fù)雜度低的優(yōu)勢,先在整炮數(shù)據(jù)上采用k-means 技術(shù)找到初至波簇,再采用DBSCAN 技術(shù)在初至波簇中拾取初始初至波,這樣的方式可以提高抗噪能力和運(yùn)行速度。微調(diào)階段通過局部線性回歸和能量比值最小化技術(shù)調(diào)整初始初至波以提高準(zhǔn)確率。通常來說,初至波具有相關(guān)性強(qiáng)、能量強(qiáng)的特點(diǎn)[13]。根據(jù)相鄰道初至波之間相關(guān)性強(qiáng)的特點(diǎn),使用局部線性回歸來補(bǔ)齊初始初至波中的缺失道;根據(jù)初至波能量強(qiáng)的特點(diǎn),使用能量比值最小化技術(shù)來調(diào)整錯(cuò)誤值。
本節(jié)詳細(xì)描述了基于k-means 和DBSCAN 技術(shù)的初至波預(yù)拾取過程。k-means 技術(shù)用于選擇初至波簇,DBSCAN 技術(shù)用于從簇中拾取初始初至波。
3.1.1 基于k-means選擇初至波簇
由于k-means 算法的時(shí)間復(fù)雜度接近于線性,它被用來選擇初至波所在的簇。對(duì)于給定的樣本集A=[aij],按照樣本之間的距離大小將樣本集劃分為e1個(gè)簇。該算法的目標(biāo)是迭代更新聚類中心ck以最小化平方誤差E:
其中:e1是聚類中心個(gè)數(shù);aij是第j道的第i個(gè)樣本能量值;ck是第k個(gè)聚類中心。ck通過如式(2)進(jìn)行迭代更新:
其中:t表示當(dāng)前迭代次數(shù);表示迭代更新后的第k個(gè)簇的中心;表示當(dāng)前迭代時(shí)屬于第k個(gè)簇的樣本索引集合。的計(jì)算公式為:
k-means 聚類結(jié)束后,選擇平均能量值最大的簇作為初至波簇,目標(biāo)函數(shù)設(shè)計(jì)如下:
其中:Sk表示迭代終止后屬于第k個(gè)簇的樣本索引集合。
算法1 描述了基于k-means 選取初至波簇的過程。第1)行將輸入得數(shù)據(jù)集歸一化到[-1,1];第2)行根據(jù)式(1)~(3)執(zhí)行k-means 聚類;第3)行根據(jù)式(4)選擇平均能量值最大的簇;第4)行根據(jù)定義1 初始化含有三個(gè)屬性(能量值、行索引和列索引)的初至波簇D;第5)行返回初至波簇D。
3.1.2 基于DBSCAN拾取初始初至波
首先,將初至波簇D中的樣本點(diǎn)依次按照下列公式標(biāo)記為核心點(diǎn)、邊界點(diǎn)或噪聲點(diǎn):
其中:Pi表示第i個(gè)點(diǎn)的類別,“1”表示核心點(diǎn),“0”表示邊界點(diǎn),“-1”表示噪聲點(diǎn);minPts表示鄰域內(nèi)的最小點(diǎn)個(gè)數(shù);Ni表示第i個(gè)點(diǎn)的鄰域。Ni的計(jì)算公式為:
其中:Eps表示鄰域半徑。
然后,循環(huán)遍歷所有的核心點(diǎn),連通與該核心點(diǎn)距離小于半徑Eps的所有核心點(diǎn)與邊界點(diǎn),并形成一個(gè)簇。
DBSCAN 聚類后,數(shù)據(jù)被分成e2個(gè)簇,對(duì)于每個(gè)簇,從每道中提取出第一個(gè)屬于該簇的樣本點(diǎn)并形成一條分割線,因此可以得到e2條分割線。
最后,選擇最完整的一條分割線作為初始初至波。目標(biāo)函數(shù)設(shè)計(jì)如下:
算法2 描述了基于DBSCAN 技術(shù)拾取初至波的過程。第1)行根據(jù)式(5)~(6)執(zhí)行DBSCAN 聚類;第2)~6)行根據(jù)定義2 計(jì)算每個(gè)簇的分割線;第7)~8)行根據(jù)式(7)選擇最完整的分割線作為初始初至波;第9)行返回找到的初始初至波F。
本節(jié)詳細(xì)描述了采用局部線性回歸和能量比值最小化技術(shù)的初至波微調(diào)階段。局部線性回歸技術(shù)用于補(bǔ)齊缺失值,能量比值最小化技術(shù)用于調(diào)整錯(cuò)誤值。
3.2.1 局部線性回歸補(bǔ)齊缺失值
由于找到的初始初至波在部分道上可能存在缺失值,在本節(jié)中,根據(jù)相鄰道初至波相關(guān)性強(qiáng)的特點(diǎn)設(shè)計(jì)了一個(gè)局部線性回歸模型來補(bǔ)齊缺失值,回歸方程定義如下:
其中:α和β是回歸系數(shù)。
該算法的目標(biāo)是最小化真實(shí)值y和預(yù)測值之間的差值,建立損失函數(shù)模型:
根據(jù)一階最優(yōu)性條件,將式(9)分別對(duì)α和β求偏導(dǎo),并令其為零,得到回歸系數(shù):
通過回歸模型,可以計(jì)算出回歸系數(shù)α和β的值,從而計(jì)算出缺失道所對(duì)應(yīng)的預(yù)測值。最后,將得到的預(yù)測值作為該道的初至波。
算法3 描述了通過局部線性回歸補(bǔ)齊缺失值的過程。第1)行初始化F′;第2)行定義了循環(huán)范圍;第3)行判斷當(dāng)前道是否存在缺失值;第4)行定義了長度為2×τ的數(shù)組x和y;第5)~9)行將缺失值前后各τ道鄰居初至波所對(duì)應(yīng)的地震道數(shù)和初至波位置分別賦值給x和y;第10)~11)行根據(jù)式(10)和式(11)計(jì)算回歸系數(shù)α和β的值;第12)行根據(jù)式(8)計(jì)算第j道的預(yù)測 值;第13)行更新第j道的值;第16)行返回補(bǔ)齊缺失值后的初至波F′。
3.2.2 能量比值最小化調(diào)整位置
由于初至波的能量值通常較大,在本節(jié)中,設(shè)定閾值ε作為能量值的約束條件,采用能量比值最小化技術(shù)對(duì)小于ε的初至波位置在初至波范圍內(nèi)進(jìn)行調(diào)整。根據(jù)初至波的特性,設(shè)定如下的優(yōu)化模型:
其中:E(low,high,j)=代表第j道的第low~high個(gè)樣本點(diǎn)能量之和;λ是初至波范圍;是第j道的初至波范圍;θ是能量比值的窗口長度;ε是能量閾值。
算法4 描述了用能量比值最小化技術(shù)對(duì)初至波能量值小于ε的位置進(jìn)行調(diào)整。第1)行定義了循環(huán)范圍;第2)~12)行是調(diào)整過程,第2)行判斷初至波能量值是否達(dá)到調(diào)整條件,小于ε才調(diào)整;第3)行初始化目標(biāo)函數(shù)的最小值;第4)行是在第j道的初至波附近滑動(dòng);第5)行根據(jù)式(12)計(jì)算能量比值;第6)行判斷目標(biāo)函數(shù)值和約束是否達(dá)到更新條件;第7)~8)行更新目標(biāo)函數(shù)值和最小能量比下標(biāo)i*;第11)行更新調(diào)整后的初至波位置;第14)行返回調(diào)整后的初至波F*。
為了驗(yàn)證算法的有效性,本文將FPCL 與改進(jìn)的能量比(Improved Modified Energy Ratio,IMER)法[8]、互相關(guān)技術(shù)(Cross Correlation Technique,CCT)[10]、基于模糊C均值聚類的微震數(shù)據(jù)自動(dòng)時(shí)間選取算法(Automatic time Picking for microseismic data based on a FuzzyC-means clustering algorithm,APF)[12]、基于兩階段優(yōu)化的初至波自動(dòng)拾取算法(First-arrival automatic Picking algorithm based on Two-stage Optimization,F(xiàn)PTO)[16]共四種算法在兩個(gè)地震數(shù)據(jù)集上進(jìn)行測試。其中,IMER 是以單道為基礎(chǔ)的能量比法,CCT 是以多道為基礎(chǔ)的互相關(guān)法,APF 是以整個(gè)數(shù)據(jù)集為基礎(chǔ)的聚類算法,F(xiàn)PTO 是兩階段優(yōu)化算法。實(shí)驗(yàn)結(jié)果表明FPCL 更準(zhǔn)確。
表1 列出了實(shí)驗(yàn)所用的兩個(gè)數(shù)據(jù)集的相關(guān)信息。
表1 地震數(shù)據(jù)集信息Tab.1 Seismic dataset information
為研究參數(shù)對(duì)算法性能的影響,本文分別測試聚類中心個(gè)數(shù)e1、鄰居通道數(shù)τ、初至波范圍λ和能量比值窗口長度θ對(duì)各數(shù)據(jù)集準(zhǔn)確率的影響。該實(shí)驗(yàn)基于單變量方法進(jìn)行測試,即固定其他參數(shù)來研究一個(gè)變量對(duì)結(jié)果的影響,實(shí)驗(yàn)結(jié)果展示在圖2 中。
圖2 參數(shù)對(duì)各數(shù)據(jù)集準(zhǔn)確率的影響Fig.2 Influence of parameters on accuracy of each dataset
根據(jù)圖2(a)的實(shí)驗(yàn)結(jié)果,聚類中心個(gè)數(shù)e1=5.00 是最佳設(shè)置;圖2(b)的結(jié)果顯示鄰居通道數(shù)τ對(duì)實(shí)驗(yàn)結(jié)果的影響不大,本文將τ設(shè)置為地震道數(shù)的5.00%;λ和θ與每道的數(shù)據(jù)規(guī)模有關(guān),根據(jù)圖2(c)和(d)所示,λ的最佳設(shè)置為數(shù)據(jù)規(guī)模的4.00%~6.00%,θ的最佳設(shè)置為數(shù)據(jù)規(guī)模的1.00%;能量閾值ε由幾個(gè)實(shí)際地震數(shù)據(jù)的初至波計(jì)算得到,選取實(shí)際地震數(shù)據(jù)的初至波,發(fā)現(xiàn)初至波的能量值通常都大于0.20,因此本文將ε=0.20 設(shè)置為初至波能量值的約束條件。
圖3 描述了在data1 上的初至波拾取過程。圖3(b)表示基于k-means 聚類后選擇的初至波簇,圖中的點(diǎn)表示簇中元素,右上角的是噪聲,和初始數(shù)據(jù)集相比,此時(shí)的數(shù)據(jù)量從300 200 減少到14 335。圖3(c)表示基于DBSCAN 在圖(b)中拾取的初始初至波,紅色的點(diǎn)表示初至波位置,從圖中可以看到,它的拾取結(jié)果對(duì)噪聲不敏感,但部分道存在缺失值和錯(cuò)誤值。圖3(d)表示用局部線性回歸和能量比值最小化技術(shù)來補(bǔ)齊缺失值和調(diào)整錯(cuò)誤值后的初至波。
圖3 在data1上的初至波拾取過程Fig.3 First-arrival picking process on dataset 1
圖4 描述了在data2 上的初至波拾取過程。圖4(b)表示基于k-means 聚類后選擇的初至波簇,圖中的點(diǎn)表示簇中元素,左上角的點(diǎn)是噪聲,和初始數(shù)據(jù)集相比,數(shù)據(jù)量大幅度減少(從600 400 減少到27 020)。圖4(c)表示DBSCAN 技術(shù)在初至波簇中拾取的初始初至波,紅色的點(diǎn)表示初至波位置,從中可看出拾取結(jié)果對(duì)噪聲不敏感,但部分道還存在缺失值和錯(cuò)誤值。圖4(d)表示用局部線性回歸和能量比值最小化來補(bǔ)齊缺失值和調(diào)整錯(cuò)誤值后的初至波。
圖4 在data2上的初至波拾取過程Fig.4 First-arrival picking process on dataset 2
將FPCL 與IMER[8]、CCT[10]、APF[12]和FPTO[16]四種算法進(jìn)行對(duì)比:IMER 用多個(gè)窗口代替了傳統(tǒng)的兩個(gè)窗口來拾取初至波;CCT 通過互相關(guān)技術(shù)來檢測初至波;APF 使用模糊C均值聚類來區(qū)分初至波和噪聲;FPTO 采用兩階段優(yōu)化拾取,先用滑動(dòng)窗口找出初至波范圍帶,然后改進(jìn)了能量比率算法在范圍帶中拾取初至波。
圖5 展示了五種算法在兩個(gè)數(shù)據(jù)集上的拾取結(jié)果,黑色方框內(nèi)為正確的初至波拾取結(jié)果。圖5(a)表示在data1 上的拾取結(jié)果,由于APF 的抗噪能力較弱,因此在強(qiáng)噪聲的情況下,拾取結(jié)果會(huì)早于初至波;除此之外,和初至波相比,CCT獲得的結(jié)果會(huì)有所延遲。圖5(b)表示在data2 上的拾取結(jié)果,在該數(shù)據(jù)集中,CCT 的偏差較大;IMER 和APF 只在較少道上存在偏差;FPTO 受噪聲干擾,拾取結(jié)果不準(zhǔn)確。
圖5 五種算法的拾取結(jié)果Fig.5 Picking results of five methods
本文在兩個(gè)地震數(shù)據(jù)集上計(jì)算了這些算法的準(zhǔn)確性,實(shí)驗(yàn)結(jié)果的準(zhǔn)確率對(duì)比如表2 所示。從表2 可看出:FPCL 與IMER 法相比,準(zhǔn)確率分別提升了4.00 個(gè)百分點(diǎn)和3.50 個(gè)百分點(diǎn);與CCT 相比,準(zhǔn)確率分別提升了38.00 個(gè)百分點(diǎn)和10.25 個(gè)百分點(diǎn);與APF 相比,準(zhǔn)確率分別提升了34.50 個(gè)百分點(diǎn)和3.50 個(gè)百分點(diǎn);與FPTO 相比,準(zhǔn)確率分別提升了5.50 個(gè)百分點(diǎn)和16.25 個(gè)百分點(diǎn)。上述實(shí)驗(yàn)結(jié)果表明FPCL的準(zhǔn)確率更高,找到的初至波位置更加準(zhǔn)確。
表2 不同方法在兩個(gè)數(shù)據(jù)集上的準(zhǔn)確率對(duì)比 單位:%Tab.2 Accuracy comparison of different methods on two data sets unit:%
目前,拾取初至波算法受到背景噪聲和復(fù)雜近地表?xiàng)l件的干擾,拾取精度會(huì)降低。本文提出了由兩階段組成的FPCL。由于DBSCAN 具有傳遞性且抗噪能力強(qiáng),適合拾取初至波,但直接在整炮數(shù)據(jù)上運(yùn)行其時(shí)間復(fù)雜度較高,因此,本文利用k-means 時(shí)間復(fù)雜度低的優(yōu)勢,預(yù)拾取階段先基于k-means 找到初至波簇,再基于DBSCAN 在初至波簇中進(jìn)行拾取,通過結(jié)合這兩種聚類算法的優(yōu)點(diǎn)來選取初始初至波,這種方式可以提高抗噪能力和運(yùn)行速度。微調(diào)階段通過局部線性回歸和能量比值最小化技術(shù)來補(bǔ)齊缺失值和調(diào)整錯(cuò)誤值,可以進(jìn)一步提高準(zhǔn)確率。
在兩個(gè)數(shù)據(jù)集上的實(shí)驗(yàn)結(jié)果表明FPCL 可以有效地拾取初至波。將來,會(huì)考慮將FPCL 運(yùn)用到其他方面,例如信號(hào)處理和地震數(shù)據(jù)去噪。