基于交疊組稀疏非凸Lp偽范數(shù)高階全變分的地震隨機(jī)噪聲壓制

梁上林 胡天躍* 崔 棟 隋京坤

(①北京大學(xué)地球與空間科學(xué)學(xué)院，北京 100871； ②中國(guó)石油勘探開(kāi)發(fā)研究院，北京 100083)

0 引言

地震記錄中常常混雜著大量的隨機(jī)噪聲，嚴(yán)重降低了資料的信噪比。有效壓制地震數(shù)據(jù)中的隨機(jī)噪聲是地震資料處理的重要環(huán)節(jié)。然而，由于其產(chǎn)生具有隨機(jī)性且無(wú)固定頻率和視速度，隨機(jī)噪聲壓制一直是勘探地球物理界的難點(diǎn)之一。目前，壓制隨機(jī)噪聲的方法可歸納為時(shí)間域?yàn)V波[1-2]、頻率域?yàn)V波[3]、空間域?yàn)V波[4]和各種變換域?yàn)V波[5-9]等。Rudin等[4]于1992年首次提出空間域的全變分(Total variation，TV)模型，利用含噪數(shù)據(jù)總變分比無(wú)噪信號(hào)總變分大的性質(zhì)，構(gòu)造能量泛函以達(dá)到去噪的目的。鑒于其有效性，該模型被廣泛應(yīng)用于噪聲壓制[10-12]、波阻抗反演[13]和全波形反演[14]等領(lǐng)域。

盡管TV模型具有保護(hù)圖像邊緣信息的優(yōu)勢(shì)，但存在嚴(yán)重的階梯效應(yīng)，導(dǎo)致信號(hào)分片光滑。針對(duì)此問(wèn)題，Selesnick等[15]率先將交疊組稀疏(Overlapping group sparsity，OGS)引入TV模型，指出信號(hào)的一階差分不僅具有稀疏特性，還具有結(jié)構(gòu)稀疏特性。之后Selesnick等[16]通過(guò)正則化技術(shù)有效識(shí)別出平滑區(qū)域內(nèi)被強(qiáng)噪聲污染的信號(hào)。為挖掘圖像一階梯度和二階梯度的結(jié)構(gòu)稀疏性，陳穎頻等[17]將OGS與TV模型相結(jié)合，構(gòu)建了一種改進(jìn)的去噪模型。

上述方法[15-17]從本質(zhì)上都是在尋找不同的正則項(xiàng)及懲罰函數(shù)，以達(dá)到最佳去噪效果。考慮到L1正則項(xiàng)是一種凸近似，往往使重構(gòu)的非零信號(hào)比真實(shí)值小而導(dǎo)致過(guò)擬合，Ahlad等[18]提出組內(nèi)使用L2范數(shù)和組間使用L1范數(shù)的加權(quán)方式。Lp偽范數(shù)本質(zhì)是一種非凸正則項(xiàng)，它具有比凸核更接近秩函數(shù)的優(yōu)點(diǎn)，能在恢復(fù)信號(hào)細(xì)節(jié)與給定殘差稀疏解之間達(dá)到平衡，被廣泛應(yīng)用于圖形修復(fù)、地震數(shù)據(jù)重建和重力反演等領(lǐng)域[19-21]。Adam等[22-23]基于高階TV模型，引入非凸正則項(xiàng)技術(shù)，在圖形去模糊方面取得一定效果； 同時(shí)指出OGS與非凸高階TV模型在去除斑塊虛假信息上是互補(bǔ)的。然而，目前尚未見(jiàn)到將Lp偽范數(shù)、OGS和TV模型三者相結(jié)合應(yīng)用于地震資料去噪的實(shí)例。

鑒于一階TV模型不能有效保護(hù)信號(hào)的邊緣及紋理特征，存在嚴(yán)重的階梯效應(yīng)，本文將高階TV模型與OGS相結(jié)合，目的是充分挖掘和利用兩者的優(yōu)點(diǎn)，在減弱階梯效應(yīng)的同時(shí)較好地保護(hù)局部信息； 同時(shí)，為解決信號(hào)重構(gòu)過(guò)程中產(chǎn)生的斑點(diǎn)和塊狀問(wèn)題，保護(hù)非零有效信號(hào)，引入非凸Lp偽范數(shù)； 最后通過(guò)模擬數(shù)據(jù)和實(shí)際資料對(duì)本文方法進(jìn)行驗(yàn)證，并與常規(guī)五種去噪方法進(jìn)行對(duì)比。

1 去噪原理

1.1 傳統(tǒng)TV去噪模型

實(shí)際地震數(shù)據(jù)可表示為

y=s+n

(1)

式中：y表示檢波器接收到的地震數(shù)據(jù)；s表示有效信號(hào)；n為隨機(jī)噪聲。TV去噪模型對(duì)應(yīng)如下數(shù)學(xué)極值問(wèn)題[4]

(2)

1.2 二維OGS正則項(xiàng)

Selesnick等[15]定義了窗尺度為K×K且中心點(diǎn)為(i,j)的二維矩陣

(3)

m1=K-12、m2=K2,

式中表示向下取整，K表示組稀疏交疊程度。二維OGS正則項(xiàng)定義為

(4)

可見(jiàn)，OGS將鄰域信息作為參考形成組合梯度，能充分挖掘信號(hào)的局部相似性。

1.3 基于OGS非凸Lp偽范數(shù)高階TV模型

考慮到OGS能有效恢復(fù)全局較大輪廓而高階TV模型可較好地保護(hù)局部弱信號(hào)，將OGS與高階TV模型結(jié)合，充分利用兩者優(yōu)點(diǎn)形成一個(gè)混合去噪模型； 同時(shí)，引入非凸Lp偽范數(shù)以更準(zhǔn)確地重構(gòu)噪聲數(shù)據(jù)中非零信號(hào)。該模型所對(duì)應(yīng)的數(shù)學(xué)極值問(wèn)題[22]表示為

(5)

改進(jìn)的方法是以信號(hào)周?chē)c(diǎn)的信息作為參考形成組合梯度，使孤立的噪聲污染點(diǎn)與鄰域的組合梯度下降； 同時(shí)保持圖像邊緣點(diǎn)與鄰域的組合梯度強(qiáng)度，通過(guò)設(shè)定合理的閾值，可在有效去除噪聲的同時(shí)，減弱階梯效應(yīng)，較好地保護(hù)有效信號(hào)。

采用分離變量法將式(5)轉(zhuǎn)化成如下受約束的極值問(wèn)題

s.t.a=s,b=2s,c=s

(6)

該式對(duì)應(yīng)的無(wú)約束增廣拉格朗日極值為

L(s,a,b,c;ξ1,ξ2,ξ3)=

(7)

為求解式(7)所示的復(fù)雜耦合問(wèn)題，本文采用交替方向乘子法(Alternating direction method of multipliers，ADMM)[25]，按照下式交替迭代更新

(8)

可見(jiàn)，式(7)所示耦合問(wèn)題被分解成四個(gè)子問(wèn)題。下面利用式(9)～式(13)求解相應(yīng)子問(wèn)題。

sk+1是一個(gè)最小平方問(wèn)題，其解析解為

sk+1=[λ1I+βT+β(2)T2+βI]-1[λ1y-Tξ1+

βTak-(2)Tξ2+β(2)Tbk-ξ3+βck]

(9)

ak+1的約束如下

(10)

式中包含了式(4)所示的OGS正則項(xiàng)，當(dāng)K=1時(shí)式(10)退化成傳統(tǒng)的高階TV模型； 當(dāng)K>1時(shí)，式(10)表示求解OGS正則項(xiàng)。鑒于最大最小值(Majorization minimization，MM)算法具有高效利用函數(shù)凹凸性搜尋原函數(shù)最值的優(yōu)點(diǎn)，當(dāng)原目標(biāo)函數(shù)難優(yōu)化時(shí)，不直接對(duì)其求解，而是求解一個(gè)易于優(yōu)化且逼近于原目標(biāo)函數(shù)的替代函數(shù)[26]。本文通過(guò)二次函數(shù)替代φ(a)實(shí)現(xiàn)MM算法。

同理，bk+1的約束為

(11)

式(11)是一個(gè)非凸優(yōu)化問(wèn)題。非凸優(yōu)化算法和閾值的選取是決定去噪效果的兩個(gè)重要因素。鑒于經(jīng)典的加權(quán)方向迭代L1算法能有效平衡信號(hào)的正負(fù)二階微分[27-28]，因此本文采用該算法求解

(12)

式中：ζ和ε都是較小量，用以避免分母為零； max表示求二者較大值； sign為符號(hào)函數(shù)。

類(lèi)似地，ck+1的約束為

(13)

另外，求解式(7)所需的三個(gè)拉格朗日乘子可表示如下

(14)

至此，四個(gè)子問(wèn)題已求解完成。

將整個(gè)算法流程總結(jié)如下：

(1)輸入二維地震數(shù)據(jù)y，給定參數(shù)λ1>0，λ2>0，K，p；

(2)初始化模型：s0=y，k=0，β=0.004，si=1,2,3=0，模型更新率εmin=1×10-6；

(3)更新式(9)中的sk+1；

(4)采用MM迭代算法更新式(10)中的ak+1；

(5)更新式(12)中的bk+1；

(6)更新式(13)中的ck+1；

2 模擬數(shù)據(jù)驗(yàn)證

2.1 評(píng)價(jià)指標(biāo)

信噪比通常用信號(hào)總能量與噪聲總能量的比值SNR(Signal to noise ratio)表征，是評(píng)價(jià)整體去噪效果的常用指標(biāo)；結(jié)構(gòu)相似性參數(shù)SSIM(Structural similarity)從亮度、對(duì)比度、結(jié)構(gòu)三方面度量?jī)煞鶊D像的相似性，常作為衡量圖像失真或噪聲壓制情況的客觀標(biāo)準(zhǔn)[29]； 同時(shí)，為表征高斯噪聲強(qiáng)弱，定義零均值噪聲的標(biāo)準(zhǔn)差為δ。三者的具體定義如下

(15)

(16)

(17)

式中：M、N分別表示二維地震數(shù)據(jù)的行數(shù)和列數(shù)；μs、μy對(duì)應(yīng)表示s、y的均值；σs、σy分別表示s、y的方差；σsy表示s與y的協(xié)方差；C1和C2是維持穩(wěn)定的兩個(gè)常量，一般分別取2.252和6.752。為了更全面地評(píng)判去噪效果，本文將SNR、SSIM和運(yùn)算耗時(shí)作為三個(gè)客觀評(píng)價(jià)指標(biāo)。

2.2 模擬數(shù)據(jù)

構(gòu)建一個(gè)水平層狀模型，通過(guò)有限差分聲波正演模擬得到圖1所示的共炮記錄。模擬采用主頻為20Hz的雷克子波，采樣間隔為4ms，采樣長(zhǎng)度為4s，接收道數(shù)為500。為了模擬地震數(shù)據(jù)中不同強(qiáng)度的隨機(jī)噪聲，向該記錄加入標(biāo)準(zhǔn)差為δ的高斯白噪聲，得到圖2所示的含噪數(shù)據(jù)。由圖2可知，隨著δ增大，背景噪聲加強(qiáng)，SNR降低，SSIM變差，弱反射信號(hào)逐漸被強(qiáng)噪聲淹沒(méi)。

圖1 正演模擬共炮記錄

圖2 不同δ的含噪數(shù)據(jù)(a)δ=10； (b)δ=20； (c)δ=30； (d)δ=40

2.3 參數(shù)優(yōu)選

OGS從一個(gè)點(diǎn)向四周延伸K×K個(gè)點(diǎn)，可有效挖掘信號(hào)的局部相似性。為考察K值對(duì)去噪效果的影響，通過(guò)不斷調(diào)整其大小，分別對(duì)圖2所示的四種不同噪聲背景下的地震數(shù)據(jù)做去噪處理。

圖3展示不同K值得到的去噪結(jié)果所對(duì)應(yīng)的三種指標(biāo)的變化曲線(xiàn)。從圖3a可見(jiàn)，隨著K值的增大，SNR先從小到大逐漸增加，達(dá)到最大值后減小。當(dāng)SNR取最大值時(shí)，K值分別為3、5、7和10，說(shuō)明K值的優(yōu)選與背景噪聲強(qiáng)度有關(guān)。當(dāng)背景噪聲較弱時(shí)，較小K值就能達(dá)到滿(mǎn)意的去噪效果；當(dāng)背景噪聲增強(qiáng)時(shí)，應(yīng)適當(dāng)增大K值以提高SNR。圖3b中的SSIM表現(xiàn)出與SNR類(lèi)似形態(tài)，區(qū)別在于當(dāng)SSIM達(dá)到最大值后緩慢減小。這說(shuō)明K值超過(guò)最優(yōu)值后，對(duì)SSIM的影響不顯著。圖3c顯示，隨著K值的增大運(yùn)算耗時(shí)也不斷增加，這是由于鄰域信息的引入，導(dǎo)致計(jì)算量增大。

圖3 去噪結(jié)果的三指標(biāo)隨K值的變化曲線(xiàn)(a)SNR； (b)SSIM； (c)耗時(shí)

因此，K值對(duì)本文方法去噪效果有顯著影響。當(dāng)K值過(guò)小時(shí)，鄰域信息不足導(dǎo)致SNR和SSIM達(dá)不到最優(yōu)，去噪效果欠佳； 而當(dāng)K值過(guò)大時(shí)，又會(huì)引入過(guò)多不相似的數(shù)據(jù)點(diǎn)，反而導(dǎo)致評(píng)價(jià)指標(biāo)減小，去噪能力下降，這不僅會(huì)產(chǎn)生平滑效應(yīng)，還會(huì)因使用過(guò)多鄰域信息增加不必要的計(jì)算量。

邊界與局部細(xì)節(jié)的恢復(fù)主要受非凸正則化階數(shù)p控制，對(duì)圖2含噪數(shù)據(jù)進(jìn)行處理并得到圖4所示三種指標(biāo)隨p的變化曲線(xiàn)。圖4a表明，當(dāng)δ=10時(shí)，SNR隨p值增大而快速增至最大值，然后快速減小；當(dāng)δ=20、30和40時(shí)，SNR曲線(xiàn)都具有平緩段、快速上升段和快速下降段。在平緩段SNR對(duì)p值不敏感，經(jīng)快速上升后達(dá)到最大值。SNR取最大值時(shí)對(duì)應(yīng)p值分別為0.23、0.46、0.62和0.75，這說(shuō)明p值的優(yōu)選與噪聲強(qiáng)度有關(guān)。圖4b表明，弱噪聲情況下SSIM曲線(xiàn)先上升，然后保持平緩；當(dāng)背景噪聲變強(qiáng)時(shí)，SSIM曲線(xiàn)呈現(xiàn)平緩段—快速上升段—平緩段分布。該指標(biāo)取最大值時(shí)對(duì)應(yīng)p值分別為0.16、0.47、0.62和0.83。圖4c表明，隨著p值增大，運(yùn)算耗時(shí)先增至最大值，然后逐漸減小，顯然該指標(biāo)與背景噪聲強(qiáng)度無(wú)關(guān)。

圖4 去噪結(jié)果的三種指標(biāo)隨p值的變化曲線(xiàn)(a)SNR； (b)SSIM； (c)耗時(shí)

綜上，p值對(duì)本文去噪方法有同樣顯著的影響。高SNR資料應(yīng)取較小p值；反之，低SNR資料應(yīng)適當(dāng)增大p值，以達(dá)到最佳去噪效果。

圖5展示K與p最優(yōu)情形下對(duì)圖2所示四種不同強(qiáng)度噪聲背景下地震數(shù)據(jù)的去噪效果。不同強(qiáng)度隨機(jī)噪聲均得到有效壓制，剖面整體干凈，反射波同相軸清晰，深層弱能量信號(hào)得到有效保護(hù)。

圖5 含有不同δ背景噪聲數(shù)據(jù)的去噪結(jié)果(a)δ=10； (b)δ=20； (c)δ=30； (d)δ=40

以上去噪試驗(yàn)表明，本文方法能壓制階梯效應(yīng)，去除不同強(qiáng)度隨機(jī)噪聲，提高信噪比，并能有效保護(hù)弱能量信號(hào)。

2.4 方法對(duì)比

為進(jìn)一步驗(yàn)證本文方法的有效性，分別與各向異性全變分(Anisotropic total variation，ATV)、交疊組稀疏全變分(Overlapping group sparsity total variation，OGSTV)、中值濾波(Median filter，MF)、三維塊匹配(Block matching 3D，BM3D)和K奇異值分解(K-singular value decomposition，KSVD)等五種方法進(jìn)行對(duì)比。其中，ATV和OGSTV與本文方法屬于同類(lèi)，后三者為常用去噪方法。

表1展示這些方法去噪后各項(xiàng)指標(biāo)對(duì)比結(jié)果。從SNR和SSIM指標(biāo)看，本文方法明顯優(yōu)于其他五種方法。從運(yùn)算耗時(shí)指標(biāo)看，本文方法用時(shí)雖然大于MF、ATV和OGSTV，但都未超過(guò)8s，在可接受范圍內(nèi)。對(duì)比后發(fā)現(xiàn)，用時(shí)明顯小于BM3D和KSVD，這是由于最后兩種方法涉及到耗時(shí)的非局部塊匹配和字典訓(xùn)練。

表1 針對(duì)四種含噪數(shù)據(jù)的六種方法去噪結(jié)果評(píng)價(jià)指標(biāo)對(duì)比

圖6展示δ=40時(shí)上述六種方法去噪結(jié)果。強(qiáng)噪聲背景下，ATV(圖6a)去噪并不理想，殘留明顯斑點(diǎn)和小塊，階梯效應(yīng)嚴(yán)重。OGSTV(圖6b)引入鄰域信息，階梯效應(yīng)得到一定壓制，但仍存部分斑點(diǎn)偽影。MF(圖6d)去噪效果最差，存在大量殘余噪聲。BM3D(圖6e)產(chǎn)生了平滑效應(yīng)，對(duì)深層弱信號(hào)的保護(hù)能力差。KSVD(圖6f)存在一定斑點(diǎn)，影響了最終去噪效果。圖6c所示剖面干凈，階梯效應(yīng)弱，深層反射波同相軸連續(xù)性好，說(shuō)明本文方法具有保護(hù)弱信號(hào)能力。

圖6 不同去噪方法對(duì)比(a)ATV； (b)OGSTV； (c)本文方法； (d)MF； (e)BM3D； (f)KSVD

從上述不同方法縱橫向?qū)Ρ瓤芍疚姆椒苡行コS機(jī)噪聲，壓制階梯效應(yīng)，更好地保護(hù)弱信號(hào)，在效率與精度上能達(dá)到良好的平衡。

3 實(shí)際資料應(yīng)用

將本文去噪方法應(yīng)用于M工區(qū)實(shí)際地震資料。所選單炮記錄(圖7a)共有240道，采樣點(diǎn)數(shù)為700，采樣間隔為2ms。由于存在大量隨機(jī)噪聲，反射波同相軸連續(xù)性差，弱能量信號(hào)被掩蓋，信噪比較低。其疊后剖面(圖7b)共有800道，采樣點(diǎn)數(shù)為800，采樣間隔為2ms。該剖面深部地層有一定起伏，且發(fā)育微斷層。隨機(jī)噪聲的存在致使剖面不清晰，同相軸錯(cuò)斷，給地質(zhì)解釋帶來(lái)諸多不便。

分別用上述六種方法進(jìn)行去噪處理并得到圖8(單炮)和圖9(剖面)所示結(jié)果。歷經(jīng)多次試驗(yàn)并優(yōu)選，針對(duì)圖7所示實(shí)際資料，本文方法參數(shù)設(shè)置為K=4、p=0.18和K=3、p=0.21。對(duì)比所得單炮記錄去噪結(jié)果可知，ATV(圖8a)和MF(圖8d)去噪效果差，信號(hào)失真嚴(yán)重，尤其是1.2～1.5s和2.1～2.5s區(qū)間的弱信號(hào)； 相對(duì)而言，OGSTV(圖8b)階梯效應(yīng)明顯減少，去噪效果有一定提升，但仍存殘余噪聲； BM3D(圖8e)能保留強(qiáng)能量信號(hào)，但弱信號(hào)恢復(fù)能力差，產(chǎn)生了嚴(yán)重的平滑效應(yīng)； KSVD(圖8f)和本文方法(圖8c)去噪效果相對(duì)較好，但KSVD在一些細(xì)節(jié)(弱信號(hào))上不如本文方法。

圖7 原始地震資料(a)單炮資料； (b)疊后剖面

圖8 不同單炮記錄去噪結(jié)果對(duì)比(a)ATV； (b)OGSTV； (c)本文方法； (d)MF； (e)BM3D； (f)KSVD

對(duì)比、分析圖9及表2可知，MF去噪效率高，但剖面(圖9d)上仍有大量噪聲殘留； OGSTV(圖9b)去噪效果有改善，但仍不理想； ATV(圖9a)去噪后產(chǎn)生的階梯效應(yīng)使部分細(xì)節(jié)丟失； 同樣地，BM3D(圖9e)出現(xiàn)平滑模糊區(qū)塊，不利于微斷裂等信息的恢復(fù)； KSVD(圖9f)和本文方法(圖9c)都有效壓制了隨機(jī)噪聲，但KSVD計(jì)算效率低。

表2 不同去噪方法的計(jì)算耗時(shí)對(duì)比

圖9 疊后資料去噪結(jié)果對(duì)比(a)ATV； (b)OGSTV； (c)本文方法； (d)MF； (e)BM3D； (f)KSVD

針對(duì)實(shí)際資料的應(yīng)用結(jié)果表明，本文方法能壓制不同強(qiáng)度隨機(jī)噪聲，減弱階梯效應(yīng)，有效保護(hù)信號(hào)細(xì)節(jié)特征，去噪后整個(gè)剖面同相軸具有更好清晰度和連續(xù)性，因此具有良好應(yīng)用潛力。

4 結(jié)論

(1)本文去噪方法適用于不同強(qiáng)度的背景噪聲，在有效壓制隨機(jī)噪聲的同時(shí)兼顧計(jì)算效率，能顯著提高地震資料的品質(zhì)，具有廣闊的應(yīng)用前景。

(2)OGS考慮了信號(hào)的鄰域信息，能充分挖掘局部相似性； 將它與高階TV模型和Lp偽范數(shù)結(jié)合，不僅能壓制階梯效應(yīng)，提高算法去噪能力，而且能有效保護(hù)圖像邊緣信息，具有較高保真度。

(3)K值決定鄰域信息的多少，若K值過(guò)小，則不能充分挖掘鄰域信息；反之，若引入非相似信息，則不僅增大了計(jì)算量，而且導(dǎo)致平滑效應(yīng)。p值關(guān)系到局部非零值信號(hào)的恢復(fù)，對(duì)于弱信號(hào)局部細(xì)節(jié)的保護(hù)有重要作用。