基于近鄰傳播聚類與曲線擬合的斷層識(shí)別

葉 濤，陳 雷

(1.青海民族大學(xué) 計(jì)算機(jī)學(xué)院，青海 西寧 810007；2.中國科學(xué)院地質(zhì)與地球物理研究所 中國科學(xué)院油氣資源研究重點(diǎn)實(shí)驗(yàn)室，北京 100029)

0 引 言

斷層識(shí)別是地震解釋中的重要環(huán)節(jié)，對(duì)油氣資源勘探具有重要意義。為克服人工識(shí)別方法較為費(fèi)時(shí)的不足，自動(dòng)識(shí)別斷層的方法獲得了廣泛的研究。

相干體算法[1](可分為三代：C1～C3)是最早出現(xiàn)的斷層自動(dòng)識(shí)別方法。它基于地震道之間的相關(guān)性對(duì)斷層進(jìn)行識(shí)別，但由于相關(guān)窗口的影響，所得斷層的均值效應(yīng)較為明顯。Wang等通過Hough變換對(duì)相干體中直線的檢測(cè)對(duì)斷層進(jìn)行了識(shí)別[2-4]，不足是對(duì)同時(shí)存在的多個(gè)斷層的識(shí)別能力較差。Wang等先將實(shí)數(shù)地震道轉(zhuǎn)換為復(fù)數(shù)形式，再基于復(fù)數(shù)形式的地震道利用相關(guān)方法識(shí)別斷層[5]，不足是耗時(shí)較高，且有較多的層位信息殘留。Di等[6]通過對(duì)地震層位在三維空間中的彎曲屬性(方向和大小等)對(duì)小斷層進(jìn)行了識(shí)別，不足是計(jì)算復(fù)雜度較高。Yan等基于蟻群算法對(duì)斷層進(jìn)行了識(shí)別[7,8]。這種方法可將微小的斷層識(shí)別出來，但所得斷層結(jié)果的可信度不高。

斷層識(shí)別方法中，人工方法最準(zhǔn)確。故，最好的斷層識(shí)別方法是用智能算法模仿人工方法的流程。模仿人工方法時(shí)層位不連續(xù)點(diǎn)的確定比較容易，基于層位不連續(xù)點(diǎn)的斷層的生成是智能算法模仿人工方法成敗的關(guān)鍵。

近鄰傳播聚類[9-11]是一種聚類算法，它可在未知聚類數(shù)目的情況下根據(jù)樣本間的相似性對(duì)待聚類樣本聚類?；诖朔椒?，本文提出了一種基于近鄰傳播聚類與曲線擬合的斷層識(shí)別方法。先采用連通區(qū)域標(biāo)注方法確定層位不連續(xù)點(diǎn)，然后用近鄰傳播聚類算法對(duì)不連續(xù)點(diǎn)聚類，最后，用曲線擬合方法對(duì)每類不連續(xù)點(diǎn)進(jìn)行擬合，模仿人工實(shí)現(xiàn)了對(duì)斷層的識(shí)別。

1 近鄰傳播聚類算法

近鄰傳播聚類算法將所有數(shù)據(jù)點(diǎn)作為潛在的聚類中心，通過消息傳遞機(jī)制迭代地對(duì)數(shù)據(jù)點(diǎn)進(jìn)行聚類并對(duì)聚類中心進(jìn)行更新，直到聚類中心和對(duì)各數(shù)據(jù)點(diǎn)的聚類不在變化為止。

1.1 相似度與參考度初始化

相似度和參考度是近鄰傳播聚類算法中兩個(gè)重要的輸入?yún)?shù)，它們?cè)谝欢ǔ潭壬蠜Q定了最終的聚類結(jié)果，其初始化方法具體如下。

1.1.1 相似度矩陣的初始化

相似度矩陣中的元素表示相應(yīng)樣本點(diǎn)之間的相似性大小，這種相似性一般用歐式距離的負(fù)值來表示。比如，待聚類樣本集Ω中的m和n兩個(gè)樣本點(diǎn)間的相似性的計(jì)算如式(1)所示

(1)

這樣，兩點(diǎn)距離越遠(yuǎn)，相似性越小，屬于同一類的可能性越?。环粗瑢儆谕活惖目赡苄栽酱?。

1.1.2 參考度的初始化

參考度決定了數(shù)據(jù)點(diǎn)可否作為聚類中心，對(duì)最終所得的聚類數(shù)目有重要影響。例如，當(dāng)對(duì)圖1(a)中所示的數(shù)據(jù)進(jìn)行近鄰傳播聚類時(shí)(用歐式距離定義相似性)，最終所得聚類數(shù)目隨參考度的變化曲線如圖1(b)所示。

圖1 參考度與聚類數(shù)目之間的關(guān)系

1.2 聚 類

采用近鄰傳播聚類算法對(duì)數(shù)據(jù)點(diǎn)進(jìn)行聚類時(shí)，初始化好相似度和參考度之后，數(shù)據(jù)點(diǎn)和聚類中心之間傳遞吸引度和歸屬度消息，以此對(duì)吸引度、歸屬度和聚類中心進(jìn)行不斷的迭代更新，同時(shí)對(duì)數(shù)據(jù)點(diǎn)聚類，直至結(jié)果穩(wěn)定。吸引度和歸屬度的傳遞方式如圖2所示。

圖2 吸引度和歸屬度傳遞

其中，r(i,k)為k點(diǎn)對(duì)i點(diǎn)的吸引度，反映了點(diǎn)k作為點(diǎn)i的聚類中心的合適程度，k′為競(jìng)爭(zhēng)的候選聚類中心。a(i,k)為i點(diǎn)對(duì)k點(diǎn)的歸屬度，表示點(diǎn)i選擇點(diǎn)k為聚類中心的合適程度，r(i′,k)表示其它點(diǎn)i′對(duì)點(diǎn)k作為聚類中心的支持程度。

迭代過程中，歸屬度、吸引度和自歸屬度的更新規(guī)則分別如式(2)～式(4)所示

s(i,k′)})+λ×r(t-1)(i,k)

(2)

max0,r(t)(i′,k)}+λ×a(t-1)(i,k)

(3)

λ×a(t-1)(k,k)

(4)

其中，λ為阻尼系數(shù)，一般在0和1之間。

式(2)右端最大函數(shù)部分表示候選聚類中心k′比聚類中心k更好的證據(jù)。吸引度更新的意義是使所有聚類中心都參與競(jìng)爭(zhēng)一個(gè)待聚類數(shù)據(jù)點(diǎn)的歸屬。

式(3)中，r(k,k)表示點(diǎn)k的自吸引度，r(k,k)為負(fù)值時(shí)表示點(diǎn)k不適合作為一個(gè)聚類中心。含最大函數(shù)的求和項(xiàng)表示候選聚類中心k從除i、k之外的點(diǎn)接收到的總的正吸引度。若某些點(diǎn)以點(diǎn)k作為聚類中心時(shí)的吸引度值為正，那么這些點(diǎn)以點(diǎn)k為聚類中心時(shí)的歸屬度值就會(huì)增大。這樣，a(i,k)也表示了作為候選聚類中心的點(diǎn)k從所有數(shù)據(jù)點(diǎn)(除點(diǎn)i)接收的吸引度之和，體現(xiàn)了點(diǎn)k是否是一個(gè)好的聚類中心。

式(4)中，a(k,k)為自歸屬度，是點(diǎn)k作為一個(gè)聚類中心的證據(jù)總和。

每次迭代完成之后，近鄰傳播聚類算法對(duì)聚類中心進(jìn)行更新并對(duì)數(shù)據(jù)點(diǎn)分類，規(guī)則如式(5)所示

(i,k) = argmax {a(i,k) +r(i,k)}

(5)

以式(5)所得結(jié)果(i,k)作為聚類依據(jù)時(shí)：如果i=k，則將點(diǎn)k作為新的聚類中心；否則，將點(diǎn)i歸屬到以點(diǎn)k為聚類中心的類別中去，實(shí)現(xiàn)對(duì)點(diǎn)i的聚類；直至聚類中心和每個(gè)點(diǎn)所歸屬的類別穩(wěn)定，近鄰傳播聚類終止。

2 本文基于近鄰傳播聚類與曲線擬合的斷層識(shí)別

本文所提斷層識(shí)別方法的具體實(shí)現(xiàn)過程如圖3所示。

圖3 本文方法流程

其中，讀入地震數(shù)據(jù)部分完成地震數(shù)據(jù)格式的轉(zhuǎn)換，將地震數(shù)據(jù)轉(zhuǎn)換成二值圖像的工作在預(yù)處理部分完成，確定層位不連續(xù)點(diǎn)的確定工作在確定層位不連續(xù)點(diǎn)部分通過連通區(qū)域標(biāo)注[12,13]方法來完成，對(duì)層位不連續(xù)點(diǎn)的聚類工作在近鄰傳播聚類部分完成，曲線擬合部分對(duì)每一類層位不連續(xù)點(diǎn)進(jìn)行曲線擬合，獲得斷層，最后，將斷層識(shí)別結(jié)果輸出。

2.1 地震層位不連續(xù)點(diǎn)確定

2.1.1 地震層位提取

地震數(shù)據(jù)中正層位與負(fù)層位交替出現(xiàn)，若只保留地震剖面中的正層位部分，斷層的數(shù)量、形狀、位置等信息并不會(huì)受到影響，所以，本文只保留正層位(下文中的層位均指正層位)。地震層位的提取過程如下：

(1)將浮點(diǎn)型二維地震數(shù)據(jù)映射成二進(jìn)值圖像；

(2)對(duì)所得二值地震剖面圖像進(jìn)行形態(tài)學(xué)濾波；

(3)對(duì)濾波后的二值地震剖面進(jìn)行連通區(qū)域標(biāo)注，將地震層位作為連通區(qū)域提取出來。

2.1.2 地震剖面中層位不連續(xù)點(diǎn)的確定

地震剖面中層位不連續(xù)的地方即為斷層，所以斷層的位置可以通過層位不連續(xù)點(diǎn)來確定。本文將層位所對(duì)應(yīng)連通區(qū)域的橫向端點(diǎn)(包含左端點(diǎn)和右端點(diǎn))作為層位不連續(xù)點(diǎn)，其確定方法如式(6)和式(7)所示

(6)

(7)

其中，xleft表示左端點(diǎn)的行坐標(biāo)，yleft表示左端點(diǎn)的列坐標(biāo)，xright表示右端點(diǎn)的行坐標(biāo)，yright表示右端點(diǎn)的列坐標(biāo)，x表示橫向端點(diǎn)的行坐標(biāo)，y表示橫向端點(diǎn)的列坐標(biāo)。M為標(biāo)記矩陣，用于存儲(chǔ)各個(gè)連通區(qū)域的標(biāo)號(hào)。y為待優(yōu)化的目標(biāo)函數(shù)。M(x,y)=b為優(yōu)化函數(shù)的約束條件，b為當(dāng)前待確定橫向端點(diǎn)的層位所對(duì)應(yīng)的連通區(qū)域的標(biāo)號(hào)。地震剖面圖像中地震層位不連續(xù)點(diǎn)的確定方法如圖4所示。

圖4 確定地震層位不連續(xù)點(diǎn)原理

其中，H表示層位(與連通區(qū)域相對(duì)應(yīng))。L和R分別表示層位的左端點(diǎn)和右端點(diǎn)，即層位不連續(xù)點(diǎn)。

2.2 層位不連續(xù)點(diǎn)聚類

要正確地識(shí)別斷層，需要對(duì)層位不連續(xù)點(diǎn)聚類以基于同類的層位不連續(xù)點(diǎn)來生成斷層。本文采用近鄰傳播聚類算法，根據(jù)確定同一條斷層的層位不連續(xù)點(diǎn)之間的相似性對(duì)其進(jìn)行聚類，用同類層位不連續(xù)點(diǎn)確定一條斷層。在聚類時(shí)，本文對(duì)相似度矩陣與參考度的初始化具體如下：

鑒于確定同一條斷層的層位不連續(xù)點(diǎn)之間水平梯度差和歐氏距離均較小，本文用層位不連續(xù)點(diǎn)之間水平梯度差和歐氏距離加權(quán)和的負(fù)值來定義相似度矩陣，如式(8)所示

ω2*Gy(i)-Gyj}

(8)

其中，s為相似度矩陣，i和j分別表示不連續(xù)點(diǎn)中的第i和第j個(gè)點(diǎn)。xi表示第i個(gè)點(diǎn)的行坐標(biāo)，xj表示第j個(gè)點(diǎn)的行坐標(biāo)，yi和yj分別相應(yīng)的列坐標(biāo)。ω1表示相似性中歐式距離的權(quán)重，ω2表示水平梯度差的權(quán)重。Gy(i)表示第i個(gè)點(diǎn)的水平梯度，Gyj表示第j個(gè)點(diǎn)的水平梯度。由式(8)可知，兩個(gè)層位不連續(xù)點(diǎn)之間相似性大小與歐式距離和水平梯度差均為反比關(guān)系。

在對(duì)參考度進(jìn)行初始化時(shí)，本文用相似度矩陣各列的均值作為參考度的初始值，以使每一個(gè)層位不連續(xù)點(diǎn)都可以成為聚類中心。參考度的初始化方法具體如式(9)所示

(9)

其中，p為參考度向量，p(k)表示用于判斷第k個(gè)數(shù)據(jù)點(diǎn)是否可以作為聚類中心的參考度，η為常數(shù)，一般取10左右，N為需要聚類的層位不連續(xù)點(diǎn)的個(gè)數(shù)，s表示相似度矩陣。近鄰傳播聚類過程中，以s(k,k)作為參考度，s(k,k)越大，第k個(gè)點(diǎn)越有可能成為聚類中心。

通過以上步驟就實(shí)現(xiàn)了地震剖面中層位不連續(xù)點(diǎn)的確定以及聚類，接下來就可以基于聚類后的不連續(xù)點(diǎn)進(jìn)行斷層的生成了。

2.3 斷層最小二乘曲線擬合

采用近鄰傳播聚類算法對(duì)層位不連續(xù)點(diǎn)聚類之后，每一類層位不連續(xù)點(diǎn)可確定一條斷層。本文提出了用最小二乘曲線擬合來確定斷層的方法。

最小二乘曲線擬合是根據(jù)已知數(shù)據(jù)找出擬合多項(xiàng)式函數(shù)的系數(shù)，使得擬合出的函數(shù)值與原離散點(diǎn)的誤差平方和最小。設(shè)N個(gè)點(diǎn)的x坐標(biāo)(自變量)分別為x1,x2,…,xN，相應(yīng)的坐標(biāo)(函數(shù)值)分別為y1,y2,…,yN,k次擬合多項(xiàng)式函數(shù)為：p=a0+a1x+a2x+…+akxk，則擬合多項(xiàng)式函數(shù)值與誤差平方和可以表示為

(10)

其中，pi為擬合多項(xiàng)式函數(shù)在第i個(gè)自變量xi處的函數(shù)值。最小二乘擬合的目標(biāo)是找出最佳的一組系數(shù)a0,a1,…,ak，使誤差平方和δ2最小，即通過下式找出p

(11)

將式(10)右邊各項(xiàng)以依次對(duì)ai,i=1,2,…,k求偏導(dǎo)數(shù)并化簡(jiǎn)，可得

(12)

式(12)又可表示為

X×A=Y

(13)

由此可得

(14)

這樣就得到了系數(shù)ai,i=1,2,…,k，也就得到了k次最佳擬合多項(xiàng)式函數(shù)的解析表達(dá)式，從而可以獲得各個(gè)擬合多項(xiàng)式函數(shù)值pi,i=1,2,…,N。

對(duì)于同一類別中的層位不連續(xù)點(diǎn)，將每個(gè)點(diǎn)的行坐標(biāo)作為自變量，列坐標(biāo)作為函數(shù)值，利用最小二乘擬合即可得出與原有離散的層位不連續(xù)點(diǎn)誤差平方和最小的連續(xù)曲線，即斷層。

3 實(shí)驗(yàn)結(jié)果和分析

為驗(yàn)證所提基于近鄰傳播聚類與曲線擬合的斷層識(shí)別方法的有效性，將其與人工斷層識(shí)別方法、經(jīng)典的C3方法和參考文獻(xiàn)[7]中所提的方法在大小為220×290×100的實(shí)際地震數(shù)據(jù)上進(jìn)行了仿真實(shí)驗(yàn)對(duì)比。計(jì)算機(jī)配置：系統(tǒng)：windows7，內(nèi)存大小：4 G；實(shí)驗(yàn)工具：Matlab2015a。

各方法在某實(shí)際地震剖面上所得的斷層識(shí)別結(jié)果的對(duì)比如圖5中所示。

圖5 本文方法與其它方法識(shí)別所得斷層結(jié)果對(duì)比

由圖5可知，人工方法識(shí)別所得的斷層最準(zhǔn)確，約與實(shí)際情況一致；由于均值效應(yīng)的存在，C3方法所得的斷層較為模糊，僅可以看出斷層的大體形狀和位置；文獻(xiàn)[7]提出的基于方向復(fù)值相關(guān)屬性的方法，在準(zhǔn)確度方面比C3方法都有了提高，斷層位置更加清晰，但有較多的層位信息殘留下來；本文所提基于近鄰傳播聚類與曲線擬合的斷層識(shí)別方法所得斷層的準(zhǔn)確度優(yōu)于C3方法及文獻(xiàn)[7]中所提的方法，從本文所提斷層識(shí)別方法的可行性得以驗(yàn)證。

為了從客觀指標(biāo)方面進(jìn)一步驗(yàn)證本文方法相對(duì)于其它方法的優(yōu)越性，本文以采用人工方法識(shí)別出的斷層為基準(zhǔn)對(duì)各方法進(jìn)行了客觀指標(biāo)對(duì)比，這些客觀指標(biāo)包括峰值信噪比(即：PSNR)、時(shí)間消耗和正確率，實(shí)驗(yàn)中采用了100幅地震剖面。對(duì)比結(jié)果如圖6～圖8所示。表1是各種方法的平均客觀指標(biāo)對(duì)比。

由圖6所示的不同方法所得識(shí)別結(jié)果的信噪比對(duì)比可知：文獻(xiàn)[7]中所提方法所得結(jié)果的信噪比最低，原因斷層結(jié)果中有較多的層位信息殘留；C3方法所得結(jié)果的信噪比也較低，與文獻(xiàn)[7]中方法相近，原因是相關(guān)性計(jì)算所帶來的均值效應(yīng)。

圖6 各方法所得識(shí)別結(jié)果的峰值信噪比對(duì)比

圖7中所示是采用不同方法對(duì)實(shí)驗(yàn)中的100幀地震剖面數(shù)據(jù)進(jìn)行斷層識(shí)別時(shí)的時(shí)間消耗對(duì)比；由圖7中所示的時(shí)間消耗曲線可知：C3方法時(shí)間消耗最多，這是因?yàn)榈卣鹌拭嬷忻恳粋€(gè)地震樣點(diǎn)都需要計(jì)算其所在地震道與立體窗口內(nèi)其它各地震道之間的相關(guān)系數(shù)，計(jì)算量較大；文獻(xiàn)[7]中方法時(shí)間消耗相對(duì)于C3有所降低，但依然較大，原因是它雖然沒有像C3那樣計(jì)算立體窗口內(nèi)各個(gè)地震道與目標(biāo)樣點(diǎn)所在地震道之間的相關(guān)系數(shù)，但仍然計(jì)算了立體窗口內(nèi)與水平方向成0°、45°、90°以及135°4個(gè)方向上各個(gè)復(fù)值地震道與目標(biāo)樣點(diǎn)所在復(fù)值地震道之間的相關(guān)系數(shù)；人工方法的時(shí)間消耗變化較大，這與人的主觀性有關(guān)；本文方法的時(shí)間消耗最小，原因是本文方法去除了較為耗時(shí)的求相關(guān)計(jì)算，時(shí)間消耗主要是對(duì)少數(shù)的層位不連續(xù)點(diǎn)進(jìn)行聚類和擬合。這就從時(shí)間消耗方面驗(yàn)證了本文所提方法的實(shí)時(shí)性。

圖7 各方法時(shí)間消耗對(duì)比

圖8 本文方法與人工方法識(shí)別所得的斷層數(shù)量比較

方法峰值信噪比斷層數(shù)量正確率平均時(shí)間消耗/s人工方法+∞100%29.9924C3方法5.7489×56.8401文獻(xiàn)[7]方法5.5783×45.6155本文方法19.886191%1.7365

圖8中所示是本文方法與手工方法在每幀剖面上所得的斷層數(shù)對(duì)比?？芍?，本文基于近鄰傳播聚類與曲線擬合的方法從每幀地震剖面中識(shí)別出的斷層數(shù)與人工方法識(shí)別所得的數(shù)目基本吻合，從而驗(yàn)證了其識(shí)別正確率。

表1是采用不同方法所得斷層結(jié)果的平均客觀指標(biāo)對(duì)比。因?yàn)镻SNR和識(shí)別正確率都是以人工方法所得斷層識(shí)別結(jié)果作為標(biāo)準(zhǔn)的，所以人工方法所得斷層結(jié)果的平均PSNR為正無窮大，平均識(shí)別正確率為100%；因?yàn)镃3方法、文獻(xiàn)[7]方法不能提供量化的斷層識(shí)別結(jié)果，所以在表1中用“×”表示它們的平均正確率。表1中各客觀指標(biāo)的對(duì)比進(jìn)一步驗(yàn)證了本文所提方法的有效性和先進(jìn)性。

4 結(jié)束語

針對(duì)現(xiàn)有斷層識(shí)別方法所存在的準(zhǔn)確性差、時(shí)間消耗大的問題，基于近鄰傳播聚類和曲線擬合，本文提出了一種新斷層識(shí)別方法。文中先是介紹和分析了傳統(tǒng)斷層識(shí)別方法的優(yōu)缺點(diǎn)；接著給出了本文所提方法的原理及實(shí)現(xiàn)步驟；最后，將本文方法與現(xiàn)在常用的人工方法以及基于方向復(fù)值相關(guān)屬性的方法在實(shí)際地震數(shù)據(jù)上進(jìn)行了實(shí)驗(yàn)對(duì)比，在所得斷層識(shí)別結(jié)果的客觀指標(biāo)等方面驗(yàn)證了本文所提方法的可行性及實(shí)時(shí)性。如何定義相似度矩陣和設(shè)置參考度的初始值以使近鄰傳播聚類算法最終所得的聚類數(shù)量與人工方法所得出的斷層條數(shù)更加一致，是本文進(jìn)一步研究的工作。