基于Jacobi算法求解結構張量的置信擴散濾波方法

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

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

0 引言

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

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

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

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

在非線性各向異性擴散濾波基礎上，采用Jacobi算法來求解結構張量D中的特征向量，并引入線狀置信度Cline來共同來構建擴散濾波系數λ，最后來求解擴散方程，實現置信擴散慮波。本文方法相比常規基于結構張量的擴散濾波方法，在構建和求解結構張量D方面及構建擴散濾波系數有所差異，最終效果較常規基于結構張量的擴散濾波方法也所提高：隨機噪聲壓制效果較好、地震同相軸連續性更強，斷層成像更加清晰準確，明顯改善資料品質。

1 常規方法

1.1 非線性各向異性擴散濾波原理

非線性各向異性擴散可模擬一個熱傳導的物理過程，它在沒有產生和損失熱量的前提下平衡溫度。將擴散應用于地震道集，用地震振幅代替熱量，在擴散過程中時間是一個固定參數。擴散進行的時間越長，地震道集濾波程度則越強，所以濾波過程相當于是提高信噪比的過程。地震振幅的梯度變化相當于熱傳導的通量，其擴散方程可以表示為[13-14]：

(1)

其中:D(|▽uσ|)為結構張量;|▽uσ|為正則化后的▽uσ，正則化過程相當于高斯濾波;t為擴散時間(也可計為迭代次數);u0(x,y)為原始地震道集；u(x,y,t)為迭代次數為t次后的濾波地震道集。

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

(2)

(3)

式中:σ為尺度函數，可以根據噪聲的均方誤差自適應調節。

1.2 最優旋轉不變算法

公式(1)可在正交坐標系下展開成如下形式：

(4)

其中結構張量可表示為式(5)。

(5)

公式(4)可表示為沿著x、y兩個方向(即地震道集道號和時間方向)上進行擴散，而最優旋轉不變算法是需要對公式(5)的結構張量D進行一定程度上的旋轉，使其沿著地震道集局部振幅值變化最快和最慢方向上擴散，如王益艷等[15]再用擴散系數來控制在不同變化方向上平滑程度，即常規基于結構張量的擴散濾波方法。

公式(4)旋轉及引入擴散系數之后可轉換成如下形式：

(6)

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

2 本文方法

2.1 結構張量的改進

為了增強同向軸連續性和提高壓制噪聲的效果，在本文方法中對結構張量的設計進行了改進，并引入線狀置信度來改進和自適應調節擴散系數。

該方法以公式(4)為基礎，并對公式(5)中的結構張量改進為如下形式：

(7)

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

λ1=α

(8)

(9)

(10)

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

(11)

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

2.2 Jacobi算法原理

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

(12)

2.3 方法的實現步驟

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

(13)

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

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

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

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

3 模型試驗及效果分析

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

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

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

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

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

3.2 去噪試驗的信噪比及誤差分析

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

3.3 去噪試驗的時頻分析

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

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

4 實際資料處理

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

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

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

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

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

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

5 結論

筆者在非線性各向異性擴散濾波基礎上，利用Jacobi算法求解結構張量的特征值，并建立線狀置信度來計算擴散濾波系數，提出了一種基于Jacobi算法求解結構張量的置信擴散濾波方法：

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

2)該方法不僅壓制隨機噪聲效果好于常規基于結構張量的擴散濾波方法，而且在同相軸的連續性增強和突出斷層位置方面都更好的效果。

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

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

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

6)該方法的參數選取較為簡單，自適應計算擴散系數，在軟件應用中不用過度注重去調節和選取參數。