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

葉 濤，陳 雷

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

0 引 言

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

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

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

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

1 近鄰傳播聚類算法

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

1.1 相似度與參考度初始化

相似度和參考度是近鄰傳播聚類算法中兩個(gè)重要的輸入?yún)?shù)，它們在一定程度上決定了最終的聚類結(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)可否作為聚類中心，對最終所得的聚類數(shù)目有重要影響。例如，當(dāng)對圖1(a)中所示的數(shù)據(jù)進(jìn)行近鄰傳播聚類時(shí)(用歐式距離定義相似性)，最終所得聚類數(shù)目隨參考度的變化曲線如圖1(b)所示。

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

1.2 聚 類

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

圖2 吸引度和歸屬度傳遞

其中，r(i,k)為k點(diǎn)對i點(diǎn)的吸引度，反映了點(diǎn)k作為點(diǎn)i的聚類中心的合適程度，k′為競爭的候選聚類中心。a(i,k)為i點(diǎn)對k點(diǎn)的歸屬度，表示點(diǎn)i選擇點(diǎn)k為聚類中心的合適程度，r(i′,k)表示其它點(diǎn)i′對點(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ù)。吸引度更新的意義是使所有聚類中心都參與競爭一個(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í)的歸屬度值就會增大。這樣，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ù)總和。

每次迭代完成之后，近鄰傳播聚類算法對聚類中心進(jìn)行更新并對數(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)對點(diǎn)i的聚類；直至聚類中心和每個(gè)點(diǎn)所歸屬的類別穩(wěn)定，近鄰傳播聚類終止。

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

本文所提斷層識別方法的具體實(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]方法來完成，對層位不連續(xù)點(diǎn)的聚類工作在近鄰傳播聚類部分完成，曲線擬合部分對每一類層位不連續(xù)點(diǎn)進(jìn)行曲線擬合，獲得斷層，最后，將斷層識別結(jié)果輸出。

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

2.1.1 地震層位提取

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

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

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

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

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

地震剖面中層位不連續(xù)的地方即為斷層，所以斷層的位置可以通過層位不連續(xù)點(diǎn)來確定。本文將層位所對應(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)記矩陣，用于存儲各個(gè)連通區(qū)域的標(biāo)號。y為待優(yōu)化的目標(biāo)函數(shù)。M(x,y)=b為優(yōu)化函數(shù)的約束條件，b為當(dāng)前待確定橫向端點(diǎn)的層位所對應(yīng)的連通區(qū)域的標(biāo)號。地震剖面圖像中地震層位不連續(xù)點(diǎn)的確定方法如圖4所示。

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

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

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

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

鑒于確定同一條斷層的層位不連續(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)系。

在對參考度進(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 斷層最小二乘曲線擬合

采用近鄰傳播聚類算法對層位不連續(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)以依次對ai,i=1,2,…,k求偏導(dǎo)數(shù)并化簡，可得

(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。

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

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

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

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

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

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

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

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

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

圖7中所示是采用不同方法對實(shí)驗(yàn)中的100幀地震剖面數(shù)據(jù)進(jìn)行斷層識別時(shí)的時(shí)間消耗對比；由圖7中所示的時(shí)間消耗曲線可知：C3方法時(shí)間消耗最多，這是因?yàn)榈卣鹌拭嬷忻恳粋€(gè)地震樣點(diǎn)都需要計(jì)算其所在地震道與立體窗口內(nèi)其它各地震道之間的相關(guān)系數(shù)，計(jì)算量較大；文獻(xiàn)[7]中方法時(shí)間消耗相對于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í)間消耗主要是對少數(shù)的層位不連續(xù)點(diǎn)進(jìn)行聚類和擬合。這就從時(shí)間消耗方面驗(yàn)證了本文所提方法的實(shí)時(shí)性。

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

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

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

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

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

4 結(jié)束語

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