基于稀疏反演的多道匹配追蹤地震信號去噪方法及其應用

任 浩,李宗杰,薛 姣,顧漢明

(1.中國地質大學(武漢)地球內部多尺度成像湖北省重點實驗室,湖北武漢430074;2.中國石油化工股份有限公司西北油田分公司勘探開發研究院,新疆烏魯木齊830011)

地震數據不可避免地包含噪聲,噪聲會影響到地震屬性提取和地震解釋等后續工作,因此壓制噪聲是地震數據處理過程中的關鍵環節[1]。

增強有效信號或減弱噪聲通常包括兩個步驟:①定義一個可區分有效信號與噪聲的特征;②研發基于該定義的分離算法[2]。考慮到相鄰地震道中有效信號是相關的,噪聲是非相關的,基于該特征進行同相位疊加可增強有效信號,并減弱噪聲。地震信號中有效信號和噪聲在時間域是混疊的,不易分離,譜分解技術將時域地震信號映射到其它域后再進行處理,取得了良好的效果。目前主流的譜分解技術包括傅里葉變換、經驗模態分解(empirical mode decomposition,EMD)、同步壓縮小波變換(synchrosqueezing wavelet transform,SWT)、小波閾值法(wavelet thresholding approach,WTA)、匹配追蹤算法(matching pursuit,MP)和反演譜分解(inverse spectral decomposition,ISD)等。傅里葉變換是將地震信號轉換到頻率域,根據噪聲頻率與有效信號頻率范圍的差異特征,設置合適的濾波器壓制噪聲,但這種方法對非平穩地震信號的處理效果不佳[3]。經驗模態分解[4]是將地震信號表示成多階固有模態函數分量之和,低階分量頻率高,高階分量頻率低,去除包含高頻噪聲的低階分量即可實現去噪,但該方法造成了部分有效信號損失,針對這一問題,相繼有學者提出了改進的方法[1,5-8]。同步壓縮小波變換[9]具有高時頻分辨率特性和信號重構能力,該方法根據有效信號的時頻分布范圍,將時頻譜中相應時頻成分重構有效信號以達到去除噪聲的目的[10-12]。小波閾值法[13]是將地震信號變換到時間-尺度域,利用噪聲能量相對較弱的特點,將小于閾值的系數設置為零,由新的小波系數經小波反變換得到去除噪聲后的信號。匹配追蹤算法[14]是將地震信號分解為一系列子波的線性組合,趙天姿等[15]人工挑選出與有效信號匹配的子波重構信號,后來有學者依據多個連續地震道中有效信號的橫向相關性,提出了多道匹配追蹤地震信號分解方法[16-17]。反演譜分解[18]與匹配追蹤算法類似,也是一種基于子波庫的分解方法,該方法預先設定誤差水平值,將地震信號表示成稀疏時頻子波的線性相加,將殘差信號當作噪聲[19-20]。

匹配追蹤分解是采用逐步搜尋子波庫中與殘差信號最佳匹配子波的方式分解地震信號,并非直接且準確地去除噪聲。搜尋最佳匹配子波是分解中最為關鍵的步驟,本文方法以子波庫為著手點來解決這一問題,采用與地震信號中有效信號匹配的子波構建的子波庫代替超完備子波庫,將地震信號表示為有效信號匹配子波的線性疊加與殘差信號之和,將殘差信號當作噪聲,即達到直接去除噪聲的目的。

1 基于稀疏反演的多道匹配追蹤分解原理

如何實現直接分離有效信號和噪聲？首先,以連續多道的中間道為標準,校正同相軸為水平后同相疊加求取疊加道;接著,稀疏反演分解疊加道,選取疊加道中有效信號的匹配子波作為連續多道的中間道匹配追蹤分解的子波庫;最后,根據過中間道上每一點的同相軸傾角大小,反校正中間道子波庫的子波中心時間,分別構建各地震道的子波庫,約束匹配追蹤分解的掃描范圍,選擇性地分離有效信號。

1.1 連續多道同相軸校正

疊后數據剖面中局部范圍內的同相軸可認為是線性的,但其方向通常并非水平,如果直接疊加會造成波形畸變,因此在疊加前需進行水平校正。HUO等[21]提出了計算連續地震道同相軸傾角的方法。該方法根據傾角的大小,以中間道為基準,校正同相軸為水平后疊加地震道,疊加道可保持連續地震道的波形特征,具體如下。

選取2W1+1個連續地震道作為一組,在中間道上坐標為(i,W1+1)各點處定義長度為2W1+1,寬度為2W2+1、傾角為p的窗口。傾角改變對應著右側采樣點的變化,p=0時窗口水平,p>0時窗口相對p=0為逆時針旋轉,反之,p<0時窗口相對p=0為順時針旋轉。設定p的取值范圍為[pmin,pmax],間隔為一個采樣點,依次計算不同傾角時窗口中連續地震道與中間道的距離,計算距離的函數定義如下。

式中:pi,k是過點(i,W1+1)的第k個傾角,當k=m時距離最小,可認為pi,m是經過該點的同相軸傾角;Xj(pi,k)表示在點(i,W1+1)處傾角為pi,k的窗口中第j道的數據向量,向量長度為2W2+1。如圖1所示,中間道第57個采樣點處的傾角為p57,m。采用該方法計算過中間道上每一個采樣點的同相軸傾角,再依據傾角大小校正同相軸,得到同相軸水平的地震道集。如圖2所示,圖2a為合成的連續地震道;經校正后得到圖2b;將圖2a和圖2b中的連續地震道疊加求平均,分別得到圖2c和圖2d中的藍色虛線。對比可知:校正后的疊加道保持了波形特征,直接疊加地震道則會造成波形畸變。

圖1 計算連續地震道同相軸傾角

圖2 合成的連續地震道同相軸校正a 合成的連續地震道; b 校正后的連續地震道; c 中間道(紅色實線)和未經校正的疊加道(藍色虛線)對比; d 中間道(紅色實線)和經校正后的疊加道(藍色虛線)對比

1.2 稀疏反演分解提取匹配子波

地震褶積模型可以描述為:一個地震信號s(t)由地震子波w(t)和反射系數r(t)的褶積加上隨機噪聲[22]:

(3)

式中:“*”表示褶積運算。(3)式無法表示信號頻率成分隨時間的變化,考慮到地震信號包含有多種頻率成分,可認為地震信號由不同頻率的子波經地下界面反射疊加而成,(3)式可表示為:

(4)

式中:wi(t)為主頻為下角標對應頻率的子波;ri(t)是與wi(t)對應的頻率依賴的偽反射系數;N為地震信號的頻率最大值。將(4)式寫成矩陣與向量乘積的形式,則地震褶積模型如下。

(5)

式中:SM×1是長度為M的地震信號向量;Wi(i=1,2,…,M×N)是長度為M的時頻子波向量;Ri(i=1,2,…,M×N)為與之對應的系數;DM×(M×N)是M行、M×N列的矩陣,即超完備子波庫;m(M×N)×1表示時頻子波的系數向量;nM×1是噪聲向量。根據地震信號傅里葉變換得到的振幅譜確定信號的頻率范圍,在該范圍內以1Hz為間隔進行采樣作為時頻子波的主頻,將這些子波沿時間軸平移建立超完備子波庫。求解(5)式可得系數向量m(M×N)×1,將其重排成M×N的矩陣,即為地震信號的時頻主頻譜。時頻主頻譜平面上的每一個坐標點均對應著一個時頻子波,所有坐標點的集合對應著超完備子波庫。不考慮噪聲的情況下(5)式可表示為(略去各個量的下角標):

(6)

(6)式是一個欠定方程,有無數個解,為了減小多解性,引入稀疏約束,即用盡量少的時頻子波表示地震信號,從而得到高分辨率的時頻主頻譜。向量的稀疏度常用L0范數表示,但含L0范數的方程不易求解,將其轉化為L1范數方程:

(7)

(7)式稱為基追蹤(basis pursuit,BP)問題。當噪聲存在時,基追蹤問題轉變為基追蹤降噪(basis pursuit de-noising,BPDN)問題:

(8)

式中:正參數σ是對數據中噪聲水平的估計。當σ=0時,即為BP問題,因此BP問題是BPDN問題的特殊情況。本文運用L1范數最小化的譜投影梯度(spectral projected-gradient for L1 minimization,SPGL1)算法求解(8)式,該算法的基本思想是將BPDN問題轉化為最小絕對值收斂和選擇算法(least absolute shrinkage and selection operator,LASSO)問題,通過迭代求得符合BPDN條件的解。LASSO問題表示為:

(9)

式中:τ為m的L1范數約束。地震信號稀疏反演分解得到的時頻主頻譜平面上有能量的坐標點對應著構成地震信號的子波,稱為匹配子波,因此通過稀疏反演可從超完備子波庫中提取與地震信號匹配的子波。

1.3 基于稀疏反演的多道匹配追蹤分解

地震信號匹配追蹤分解受到噪聲的強烈影響[16],結果具有非唯一性,本文提出的基于稀疏反演的多道匹配追蹤算法能夠減小噪聲干擾,直接分離有效信號,具體實現步驟如下。

選出2W1+1個連續地震道{S1,…,SW1+1,…,S2W1+1}作為一組,經校正后疊加表示為:

(10)

(11)

(12)

(13)

因此,每次迭代更新后的殘差信號為:

(14)

(15)

2 模型試算

2.1 含噪單道信號試算

圖3a為合成的單道信號,由不同主頻、時移和相位的Ricker子波組成,其子波主頻從上到下依次為60,40,20,35,35,18,56Hz,加入高斯白噪聲后如圖3b所示。圖3c為有效信號的理論時頻主頻譜,圖中每一個具有能量的坐標點均代表著一個時頻子波,表明信號中含有該時頻子波成分,能量越強,與信號越匹配。建立超完備復Ricker子波庫,設定σ值,圖3d為合成的含噪信號稀疏反演分解得到的時頻主頻譜,與圖3c對比發現,雖然能量分布范圍有所擴大,但準確地定位出了時頻主頻譜的實際分布位置。稀疏反演分解方法具有較強的抗噪能力,能夠精確地從超完備子波庫中選取出與有效信號匹配的子波。

圖3 合成的含噪單道信號時頻主頻譜a 合成的單道信號; b 加噪單道信號; c 有效信號的理論時頻主頻譜; d 合成的含噪信號稀疏反演分解的時頻主頻譜

圖4a是根據圖3b所示的信號單道匹配追蹤分解得到的重構信號(紅色虛線),由于分解過程中無法區分有效信號(綠色實線)和噪聲,重構信號時引入了部分噪聲。將理論時頻主頻譜(圖3c)中有能量的坐標點對應的時頻子波(有效信號的匹配子波)設定為子波庫,采用與單道匹配追蹤分解相同的迭代次數,得到圖4b所示的重構信號(紅色虛線),可以發現重構信號和有效信號基本重合,說明將地震信號中有效信號的最佳匹配子波作為子波庫,可分離有效信號和噪聲,準確重構有效信號。稀疏反演方法能較精確地選出有效信號的匹配子波,由此與匹配追蹤分解結合也可實現有效信號的準確重構。

2.2 含噪連續多道信號試算

圖5a為合成的地震數據剖面,包括兩個斜線型和1個拋物線型同相軸;加入噪聲后的剖面如圖5b所示,一共有50道,分成10組,每組5道;水平校正每1組地震道后的剖面如圖5c所示;圖5d為所有組的疊加道,疊加道信噪比相較于中間道信噪比更高。第1組中間道的信噪比RSN=0.77,疊加道信噪比RSN=1.67(RSN=Jsignal/Jnoise,其中Jsignal為有效信號能量,Jnoise為噪聲能量),高信噪比有利于稀疏反演提取出更準確匹配子波。每道有500個時間采樣點,間隔為1ms,經傅里葉變換分析,所有道的頻率范圍為0～100Hz,以1Hz為主頻率間隔采樣,即有100個主頻采樣點,可建立矩陣大小為500×50000的超完備復Ricker子波庫。圖6a和圖6b分別為采用本文方法去噪后的疊加剖面和重構誤差剖面;圖6c和圖6d 分別為采用單道匹配追蹤方法去噪后的疊加剖面和重構誤差剖面,迭代次數均為20次。對比發現,兩種方法均可有效壓制數據中噪聲,但采用本文方法得到的重構結果與原地震數據誤差小,去噪效果佳。

圖4 有效信號(綠色實線)的重構a 單道匹配追蹤分解的重構信號(紅色虛線); b 有效信號匹配子波作為子波庫的匹配追蹤分解的重構信號(紅色虛線)

圖5 合成的地震數據剖面及同相軸校正a 合成的地震數據剖面; b 加入噪聲后的地震數據剖面; c 同相軸校正后剖面; d 所有組的疊加道(共10組)

圖6 不同方法的去噪效果a 本文方法去噪后的疊加剖面; b 本文方法去噪后的重構誤差剖面; c 單道匹配追蹤分解去噪后的疊加剖面; d 單道匹配追蹤分解的重構誤差剖面

3 實際數據處理

圖7a為某探區實際疊后數據剖面,共有50個地震道,每道400個采樣點,采樣間隔為4ms。采用本文方法對實際數據進行處理:設置7個連續地震道為1組,為了提高去噪后剖面的橫向連續性,相鄰兩組共享3個地震道,除最后1組對該組內全部道進行分解外,其它組只分解前4道。采用單道匹配追蹤方法和本文方法分解每一個地震道的迭代次數均為100次,圖7b和圖7c分別為采用本文方法去噪后得到的疊加剖面和殘差剖面,處理結果表明,采用本文方法可較好地壓制隨機噪聲,同時保持地震數據的橫向連續性。這是由于子波庫中的子波與地震道中有效信號是匹配的,在匹配追蹤分解的過程中可有選擇性地提取有效信號成分,避免了將噪聲引入重構信號,同時,疊加道很好地保留了連續地震道中橫向連續的有效信號成分,利用疊加道稀疏反演分解提取的子波建立子波庫,匹配追蹤分解識別了橫向連續的有效信號,因此重構信號具有清晰的同相軸。圖7d和圖7e分別為采用單道匹配追蹤方法去噪后得到的疊加剖面和殘差剖面,相較于采用本文方法得到的處理結果,該方法得到結果中仍有部分噪聲未去除,且同相軸不清晰。

圖7 實際地震數據去噪效果a 實際疊后數據剖面; b 采用本文方法去噪后得到的疊加剖面; c 采用本文方法去噪后得到的殘差剖面; d 單道匹配追蹤分解去噪后得到的疊加剖面; e 單道匹配追蹤分解去噪后得到的殘差剖面

4 結論

本文提出的基于稀疏反演的多道匹配追蹤去噪方法,利用疊加道的高信噪比、稀疏反演分解地震信號的抗噪性和解的稀疏表示特征,提取與有效信號匹配的子波作為子波庫,在匹配追蹤分解過程中分離有效信號和噪聲,實現了信號的精確重構。擬合算例和實際地震數據處理結果表明,本文方法相較于單道匹配追蹤方法具有以下優點:①去除噪聲的效果佳;②去噪后數據剖面上同相軸橫向連續性好;③合理地增加每一組中地震道數,運算效率高。但本文方法也存在著不足之處,當地震剖面中存在斷層或者較多跳點時,連續地震道同相軸校正不準確,造成疊加道失真,構建的子波庫中子波與有效信號不匹配,導致匹配追蹤分解無法分離噪聲和有效信號,這些問題有待于今后進一步的研究。