地震數(shù)據(jù)譜反演壓縮感知算法實現(xiàn)及應(yīng)用

夏紅敏 劉蘭鋒 張顯輝 陳雙全*

(①中國石化石油勘探開發(fā)研究院，北京100083；②中國石油大學(xué)(北京)油氣資源與探測國家重點實驗室，北京102249；③中國石油大學(xué)(北京)CNPC物探重點實驗室，北京102249)

0 引言

為了解決油氣勘探中薄儲層的識別問題，地球物理學(xué)家們嘗試了許多提高地震資料分辨率的處理方法。Portniaguine等[1-2]和Chopra等[3-5]最早基于地層反射系數(shù)奇、偶分解原理提出了譜反演方法；Puryear等[6]將反射系數(shù)對分解成奇、偶兩個分量，根據(jù)奇、偶分量對薄層的分辨能力不同，首先構(gòu)建了地震數(shù)據(jù)譜反演方法模型，合成數(shù)據(jù)測試和實際地震數(shù)據(jù)反演驗證了方法的正確性，為提高地震數(shù)據(jù)分辨率提供了一種新方法。在隨后的幾年，譜反演方法得到了快速的發(fā)展和完善[7-8]。

Zhang等[9]介紹了一種反射系數(shù)基追蹤反演方法，該方法通過建立反映反射模式的函數(shù)字典，將地震道數(shù)據(jù)看成是這些函數(shù)的疊加，并在反演過程中引入先驗信息進(jìn)行反演約束。柴新濤等[10]提出基于最小二乘奇異值分解(LSQR)算法與模擬退火算法相結(jié)合的方法，完成了稀疏譜與非稀疏譜兩種假設(shè)條件下的反演，同時分析了譜反演結(jié)果的影響因素。Yuan等[11]將反射系數(shù)反演與稀疏貝葉斯學(xué)習(xí)相結(jié)合實現(xiàn)了譜反演，通過添加、刪除或重新估計等方式檢索稀疏反射系數(shù)序列，無需預(yù)先設(shè)置非零反射系數(shù)的個數(shù)，并在合成數(shù)據(jù)、物理模型數(shù)據(jù)和實際數(shù)據(jù)中驗證了該方法，結(jié)果表明譜反演方法能較好地識別地下薄層。劉萬金等[12]利用譜反演法有效地提高了地震數(shù)據(jù)的分辨率，并將其應(yīng)用于地震屬性分析。田立新等[13]基于貝葉斯理論框架，將測井等信息轉(zhuǎn)換為地質(zhì)統(tǒng)計先驗信息，引入到譜反演過程，提出地質(zhì)統(tǒng)計學(xué)譜反演法。張華等[14]及劉潔等[15]在實際地震資料處理中應(yīng)用譜反演，均取得不錯的效果。張軍華等[16]實現(xiàn)了基于Moore-Penrose算法的譜反演；陳祖慶等[17]實現(xiàn)了基于壓縮感知的譜反演算法。

地震數(shù)據(jù)譜反演的假設(shè)條件是反射系數(shù)為稀疏序列，而稀疏性也是信號分析中壓縮感知算法的條件之一。因此，將壓縮感知算法引入地震數(shù)據(jù)譜反演，可以很好地實現(xiàn)地震數(shù)據(jù)譜反演方法的求解。另外，針對實際地震數(shù)據(jù)反演過程中地震子波譜難以估算的問題，唐文博等[18]提出利用地震數(shù)據(jù)二次譜估算子波譜的方法，直接將頻率域的地震子波譜用于譜反演。

本文提出聯(lián)合地震數(shù)據(jù)二次譜計算方法及基于壓縮感知算法的譜反演方法，試算結(jié)果表明，可以很好地解決實際數(shù)據(jù)譜反演中的穩(wěn)定性問題。

1 方法原理

1.1 譜反演方法

在地震數(shù)據(jù)時間域，利用脈沖函數(shù)的性質(zhì)，將時間間隔為T的地層頂、底反射系數(shù)對(r1，r2)表示為[7]

g(t)=r1δ(t-t1)+r2δ(t-t1-T)

(1)

式中：r1、r2分別是地層頂、底界面的反射系數(shù)；t為時間序列采樣；t1、t2分別是頂、底界面的反射時間，T=t2-t1；δ(·)為脈沖函數(shù)。如果零時間設(shè)置到地層中點，則有

(2)

對式(2)進(jìn)行傅里葉變換得到

(3)

根據(jù)反射系數(shù)對的奇偶分解原理[6]可得

g(t,f)=2recos(πfT)+j2rosin(πfT)

(4)

其中re、ro為地層反射系數(shù)對(r1，r2)的偶分量和奇分量，且2re=r1+r2，2ro=r1-r2。

頻率域反射系數(shù)對的實部Re[·]和虛部Im[·]為

Re[g(t,f)]=2recos(πfT)

(5)

Im[g(t,f)]=2rosin(πfT)

(6)

則單一層的頂、底反射系數(shù)對 Δt=t1+T/2的時延譜為

米多想到自己莫名丟失的錄音筆，正想問鮑澤，對方卻神色一變，突然擺出一副不耐煩的樣子。隨即，米多的耳旁傳來了一個熟悉的聲音。

Re[e-j2πfΔtg(t,f)]=2recos(πfT)cos(2πfΔt)+

2rosin(πfT)sin(2πfΔt)

(7)

Im[e-j2πfΔtg(t,f)]=2rosin(πfT)cos(2πfΔt)-

2recos(πfT)sin(2πfΔt)

(8)

根據(jù)以上推導(dǎo)可得到N個稀疏反射系數(shù)序列r(ti)，i=1，2，…，N，其頻率域反射系數(shù)的奇、偶分量形式如下

2ro(i，i+1)sin(πfTi)sin(2πfΔti)]

(9)

2re(i，i+1)cos(πfTi)sin(2πfΔti)]

(10)

式中：2re(i，i+1)=ri+ri+1；2ro(i，i+1)=ri-ri+1。

式(9)和式(10)表明，可以利用地震記錄和地震子波的部分譜信息進(jìn)行反演進(jìn)而得到反射系數(shù)序列，該方法即地震數(shù)據(jù)譜反演[7]。因此，目標(biāo)函數(shù)式為

(11)

譜反演目標(biāo)函數(shù)求解過程是利用部分地震信息求解整個時間域內(nèi)稀疏反射系數(shù)的過程，壓縮感知的基本思想是在滿足有效等距條件(RIP)下，利用少量采樣信號重構(gòu)出完整的原信號。二者思路不謀而合，因此本文選用壓縮感知算法解決譜反演問題。假設(shè)反射系數(shù)為稀疏信號，可將譜反演目標(biāo)函數(shù)重寫為最優(yōu)化目標(biāo)函數(shù)公式

(12)

式中：r是待求解的反射系數(shù)；b是地震記錄頻譜與地震子波頻譜之比；A是構(gòu)建的測量矩陣；τ表示反射系數(shù)的最大幅值，一般取0.2。

可將式(12)的約束優(yōu)化條件改寫為無約束優(yōu)化條件

(13)

其中λ是用于平衡二次項和稀疏項結(jié)果權(quán)重的正則化參數(shù)。上述最優(yōu)化過程為帶邊界約束的二次規(guī)劃問題，正則化參數(shù)用于控制稀疏解。為了將非凸優(yōu)化問題轉(zhuǎn)化為凸優(yōu)化問題，采用梯度投影稀疏重構(gòu)(Gradient Projection for Sparse Reconstruction， GPSR)算法[19]進(jìn)行求解。

1.2 地震子波譜的求取

式(11)表明，地震數(shù)據(jù)譜反演目標(biāo)函數(shù)與地震子波譜相關(guān)。因此，地震子波譜估算的準(zhǔn)確度直接影響譜反演結(jié)果的準(zhǔn)確性。本文通過地震數(shù)據(jù)二次譜求取地震子波譜，具體方法如下。

首先，根據(jù)地震數(shù)據(jù)褶積模型，得到頻率域的地震數(shù)據(jù)振幅譜關(guān)系式

AS(f)=AW(f)Ar(f)

(14)

式中AS(f)、AW(f)和Ar(f)分別為地震記錄、地震子波和反射系數(shù)的振幅譜。

然后，設(shè)法提取地震子波振幅譜。在頻率域，地震信號中同時包含反射系數(shù)與地震子波的頻譜信息，無法將其有效分離。而在譜模擬假設(shè)條件下，地震子波的振幅譜主要為地震記錄頻譜中光滑的低頻分量，因此，對式(14)兩邊同時進(jìn)行傅里葉變換，可得到地震數(shù)據(jù)的二次譜，即地震子波二次譜與反射系數(shù)二次譜的褶積，表達(dá)式為

(15)

2 合成數(shù)據(jù)測試

進(jìn)行合成數(shù)據(jù)測試時，首先利用合成地震道數(shù)據(jù)驗證由地震記錄二次譜求取地震子波譜方法的準(zhǔn)確性。圖1顯示了30Hz Ricker子波和三種不同信噪比的合成地震記錄的二次譜曲線，可以看出：在0～45Hz頻率范圍內(nèi)，合成地震記錄的二次譜可以很好地模擬地震子波的振幅譜。圖2顯示合成地震記錄振幅譜(黑線)、利用二次譜得到的子波譜(紅線)和已知的地震子波譜(藍(lán)線)。對比表明，利用地震記錄二次譜可以計算得到最吻合的地震子波頻譜，為實際地震數(shù)據(jù)譜反演中的地震子波譜估算提供了實用且有效的方法。

圖1 地震子波與合成地震記錄二次譜對比

利用地震記錄二次譜估算得到較準(zhǔn)確的地震子波譜后，通過GPSR算法對合成地震記錄進(jìn)行譜反演試算。圖3分別展示了30Hz地震Ricker子波(圖3a)、無噪聲的合成地震記錄(圖3b)以及信噪比SNR分別為10、5、2的合成地震記錄(圖3c～圖3e)。譜反演結(jié)果如圖4所示。對比表明，對于信噪比較高的地震數(shù)據(jù)，通過譜反演可以得到準(zhǔn)確的反射系數(shù)。但是當(dāng)SNR < 5時，反演結(jié)果受信噪比影響較大，低信噪比(如SNR=2)地震數(shù)據(jù)反演結(jié)果與實際反射系數(shù)的偏差較大。因此，譜反演輸入地震數(shù)據(jù)的信噪比不能太低。

圖2 合成地震記錄振幅譜(黑線)、地震子波振幅譜(藍(lán)線)和計算得到的地震子波振幅譜(紅線)對比(a)無噪聲； (b)SNR=10； (c)SNR=5

圖3 地震子波及合成地震記錄對比(a)Ricker子波； (b)無噪； (c)SNR=10； (d)SNR=5； (e)SNR=2

圖4 原始反射系數(shù)與譜反演結(jié)果對比(a)原始反射系數(shù)序列； (b)無噪； (c)SNR=10； (d)SNR=5； (e)SNR=2

3 實際數(shù)據(jù)應(yīng)用

地震數(shù)據(jù)譜反演最終得到的是反射系數(shù)序列，之后，可以進(jìn)一步反演波阻抗，也可以將反射系數(shù)與寬頻地震子波合成高分辨率地震記錄。本文利用基于壓縮感知算法的譜反演結(jié)果進(jìn)行實際地震數(shù)據(jù)提高分辨率處理的應(yīng)用，以驗證該反演算法的穩(wěn)定性和有效性。處理流程包括：

(1)對地震數(shù)據(jù)進(jìn)行頻譜分析，確定反演的頻率范圍；

(2)計算地震記錄的二次譜，進(jìn)行低通濾波處理；

(3)求取地震子波譜；

(4)利用式(11)求解反射系數(shù)序列；

(5)利用寬頻帶子波合成地震數(shù)據(jù)，得到提高分辨率后的處理結(jié)果。

所選實際地震工區(qū)的目的層為奧陶系風(fēng)化殼氣藏儲層，由于受上覆煤層影響，目的層的地震資料分辨率較低，后續(xù)進(jìn)行高精度儲層預(yù)測難度很大。圖5為工區(qū)內(nèi)一條地震剖面，圖中的強反射同相軸為煤層反射信號。

圖5 實際地震數(shù)據(jù)剖面

根據(jù)上述處理流程，首先通過地震數(shù)據(jù)頻譜分析選取合適的反演頻率范圍，結(jié)合反演測試，最終選取的頻率范圍為10～70Hz。其次，利用地震資料的二次譜進(jìn)行地震子波譜估計。針對該實際數(shù)據(jù)，采用多道平均譜估算地震子波譜。圖6為地震數(shù)據(jù)的二次譜(藍(lán)色)與低通濾波后的二次譜(紅色)對比，圖7為計算得到的地震子波譜。利用本文提出的譜反演方法對地震數(shù)據(jù)進(jìn)行反演，并將反演得到的反射系數(shù)序列與寬頻子波合成，輸出提高分辨率處理后的結(jié)果(圖8a)。

為了驗證本文方法針對實際地震數(shù)據(jù)的處理效果，同時對地震數(shù)據(jù)進(jìn)行了譜白化處理，結(jié)果如圖8b所示。對比兩種方法處理后的地震剖面可以看到：本文方法處理后的地震數(shù)據(jù)信噪比更高，剖面信息更加豐富，視分辨率得到明顯提高；而利用譜白化處理得到的地震剖面，分辨率提高不明顯，并且信噪比偏低。

圖6 實際地震數(shù)據(jù)二次譜(藍(lán)色)與濾波后的二次譜(紅色)對比

圖7 計算得到的地震子波振幅譜

圖8 實際數(shù)據(jù)反演提高分辨率處理(a)及譜白化處理(b)結(jié)果對比

對原始地震數(shù)據(jù)及處理后數(shù)據(jù)進(jìn)行分頻掃描，可更好地分析不同方法處理效果。圖9對比了原始數(shù)據(jù)、譜白化處理及本文方法處理結(jié)果的分頻掃描結(jié)果，掃描頻率分別為10～40Hz、40～80Hz和80～120Hz?？梢钥吹?，在10～40Hz掃描頻段，譜白化處理(圖9d)及本文方法處理(圖9g)結(jié)果均能很好地保持原始數(shù)據(jù)的信息；40～80Hz掃描時譜白化處理結(jié)果信噪比明顯降低(圖9e)；在80～120Hz的掃描結(jié)果中，原始數(shù)據(jù)中80Hz以上的信息不足(圖9c)，譜白化處理結(jié)果的信噪比更低(圖9f)，只有本文方法處理結(jié)果恢復(fù)得到了高頻段的有效信號(圖9i)。三個頻段的掃描結(jié)果均表明，本文方法處理后的地震剖面的信噪比最高。

圖9 數(shù)據(jù)分頻掃描剖面對比(a～c)原始數(shù)據(jù)； (d～f)譜白化處理結(jié)果； (g～i)本文方法處理結(jié)果

4 結(jié)論

本文結(jié)合地震記錄的二次譜估算地震子波譜的方法，提出了一種基于壓縮感知算法的地震數(shù)據(jù)譜反演方法，形成了利用譜反演結(jié)果提高地震資料分辨率的數(shù)據(jù)處理流程。實際地震資料應(yīng)用結(jié)果表明，該處理流程計算步驟和參數(shù)選取均較簡單，在實際數(shù)據(jù)處理中可以很好地恢復(fù)高頻段有效信號，而且經(jīng)過處理后的地震數(shù)據(jù)信噪比高、橫向連續(xù)性好。合成數(shù)據(jù)測試表明，該方法受原始地震數(shù)據(jù)信噪比的影響較大，因此對于低信噪比資料的數(shù)據(jù)反演還需要進(jìn)一步完善。