林紅波, 張丹丹

(吉林大學(xué) 通信工程學(xué)院, 長(zhǎng)春 130012)

0 引 言

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

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

1 線性復(fù)擴(kuò)散與非線性復(fù)擴(kuò)散

1.1 線性復(fù)擴(kuò)散

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

(1)

線性擴(kuò)散方程定義為

(2)

其中c是常擴(kuò)散系數(shù), 初始條件為u(x)|t=0=u0(x),u0(x)為初始的含噪數(shù)據(jù)。基于自由薛定諤方程, 將擴(kuò)散方程推廣到復(fù)數(shù)域, 得到線性復(fù)擴(kuò)散方程

(3)

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

(4)

1.2 非線性復(fù)擴(kuò)散

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

(5)

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

(6)

(7)

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

2 自適應(yīng)結(jié)構(gòu)方向復(fù)擴(kuò)散原理

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

在離散NCD模型的基礎(chǔ)上, 引入基于結(jié)構(gòu)方向的方向?qū)?shù), 得到自適應(yīng)結(jié)構(gòu)方向復(fù)擴(kuò)散模型

(8)

圖1 掃描示意圖

其中cT表示結(jié)構(gòu)方向的擴(kuò)散系數(shù),cN表示為垂直于結(jié)構(gòu)方向的擴(kuò)散系數(shù),DR(un)表示結(jié)構(gòu)方向的方向?qū)?shù)估計(jì),DL(un)表示與結(jié)構(gòu)方向相反方向的方向?qū)?shù)估計(jì)。如圖1所示, 方向?qū)?shù)估計(jì)為沿同相軸方向求差分, 定義為

(9)

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

2.1 自適應(yīng)復(fù)擴(kuò)散系數(shù)

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

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

(10)

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

(11)

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

(12)

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

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

(13)

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

τ=σuσn

(14)

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

(15)

2.2 結(jié)構(gòu)方向估計(jì)

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

(16)

(17)

3 仿真實(shí)驗(yàn)

3.1 模擬記錄的處理結(jié)果

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

a 純凈信號(hào) b 含噪信號(hào) c NCD去噪結(jié)果

d NCD濾除的噪聲 e ASOCD去噪結(jié)果 f ASOCD濾除的噪聲

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

(18)

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

表1 不同沙漠噪聲水平下的信噪比對(duì)比

3.2 實(shí)際記錄的處理結(jié)果

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

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

4 結(jié) 語(yǔ)

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