基于自適應(yīng)秩選方法的地震降秩去噪

楊志鵬, 陳秀清, 李 勇, 陳力鑫, 龐 謙, 馬澤川, 王登偉

(1.四川省地震局 西昌地震中心站，四川 西昌 615022；2.成都理工大學(xué) 地球物理學(xué)院，成都 610059；3.四川省地震局 地震監(jiān)測中心，成都 610041)

人工地震勘探是資源勘查、活斷層探測、地質(zhì)構(gòu)造分析的重要手段[1-3]?？碧江h(huán)境的復(fù)雜化和勘探深度的增加致使采集得到的地震數(shù)據(jù)中有效反射信號較弱，且往往被強隨機噪聲覆蓋[4]。地震隨機噪聲一般具有分布頻帶寬、振動無規(guī)律、隨機分布于整個地震道等特征，由地下復(fù)雜介質(zhì)構(gòu)造和其他不確定因素共同作用產(chǎn)生。地震隨機噪聲會對后續(xù)的諸如多次波衰減、地震波反演成像、地震道重建、地質(zhì)構(gòu)造分析等數(shù)據(jù)處理與解釋流程產(chǎn)生嚴(yán)重影響[5-6]，因此壓制地震隨機噪聲是地震資料弱信號提取與增強的重要步驟。目前，主流的地震隨機噪聲壓制方法可以分為5個大類[7-8]：基于域變換和濾波理論方法[9-10]、基于稀疏變換理論方法[11-13]、基于信號分解理論方法[14-16]、基于機器學(xué)習(xí)理論方法[17-18]、基于降秩理論方法[19-25]。

基于降秩理論的去噪方法主要是對地震數(shù)據(jù)構(gòu)建特殊數(shù)學(xué)矩陣，再利用低秩逼近(秩約化)方法對其進行降秩處理，達到去除隨機干擾的目的[26]。對頻域地震數(shù)據(jù)構(gòu)造Hankel矩陣是常見做法之一，并相繼發(fā)展出利用截斷奇異值分解(TSVD)方法[19]、阻尼奇異值分解(DSVD)方法[20]、阻尼最優(yōu)異值收縮(DOptShrink)方法[24]，低秩-稀疏混合約束(HRSC)方法[27]、結(jié)構(gòu)低秩約束(HLR)方法[22]等秩約化方法對構(gòu)造的Hankel矩陣進行降秩去噪，其中DSVD方法和DOptShrink方法還被拓展應(yīng)用于五維地震資料的同步去噪與重建[21,28]。

基于Hankel矩陣的降秩去噪結(jié)果除了受秩約化方法限制外，另一關(guān)鍵因素為秩選參數(shù)，即純凈數(shù)據(jù)的秩值。尤其是在處理復(fù)雜構(gòu)造的含噪地震資料及其分塊計算時，由于純凈數(shù)據(jù)未知，秩約化方法的秩選參數(shù)難以準(zhǔn)確測定，容易導(dǎo)致去噪不足或去噪過度問題[29]。為解決傳統(tǒng)的固定秩選參數(shù)難以確定的問題，本文引入了一種基于奇異值最優(yōu)硬閾值的自適應(yīng)秩選方法[30-31]，通過最優(yōu)硬閾值系數(shù)和估計數(shù)據(jù)噪聲水平計算出信噪分界奇異值大小，從而自適應(yīng)地確定秩選參數(shù)。本文將其與TSVD方法、DSVD方法和DOptShrink方法相結(jié)合，提出自適應(yīng)降秩去噪方法(ATSVD、ADSVD、ADOptShrink)，利用理論模型和實際資料證明所提自適應(yīng)降秩去噪方法從壓噪性和保幅性上均較傳統(tǒng)的基于固定秩選的去噪方法有所提高，且從計算效率等方面對新方法進行了探討。

1 基于自適應(yīng)秩選方法的地震數(shù)據(jù)去噪

1.1 Hankel矩陣降秩去噪

Hankel矩陣為沿反對角線元素相等的特殊數(shù)學(xué)矩陣，被廣泛應(yīng)用于矩陣降秩去噪、重建等領(lǐng)域?？紤]時域三維地震數(shù)據(jù)可以表示為

Dtime(x,y,t)∈RNx×Ny×Nt

其中：Nx為Inline測線道數(shù)；Ny為Xline測線道數(shù)；Nt為時間軸采樣點數(shù)。首先通過傅里葉變換將Dtime(x，y，t)由時域轉(zhuǎn)換到頻域

(1)

式中：Dfreq(x,y,ω)∈RNx×Ny×Nω表示頻域數(shù)據(jù)；Nω為頻率切片數(shù)。對于每個特定的頻率切片Dfreq(jΔω)∈RNx×Ny(1≤j≤Nω)可以表示為矩陣

(2)

對該頻率切片分量Dfreq(jΔω)中每行數(shù)據(jù)構(gòu)造Hankel矩陣，如第i(1≤i≤Nx)行數(shù)據(jù)構(gòu)建的Hankel矩陣可以表示為[30]

(3)

(4)

理論證明，對于由線性同相軸構(gòu)成的無噪地震數(shù)據(jù)，其每個頻率切片構(gòu)造的塊Hankel矩陣H具有低秩性質(zhì)，且秩的大小等于剖面中不同傾角線性同相軸的個數(shù)[30,32-33]，數(shù)據(jù)的缺失或隨機噪聲的存在都會顯著增加矩陣H的秩值大小(圖2)。因此通過對塊Hankel矩陣H的進行秩約化處理，可以壓制隨機噪聲，恢復(fù)有效信號[25]。

1.2 秩約化方法

1.2.1 TSVD方法

考慮觀測矩陣H∈Ra×b由信號分量X∈Ra×b(rank(X) =r

H=X+N

(5)

通常，低秩矩陣去噪問題在最小二乘意義下的一般形式可以表示為

(6)

(7)

(8)

(7)式為觀測矩陣H的奇異值分解形式，xi為H的奇異值，ui和νi分別為H的左、右奇異向量。

1.2.2 DSVD方法

(9)

(10)

式中：I為單位陣；k為控制Ti衰減程度的阻尼參數(shù)。

1.2.3 DOptShrink方法

為了克服直接使用核范數(shù)對秩函數(shù)進行凸松弛近似無法得到最優(yōu)解的問題，R.R.Nadakuditi[34]提出利用觀測矩陣低秩特性先驗，通過TSVD逼近觀測矩陣非凸秩函數(shù)的優(yōu)化方法OptShrink，利用門限閾值對奇異值進行收縮處理?？紤]建立非凸優(yōu)化去噪問題，將(6)式轉(zhuǎn)化為一般形式

(11)

(12)

(13)

式中：xi為H的第i(1≤i≤r)個奇異值;Σ為由奇異值序列xr+1,xr+2, …,xa-1,xa構(gòu)成的對角陣;D(·)和D′(·)為互逆的特定變換算子，分別定義為

(14)

(15)

利用(9)式的阻尼項進一步約束OptShrink，結(jié)果可得到DOptShrink方法降秩矩陣[24,28]

(16)

1.3 自適應(yīng)秩選方法

為解決上述問題，本文在秩約化方法中引入了基于奇異值最優(yōu)硬閾值的自適應(yīng)秩參數(shù)選擇方法[35-36]。該方法在漸近均方誤差(asymptotic mean squared error, AMSE)框架下，利用大維隨機矩陣?yán)碚摵蜐u近矩陣?yán)碚撛O(shè)計奇異值最優(yōu)硬閾值的選取，進而自適應(yīng)地選取秩參數(shù)。假設(shè)Z∈Ra×b為標(biāo)準(zhǔn)正態(tài)分布的高斯白噪聲，σ為噪聲的標(biāo)準(zhǔn)差(噪聲水平)，則(5)式模型可以簡化為

H=X+σZ

(17)

(18)

令β=a/b(0<β≤1)，矩陣維數(shù)a和b趨于無窮大，而X的秩和非零奇異值固定，即

(19)

基于(18)式、(19)式的漸近模型，則此時奇異值硬閾值系數(shù)λ的AMSE框架可定義為

(20)

(21)

(22)

(23)

(24)

但當(dāng)噪聲水平σ未知條件下，在應(yīng)用(24)式時還應(yīng)首先對σ進行噪聲水平估計

(25)

式中:xmed是觀測矩陣H的奇異值序列的中值;μβ是β基于隨機矩陣?yán)碚摰腗arcenko-Pastur分布定律的中值解，即求解

(26)

聯(lián)立(23)式、(24)式、(25)式、(26)式可以得到在噪聲水平σ未知情況下奇異值最優(yōu)硬閾值τ*為

(27)

最后，通過判斷奇異值xi是否大于最優(yōu)硬閾值τ*來判斷自適應(yīng)秩選參數(shù)r*的大小

r*=#{i,xi≥τ*}

(28)

1.4 自適應(yīng)降秩去噪方法

綜上分析，將自適應(yīng)秩選方法分別與常規(guī)的TSVD、DSVD、DOptShrink等秩約化方法相結(jié)合，分別提出了自適應(yīng)截斷奇異值分解(ATSVD)、自適應(yīng)阻尼奇異值分解(ADSVD)、自適應(yīng)阻尼最優(yōu)奇異值收縮(ADOptShrink)等自適應(yīng)降秩去噪方法。

自適應(yīng)降秩去噪方法的去噪流程如下：

(1)首先將三維地震數(shù)據(jù)體Dtime(x,y,t)進行分塊預(yù)處理，時空窗滑動，讀取第i個時空窗內(nèi)子數(shù)據(jù)塊Dtime(x,y,t)i。

(2)利用傅里葉變換將Dtime(x,y,t)i從時域轉(zhuǎn)換到頻域Dfreq(x,y,ω)i，對其每個頻率切片Dfreq(x,y,jΔω)i構(gòu)造塊Hankel矩陣H。

(5)時空窗繼續(xù)滑動，讀取第i+1個時空窗內(nèi)子數(shù)據(jù)塊Dtime(x,y,t)i+1，重復(fù)第(2)-第(4)步操作，直至分塊子數(shù)據(jù)全部處理完畢。

2 模型與算例

為定量評估上述方法的去噪性能，本文采用信噪比(RSN)和局部相似屬性(local similarity, LS)作為衡量去噪效果的標(biāo)準(zhǔn)[37-38]。

(29)

局部相似屬性是基于隨機噪聲與有效信號局部正交的先驗，可用于無參考場景的壓噪效果評價及檢測是否存在有效信號泄漏，其數(shù)值在區(qū)間[0, 1]變化，值越大的區(qū)域說明局部相似程度越高，有效信號的能量泄漏問題就越嚴(yán)重，即方法保幅能力越差[7,30,38]。同維矩陣A和B之間的局部相似屬性LS定義為

(30)

(31)

(32)

其中：a=diag(A),b=diag(B),diag(·)為矩陣取主對角線上元素的算子。(30)式和(31)式中的最小平方目標(biāo)函數(shù)優(yōu)化問題可由帶局部光滑約束的整形正則化方法求解

(33)

(34)

其中λ是控制物理維度和收斂速度的參數(shù)，Γ為整形算子。

2.1 線性模型測試

首先將本文的方法應(yīng)用于雷克子波構(gòu)建的三維線性模型[39](圖7-A),該模型包含4條不同主頻的線性同相軸,數(shù)據(jù)大小為300×40×40(時間采樣點數(shù)×Inline道數(shù)×Xline道數(shù)),采樣間隔為4 ms。在模擬信號基礎(chǔ)上加入50%的高斯白噪聲,使得該線性模型的信噪比為 -10.780 2 dB，得到含噪數(shù)據(jù)(圖7-B)，可以看出有效信號被強噪聲所覆蓋，同相軸細(xì)節(jié)變得十分模糊，很難識別出有效信號。由于該模型嚴(yán)格滿足同相軸為線性的假設(shè)條件，因此可對該模型數(shù)據(jù)整體直接應(yīng)用降秩去噪方法進行處理，選取固定秩選參數(shù)r為4，阻尼參數(shù)k為5。實驗中，我們使用傳統(tǒng)的基于固定秩選參數(shù)的降秩去噪方法TSVD、DSVD、DOptShrink和本文所提基于自適應(yīng)秩選方法的降秩去噪方法ATSVD、ADSVD、ADOptShrink作為對比研究，結(jié)果如圖7-C、E、G、D、F和H所示。從圖7的各種方法縱向?qū)Ρ葋砜?，基于本文改進方法的去噪結(jié)果背景更為干凈，殘余噪聲更少；橫向?qū)Ρ葋砜?，DSVD和DOptshrink去噪結(jié)果明顯優(yōu)于TSVD的去噪結(jié)果，ADSVD和ADOptShrink去噪結(jié)果也優(yōu)于ATSVD。

表1給出了6種方法處理后的信噪比，也可以看出本文改進方法相較于原方法在信噪比上均有所提高，其中ADSVD方法和ADOptShrink方法信噪比分別提高到 13.211 5 dB和 13.304 9 dB。從圖8各種方法的局部相似屬性剖面可以看出，全局均接近于零值，指示各種方法在處理由線性同相軸構(gòu)成的地震數(shù)據(jù)時，均具有較好的保幅性。

表1 三維線性模型各種方法去噪前后的信噪比Table 1 Signal to noise ratio (SNR) of de-noising results by various methods of 3D linear model

圖9為模擬信號在加入不同強度高斯白噪聲后，各種方法去噪后的信噪比，可以看到DOptShrink方法、ADSVD方法、ADOptShrink方法去噪后的信噪比最高，且十分接近。圖10為各種方法在不同秩選參數(shù)下，對圖7-B含噪數(shù)據(jù)處理后的信噪比?？梢钥闯鰝鹘y(tǒng)降秩去噪方法需在秩選參數(shù)r嚴(yán)格取為同相軸的傾角數(shù)時，才能取得最高信噪比；此后隨著r的增大，更多的噪聲成分會被引入去噪結(jié)果，導(dǎo)致信噪比下降。但相對而言DOptShrink方法下降趨勢更緩，因為該方法對秩選參數(shù)并不敏感[28]，而本文方法結(jié)果則完全不受固定值r影響。圖11為在ubuntu環(huán)境下利用Matlab2016b計算平臺，配置為Intel Core i5-8400CPU@2.80 GHz、4GB RAM時，統(tǒng)計的各種方法在不同秩選參數(shù)下的運算時間。由于降秩去噪方法的時間復(fù)雜度集中在奇異值分解步驟上，因此引入自適應(yīng)秩選方法并不會顯著增加運算時間。另一方面，可以看到DOptShrink方法結(jié)果雖然對秩選參數(shù)r不敏感，但其運算時間會隨r的增大而顯著增加，因此為獲得較好的去噪效果而對其選擇一個較大的秩值并不理想。

2.2 雙曲模型測試

為了測試本文的方法，在處理由非線性同相軸構(gòu)成的地震數(shù)據(jù)時的去噪表現(xiàn)，引入一個合成的三維雙曲模型[40](圖12-A)。該三維合成數(shù)據(jù)中包含5個不同曲率的雙曲線同相軸，數(shù)據(jù)大小為128×101×101，同樣在模擬信號基礎(chǔ)上加入50%高斯白噪聲，得到含噪數(shù)據(jù)(圖12-B)，由于噪音的加入，有效波的同相軸能量被嚴(yán)重干擾，只剩下大致輪廓可見。為了滿足線性同相軸假設(shè)，提高處理效果，將該模型分割為50×50×50的計算子窗口。由于受秩非一致性問題影響，對各個不同的處理子窗口，只能經(jīng)驗性地選取統(tǒng)一的固定秩選參數(shù)r為15，阻尼參數(shù)k為5。各種方法處理結(jié)果如圖12-C、D、E、D、F、G和H所示。從圖中可見TSVD方法去噪結(jié)果仍存在有大量的殘余噪聲，去噪效果最差；而DSVD方法和OptShrink方法的去噪結(jié)果雖然較好地去除了背景噪聲，但受計算窗口內(nèi)不當(dāng)?shù)闹冗x參數(shù)影響，同相軸在曲部發(fā)生了扭曲和斷裂，且信號邊緣存在偽影現(xiàn)象。而本文所提的ATSVD、ADSVD、ADOptShrink這3種方法去噪結(jié)果的同相軸連續(xù)、清晰，且背景無明顯噪點。這是由于新方法可以通過自適應(yīng)秩選方法調(diào)整秩值大小，在每個計算子窗口內(nèi)數(shù)據(jù)驅(qū)動地獲得理想秩值，因此可以得到優(yōu)化后的去噪結(jié)果。從圖13所示的各種方法局部相似剖面也可以看到傳統(tǒng)的基于固定秩選的降秩去噪方法在局部有明顯的信號泄漏現(xiàn)象，相對的基于自適應(yīng)秩選方法的降秩去噪方法則具有較好的保幅性能。

觀察該雙曲模型在不同強度噪聲下各種方法去噪后的信噪比(圖14)、各種方法在不同秩選參數(shù)下對含噪50%的雙曲模型處理后的信噪比(圖15)，以及各種方法在不同秩選參數(shù)下對該模型的運算時間(圖16)，可以看出基于自適應(yīng)秩選方法的降秩去噪結(jié)果去噪性能相對較好，且不受秩選參數(shù)取值影響，也不會顯著增加時間復(fù)雜度，尤其是對于DOptShrink方法，結(jié)合自適應(yīng)秩選方法后，還會顯著提高運算效率。

2.3 實際資料測試

為驗證本文方法的實際處理效果，將本文的方法應(yīng)用于某工區(qū)陸上三維疊后實際地震記錄[41](圖17-A)，該數(shù)據(jù)大小為300×200×200，時間采樣間隔4 ms， 從圖中可見隨機噪聲分布較多，同相軸模糊、錯斷、不連續(xù)，不利于構(gòu)造解釋工作。將該實際數(shù)據(jù)分割為50×50×50的計算子窗口，選取固定秩選參數(shù)r為15，阻尼參數(shù)k為5進行處理。為更清晰地展示資料處理前后的地震反射差異細(xì)節(jié)，抽取Inline為10的單道剖面進行對比(圖17-B)。

圖18為傳統(tǒng)的基于固定秩選方法的去噪結(jié)果及其噪聲差剖面，可以看到傳統(tǒng)方法在箭頭所指區(qū)域丟失了大量有效波信號，且同相軸邊緣信息被處理得過于平滑，斷層構(gòu)造被涂抹。圖19為基于自適應(yīng)秩選方法的降秩去噪結(jié)果及其噪聲差剖面，圖中可見去噪結(jié)果視覺表現(xiàn)較好，不僅充分壓制了隨機噪聲，使得整個剖面同相軸細(xì)節(jié)、連續(xù)性、斷裂構(gòu)造走向表現(xiàn)清晰；從噪聲差剖面和圖20所示局部相似屬性剖面對比來看，有效信號也得到了充分保留，沒有明顯的同相軸損傷和丟失現(xiàn)象。

3 結(jié)論和建議

本文引入了一種基于奇異值最優(yōu)硬閾值的自適應(yīng)秩選方法替換秩約化方法中固定秩選參數(shù)，通過計算最優(yōu)硬閾值系數(shù)和估計數(shù)據(jù)噪聲水平得出潛在有效信號的秩值，從而克服了傳統(tǒng)的基于固定秩值的降秩去噪方法在實際應(yīng)用中的秩非一致性問題，最大限度地衰減隨機噪聲和保留有效信號。三維理論模型和實際數(shù)據(jù)試驗證明，與傳統(tǒng)的固定秩值降秩去噪方法相比，本文方法去噪后的信噪比相對更高，信號保持能力也相對最好，在計算效率上也有一定提升，因此建議將ADSVD方法和ADOptShrink方法推廣應(yīng)用于五維地震資料的自適應(yīng)同步去噪與插值處理中。

浙江大學(xué)地球科學(xué)學(xué)院陳陽康研究員提供了數(shù)據(jù)并進行了討論，作者借此表達感謝。