基于Jacobi算法求解結(jié)構(gòu)張量的置信擴(kuò)散濾波方法

饒 溯， 李錄明， 劉力輝， 胡 濱， 馮 鑫

(1.中海石油國(guó)際能源服務(wù)(北京)有限公司，北京 100028；2.油氣藏地質(zhì)及開(kāi)發(fā)工程國(guó)家重點(diǎn)實(shí)驗(yàn)室(成都理工大學(xué))，成都 610059；3.北京諾克斯達(dá)石油科技有限公司，北京 100192)

0 引言

地震資料去噪在疊前和疊后資料應(yīng)用中均非常重要。壓制疊前角道集的噪聲一直是地震數(shù)據(jù)疊前處理中的重要問(wèn)題，其對(duì)疊前振幅屬性分析及應(yīng)用有重要影響(如噪聲壓制的好壞將直接影響疊前地震資料振幅的應(yīng)用(如疊前AVA反演)效果)。而壓制疊后地震資料的噪聲也將更有利于斷層斷點(diǎn)準(zhǔn)確位置的解釋、層位的自動(dòng)追蹤及構(gòu)造圈閉的準(zhǔn)確落實(shí)。

地震勘探中噪聲一般包括隨機(jī)噪聲和相干噪聲。通常相干噪聲屬低頻信號(hào)，而隨機(jī)噪聲屬高頻信號(hào)，并具有統(tǒng)計(jì)特性。噪聲的壓制通常采用濾波的方法，針對(duì)噪聲的不同特征及疊前、疊后處理，采用多種濾波方法，其中具有代表性方法包括奇異值分解(SVD)[1]、F-xy域圖像特征值濾波[2]、F-k域?yàn)V波[3]、τ-p域?yàn)V波[4]、神經(jīng)網(wǎng)絡(luò)等一些模糊濾波去噪[5]和擴(kuò)散濾波[6-9]等。

Koenderink[10]發(fā)現(xiàn)擴(kuò)散濾波等價(jià)于一個(gè)熱傳導(dǎo)的物理過(guò)程。而非線性各向異性擴(kuò)散方程(P-M方程)是由Perona等[11]把它作為圖像處理中有選擇性地保留邊界信息的工具而首先提出。P-M方程有可能是病態(tài)的， You[12]通過(guò)正則化的方法消除了這種病態(tài)性。

Weickert[13-14]對(duì)P-M方程進(jìn)行了改進(jìn)，引進(jìn)了結(jié)構(gòu)張量D，并提出最優(yōu)旋轉(zhuǎn)不變性的新型算法;孫夕平等[6]對(duì)Weickert的新型算法公式采用有限差分優(yōu)化求解并設(shè)計(jì)近似旋轉(zhuǎn)不變性模板代替常規(guī)的平滑算子；王益艷等[15]在非線性各向異性擴(kuò)散濾波基礎(chǔ)上構(gòu)建梯度和切線方向上的擴(kuò)散濾波系數(shù)，并加入了數(shù)據(jù)保真項(xiàng)。他們所提出的方法，筆者都將其歸為常規(guī)基于結(jié)構(gòu)張量的擴(kuò)散濾波方法，而其中結(jié)構(gòu)張量的設(shè)計(jì)和求解方面還存在一定的改進(jìn)空間，筆者是在常規(guī)方法的基礎(chǔ)上進(jìn)行改進(jìn)。

在非線性各向異性擴(kuò)散濾波基礎(chǔ)上，采用Jacobi算法來(lái)求解結(jié)構(gòu)張量D中的特征向量，并引入線狀置信度Cline來(lái)共同來(lái)構(gòu)建擴(kuò)散濾波系數(shù)λ，最后來(lái)求解擴(kuò)散方程，實(shí)現(xiàn)置信擴(kuò)散慮波。本文方法相比常規(guī)基于結(jié)構(gòu)張量的擴(kuò)散濾波方法，在構(gòu)建和求解結(jié)構(gòu)張量D方面及構(gòu)建擴(kuò)散濾波系數(shù)有所差異，最終效果較常規(guī)基于結(jié)構(gòu)張量的擴(kuò)散濾波方法也所提高：隨機(jī)噪聲壓制效果較好、地震同相軸連續(xù)性更強(qiáng)，斷層成像更加清晰準(zhǔn)確，明顯改善資料品質(zhì)。

1 常規(guī)方法

1.1 非線性各向異性擴(kuò)散濾波原理

非線性各向異性擴(kuò)散可模擬一個(gè)熱傳導(dǎo)的物理過(guò)程，它在沒(méi)有產(chǎn)生和損失熱量的前提下平衡溫度。將擴(kuò)散應(yīng)用于地震道集，用地震振幅代替熱量，在擴(kuò)散過(guò)程中時(shí)間是一個(gè)固定參數(shù)。擴(kuò)散進(jìn)行的時(shí)間越長(zhǎng)，地震道集濾波程度則越強(qiáng)，所以濾波過(guò)程相當(dāng)于是提高信噪比的過(guò)程。地震振幅的梯度變化相當(dāng)于熱傳導(dǎo)的通量，其擴(kuò)散方程可以表示為[13-14]：

(1)

其中:D(|▽uσ|)為結(jié)構(gòu)張量;|▽uσ|為正則化后的▽uσ，正則化過(guò)程相當(dāng)于高斯濾波;t為擴(kuò)散時(shí)間(也可計(jì)為迭代次數(shù));u0(x,y)為原始地震道集；u(x,y,t)為迭代次數(shù)為t次后的濾波地震道集。

▽uσ=▽(Kσ*u(x,y,t))

(2)

(3)

式中:σ為尺度函數(shù)，可以根據(jù)噪聲的均方誤差自適應(yīng)調(diào)節(jié)。

1.2 最優(yōu)旋轉(zhuǎn)不變算法

公式(1)可在正交坐標(biāo)系下展開(kāi)成如下形式：

(4)

其中結(jié)構(gòu)張量可表示為式(5)。

(5)

公式(4)可表示為沿著x、y兩個(gè)方向(即地震道集道號(hào)和時(shí)間方向)上進(jìn)行擴(kuò)散，而最優(yōu)旋轉(zhuǎn)不變算法是需要對(duì)公式(5)的結(jié)構(gòu)張量D進(jìn)行一定程度上的旋轉(zhuǎn)，使其沿著地震道集局部振幅值變化最快和最慢方向上擴(kuò)散，如王益艷等[15]再用擴(kuò)散系數(shù)來(lái)控制在不同變化方向上平滑程度，即常規(guī)基于結(jié)構(gòu)張量的擴(kuò)散濾波方法。

公式(4)旋轉(zhuǎn)及引入擴(kuò)散系數(shù)之后可轉(zhuǎn)換成如下形式：

(6)

式中:λ為擴(kuò)散系數(shù)，擴(kuò)散系數(shù)小則表示擴(kuò)散較慢，反之?dāng)U散系數(shù)大則表示擴(kuò)散較快；η表示地震道集局部振幅值變化較快和變化劇烈處，其擴(kuò)散程度較慢；ξ表示地震道集局部振幅值變化較慢和變化均勻處，其擴(kuò)散程度較快。

2 本文方法

2.1 結(jié)構(gòu)張量的改進(jìn)

為了增強(qiáng)同向軸連續(xù)性和提高壓制噪聲的效果，在本文方法中對(duì)結(jié)構(gòu)張量的設(shè)計(jì)進(jìn)行了改進(jìn)，并引入線狀置信度來(lái)改進(jìn)和自適應(yīng)調(diào)節(jié)擴(kuò)散系數(shù)。

該方法以公式(4)為基礎(chǔ)，并對(duì)公式(5)中的結(jié)構(gòu)張量改進(jìn)為如下形式：

(7)

指示地震剖面局部振幅值變化劇烈位置處擴(kuò)散系數(shù)較小，記為λ1，指示地震剖面局部振幅值變化均勻段擴(kuò)散系數(shù)較大，記為λ2，有如下形式:

λ1=α

(8)

(9)

(10)

最后根據(jù)特征值u、擴(kuò)散系數(shù)λ和線狀置信度Cline可以將公式(4)轉(zhuǎn)換成形式如下：

(11)

當(dāng)η值和ξ值比較接近時(shí)，從公式(11)可以看出兩個(gè)方向上擴(kuò)散速率較為一致。公式(11)中振幅劇烈變化處位置值η值必然大于振幅平穩(wěn)變化處位置值ξ，而當(dāng)ξ值遠(yuǎn)小于η值時(shí)即代表振幅平穩(wěn)位置有很快的擴(kuò)散速率。

2.2 Jacobi算法原理

考慮到結(jié)構(gòu)張量D矩陣是一個(gè)2×2階的對(duì)稱矩陣，采用2×2的對(duì)稱Schur分解的Jacobi算法來(lái)求解效率較快。Jacobi算法的思想是逐步地減小特征分量矩陣v中非對(duì)角位置上值的大小[17]，也可稱為Givens旋轉(zhuǎn)變換：

(12)

2.3 方法的實(shí)現(xiàn)步驟

公式(11)可化為離散差分形式如下：

(13)

本文方法即是對(duì)公式(13)進(jìn)行求解，完成對(duì)輸入原始地震道集u0的t次迭代運(yùn)算。主要可分為以下幾步：

1)先對(duì)原始地震道集u0進(jìn)行高斯濾波，這個(gè)高斯濾波相當(dāng)于是正則化，然后分別求道號(hào)和時(shí)間方向上的二階導(dǎo)數(shù)和混合偏導(dǎo)數(shù)v。

2)根據(jù)公式(7)中構(gòu)建特征分量v，并用Jacobi算法來(lái)求解其特征值u，再根據(jù)公式(10)計(jì)算線狀置信度Cline，將Cline代入擴(kuò)散系數(shù)公式(9)中計(jì)算出λ1和λ2。

3)將原始地震道集u0、及其在道號(hào)和時(shí)間方向上的二階導(dǎo)數(shù)λ1和λ2代入到公式(13)，即完成一次擴(kuò)散濾波計(jì)算過(guò)程，有擴(kuò)散濾波計(jì)算結(jié)果ut。

3 模型試驗(yàn)及效果分析

3.1 理論模型正演及去噪試驗(yàn)

表1給出了單層滿足四類AVA曲線的各向同性介質(zhì)理論模型，圖1(a)～圖1(c)為由表1中給出的模型組合在一個(gè)道集上的AVA正演角道集，入射角為1°～45°,AVA正演角道集采用合成記錄正演。從圖1(c)擴(kuò)散濾波方法的結(jié)果中仍然可見(jiàn)四類AVA曲線地震振幅特征，即該方法有一定程度的保幅效果，且地震同相軸也較為連續(xù)。

圖1 含四類AVA曲線的正演模型角道集Fig.1 Forward modeling of angle gather with four kinds of AVA curves(a)原始正演模型；(b)加20%隨機(jī)噪聲模型；(c)20%隨機(jī)噪聲模型去噪后

表1 四類AVA曲線的各向同性介質(zhì)理論模型

圖2 SNR和MSE隨迭代次數(shù)變化的曲線Fig.2 Curves of SNR and MSE with iteration times

3.2 去噪試驗(yàn)的信噪比及誤差分析

圖像濾波和信號(hào)降噪通常采用SNR(信噪比)和MSE(均方誤差)兩項(xiàng)指標(biāo)來(lái)進(jìn)行衡量。本文分別以含10%和20%隨機(jī)噪聲的模型去噪試驗(yàn)為例，對(duì)比分析了SNR和MSE曲線的變化。如圖2所示，隨著迭代次數(shù)的增加，SNR值合理的變大，MSE值有效的減小，在滿足一定迭代次數(shù)的時(shí)候，10%和20%噪聲下的SNR值和MSE值分別達(dá)到極值，即擴(kuò)散到一定程度，擴(kuò)散過(guò)程停止，地震振幅濾波效果將不再發(fā)生變化。

3.3 去噪試驗(yàn)的時(shí)頻分析

去噪試驗(yàn)壓制噪聲效果的好壞也可通過(guò)時(shí)頻域的響應(yīng)特性來(lái)觀察。如圖3(a)～圖3(c)所示，對(duì)比了加噪正演模型帶通濾波和本文方法去噪后的時(shí)頻圖，其中時(shí)頻圖采用廣義S變換方法計(jì)算[18]。通過(guò)對(duì)比圖3(b)和圖3(c)的壓制噪聲效果好壞，可觀察到圖3(c)中紅色方框內(nèi)時(shí)頻能量團(tuán)更為聚焦，其對(duì)應(yīng)的時(shí)間刻度更為準(zhǔn)確。圖3(b)為帶通濾波方法的結(jié)果，它僅僅壓制了頻帶范圍外的噪聲，對(duì)于有效信號(hào)頻帶范圍內(nèi)噪音沒(méi)有起到壓制作用。而本文方法可以消除部分有效信號(hào)頻帶范圍內(nèi)的噪聲，但仍有部分噪聲信號(hào)未被消除。對(duì)比圖3(b)和圖3(c)的淺層黑色方框位置，可見(jiàn)帶通濾波方法仍可見(jiàn)時(shí)頻能量團(tuán)的信息，而本文方法中弱的時(shí)頻能量已經(jīng)掩蓋在噪聲中。這原因可能歸結(jié)于擴(kuò)散濾波方法僅是在時(shí)間域上進(jìn)行相干增強(qiáng)濾波，有效的強(qiáng)地震同向軸可明顯增強(qiáng)，但要保證弱地震同向軸信息可嘗試變換域內(nèi)去除噪聲，可參考小波域擴(kuò)散濾波方法，但其計(jì)算速度較慢[19-20]。

圖3 含四類AVA曲線的正演模型的時(shí)頻圖Fig.3 Forward modeling of the time-frequency diagrams with four kinds of AVA curves(a)原始模型；(b)加20%隨機(jī)噪聲模型帶通濾波結(jié)果；(c)加20%隨機(jī)噪聲模型擴(kuò)散濾波結(jié)果

4 實(shí)際資料處理

為了驗(yàn)證基于Jacobi算法求解結(jié)構(gòu)張量的置信擴(kuò)散濾波方法，在實(shí)際地震資料處理中的效果，將該方法應(yīng)用在Geoscope軟件平臺(tái)上進(jìn)行疊后和疊前地震數(shù)據(jù)去噪處理[21]。本文方法源自于圖像去噪處理，對(duì)圖像質(zhì)量沒(méi)有特定要求。由方法原理可知，該方法只需輸入地震數(shù)據(jù)，設(shè)置控制迭代次數(shù)的擴(kuò)散變換率閾值δ0、與擴(kuò)散系數(shù)有關(guān)的系數(shù)α、高斯窗尺度σ的值，其對(duì)其他條件并無(wú)特別要求，故對(duì)低、高信噪比的二維和三維地震資料均可適用。根據(jù)地震資料情況，僅去噪效果有一定的差別；對(duì)于信噪比極低的地震資料去噪效果本身就較差，而對(duì)于信噪比極高的地震資料去噪效果可能會(huì)不太顯著。三項(xiàng)關(guān)鍵參數(shù)δ0、α、σ值的選取參考方法設(shè)計(jì)原理和試驗(yàn)測(cè)試的結(jié)果，δ0和α取值都不宜過(guò)大且一般都取接近于零的極小值，略微調(diào)節(jié)該兩參數(shù)的大小對(duì)實(shí)際地震資料處理的影響不大，σ高斯窗尺度選取3或5。

圖4(a)～4(b)展示了鄰近幾道的疊前道集去噪對(duì)比圖。Geoscope軟件中采用作業(yè)流形式對(duì)逐個(gè)道集進(jìn)行處理。從圖4(a)可以看出，這幾個(gè)道集含有較強(qiáng)的隨機(jī)噪聲，信噪比不高。圖4(c)為去噪前后幾個(gè)道集的差值剖面，與隨機(jī)噪音比較符合。說(shuō)明經(jīng)過(guò)本文方法去噪后，該道集中的隨機(jī)噪聲得到了有效地壓制，同相軸也變得更加清晰，但是道集上仍有部分有規(guī)則形態(tài)的斜干擾噪聲不能很好去除。圖5(a)～圖5(c)展示了疊后地震資料的常規(guī)方法和本文方法去噪對(duì)比剖面圖。從中可知，本文方法壓制噪聲效果明顯優(yōu)于常規(guī)方法的效果：更加有效地壓制了隨機(jī)噪聲，更加增強(qiáng)了同相軸連續(xù)性，更明顯地突出了斷層位置，削弱了同相軸的蚯蚓狀現(xiàn)象。

實(shí)際資料處理表明，本文方法不僅在疊前去噪處理中為后續(xù)的AVO屬性技術(shù)和疊加處理提供了更有力地保證，而且在疊后地震數(shù)據(jù)體處理將更有利于層位自動(dòng)追蹤和斷層位置精確解釋。

此外，本文方法和常規(guī)方法的中δ0、σ值可選同一值，僅α值為本文方法中特有取值。當(dāng)α和σ值恒定時(shí)，本文方法同其常規(guī)方法一致，其計(jì)算速度和效率僅與控制迭代次數(shù)的擴(kuò)散變換率閾值δ0和待處理的地震數(shù)據(jù)體大小有關(guān)，δ0越小則迭代次數(shù)越多、地震數(shù)據(jù)體越大則計(jì)算量增加，導(dǎo)致計(jì)算效率有所降低。當(dāng)δ0一定時(shí)，本文方法僅與地震數(shù)據(jù)體大小有關(guān)。在Intel至強(qiáng)E3系列CPU處理器、64G內(nèi)存和64位操作系統(tǒng)的個(gè)人工作電腦上使用Geo-scope軟件平臺(tái)windows版本，計(jì)算1.7GB的疊后三維地震數(shù)據(jù)體的去噪結(jié)果，本文方法計(jì)算用時(shí)約為38 min，常規(guī)方法計(jì)算用時(shí)約為5 min，一定程度上認(rèn)為本文方法的計(jì)算用時(shí)在可接收范圍，而對(duì)較為龐大的疊前地震數(shù)據(jù)處理情況，可嘗試改進(jìn)算法及采用GPU計(jì)算來(lái)達(dá)到提高計(jì)算效率的目的。

圖4 擴(kuò)散濾波在疊前地震資料中的應(yīng)用效果Fig.4 Application results of diffusion filtering in pre-stack seismic data(a)原始疊前道集；(b)本文方法去噪后疊前道集；(c)去噪前后的差值剖面

圖5 擴(kuò)散濾波在疊后地震資料中的應(yīng)用效果Fig.5 Application results of diffusion filtering in post-stack seismic data(a)原始疊后剖面圖；(b)常規(guī)基于結(jié)構(gòu)張量的擴(kuò)散濾波方法去噪后剖面圖；(c)本文方法去噪后剖面圖

5 結(jié)論

筆者在非線性各向異性擴(kuò)散濾波基礎(chǔ)上，利用Jacobi算法求解結(jié)構(gòu)張量的特征值，并建立線狀置信度來(lái)計(jì)算擴(kuò)散濾波系數(shù)，提出了一種基于Jacobi算法求解結(jié)構(gòu)張量的置信擴(kuò)散濾波方法：

1)該方法對(duì)疊前道集的AVA曲線特征有一定保幅性，保證了后續(xù)AVO屬性分析和疊加處理過(guò)程的順利進(jìn)行。

2)該方法不僅壓制隨機(jī)噪聲效果好于常規(guī)基于結(jié)構(gòu)張量的擴(kuò)散濾波方法，而且在同相軸的連續(xù)性增強(qiáng)和突出斷層位置方面都更好的效果。

3)該方法可為層位自動(dòng)追蹤和斷層位置精確解釋，提供高信噪比資料。

4)該方法隨迭代次數(shù)的增加，SNR合理的變大，MSE有效的減小，最終達(dá)到一定迭代次數(shù)，擴(kuò)散變化率極低，擴(kuò)散過(guò)程停止，該方法將不再有濾波效果。

5)該方法相比帯通濾波方法，部分有效信號(hào)頻帶范圍內(nèi)的噪聲能被有效壓制，強(qiáng)地震同向軸位置的時(shí)頻能量團(tuán)更為聚焦，而弱地震同向軸的時(shí)頻能量團(tuán)掩蓋在噪聲中。

6)該方法的參數(shù)選取較為簡(jiǎn)單，自適應(yīng)計(jì)算擴(kuò)散系數(shù)，在軟件應(yīng)用中不用過(guò)度注重去調(diào)節(jié)和選取參數(shù)。