基于聚類和局部線性回歸的初至波自動拾取算法

高 磊，羅關(guān)鳳，劉 蕩，閔 帆，2*

（1.西南石油大學(xué)計算機科學(xué)學(xué)院，成都 610500；2.西南石油大學(xué)人工智能研究所，成都 610500）

0 引言

在地震數(shù)據(jù)處理過程中，確定初至波的到達時間是非常重要的，它將直接影響后續(xù)步驟的精度，例如，動校正［1］、靜校正［2］和速度分析［3］等。除此之外，錯誤的初至波拾取結(jié)果也會對后續(xù)子波的分析造成干擾，使專家不能正確地進行事件識別、震源反演，以及震源判斷等。然而，在強烈的背景噪聲和復(fù)雜近地表條件地干擾下，現(xiàn)有算法的準(zhǔn)確率會降低，如何準(zhǔn)確地拾取初至波是本領(lǐng)域目前亟待解決的問題。

早期的初至波拾取由專家根據(jù)經(jīng)驗手動完成，這種方式過分依賴于專家經(jīng)驗，耗時且主觀［4］。后來發(fā)展為半自動拾取，由專家先拾取一部分初至波，然后通過軟件確定剩下的初至波，這種方式雖然在一定程度上提升了效率，但仍不能滿足實際需求。隨著自動化拾取算法開始興起，出現(xiàn)了能量比法、相關(guān)法、聚類法等各種初至波自動拾取算法，其中聚類法是目前較為流行的初至波拾取算法。本文通過分析各種聚類算法的優(yōu)缺點，選取了k均值（k-means）聚類和基于密度的噪聲應(yīng)用空間聚類（Density-Based Spatial Clustering of Applications with Noise，DBSCAN）這兩種算法來拾取初至波。

由于DBSCAN 對噪聲不敏感，聚類效果好，但對于龐大的地震數(shù)據(jù)，它的時間復(fù)雜度高；而k-means 對噪聲比較敏感，但其時間復(fù)雜度低。因此，本文提出將這兩種聚類相結(jié)合，先用k-means 技術(shù)找到初至波簇；再用DBSCAN 技術(shù)在初至波簇中拾取初始初至波。這種方式充分利用了兩種聚類算法的優(yōu)點，DBSCAN 不直接在整炮數(shù)據(jù)上運行，而是在初至波簇中進行拾取，這種方式可以提高運行速度，準(zhǔn)確率也得到了很大的提升。

本文將DBSCAN 用于初至波拾取，取得了較好的效果，主要工作包括以下幾個方面：

1）基于k-means 技術(shù)選取初至波簇；

2）基于DBSCAN 技術(shù)在初至波簇中拾取初始初至波；

3）通過局部線性回歸來補齊初始初至波中的缺失值；

4）通過能量比值對初始初至波中的錯誤值進行調(diào)整。

1 相關(guān)工作

初至波拾取在地震勘探領(lǐng)域中發(fā)揮著非常重要的作用，是地震數(shù)據(jù)分析的預(yù)處理階段，被廣泛地應(yīng)用到各領(lǐng)域中，具有重要的研究意義。接下來本文將對初至波拾取的發(fā)展歷程進行一個簡單的介紹。

目前，國內(nèi)外學(xué)者進行了很多研究，提出了多種初至波拾取方法，本文主要從三個方面進行介紹：

1）以單道為基礎(chǔ)的能量比法，主要以Coppens 算法［5］及其變體為代表。Coppens 使用兩個窗口內(nèi)的能量比來處理數(shù)據(jù)；Sabbione 和Velis［6］改進了Coppens 算法，以熵和分形維數(shù)算法來拾取初至波；在此基礎(chǔ)上，Allen［7］利用短期平均值與長期平均值之比（Short-Term Average/Long-Term Average，STA/LTA）來提高準(zhǔn)確率；Lee 等［8］用多個窗口代替了傳統(tǒng)的兩個窗口以獲得更可靠的結(jié)果。

2）以多道為基礎(chǔ)的互相關(guān)法。Peraldi 等［9］考慮了相鄰道之間的能量，提出了互相關(guān)算法；Senkaya 等［10］進一步考慮了每一道數(shù)據(jù)和源小波之間的互相關(guān)；Molyneux 等［11］提出使用最大互相關(guān)值作為評價準(zhǔn)則的直接相關(guān)算法來拾取初至波。

3）以整個數(shù)據(jù)集為基礎(chǔ)的聚類算法。Zhu 等［12］提出了一種基于聚類的初至波拾取算法；周竹生等［13］提出了基于加權(quán)k均值聚類的多屬性初至拾取算法；Gao 等［14］將滑動窗口與模糊c均值聚類結(jié)合，采用滑動窗口獲得范圍帶，模糊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)確率有所下降；以整個數(shù)據(jù)集為基礎(chǔ)的大多數(shù)聚類算法易受到噪聲干擾，聚類效果差，但這種算法考慮了整個數(shù)據(jù)的分布，并且可以根據(jù)初至波的特性將其最大限度地劃分為一簇。因此，本文將繼續(xù)沿用基于聚類的初至波拾取方法，提出了基于聚類和局部線性回歸的初至波自動拾取算法（First-arrival automatic Picking algorithm based on Clustering and Local linear regression，F(xiàn)PCL），取得了較好的效果。

2 問題及定義

在本章中，將建立一個數(shù)據(jù)模型，并介紹初至波拾取問題及其相關(guān)定義。

2.1 數(shù)據(jù)模型

在地震勘探中，野外數(shù)據(jù)的采集主要是通過炮點激發(fā)，并由不同的檢波器接收。其中，每個檢波器接收到的波形記錄對應(yīng)一道地震數(shù)據(jù)，所有地震道的數(shù)據(jù)集合組成了共炮點道集。本文采用共炮點道集作為數(shù)據(jù)集來進行初至波拾取。

圖1 是輸入的原始數(shù)據(jù)和拾取結(jié)果。它表示一個包含了210 道數(shù)據(jù)，每道數(shù)據(jù)含有700 個樣本點的共炮點道集，水平坐標(biāo)表示地震道的數(shù)量，垂直坐標(biāo)表示每道的樣本數(shù)量。對于每一道數(shù)據(jù)，都能找到一個位于波峰的初至波，其中，白點是每道的初至波位置。

圖1 共炮點道集和初至波Fig.1 Common shot gather and first-arrivals

2.2 初至波拾取問題及其相關(guān)定義

在地震勘探中，由炮點激發(fā)并最先到達檢波器的地震波稱為初至波［16］。初至波拾取是地震數(shù)據(jù)分析的重要預(yù)處理階段，它可以定義如下：

問題1 初至波拾取問題。

輸入 共炮點道集A=[aij]是一個m×n的矩陣，n是地震道數(shù)，m是每道的樣本數(shù)，aij是第j道的第i個樣本能量值。

輸出 初至波位置數(shù)組F=[f1，f2，…，fn]，fj為第j道的初至波位置。

定義1初至波簇是一個l×3 的矩陣D=[di′j′]，其中，l=|Sk*|代表初至波簇的樣本點個數(shù)，di′1=aij，(i，j)∈Sk*代表能量值，di′2存儲行下標(biāo)i，di′3存儲列下標(biāo)j。

定義2分割線Zk=，其中Zk代表第k簇的分割線，

3 本文算法

FPCL 算法由預(yù)拾取和微調(diào)兩個階段來實現(xiàn)。預(yù)拾取階段通過k-means 和DBSCAN 技術(shù)來找到初始初至波。由于DBSCAN 具有傳遞性且抗噪能力強，適合用于拾取初至波，但它直接在整炮數(shù)據(jù)上運行時間復(fù)雜度較高，因此，本文利用k-means 時間復(fù)雜度低的優(yōu)勢，先在整炮數(shù)據(jù)上采用k-means 技術(shù)找到初至波簇，再采用DBSCAN 技術(shù)在初至波簇中拾取初始初至波，這樣的方式可以提高抗噪能力和運行速度。微調(diào)階段通過局部線性回歸和能量比值最小化技術(shù)調(diào)整初始初至波以提高準(zhǔn)確率。通常來說，初至波具有相關(guān)性強、能量強的特點［13］。根據(jù)相鄰道初至波之間相關(guān)性強的特點，使用局部線性回歸來補齊初始初至波中的缺失道；根據(jù)初至波能量強的特點，使用能量比值最小化技術(shù)來調(diào)整錯誤值。

3.1 預(yù)拾取階段

本節(jié)詳細描述了基于k-means 和DBSCAN 技術(shù)的初至波預(yù)拾取過程。k-means 技術(shù)用于選擇初至波簇，DBSCAN 技術(shù)用于從簇中拾取初始初至波。

3.1.1 基于k-means選擇初至波簇

由于k-means 算法的時間復(fù)雜度接近于線性，它被用來選擇初至波所在的簇。對于給定的樣本集A=[aij]，按照樣本之間的距離大小將樣本集劃分為e1個簇。該算法的目標(biāo)是迭代更新聚類中心ck以最小化平方誤差E：

其中：e1是聚類中心個數(shù)；aij是第j道的第i個樣本能量值；ck是第k個聚類中心。ck通過如式（2）進行迭代更新：

其中：t表示當(dāng)前迭代次數(shù)；表示迭代更新后的第k個簇的中心；表示當(dāng)前迭代時屬于第k個簇的樣本索引集合。的計算公式為：

k-means 聚類結(jié)束后，選擇平均能量值最大的簇作為初至波簇，目標(biāo)函數(shù)設(shè)計如下：

其中：Sk表示迭代終止后屬于第k個簇的樣本索引集合。

算法1 描述了基于k-means 選取初至波簇的過程。第1）行將輸入得數(shù)據(jù)集歸一化到［-1，1］；第2）行根據(jù)式（1）～（3）執(zhí)行k-means 聚類；第3）行根據(jù)式（4）選擇平均能量值最大的簇；第4）行根據(jù)定義1 初始化含有三個屬性（能量值、行索引和列索引）的初至波簇D；第5）行返回初至波簇D。

3.1.2 基于DBSCAN拾取初始初至波

首先，將初至波簇D中的樣本點依次按照下列公式標(biāo)記為核心點、邊界點或噪聲點：

其中：Pi表示第i個點的類別，“1”表示核心點，“0”表示邊界點，“-1”表示噪聲點；minPts表示鄰域內(nèi)的最小點個數(shù)；Ni表示第i個點的鄰域。Ni的計算公式為：

其中：Eps表示鄰域半徑。

然后，循環(huán)遍歷所有的核心點，連通與該核心點距離小于半徑Eps的所有核心點與邊界點，并形成一個簇。

DBSCAN 聚類后，數(shù)據(jù)被分成e2個簇，對于每個簇，從每道中提取出第一個屬于該簇的樣本點并形成一條分割線，因此可以得到e2條分割線。

最后，選擇最完整的一條分割線作為初始初至波。目標(biāo)函數(shù)設(shè)計如下：

算法2 描述了基于DBSCAN 技術(shù)拾取初至波的過程。第1）行根據(jù)式（5）～（6）執(zhí)行DBSCAN 聚類；第2）～6）行根據(jù)定義2 計算每個簇的分割線；第7）～8）行根據(jù)式（7）選擇最完整的分割線作為初始初至波；第9）行返回找到的初始初至波F。

3.2 微調(diào)階段

本節(jié)詳細描述了采用局部線性回歸和能量比值最小化技術(shù)的初至波微調(diào)階段。局部線性回歸技術(shù)用于補齊缺失值，能量比值最小化技術(shù)用于調(diào)整錯誤值。

3.2.1 局部線性回歸補齊缺失值

由于找到的初始初至波在部分道上可能存在缺失值，在本節(jié)中，根據(jù)相鄰道初至波相關(guān)性強的特點設(shè)計了一個局部線性回歸模型來補齊缺失值，回歸方程定義如下：

其中：α和β是回歸系數(shù)。

該算法的目標(biāo)是最小化真實值y和預(yù)測值之間的差值，建立損失函數(shù)模型：

根據(jù)一階最優(yōu)性條件，將式（9）分別對α和β求偏導(dǎo)，并令其為零，得到回歸系數(shù)：

通過回歸模型，可以計算出回歸系數(shù)α和β的值，從而計算出缺失道所對應(yīng)的預(yù)測值。最后，將得到的預(yù)測值作為該道的初至波。

算法3 描述了通過局部線性回歸補齊缺失值的過程。第1）行初始化F′；第2）行定義了循環(huán)范圍；第3）行判斷當(dāng)前道是否存在缺失值；第4）行定義了長度為2×τ的數(shù)組x和y；第5）～9）行將缺失值前后各τ道鄰居初至波所對應(yīng)的地震道數(shù)和初至波位置分別賦值給x和y；第10）～11）行根據(jù)式（10）和式（11）計算回歸系數(shù)α和β的值；第12）行根據(jù)式（8）計算第j道的預(yù)測 值；第13）行更新第j道的值；第16）行返回補齊缺失值后的初至波F′。

3.2.2 能量比值最小化調(diào)整位置

由于初至波的能量值通常較大，在本節(jié)中，設(shè)定閾值ε作為能量值的約束條件，采用能量比值最小化技術(shù)對小于ε的初至波位置在初至波范圍內(nèi)進行調(diào)整。根據(jù)初至波的特性，設(shè)定如下的優(yōu)化模型：

其中：E(low，high，j)=代表第j道的第low～high個樣本點能量之和；λ是初至波范圍；是第j道的初至波范圍；θ是能量比值的窗口長度；ε是能量閾值。

算法4 描述了用能量比值最小化技術(shù)對初至波能量值小于ε的位置進行調(diào)整。第1）行定義了循環(huán)范圍；第2）～12）行是調(diào)整過程，第2）行判斷初至波能量值是否達到調(diào)整條件，小于ε才調(diào)整；第3）行初始化目標(biāo)函數(shù)的最小值；第4）行是在第j道的初至波附近滑動；第5）行根據(jù)式（12）計算能量比值；第6）行判斷目標(biāo)函數(shù)值和約束是否達到更新條件；第7）～8）行更新目標(biāo)函數(shù)值和最小能量比下標(biāo)i*；第11）行更新調(diào)整后的初至波位置；第14）行返回調(diào)整后的初至波F*。

4 實驗與結(jié)果分析

為了驗證算法的有效性，本文將FPCL 與改進的能量比（Improved Modified Energy Ratio，IMER）法［8］、互相關(guān)技術(shù)（Cross Correlation Technique，CCT）［10］、基于模糊C均值聚類的微震數(shù)據(jù)自動時間選取算法（Automatic time Picking for microseismic data based on a FuzzyC-means clustering algorithm，APF）［12］、基于兩階段優(yōu)化的初至波自動拾取算法（First-arrival automatic Picking algorithm based on Two-stage Optimization，F(xiàn)PTO）［16］共四種算法在兩個地震數(shù)據(jù)集上進行測試。其中，IMER 是以單道為基礎(chǔ)的能量比法，CCT 是以多道為基礎(chǔ)的互相關(guān)法，APF 是以整個數(shù)據(jù)集為基礎(chǔ)的聚類算法，F(xiàn)PTO 是兩階段優(yōu)化算法。實驗結(jié)果表明FPCL 更準(zhǔn)確。

4.1 數(shù)據(jù)集

表1 列出了實驗所用的兩個數(shù)據(jù)集的相關(guān)信息。

表1 地震數(shù)據(jù)集信息Tab.1 Seismic dataset information

4.2 參數(shù)敏感性分析

為研究參數(shù)對算法性能的影響，本文分別測試聚類中心個數(shù)e1、鄰居通道數(shù)τ、初至波范圍λ和能量比值窗口長度θ對各數(shù)據(jù)集準(zhǔn)確率的影響。該實驗基于單變量方法進行測試，即固定其他參數(shù)來研究一個變量對結(jié)果的影響，實驗結(jié)果展示在圖2 中。

圖2 參數(shù)對各數(shù)據(jù)集準(zhǔn)確率的影響Fig.2 Influence of parameters on accuracy of each dataset

4.3 實驗參數(shù)設(shè)置

根據(jù)圖2（a）的實驗結(jié)果，聚類中心個數(shù)e1=5.00 是最佳設(shè)置；圖2（b）的結(jié)果顯示鄰居通道數(shù)τ對實驗結(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%；能量閾值ε由幾個實際地震數(shù)據(jù)的初至波計算得到，選取實際地震數(shù)據(jù)的初至波，發(fā)現(xiàn)初至波的能量值通常都大于0.20，因此本文將ε=0.20 設(shè)置為初至波能量值的約束條件。

4.4 實驗結(jié)果描述

圖3 描述了在data1 上的初至波拾取過程。圖3（b）表示基于k-means 聚類后選擇的初至波簇，圖中的點表示簇中元素，右上角的是噪聲，和初始數(shù)據(jù)集相比，此時的數(shù)據(jù)量從300 200 減少到14 335。圖3（c）表示基于DBSCAN 在圖（b）中拾取的初始初至波，紅色的點表示初至波位置，從圖中可以看到，它的拾取結(jié)果對噪聲不敏感，但部分道存在缺失值和錯誤值。圖3（d）表示用局部線性回歸和能量比值最小化技術(shù)來補齊缺失值和調(diào)整錯誤值后的初至波。

圖3 在data1上的初至波拾取過程Fig.3 First-arrival picking process on dataset 1

圖4 描述了在data2 上的初至波拾取過程。圖4（b）表示基于k-means 聚類后選擇的初至波簇，圖中的點表示簇中元素，左上角的點是噪聲，和初始數(shù)據(jù)集相比，數(shù)據(jù)量大幅度減少（從600 400 減少到27 020）。圖4（c）表示DBSCAN 技術(shù)在初至波簇中拾取的初始初至波，紅色的點表示初至波位置，從中可看出拾取結(jié)果對噪聲不敏感，但部分道還存在缺失值和錯誤值。圖4（d）表示用局部線性回歸和能量比值最小化來補齊缺失值和調(diào)整錯誤值后的初至波。

圖4 在data2上的初至波拾取過程Fig.4 First-arrival picking process on dataset 2

4.5 實驗結(jié)果與分析

將FPCL 與IMER［8］、CCT［10］、APF［12］和FPTO［16］四種算法進行對比：IMER 用多個窗口代替了傳統(tǒng)的兩個窗口來拾取初至波；CCT 通過互相關(guān)技術(shù)來檢測初至波；APF 使用模糊C均值聚類來區(qū)分初至波和噪聲；FPTO 采用兩階段優(yōu)化拾取，先用滑動窗口找出初至波范圍帶，然后改進了能量比率算法在范圍帶中拾取初至波。

圖5 展示了五種算法在兩個數(shù)據(jù)集上的拾取結(jié)果，黑色方框內(nèi)為正確的初至波拾取結(jié)果。圖5（a）表示在data1 上的拾取結(jié)果，由于APF 的抗噪能力較弱，因此在強噪聲的情況下，拾取結(jié)果會早于初至波；除此之外，和初至波相比，CCT獲得的結(jié)果會有所延遲。圖5（b）表示在data2 上的拾取結(jié)果，在該數(shù)據(jù)集中，CCT 的偏差較大；IMER 和APF 只在較少道上存在偏差；FPTO 受噪聲干擾，拾取結(jié)果不準(zhǔn)確。

圖5 五種算法的拾取結(jié)果Fig.5 Picking results of five methods

本文在兩個地震數(shù)據(jù)集上計算了這些算法的準(zhǔn)確性，實驗結(jié)果的準(zhǔn)確率對比如表2 所示。從表2 可看出：FPCL 與IMER 法相比，準(zhǔn)確率分別提升了4.00 個百分點和3.50 個百分點；與CCT 相比，準(zhǔn)確率分別提升了38.00 個百分點和10.25 個百分點；與APF 相比，準(zhǔn)確率分別提升了34.50 個百分點和3.50 個百分點；與FPTO 相比，準(zhǔn)確率分別提升了5.50 個百分點和16.25 個百分點。上述實驗結(jié)果表明FPCL的準(zhǔn)確率更高，找到的初至波位置更加準(zhǔn)確。

表2 不同方法在兩個數(shù)據(jù)集上的準(zhǔn)確率對比 單位：%Tab.2 Accuracy comparison of different methods on two data sets unit：%

5 結(jié)語

目前，拾取初至波算法受到背景噪聲和復(fù)雜近地表條件的干擾，拾取精度會降低。本文提出了由兩階段組成的FPCL。由于DBSCAN 具有傳遞性且抗噪能力強，適合拾取初至波，但直接在整炮數(shù)據(jù)上運行其時間復(fù)雜度較高，因此，本文利用k-means 時間復(fù)雜度低的優(yōu)勢，預(yù)拾取階段先基于k-means 找到初至波簇，再基于DBSCAN 在初至波簇中進行拾取，通過結(jié)合這兩種聚類算法的優(yōu)點來選取初始初至波，這種方式可以提高抗噪能力和運行速度。微調(diào)階段通過局部線性回歸和能量比值最小化技術(shù)來補齊缺失值和調(diào)整錯誤值，可以進一步提高準(zhǔn)確率。

在兩個數(shù)據(jù)集上的實驗結(jié)果表明FPCL 可以有效地拾取初至波。將來，會考慮將FPCL 運用到其他方面，例如信號處理和地震數(shù)據(jù)去噪。