一種自動判定保留的奇異值個數(shù)的地震隨機(jī)噪聲壓制算法

朱躍飛 曹靜杰* 殷晗鈞

(①自然資源部京津冀城市群地下空間智能探測與裝備重點(diǎn)實(shí)驗(yàn)室，河北石家莊 050031； ②河北省戰(zhàn)略性關(guān)鍵礦產(chǎn)資源重點(diǎn)實(shí)驗(yàn)室，河北石家莊 050031)

0 引言

受人類活動、儀器、環(huán)境、天氣等多種因素影響，野外采集的地震數(shù)據(jù)往往含有各種噪聲，嚴(yán)重影響速度分析和靜校正、速度建模及偏移成像等處理的效果。因此，消除噪聲以獲取高信噪比地震數(shù)據(jù)一直是地震勘探面臨的難題[1]。地震數(shù)據(jù)去噪一般依賴信號和噪聲在頻率、統(tǒng)計(jì)規(guī)律、振幅等方面的差異分離信號和噪聲。地震噪聲分為隨機(jī)噪聲和相干噪聲。地震隨機(jī)噪聲的消除方法眾多，大體上分為濾波類方法、基于變換的方法、降秩方法和深度學(xué)習(xí)方法等。

濾波類方法基于地震數(shù)據(jù)時間域分布特點(diǎn)構(gòu)建濾波函數(shù)去除噪聲，主要方法有中值濾波[2-3]、各向異性擴(kuò)散濾波[4]等。基于變換的方法假設(shè)地震數(shù)據(jù)經(jīng)過某個變換后的系數(shù)具有稀疏特征，選取較大的系數(shù)，通過閾值運(yùn)算去掉小的系數(shù)，最后反變換到時間域?qū)崿F(xiàn)去噪[5]，常用的變換有傅里葉變換[6]、Radon變換[7]、Wavelet變換[8-10]、S變換[11]、曲波變換[12]等。

深度學(xué)習(xí)去噪方法是目前的研究熱點(diǎn)，基本原理是利用大量的樣本數(shù)據(jù)的特征，通過多層卷積的方式提取數(shù)據(jù)的時域特征，然后采用深度學(xué)習(xí)的非線性逼近能力調(diào)整網(wǎng)絡(luò)參數(shù)，從而建立一個復(fù)雜的去噪模型實(shí)現(xiàn)去噪。目前卷積神經(jīng)網(wǎng)絡(luò)[13]、殘差學(xué)習(xí)[14-15]、生成對抗網(wǎng)絡(luò)[16]、降噪自編碼[17]等深度學(xué)習(xí)網(wǎng)絡(luò)被用于地震數(shù)據(jù)去噪。深度學(xué)習(xí)方法需對不同的數(shù)據(jù)大量訓(xùn)練，因此計(jì)算量大。

多道奇異譜分析(MSSA)是一種基于奇異值分解的降秩去噪方法，通過奇異值分解將原始數(shù)據(jù)分解為信號子空間和噪聲子空間，然后將噪聲子空間的能量置為零(截?cái)?，再通過反變換去噪[18]。MSSA用于多道時間序列分析，是單道奇異譜分析(SSA)的推廣[19-20]。Read[21]率先將SSA拓展到多變量MSSA方法研究，基于線性同相軸的假設(shè)，利用相鄰地震道的頻譜相似性與可預(yù)測性組成低秩的Hankel矩陣[22]，噪聲破壞了數(shù)據(jù)頻率切片Hankel矩陣的低秩結(jié)構(gòu)[23]，常用截?cái)嗥娈愔捣纸夥椒ń鉀Q低秩近似問題。在地震信號處理領(lǐng)域，MSSA和Cadzow濾波是等效的但卻來自不同的領(lǐng)域[24]，即Cadzow濾波源于信號和圖像去噪，MSSA則源于分析由動力系統(tǒng)引起的時間序列，本文采用MSSA的名稱表示這類方法。Oropeza等[25]利用MSSA同時對疊前三維數(shù)據(jù)去噪和重建，數(shù)值實(shí)驗(yàn)表明無法完全消除隨機(jī)噪聲，其去噪效果有很大的提升空間。Huang等[26]將阻尼算子引入傳統(tǒng)MSSA中，提出了阻尼多道奇異譜分析(DMSSA)算法。通過融合軟閾值移動平均算子和阻尼算子的優(yōu)點(diǎn)，Oboue等[27]利用魯棒阻尼降秩方法提高地震數(shù)據(jù)的信噪比。阻尼降秩方法已成為一種有效的去噪方法，可以從含噪和不完備的觀測數(shù)據(jù)中恢復(fù)有效信號。

對于海量地震數(shù)據(jù)來說，基于降秩的方法需要將地震數(shù)據(jù)分成不同的塊，然而每個塊對應(yīng)的奇異值個數(shù)不同，目前需要人工估計(jì)每個塊的有效奇異值個數(shù)，計(jì)算效率低，無法實(shí)現(xiàn)產(chǎn)業(yè)化。為此，本文利用Akaike信息準(zhǔn)則自動地確定地震信號的奇異值個數(shù)，然后基于DMSSA方法去噪。首先介紹了MSSA方法的去噪原理，然后給出確定有效奇異值個數(shù)的Akaike信息準(zhǔn)則和經(jīng)驗(yàn)方法，經(jīng)驗(yàn)方法可以驗(yàn)證Akaike信息準(zhǔn)則的有效性。模擬和實(shí)際數(shù)據(jù)實(shí)驗(yàn)表明，Akaike信息準(zhǔn)則能夠自動確定有效奇異值個數(shù)，避免了人工操作，有利于實(shí)現(xiàn)產(chǎn)業(yè)化。

1 算法原理

1.1 DMSSA方法[18]

假設(shè)一個含噪三維地震數(shù)據(jù)為D(x,y,t)，其中x=(x1,…,xm)、y=(y1,…,yn)表示空間坐標(biāo)，t=(t1,…,ts)表示時間坐標(biāo)，m、n為道數(shù)，s為采樣點(diǎn)數(shù)。根據(jù)DMSSA理論，使用以下步驟去噪。

(1)通過離散傅里葉變換將D(x,y,t)從時間域變換為頻率域數(shù)據(jù)F(x,y,ω)，其中ω=(ω1,…,ωj)為離散的頻率序列，j為頻率切片個數(shù)。

(2)在給定的頻率范圍內(nèi)將不同頻率切片數(shù)據(jù)排列成塊Hankel矩陣。當(dāng)頻率為ωi(i=1,…,j)時，有

(1)

首先，將F(x,y,ωi)的每一行構(gòu)造成Hankel矩陣

(2)

Rk表示由F(x,y,ωi)的第k行構(gòu)造的Hankel矩陣，大小為v×h，v=

n/2

+1，h=n-v+1，

·

表示向下取整，其中2

(3)

H為(v×l)×(h×f)階塊Hankel矩陣，l=

m/2

+1，f=m-l+1，2

(3)對H進(jìn)行奇異值分解，并且選擇和截?cái)嗥娈愔?，是MSSA類方法的關(guān)鍵。如果有效信號對應(yīng)的奇異值個數(shù)為N，則奇異值對角矩陣僅保留前N個奇異值，而其他所有奇異值均設(shè)置為零。對H進(jìn)行奇異值分解，得到

(4)

其中

(5)

式中：U為H的左奇異值向量組成的(v×l)×(v×l)階正交矩陣；VT為H的右奇異值向量組成的 (h×f)×(h×f)階正交矩陣；Σ為按奇異值遞減順序σ1≥σ2≥…≥σd組成的對角矩陣，非零奇異值的個數(shù)d等于H的秩。

(4)基于截?cái)嗟钠娈愔涤?jì)算去噪結(jié)果。通過將Hankel矩陣反變換到頻率域，再通過離散傅里葉逆變換得到時間—空間域去噪地震數(shù)據(jù)，即

Σ1=diag(σ1,σ2,…,σN) 0≤N

(6)

(7)

對所有的頻率域數(shù)據(jù)進(jìn)行上述操作即可得到去噪地震數(shù)據(jù)。

若采集的地震數(shù)據(jù)中不含噪聲，則Σ僅包含與有效信號相關(guān)的非零σ。若采集的地震數(shù)據(jù)中含有噪聲，所有σ都會發(fā)生改變，非零σ的個數(shù)也將增加。原始MSSA方法僅保留了N個σ，對σ的大小并沒有影響，因此去噪結(jié)果有很大的改進(jìn)空間。Huang等[18]提出的DMSSA方法可以減小σ，因此去噪效果更好。

DMSSA方法通過添加阻尼因子減弱由噪聲引起的σ增量，即

(8)

(9)

式中：T為阻尼算子；I為單位矩陣；D為阻尼因子，其值越小，阻尼效果越強(qiáng)，反之亦然。DMSSA去噪的本質(zhì)就是利用D對第N+1個σ放大或縮小，然后使用前N個σ與其求差，并對第N+1之后的σ置零，以達(dá)到壓制噪聲的目的。

確定N是DMSSA去噪最關(guān)鍵的一步，將影響噪聲抑制效果和有效信號的保護(hù)程度。式(4)和式(9)是假定N已知的情況得到的，如果選擇N太小，將損壞有效信號； 如果選擇N太大，將降低噪聲壓制效果。對于MSSA類方法，自動確定N是關(guān)鍵。面對復(fù)雜多變的實(shí)際地震數(shù)據(jù)，在沒有充足的地質(zhì)資料時，確定數(shù)據(jù)塊中需要保留的N是一個值得研究的問題。實(shí)際數(shù)據(jù)去噪時要對數(shù)據(jù)分塊，需要人工估計(jì)每個塊的有效N，計(jì)算效率低，無法實(shí)現(xiàn)產(chǎn)業(yè)化。為此，首先引入一種確定有效N的經(jīng)驗(yàn)估計(jì)方法，然后給出了一種自動確定N的方法，該方法基于Akaike信息準(zhǔn)則自動確定有效N，有利于MSSA產(chǎn)業(yè)化。

1.2 經(jīng)驗(yàn)公式法

圖1為模擬地震數(shù)據(jù)A三維視圖。模擬數(shù)據(jù)使用主頻為40Hz的雷克子波作為震源，信噪比(SNR)定義為

圖1 模擬地震數(shù)據(jù)A三維視圖(a)無噪數(shù)據(jù)； (b)加入10%隨機(jī)噪聲數(shù)據(jù)(信噪比為-1.322dB)數(shù)據(jù)含有5個不同傾角的地震同相軸，時間方向有300個采樣點(diǎn)，采樣率為2ms，主測線和聯(lián)絡(luò)測線方向均為60個采樣點(diǎn)

(10)

式中：d為不含噪數(shù)據(jù)；r為去噪后數(shù)據(jù)。

圖2為圖1數(shù)據(jù)在ω12處的σ曲線。由圖可見，地震數(shù)據(jù)在無噪聲時僅出現(xiàn)少數(shù)較大的σ且個數(shù)為N，其他σ均較小，因此σ曲線出現(xiàn)明顯的彎折現(xiàn)象(圖2b紅線)。從理論上來說，無噪時N以后的σ應(yīng)該全部為零。數(shù)據(jù)加入噪聲后，σ發(fā)生改變的同時也出現(xiàn)大量非零σ，因此加噪后的σ曲線下降相對平緩，但在第N和第N+1項(xiàng)之間，依然存在巨大的落差。

通過區(qū)分?jǐn)?shù)據(jù)含噪和不含噪情況的σ，從而確定數(shù)據(jù)塊中要保留的N。基于錯位相除的思路，提出一種確定數(shù)據(jù)塊中有效N的經(jīng)驗(yàn)估計(jì)方法。首先將Σ的σ排列成新的向量

Q=(σ1,σ2,…,σd)

(11)

對Q錯位相除，得

(12)

定義新的向量

Q′=(q1,q2,…,qd-1)

(13)

式中Q′是基于Q的元素錯位相除構(gòu)造的新向量，其第i個元素等于Q的第i+1個元素除以第i個元素。若σi和σi+1的值較接近，則qi的值趨于1，否則將出現(xiàn)一個較小的值。對Q′取極小，并獲取極小值對應(yīng)的索引，有

N=arg minQ′=(q1,q2,…,qd-1)

(14)

在實(shí)際處理中，只需截取Q′的前若干項(xiàng)即可。圖3為錯位相除向量曲線。由圖可見，Q′的第5個元素值明顯小于第4個元素和第6個元素，說明在該頻率下應(yīng)該保留5個奇異值(圖3b)。

為了排除偶然性，使用上述方法對數(shù)據(jù)塊的主頻數(shù)據(jù)進(jìn)行相同的計(jì)算，并得到新向量

flag1=(arg minQ′ωL,…,arg minQ′ωH)

(15)

其中ωL～ωH為信號的主要頻率范圍。通過統(tǒng)計(jì)flag1中不同的N出現(xiàn)的占比，根據(jù)經(jīng)驗(yàn)選擇合適的N。圖4為由經(jīng)驗(yàn)公式法確定的圖1a的σ分布。由圖可見：不論是否含有噪聲，在主頻范圍內(nèi)N均為5(圖4a)； 含噪聲數(shù)據(jù)和不含噪聲數(shù)據(jù)N=5的占比分別為0.5、0.35，故確定N為5(圖4b)。因此，無論數(shù)據(jù)中是否包含噪聲，利用經(jīng)驗(yàn)公式法均可準(zhǔn)確估計(jì)N。經(jīng)驗(yàn)公式法所用的錯位相除策略和一階差分具有異曲同工的效果，可以驗(yàn)證其他方法的效果，為確定N提供了有效工具。

圖2 圖1數(shù)據(jù)在ω12處的σ曲線(a)全部σ分布； (b)圖a前30個σ分布局部放大

圖3 錯位相除向量曲線(a)完整曲線； (b)圖a前10項(xiàng)元素的局部放大

圖4 由經(jīng)驗(yàn)公式法確定的圖1a的σ分布(a)全頻率掃描結(jié)果； (b)統(tǒng)計(jì)結(jié)果

1.3 ADMSSA算法

處理海量地震數(shù)據(jù)時需要分塊地震數(shù)據(jù)，然后對每個數(shù)據(jù)塊去噪。但是每個數(shù)據(jù)塊對應(yīng)的N并不相等，人工預(yù)估方法不利于算法實(shí)施，因此需要研究自適應(yīng)算法。

對σ序列Q可依據(jù)第N和第N+1個σ之間的巨大落差并伴隨嚴(yán)重的彎折現(xiàn)象確定有效N。事實(shí)上，N值的選擇就是檢測σ序列中的拐點(diǎn)位置。本文利用Akaike信息準(zhǔn)則[28-29]自動判定保留的N。首先對Q作如下變換

fωi(σμ)=σμ+1-2σμ+σμ-11≤μ≤d

(16)

式(16)實(shí)際上是求σ序列曲線的二階導(dǎo)數(shù)，fωi(σμ)描述σ曲線斜率的變化率。確定頻率為ωi、第R點(diǎn)的有效N值的Akaike信息準(zhǔn)則為

AICωi(R)=Rlg[var(fωi[σi,σR])]+

(d-R-1)lg[var(fωi[σR+1,σd])]

(17)

式中：var表示數(shù)據(jù)序列的方差； AICωi(R)是長度為d的序列，其全局最小值對應(yīng)的位置即為拐點(diǎn)，按

flag2=(arg minAICωL,…,arg min AICωH)

(18)

求出所有頻率中的最小值。式(18)中元素最小值即是整個數(shù)據(jù)塊中需要保留的N。圖5為由ADMSSA算法確定的圖1a的σ分布。由圖可見：①在N=8時AICωi(R)曲線取極小值(圖5a)。由于主頻范圍內(nèi)數(shù)值結(jié)果較穩(wěn)定，其他范圍則經(jīng)常出現(xiàn)異常值，為了提高精度將頻率控制在10～90Hz內(nèi)。②ADMSSA算法確定的N為8(圖5b)，經(jīng)驗(yàn)公式法確定的N為5(圖4b)，證明利用基于Akaike信息準(zhǔn)則的方法去噪，可以自動地估計(jì)出與真實(shí)值較接近的N。

在使用上述方法確定N后，去噪過程采用DMSSA方法的框架，在ADMSSA算法中僅需要確定信號的主頻范圍，就可以自動地去噪。

圖5 由ADMSSA算法確定的圖1a的σ分布(a)信息準(zhǔn)則曲線； (b)信息準(zhǔn)則確定的結(jié)果

2 數(shù)值實(shí)驗(yàn)

2.1 模擬數(shù)據(jù)實(shí)驗(yàn)

分別使用ADMSSA、DMSSA方法對圖1b去噪，結(jié)果(圖6)表明：由ADMSSA、DMSSA方法確定的N分別為8、5，令阻尼因子D=3，去噪結(jié)果的信噪比分別為22.110(圖6a)、22.438dB(圖6c)，兩者的去噪效果較接近，信噪比較高，同相軸清晰連貫，局部細(xì)節(jié)得以保留。

圖6 圖1b的去噪效果對比(a)ADMSSA方法去噪結(jié)果； (b)ADMSSA方法去除的噪聲； (c) DMSSA方法去噪結(jié)果； (d) DMSSA方法去除的噪聲

圖7 模擬地震數(shù)據(jù)B三維視圖(a)無噪數(shù)據(jù)； (b)加入15%隨機(jī)噪聲數(shù)據(jù)(信噪比為-4.659dB)數(shù)據(jù)含有4個不同傾角的地震同相軸，時間方向有300個采樣點(diǎn)，主測線和聯(lián)絡(luò)測線方向均為60個采樣點(diǎn)，采用主頻為40Hz的雷克子波正演

能量較強(qiáng)的噪聲經(jīng)常使地震信號發(fā)生嚴(yán)重畸變，導(dǎo)致塊Hankel矩陣σ變化復(fù)雜。為了驗(yàn)證ADMSSA方法的有效性和對噪聲的敏感程度，對含噪地震數(shù)據(jù)(圖7b)去噪，結(jié)果表明，DMSSA方法確定的N為3(圖8紅色實(shí)線)，ADMSSA方法確定的N為7(圖8藍(lán)色實(shí)線)。

圖9為圖7b的去噪效果對比。由圖可見，令阻尼因子D=3，ADMSSA、DMSSA方法去噪結(jié)果的信噪比分別為21.263(圖9a)、21.778dB(圖9c)，即前者的信噪比略低，但去噪結(jié)果的同相軸清晰，噪聲殘留較少(圖9b)，說明ADMSSA方法高效、精確。

2.2 實(shí)際地震數(shù)據(jù)實(shí)例

為了證明ADMSSA算法對實(shí)際地震數(shù)據(jù)的去噪效果，分別使用二維和三維疊后地震數(shù)據(jù)(圖10)驗(yàn)證。二維地震數(shù)據(jù)(圖10a)信噪比低，地震同相軸連續(xù)性差，隨機(jī)噪聲能量強(qiáng)，剖面中間部分以及下部存在斷層，尤其是下部存在多處斷裂構(gòu)造。三維地震數(shù)據(jù)(圖10b)信噪比低，噪聲能量較強(qiáng)，有效信號被噪聲嚴(yán)重污染，中間部分同相軸出現(xiàn)彎曲、斷裂現(xiàn)象。

圖8 兩種方法對圖7b數(shù)據(jù)確定的N

圖9 圖7b的去噪效果對比(a)ADMSSA方法去噪結(jié)果； (b)ADMSSA方法去除的噪聲； (c)DMSSA方法去噪結(jié)果； (d)DMSSA方法去除的噪聲

圖11為圖10a的σ分布。由圖可見：ADMSSA方法確定的N為8(圖11a藍(lán)色實(shí)線)； DMSSA方法統(tǒng)計(jì)N出現(xiàn)的百分比(圖11b)確定的N為10。圖12為二維地震數(shù)據(jù)去噪效果對比。由圖可見，令阻尼因子D=5，地震同相軸邊緣刻畫清晰，噪聲去除徹底，對構(gòu)造細(xì)節(jié)保護(hù)較好，去噪效果均較好。圖13為圖10b的σ分布。由圖可見：ADMSSA方法確定的N為5(圖13a藍(lán)色實(shí)線)； DMSSA方法統(tǒng)計(jì)N出現(xiàn)的百分比確定的N為4(圖13b)。

圖14為三維地震數(shù)據(jù)去噪效果對比。由圖可見，兩種方法去噪結(jié)果基本相同，無論是在地下結(jié)構(gòu)較穩(wěn)定區(qū)域還是在斷點(diǎn)附近，ADMSSA方法的去噪結(jié)果很好地保護(hù)了構(gòu)造細(xì)節(jié)(圖14a)，同相軸的輪廓清晰，去噪效果很好。

圖10 二維(a)、三維(b)地震數(shù)據(jù)

圖11 圖10a的σ分布(a)DMSSA、ADMSSA方法確定的N； (b)DMSSA方法統(tǒng)計(jì)結(jié)果

圖12 二維地震數(shù)據(jù)去噪效果對比(a)ADMSSA方法去噪結(jié)果； (b)ADMSSA方法去除的噪聲； (c)DMSSA方法去噪結(jié)果； (d)DMSSA方法去除的噪聲

圖13 圖10b的σ分布(a) DMSSA、ADMSSA方法確定的N； (b)DMSSA方法統(tǒng)計(jì)結(jié)果

圖14 三維地震數(shù)據(jù)去噪效果對比(a)ADMSSA方法去噪結(jié)果； (b)ADMSSA方法去除的噪聲； (c)DMSSA方法去噪結(jié)果； (d)DMSSA方法去除的噪聲

3 討論

對于MSSA類方法來說，劃分的數(shù)據(jù)塊的尺度對計(jì)算時間和去噪效果均有影響。ADMSSA方法同樣受到劃分?jǐn)?shù)據(jù)塊尺度的影響，因此在劃分?jǐn)?shù)據(jù)塊時應(yīng)先做試驗(yàn)，再確定劃分尺度。ADMSSA方法需要選擇一個合適的頻率范圍，一般取信號有效頻率范圍即可取得較好去噪結(jié)果。DMSSA方法受阻尼因子D的影響，當(dāng)噪聲能量較強(qiáng)時，選擇D約為3，當(dāng)噪聲信號較弱時，應(yīng)選擇較大的D值。ADMSSA方法對強(qiáng)脈沖噪聲具有一定壓制作用，但是去噪結(jié)果中仍殘存一些強(qiáng)脈沖噪聲。MSSA方法對奇異值個數(shù)N依賴很強(qiáng)，當(dāng)N較小時，會影響斷層的識別。對于曲率較大的彎曲同相軸，需要合理劃分?jǐn)?shù)據(jù)塊尺度。MSSA方法的計(jì)算量主要為奇異值分解，為了提升算法的計(jì)算效率，可以采用隨機(jī)奇異值分解等方法。

DMSSA方法和ADMSSA方法確定的N存在差別，這是由于前者依靠經(jīng)驗(yàn)公式法估計(jì)N，存在一定誤差，而ADMSSA方法依靠數(shù)據(jù)信息通過概率分析提取信號特征。

降秩類去噪算法在消除隨機(jī)噪聲的同時會對有效信號造成一定損傷，其原因是該類方法基于線性同相軸假設(shè)，實(shí)際數(shù)據(jù)很難滿足條件。因此，盡管降秩類去噪算法獲得了較好效果，但是去噪效果尚有很大的提升空間。

4 結(jié)論

對于MSSA類去噪方法來說，確定有效的奇異值個數(shù)是關(guān)鍵。目前都是依靠人工經(jīng)驗(yàn)估計(jì)奇異值個數(shù)，不利于該類方法的產(chǎn)業(yè)化。本文提出了一種確定有效奇異值個數(shù)的方法，該方法基于Akaike信息準(zhǔn)則自動區(qū)分有效信號對應(yīng)的奇異值與噪聲相關(guān)的奇異值，克服了人工選擇奇異值個數(shù)的問題，有利于海量地震數(shù)據(jù)去噪。此外，本文還提出利用經(jīng)驗(yàn)公式法驗(yàn)證ADMSSA方法的可靠性。數(shù)值實(shí)驗(yàn)證明，ADMSSA方法能夠自動地確定可靠的奇異值個數(shù)，并且獲得高信噪比的去噪結(jié)果，該算法在工業(yè)化應(yīng)用中具有巨大潛力。