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

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

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

0 引言

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

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

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

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

在非線性各向異性擴(kuò)散濾波基礎(chǔ)上，采用Jacobi算法來求解結(jié)構(gòu)張量D中的特征向量，并引入線狀置信度Cline來共同來構(gòu)建擴(kuò)散濾波系數(shù)λ，最后來求解擴(kuò)散方程，實現(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ò)散可模擬一個熱傳導(dǎo)的物理過程，它在沒有產(chǎn)生和損失熱量的前提下平衡溫度。將擴(kuò)散應(yīng)用于地震道集，用地震振幅代替熱量，在擴(kuò)散過程中時間是一個固定參數(shù)。擴(kuò)散進(jìn)行的時間越長，地震道集濾波程度則越強(qiáng)，所以濾波過程相當(dāng)于是提高信噪比的過程。地震振幅的梯度變化相當(dāng)于熱傳導(dǎo)的通量，其擴(kuò)散方程可以表示為[13-14]：

(1)

其中:D(|▽uσ|)為結(jié)構(gòu)張量;|▽uσ|為正則化后的▽uσ，正則化過程相當(dāng)于高斯濾波;t為擴(kuò)散時間(也可計為迭代次數(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)系下展開成如下形式：

(4)

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

(5)

公式(4)可表示為沿著x、y兩個方向(即地震道集道號和時間方向)上進(jìn)行擴(kuò)散，而最優(yōu)旋轉(zhuǎn)不變算法是需要對公式(5)的結(jié)構(gòu)張量D進(jìn)行一定程度上的旋轉(zhuǎn)，使其沿著地震道集局部振幅值變化最快和最慢方向上擴(kuò)散，如王益艷等[15]再用擴(kuò)散系數(shù)來控制在不同變化方向上平滑程度，即常規(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ù)性和提高壓制噪聲的效果，在本文方法中對結(jié)構(gòu)張量的設(shè)計進(jìn)行了改進(jìn)，并引入線狀置信度來改進(jìn)和自適應(yīng)調(diào)節(jié)擴(kuò)散系數(shù)。

該方法以公式(4)為基礎(chǔ)，并對公式(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)η值和ξ值比較接近時，從公式(11)可以看出兩個方向上擴(kuò)散速率較為一致。公式(11)中振幅劇烈變化處位置值η值必然大于振幅平穩(wěn)變化處位置值ξ，而當(dāng)ξ值遠(yuǎn)小于η值時即代表振幅平穩(wěn)位置有很快的擴(kuò)散速率。

2.2 Jacobi算法原理

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

(12)

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

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

(13)

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

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

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

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

3 模型試驗及效果分析

3.1 理論模型正演及去噪試驗

表1給出了單層滿足四類AVA曲線的各向同性介質(zhì)理論模型，圖1(a)～圖1(c)為由表1中給出的模型組合在一個道集上的AVA正演角道集，入射角為1°～45°,AVA正演角道集采用合成記錄正演。從圖1(c)擴(kuò)散濾波方法的結(jié)果中仍然可見四類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 去噪試驗的信噪比及誤差分析

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

3.3 去噪試驗的時頻分析

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

圖3 含四類AVA曲線的正演模型的時頻圖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 實際資料處理

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

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

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

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

圖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)張量的特征值，并建立線狀置信度來計算擴(kuò)散濾波系數(shù)，提出了一種基于Jacobi算法求解結(jié)構(gòu)張量的置信擴(kuò)散濾波方法：

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

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

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

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

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

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