自適應結構方向復擴散壓制沙漠地震隨機噪聲

林紅波, 張丹丹

(吉林大學 通信工程學院, 長春 130012)

0 引 言

在地震勘探中, 所獲得的地震勘探數(shù)據(jù)往往伴隨著大量的隨機噪聲, 嚴重影響地震數(shù)據(jù)處理的精度, 給進一步的地質(zhì)解釋工作帶來很多困難。因此, 壓制地震勘探隨機噪聲, 提高地震勘探數(shù)據(jù)質(zhì)量是地震勘探數(shù)據(jù)處理中重要的基礎環(huán)節(jié)。近些年來, 學者們提出很多地震勘探隨機噪聲壓制方法, 如維納濾波[1]、經(jīng)驗模式分解[2]、奇異值分解[3]、小波變換[4]和曲波變換[5]等。這些方法在地震勘探隨機噪聲壓制方面取得了很好的效果。但在低信噪比的情況下, 有效信號的恢復不夠理想, 尤其在沙漠地震勘探數(shù)據(jù)中, 沙漠隨機噪聲頻率低, 能量強, 與有效信號具有相似的時空特性和頻率特性, 使辨識和提取微弱信號變得更加困難。

基于PDE(Partial Differential Equation)的多尺度分析方法在圖像去噪方面效果顯著, 最早的PDE模型是實數(shù)域線性擴散方程[6], 擴散方程的解等同于含噪圖像和時變標準偏差的高斯核的卷積, 擴散系數(shù)為常數(shù)。該方法能有效平滑隨機噪聲, 同時也模糊了一些邊緣和細節(jié)。在此基礎上, 非線性擴散模型構建了基于梯度變化的擴散函數(shù), 有效改善了圖像結構信息的保留效果[7-10]。非線性擴散濾波方法不依賴先驗知識, 如頻率、線性、低秩以及噪聲統(tǒng)計特性等, 因此在復雜隨機噪聲壓制方面具有明顯的優(yōu)勢。然而, 這種實數(shù)域擴散方法利用梯度辨識邊緣并調(diào)整擴散強度, 在低信噪比情況下易混淆邊緣和噪聲?；谶@個問題, Gilboa等[11]提出復數(shù)域擴散方法。非線性復擴散(NCD: Nonlinear Complex Diffusion)的虛部為圖像平滑的二階導, 能抑制噪聲對邊緣檢測的影響, 擴散系數(shù)是基于虛部的單調(diào)遞減函數(shù), 相比于基于梯度變化的擴散系數(shù)能更好適應復雜的結構特征?；贜CD的算法適應低信噪比數(shù)據(jù)和復雜結構, 成功應用于地震勘探數(shù)據(jù)去噪[12]。然而, 該擴散方法仍然利用梯度描述同相軸的局部特征, 無法準確刻畫沙漠地震勘探數(shù)據(jù)中具有復雜結構同相軸的方向特征, 易導致同相軸衰減。筆者提出了自適應結構方向復擴散算法(ASOCD：Adaptive Structure-Oriented Complex Diffusion), 通過在復擴散方程中引入基于結構方向的方向導數(shù), 使擴散過程沿同相軸方向演化, 更好地保留同相軸；此外, 構建了基于局部協(xié)方差矩陣的自適應擴散閾值, 有效地解決了沙漠隨機噪聲壓制不徹底問題。合成數(shù)據(jù)和實際數(shù)據(jù)處理驗證了ASOCD方法的有效性。

1 線性復擴散與非線性復擴散

1.1 線性復擴散

Gilboa等[11]受到量子力學啟發(fā), 將實數(shù)域的線性擴散和自由薛定諤方程結合得到線性復擴散。自由薛定諤方程定義為

(1)

線性擴散方程定義為

(2)

其中c是常擴散系數(shù), 初始條件為u(x)|t=0=u0(x),u0(x)為初始的含噪數(shù)據(jù)?；谧杂裳Χㄖ@方程, 將擴散方程推廣到復數(shù)域, 得到線性復擴散方程

(3)

其中c=rejθ,r是c的模。當θ趨近于0時, 式(3)的基本解可以表示為

(4)

1.2 非線性復擴散

線性復擴散的常擴散系數(shù)會導致在去除噪聲的同時平滑掉邊緣。為解決這個問題, Gilboa等[11]引入以虛部為變量的擴散系數(shù), 提出了非線性復擴散NCD。對于含噪數(shù)據(jù)u∈R2, 非線性復擴散定義為

(5)

其中擴散系數(shù)c是關于虛部單調(diào)遞減的函數(shù), 擴散閾值k和虛部共同作用調(diào)節(jié)擴散強度。為了應用復擴散算法, 需要在離散形式下求解, 用u(p,q)表示地震勘探數(shù)據(jù)u的第q道的第p個采樣點。對式(5)的離散化, 筆者采用有限差分離散化方法, 得NCD方程離散化公式為

(6)

(7)

從式(7)中可看出, NCD方法利用相鄰點表征式(5)的中的u, 從而實現(xiàn)擴散沿著含噪數(shù)據(jù)的邊緣方向和梯度方向, 并結合擴散強度實現(xiàn)噪聲的抑制以及邊緣的保留。

2 自適應結構方向復擴散原理

由于接收沙漠地震數(shù)據(jù)的檢波器間隔較遠, 相鄰檢波器接收到反射子波的時間延遲比較長, 因此同相軸往往具有較大的傾角。而NCD方法僅利用相鄰點表征結構信息, 所描述的邊緣方向和沙漠地震數(shù)據(jù)的同相軸方向不能吻合, 會導致同相軸的失真。另外, 部分沙漠地震隨機噪聲和有效信號具有相似的虛部, 對信號和噪聲采用相同的擴散閾值, 無法在保留有效信號的同時抑制隨機噪聲。針對這些問題, 筆者提出自適應結構方向復擴散(ASOCD)方法, 構造基于結構方向的方向導數(shù), 使擴散能沿同相軸方向演化, 并根據(jù)區(qū)域特性自適應調(diào)節(jié)擴散閾值, 從而在抑制噪聲的同時增強有效信號。

在離散NCD模型的基礎上, 引入基于結構方向的方向導數(shù), 得到自適應結構方向復擴散模型

(8)

圖1 掃描示意圖

其中cT表示結構方向的擴散系數(shù),cN表示為垂直于結構方向的擴散系數(shù),DR(un)表示結構方向的方向導數(shù)估計,DL(un)表示與結構方向相反方向的方向導數(shù)估計。如圖1所示, 方向導數(shù)估計為沿同相軸方向求差分, 定義為

(9)

其中sp為結構方向的斜率。式(8)表明, 自適應結構方向復擴散能沿著同相軸方向定向擴散, 擴散強度由擴散系數(shù)調(diào)節(jié)。為了實現(xiàn)對隨機噪聲的抑制, 需要增強噪聲區(qū)域的擴散系數(shù)cT和cN；而為了更好地恢復同相軸, 則信號區(qū)域的擴散系數(shù)cT要比噪聲區(qū)域小, 并盡可能減少擴散系數(shù)cN。為此, 筆者構建了自適應復擴散系數(shù)。

2.1 自適應復擴散系數(shù)

基于地震勘探數(shù)據(jù)的區(qū)域特征, 自適應地調(diào)整ASOCD模型的擴散系數(shù)。結構方向的擴散系數(shù)cT定義為

cT(Im(u(p,q)))=ejθ/(1+(Im(u(p,q))/k′(p,q)θ)2)

(10)

其中閾值k′(p,q)能在信號區(qū)域Ω和噪聲區(qū)域φ自適應調(diào)整, 定義為

(11)

滿足Kn>Ks, 垂直與結構方向的擴散系數(shù)cN定義為

(12)

由式(10)和式(12)可知, 結構方向復擴散模型在噪聲區(qū)域和信號區(qū)域具有不同的擴散強度。沙漠隨機噪聲具有弱相似性, 和有效信號具有相似虛部, 此時利用較大的閾值Kn可以增強擴散強度, 能解決常閾值擴散濾波對噪聲壓制不徹底的問題。在信號區(qū)域Ω, 結構方向復擴散只進行結構方向的定向擴散, 在同相軸方向采用較小閾值控制擴散強度, 在垂直同相軸方向的擴散系數(shù)為0, 避免同相軸模糊。

筆者利用局部協(xié)方差矩陣辨識信號區(qū)域Ω和噪聲區(qū)域φ。由于有效信號存在較強的空間相關性, 沙漠地震勘探隨機噪聲沒有或存在較弱的空間相關性, 而協(xié)方差矩陣能反映矩陣中任意兩個變量之間的相關性[13], 則能反映在一定窗長內(nèi)的任意兩道地震數(shù)據(jù)之間的相關性, 因此能識別信號區(qū)域和噪聲區(qū)域。對于點u(p,q)的局部協(xié)方差矩陣Fp,q的計算, 以該點為中心, 計算大小為A×B的鄰域L的協(xié)方差矩陣, 則局部協(xié)方差矩陣Fp,q中任意一點Fp,q(p′,q′)可表示為

(13)

其中up′表示矩陣L中的第p′列,μp′表示矩陣L中p′列的均值, T代表轉置。局部協(xié)方差矩陣最大的值在對角線上取得, 對于噪聲區(qū)域, 局部協(xié)方差矩陣對角線上的值即是噪聲方差, 對于信號區(qū)域, 局部協(xié)方差矩陣對角線上的值即是含噪信號的方差?？梢栽O置區(qū)分信號和噪聲的閾值τ, 表示為

τ=σuσn

(14)

其中σn是利用Chen等[14]提出的噪聲估計方法估計出的噪聲標準差,σu表示含噪數(shù)據(jù)的標準差。由于含噪信號方差大于噪聲方差, 則對于信號區(qū)域, 局部協(xié)方差矩陣對角線上的值要比閾值大, 對于噪聲區(qū)域, 局部協(xié)方差矩陣對角線上的值要比閾值小。則識別的信號區(qū)域Ω和噪聲區(qū)域φ表示為

(15)

2.2 結構方向估計

結構方向復擴散模型基于方向導數(shù)沿地震勘探數(shù)據(jù)的結構方向自適應擴散, 其中結構方向的估計至關重要。基于地震勘探同相軸的空間相關性, 筆者采用非局部統(tǒng)計特性掃描方法估計地震勘探數(shù)據(jù)的結構方向。由于一部分沙漠隨機噪聲具有弱空間相似性, 若對隨機噪聲也進行方向擴散, 易產(chǎn)生具有一定空間相似性的噪聲殘留。為更好地壓制隨機噪聲, 筆者在估計結構方向時增加約束條件, 使噪聲區(qū)域的結構方向為水平方向, 此時結構方向擴散退化為非線性復擴散, 并通過自適應調(diào)節(jié)擴散系數(shù)增強對沙漠隨機噪聲抑制能力。

(16)

(17)

3 仿真實驗

3.1 模擬記錄的處理結果

為了驗證ASOCD方法的有效性, 對合成的沙漠地震記錄進行處理。合成的沙漠地震記錄包含由Ricker構成的3條同相軸, 主頻率從上到下依次為20 Hz, 18 Hz和16 Hz。采樣頻率是1 000 Hz, 如圖2a所示。添加真實的沙漠隨機噪聲, 使信噪比達到-5.46 dB(見圖2b)。將ASOCD濾波方法應用于含噪數(shù)據(jù)處理并和NCD方法進行對比, 去噪結果以及濾除的噪聲如圖2c～圖2f所示。可以看出, 這兩種方法都能抑制噪聲, 但ASOCD濾波方法能獲得更干凈的背景。在信號保留方面, NCD方法濾除的噪聲中可以看出同相軸, 則NCD方法在抑制噪聲的同時損失了有效信號, 而ASOCD算法去除的噪聲中幾乎不存在有效信號, 因此ASOCD算法能夠抑制噪聲的同時保留有效信號。

a 純凈信號 b 含噪信號 c NCD去噪結果

d NCD濾除的噪聲 e ASOCD去噪結果 f ASOCD濾除的噪聲

對去噪結果進行量化分析, 表1給出了不同沙漠噪聲水平下, 去噪前后數(shù)據(jù)的信噪比比較, 其中使用信噪比公式為

(18)

其中us和ud分別為原始無噪數(shù)據(jù)和去噪后的數(shù)據(jù), 數(shù)據(jù)結果如表1所示。從表1中可看出, 不同噪聲水平下, ASOCD算法去噪后的信噪比相比于傳統(tǒng)NCD算法平均提高了6 dB。

表1 不同沙漠噪聲水平下的信噪比對比

3.2 實際記錄的處理結果

為了進一步驗證提出方法的有效性, 對真實的沙漠地震數(shù)據(jù)進行處理。圖3a為中國沙漠地區(qū)包含107道和1 400個采樣點的真實沙漠地震數(shù)據(jù)。顯然有效信號被噪聲嚴重損壞了, 某些噪聲區(qū)域的幅值甚至超過了有效信號。對該地震數(shù)據(jù)分別應用NCD方法和ASOCD方法去噪, 對比去噪后的結果, ASOCD方法(見圖3c)相比NCD方法(見圖3b)可以獲得更加干凈的背景。尤其在黑色框所示區(qū)域, ASOCD方法能獲得更加清晰更加連續(xù)的同相軸, 噪聲也能得到更好的抑制。

a 原始地震數(shù)據(jù) b ASOCD去噪結果 c NCD去噪結果

4 結 語

筆者提出ASOCD算法對沙漠地震數(shù)據(jù)中的隨機噪聲進行壓制, 利用非局部特性掃描方法計算同相軸的方向, 并構造基于結構方向的方向導數(shù), 使擴散能沿著同相軸方向進行。此外, 利用局部協(xié)方差矩陣提高了信號區(qū)域和噪聲區(qū)域的識別能力, 能更有效區(qū)分信號和噪聲在低信噪比情況下, 使得擴散閾值能自適應區(qū)域特性進行調(diào)整。ASOCD算法實現(xiàn)了在噪聲區(qū)域進行擴散閾值較大的非線性復擴散, 有效地抑制沙漠隨機噪聲, 在信號區(qū)域沿著同相軸方向進行擴散閾值較小的定向復擴散, 有效地保留同相軸。合成的沙漠地震數(shù)據(jù)以及真實的沙漠地震數(shù)據(jù)驗證了所提算法的有效性。同NCD算法相比, ASOCD算法具有更好的信號恢復能力, 以及更優(yōu)越的噪聲抑制能力。