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

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

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

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

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

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

1.1 Hankel矩陣降秩去噪

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

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

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

(1)

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

(2)

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

(3)

(4)

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

(11)

(12)

(13)

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

(14)

(15)

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

(16)

1.3 自適應秩選方法

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

H=X+σZ

(17)

(18)

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

(19)

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

(20)

(21)

(22)

(23)

(24)

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

(25)

式中:xmed是觀測矩陣H的奇異值序列的中值;μβ是β基于隨機矩陣理論的Marcenko-Pastur分布定律的中值解，即求解

(26)

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

(27)

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

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

(28)

1.4 自適應降秩去噪方法

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

自適應降秩去噪方法的去噪流程如下：

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

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

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

2 模型與算例

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

(29)

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

(30)

(31)

(32)

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

(33)

(34)

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

2.1 線性模型測試

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

表1給出了6種方法處理后的信噪比，也可以看出本文改進方法相較于原方法在信噪比上均有所提高，其中ADSVD方法和ADOptShrink方法信噪比分別提高到 13.211 5 dB和 13.304 9 dB。從圖8各種方法的局部相似屬性剖面可以看出，全局均接近于零值，指示各種方法在處理由線性同相軸構成的地震數(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嚴格取為同相軸的傾角數(shù)時，才能取得最高信噪比；此后隨著r的增大，更多的噪聲成分會被引入去噪結果，導致信噪比下降。但相對而言DOptShrink方法下降趨勢更緩，因為該方法對秩選參數(shù)并不敏感[28]，而本文方法結果則完全不受固定值r影響。圖11為在ubuntu環(huán)境下利用Matlab2016b計算平臺，配置為Intel Core i5-8400CPU@2.80 GHz、4GB RAM時，統(tǒng)計的各種方法在不同秩選參數(shù)下的運算時間。由于降秩去噪方法的時間復雜度集中在奇異值分解步驟上，因此引入自適應秩選方法并不會顯著增加運算時間。另一方面，可以看到DOptShrink方法結果雖然對秩選參數(shù)r不敏感，但其運算時間會隨r的增大而顯著增加，因此為獲得較好的去噪效果而對其選擇一個較大的秩值并不理想。

2.2 雙曲模型測試

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

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

2.3 實際資料測試

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

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

3 結論和建議

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

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