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

朱躍飛 曹靜杰* 殷晗鈞

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

0 引言

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

濾波類方法基于地震數據時間域分布特點構建濾波函數去除噪聲，主要方法有中值濾波[2-3]、各向異性擴散濾波[4]等。基于變換的方法假設地震數據經過某個變換后的系數具有稀疏特征，選取較大的系數，通過閾值運算去掉小的系數，最后反變換到時間域實現去噪[5]，常用的變換有傅里葉變換[6]、Radon變換[7]、Wavelet變換[8-10]、S變換[11]、曲波變換[12]等。

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

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

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

1 算法原理

1.1 DMSSA方法[18]

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

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

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

(1)

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

(2)

Rk表示由F(x,y,ωi)的第k行構造的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進行奇異值分解，并且選擇和截斷奇異值，是MSSA類方法的關鍵。如果有效信號對應的奇異值個數為N，則奇異值對角矩陣僅保留前N個奇異值，而其他所有奇異值均設置為零。對H進行奇異值分解，得到

(4)

其中

(5)

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

(4)基于截斷的奇異值計算去噪結果。通過將Hankel矩陣反變換到頻率域，再通過離散傅里葉逆變換得到時間—空間域去噪地震數據，即

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

(6)

(7)

對所有的頻率域數據進行上述操作即可得到去噪地震數據。

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

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

(8)

(9)

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

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

1.2 經驗公式法

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

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

(10)

式中：d為不含噪數據；r為去噪后數據。

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

通過區分數據含噪和不含噪情況的σ，從而確定數據塊中要保留的N。基于錯位相除的思路，提出一種確定數據塊中有效N的經驗估計方法。首先將Σ的σ排列成新的向量

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

(11)

對Q錯位相除，得

(12)

定義新的向量

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

(13)

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

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

(14)

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

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

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

(15)

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

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

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

圖4 由經驗公式法確定的圖1a的σ分布(a)全頻率掃描結果； (b)統計結果

1.3 ADMSSA算法

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

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

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

(16)

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

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

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

(17)

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

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

(18)

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

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

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

2 數值實驗

2.1 模擬數據實驗

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

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

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

能量較強的噪聲經常使地震信號發生嚴重畸變，導致塊Hankel矩陣σ變化復雜。為了驗證ADMSSA方法的有效性和對噪聲的敏感程度，對含噪地震數據(圖7b)去噪，結果表明，DMSSA方法確定的N為3(圖8紅色實線)，ADMSSA方法確定的N為7(圖8藍色實線)。

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

2.2 實際地震數據實例

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

圖8 兩種方法對圖7b數據確定的N

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

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

圖14為三維地震數據去噪效果對比。由圖可見，兩種方法去噪結果基本相同，無論是在地下結構較穩定區域還是在斷點附近，ADMSSA方法的去噪結果很好地保護了構造細節(圖14a)，同相軸的輪廓清晰，去噪效果很好。

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

圖11 圖10a的σ分布(a)DMSSA、ADMSSA方法確定的N； (b)DMSSA方法統計結果

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

圖13 圖10b的σ分布(a) DMSSA、ADMSSA方法確定的N； (b)DMSSA方法統計結果

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

3 討論

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

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

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

4 結論

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