許李囡 高靜懷* 楊 陽 高照奇 王 前

(①西安交通大學(xué)信息與通信工程學(xué)院，陜西西安 710049； ②湖北文理學(xué)院數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，湖北襄陽 441053)

0 引言

地震波在地層傳播過程中，會受到介質(zhì)的吸收和散射效應(yīng)，導(dǎo)致地震波的主頻降低、頻寬變窄，從而降低中、深層地震資料分辨率，給油氣勘探帶來諸多困難[1]。地震波在地下介質(zhì)中傳播時產(chǎn)生的衰減效應(yīng)，可以用品質(zhì)因子Q定量表示[2]。實(shí)驗(yàn)室測量表明，地層對地震波的吸收作用主要取決于巖石物性、孔隙度、孔隙流體成分及飽和度等[3],因此，Q值估計(jì)是預(yù)測和直接尋找油氣儲層的有效途徑之一[4]。準(zhǔn)確的Q模型還可用于反Q濾波[5-6]以提高地震資料的分辨率，為后續(xù)精細(xì)解釋和準(zhǔn)確成像提供高分辨率的地震資料。

近年來，學(xué)者們已經(jīng)提出了許多種Q估計(jì)方法。常見的Q估計(jì)方法一般可分為時間域、頻率域和時頻域三大類[7]。時間域方法包括脈沖上升時間法、子波模擬法等； 頻率域方法包括對數(shù)譜比法(Logarithmic Spectral Ratio， LSR)[8]、中心頻率偏移法(Center Frequency Shift， CFS)[9]、峰值頻率偏移法(Peak Frequency Shift, PFS)[9]等； 時頻域方法包括子波包絡(luò)峰值頻率偏移法[10]、基于Gabor譜的常Q估計(jì)方法(Constant-QEstimation via Gabor Analysis, CVG)[11]等。其中，CVG方法是一種基于Gabor變換的Q值估計(jì)方法，它不需對地震子波進(jìn)行假設(shè)，已廣泛應(yīng)用于提取地震數(shù)據(jù)的品質(zhì)因子Q。但是，由于Gabor變換的窗函數(shù)固定，該方法對窗函數(shù)的參數(shù)比較敏感； 同時，CVG方法對含噪地震數(shù)據(jù)的估計(jì)精度并不理想。

為了克服上述缺點(diǎn)，本文提出了一種新的基于S變換和變分法的Q估計(jì)方法(QEstimation Based on the Calculus of Variations and S Transform, CVS)。與Gabor變換相比，S變換[12]具有多尺度分辨率和自適應(yīng)調(diào)整窗長度等優(yōu)點(diǎn)，可有效避免CVG方法對窗函數(shù)比較敏感的問題。首先，推導(dǎo)了非平穩(wěn)地震數(shù)據(jù)在S變換域中的近似表示。然后，根據(jù)該近似表示構(gòu)建品質(zhì)因子Q和子波之間的優(yōu)化函數(shù)，利用微分和變分法對該優(yōu)化函數(shù)進(jìn)行求解，進(jìn)而得到品質(zhì)因子Q的計(jì)算表達(dá)式； 與CVG方法相比，提出的Q估計(jì)方法無需假設(shè)地震子波的類型，也不必人為調(diào)節(jié)窗函數(shù)的長度。最后，為提高CVS方法的準(zhǔn)確性和抗噪性，設(shè)計(jì)了一種自適應(yīng)參數(shù)選擇的方案，對計(jì)算過程中的積分區(qū)間進(jìn)行自動選取。

綜上所述，CVS方法有效解決了CVG方法中存在的需要固定窗長和抗噪性較差的問題，顯著提高了Q值計(jì)算的準(zhǔn)確性和穩(wěn)定性。模型數(shù)據(jù)和實(shí)際數(shù)據(jù)算例均驗(yàn)證了該方法的有效性。

1 基本原理

1.1 非平穩(wěn)地震數(shù)據(jù)的S域近似表示

常規(guī)的地震資料高分辨率處理一般建立在傳統(tǒng)的平穩(wěn)褶積模型[13]上?？紤]到地震波在傳播過程中存在衰減效應(yīng)，需在平穩(wěn)褶積模型中引入衰減，得到更符合實(shí)際地下情況的非平穩(wěn)褶積模型[14]。令x(t)為非平穩(wěn)反射地震數(shù)據(jù)，其頻域形式在非平穩(wěn)褶積模型中可表示為

(1)

(2)

式中fR表示參考頻率。

對式(1)進(jìn)行傅里葉逆變換，x(t)可表示為

(3)

式中u表示與時間和頻率相關(guān)的衰減模型中的時間變量。

由于S變換的窗函數(shù)隨著頻率變化，所以低頻與高頻的時頻分辨率具有差異，這有利于刻畫地震數(shù)據(jù)的非平穩(wěn)性。引入S變換分析非平穩(wěn)地震數(shù)據(jù)。地震數(shù)據(jù)x(t)的S變換[7]定義為

(4)

式中：τ和ν分別是時間變量和頻率變量；g是高斯窗函數(shù)。由式(4)能看出，與Gabor變換的固定高斯窗函數(shù)相比，S變換的窗函數(shù)可隨頻率進(jìn)行調(diào)節(jié)，一定程度上解決了Gabor變換中固定窗的問題。

為更清楚說明地震信號的S變換多尺度分辨率的優(yōu)勢，選擇主頻為40Hz的Ricker子波分別進(jìn)行Gabor變換和S變換(圖1)。由圖可見，S變換中高頻分量具有更好的時間分辨率，低頻分量具有更好的頻率分辨率，因此S變換比Gabor變換能夠更好地描述非平穩(wěn)地震數(shù)據(jù)。

圖1 地震數(shù)據(jù)Gabor變換結(jié)果(a)和S變換結(jié)果(b)

將式(3)代入式(4)，式(4)可表示為

(5)

(6)

因此，非平穩(wěn)地震數(shù)據(jù)的S域可近似表示為

(7)

式中Sr(τ,ν)表示反射系數(shù)的S變換結(jié)果。

1.2 品質(zhì)因子Q的計(jì)算公式

在近似表示的基礎(chǔ)上，僅考慮地震記錄的振幅譜，式(7)可表示為

(8)

假設(shè)反射系數(shù)序列的譜是白的，在非平穩(wěn)地震數(shù)據(jù)變換到S域后，式(8)可以忽略|Sr(τ,ν)|項(xiàng)，進(jìn)一步簡化并兩邊取對數(shù)可得

(9)

式中Q是待求的參數(shù)。一個魯棒的思路是優(yōu)化關(guān)于Q的目標(biāo)函數(shù)，從而得到一個穩(wěn)定解。通過式(9)，能夠在最小二乘的意義下，基于S變換建立如下的目標(biāo)函數(shù)

(10)

式中Ω表示有效能量的積分區(qū)間。由于地震數(shù)據(jù)是帶限的，計(jì)算Q過程中的積分范圍也是有限的，實(shí)際積分中只考慮S變換譜中顯著的數(shù)值區(qū)域Ω，為此，定義一個加權(quán)函數(shù)作為被積函數(shù)的因子，即

(11)

式中：h(τ,ν)為被積函數(shù)；FΩ(τ,ν)為加權(quán)函數(shù)；τmin和τmax分別表示地震信號的傳播開始時間和結(jié)束時間；νmin和νmax分別表示積分區(qū)間的上、下頻率邊界。

(12)

(13)

(14)

(15)

式中

(16)

(17)

將式(15)代入式(12)，消除子波變量，得到最終的Q值計(jì)算公式

(18)

1.3 自適應(yīng)積分區(qū)間的確定方法

上述積分區(qū)域Ω對Q值計(jì)算精度起著重要的作用。對于含有強(qiáng)烈噪聲的地震數(shù)據(jù)，傳統(tǒng)的Q估計(jì)方法很難估計(jì)高精度的Q值，如CFS、CVG和LSR。為了解決這一問題，通常采用兩步估計(jì)方法：第一步是應(yīng)用某種去噪算法，對含噪地震數(shù)據(jù)進(jìn)行去噪處理，提高地震數(shù)據(jù)的質(zhì)量； 第二步是從去噪后的地震數(shù)據(jù)中提取品質(zhì)因子Q。顯而易見，兩步法需要較高的時間成本以及依賴去噪算法的效果。因此，本文提出了一個自適應(yīng)的積分區(qū)間選擇算法，這種方法可以在地震數(shù)據(jù)中出現(xiàn)較高噪聲的情況下得到相對準(zhǔn)確的Q值。

分析式(18)，估計(jì)Q值時需要對時頻譜的時間變量τ和頻率變量ν進(jìn)行二重積分，積分區(qū)間Ω可由τmin、τmax、νmin、νmax四個參數(shù)確定。顯而易見，τmin和τmax能夠被地震信號的先驗(yàn)信息確定。根據(jù)地震信號的傳播時間，本文在合成數(shù)據(jù)算例中設(shè)置τmin=0和τmax=1s。這樣就僅需要考慮頻率參數(shù)νmin和νmax如何設(shè)置。

小波變換中常用VisuShrink算法[17]設(shè)置一個全局閾值[18]，從而確定有效能量空間。全局閾值的計(jì)算公式為

(19)

式中：σ為高斯噪聲標(biāo)準(zhǔn)方差；N為地震信號采樣點(diǎn)數(shù)。

受VisuShrink去噪算法的啟發(fā)，本文提出一種時變閾值的自適應(yīng)頻率參數(shù)選擇方案，其具體步驟如下：

(1)利用S變換得到給定地震數(shù)據(jù)的時頻譜|Sx(τ,ν)|。

(2)由于地震資料的衰減是隨波的傳播而累積的，所以對每個時刻τi∈[τmin,τmax]計(jì)算當(dāng)前的閾值ρτi，公式如下

(20)

式中：k是與信噪比相關(guān)的調(diào)整參數(shù)；Mτi表示時刻τi對應(yīng)的S變換譜系數(shù)的均值。

(3)對于每個時刻的τi，當(dāng)滿足|Sx(τi,νmin,τi)|=ρτi和|Sx(τi,νmax,τi)|=ρτi時，確定當(dāng)前低頻νmin,τi和高頻νmax,τi。

(4)所有時刻計(jì)算的上、下限頻率構(gòu)成了不規(guī)則的積分區(qū)間Ω。其中：νmin=min(νmin,τ1,νmin,τ2，…，νmin,τN)；νmax=max(νmax,τ1,νmax,τ2,…,νmax,τN)。

CVS方法將上述方案確定的積分區(qū)間應(yīng)用到式(18)，并計(jì)算非平穩(wěn)地震數(shù)據(jù)中的品質(zhì)因子Q。

2 模型測試

在上述自適應(yīng)積分區(qū)間選擇方法的基礎(chǔ)上，首先研究不同參數(shù)k對CVS方法的Q值估計(jì)結(jié)果的影響，選擇相對誤差最小時對應(yīng)的k值，并代入時變閾值計(jì)算公式，從而確定合適的積分區(qū)間； 然后，通過合成的無噪地震數(shù)據(jù)和不同信噪比的含噪地震數(shù)據(jù)驗(yàn)證本文所提出的CVS方法，并與CVG、CFS和LSR三種方法進(jìn)行對比。

2.1 自適應(yīng)積分區(qū)間中參數(shù)k的選擇

在上述提出的自適應(yīng)積分區(qū)間的確定方法中，參數(shù)k需要人為設(shè)置。為驗(yàn)證不同參數(shù)k對CVS方法所估計(jì)Q值的影響，實(shí)驗(yàn)采用主頻為30Hz的Ricker子波、理論Q為50以及如圖2所示的反射系數(shù)合成無噪地震數(shù)據(jù)。向無噪地震數(shù)據(jù)中加入不同信噪比的高斯白噪聲后，得到信噪比為5、12和20dB的含噪地震數(shù)據(jù)。針對無噪與含噪地震數(shù)據(jù)，分別選取不同范圍的k進(jìn)行Q值估計(jì)。

根據(jù)實(shí)驗(yàn)結(jié)果可得，選擇合理的k值有利于CVS方法最終估計(jì)出高精度的Q值。如圖3所示，對于無噪地震數(shù)據(jù)，應(yīng)該設(shè)置k=0.07； 對于信噪比為20、12和5dB的含噪地震數(shù)據(jù)，分別設(shè)置k的值為0.14、0.18、0.25。由圖可見，參數(shù)k與信噪比相關(guān)。當(dāng)k值設(shè)置合理時，k和Mτi的乘積能夠作為高斯標(biāo)準(zhǔn)方差的近似。無噪地震數(shù)據(jù)應(yīng)選擇較小的k，以避免S變換中幅值過小的系數(shù)，從而避免Q值估計(jì)過程中不穩(wěn)定的情況出現(xiàn)。而對于含噪地震數(shù)據(jù)還需在一定程度上抑制噪聲，所以k值隨著信噪比的增加而減小。

2.2 無噪合成地震數(shù)據(jù)算例

在本實(shí)驗(yàn)中，非平穩(wěn)地震數(shù)據(jù)由圖2所示白的偽隨機(jī)反射系數(shù)序列和主頻為30Hz的Ricker子波生成。圖4分別為Q=25、50、75、100和125的5道非平穩(wěn)地震信號。對于CVS方法，式(18)中的參數(shù)k設(shè)置為0.07。

圖4 5種理論Q合成的非平穩(wěn)地震數(shù)據(jù)

上述四種方法的Q值估計(jì)結(jié)果如圖5所示，為了更明顯地展示不同方法的估計(jì)精度，在每列柱狀圖上顯示了相應(yīng)實(shí)驗(yàn)條件下的估計(jì)結(jié)果。與理論Q值(紅色)相比，CFS(綠色)的估計(jì)結(jié)果偏大，這是因?yàn)镃FS需要假設(shè)源子波的譜是高斯譜[9]，而實(shí)驗(yàn)中采用的Ricker子波并不滿足這一假設(shè)。LSR方法(紫色)由于沒有假設(shè)地震子波，所以其估計(jì)精度高于CFS，但是低于CVG(黃色)和CVS(藍(lán)色)。CVG方法和所提出的CVS方法都有較高的精度。具體來說，當(dāng)Q=25和50時，CVG和CVS方法的估計(jì)結(jié)果比較接近；但當(dāng)Q高于75時，CVS方法估計(jì)結(jié)果更接近理論Q值。

圖5 地震子波為Ricker時，CVS、CVG、CFS和LSR方法Q值估計(jì)結(jié)果

為了進(jìn)一步對比本文所提的方法和CVG方法，本實(shí)驗(yàn)測試了Gabor變換中選取不同長度的窗函數(shù)對CVG方法所估計(jì)Q值的影響。仍對圖4中的合成地震數(shù)據(jù)進(jìn)行分析，其中Gabor變換中高斯窗的長度twin分別設(shè)置為20、40、80、120ms。圖6顯示Q=25、50、75、100和125時，采用不同窗長度的CVG方法和CVS方法的估計(jì)結(jié)果。顯然，與各種窗長的CVG方法相比，所提出的CVS方法在所有情況下都更接近理論Q值。與CVS方法相比，CVG方法對窗函數(shù)的長度很敏感。在四種時窗中，twin為80ms的窗長是最合適的，其對應(yīng)的CVG估計(jì)精度最高。twin為40ms時，CVG的估計(jì)精度略低于80ms的估值精度。選擇twin為20ms和120ms的窗函數(shù)時，CVG方法的估計(jì)誤差較大。這是因?yàn)閠win為80ms的時間窗與完整的震源子波的持續(xù)時間接近，過長或過短的時間窗都會影響時頻分辨率，從而影響Q估計(jì)的精度。因此，所提CVS方法的一個重要優(yōu)勢為S變換的窗函數(shù)和頻率相關(guān)，無需設(shè)置窗函數(shù)的長度。

圖6 地震子波是Ricker時，四種窗函數(shù)長度的CVG和CVS的Q值估計(jì)結(jié)果

為了揭示四種Q值估計(jì)方法對地震子波類型的依賴性，實(shí)驗(yàn)選擇不同類型的地震子波與圖2中的反射系數(shù)合成地震記錄，包括Ricker子波、Gauss子波和廣義地震子波(GSW)[19]。將這三個地震子波的主頻均設(shè)置為30Hz。實(shí)驗(yàn)測試了三種子波類型下，CVS、CVG、CFS和LSR方法估計(jì)Q值的相對誤差，結(jié)果如圖7所示。由圖可見，CVG方法、CVS方法和LSR方法對子波的類型并不敏感，在三種子波情況下，LSR方法的估計(jì)結(jié)果在5%附近波動； 而對于CVG和CVS方法，所有估計(jì)結(jié)果的相對誤差均保持在5%以下。進(jìn)一步比較這兩種方法，所提出的CVS方法在大多數(shù)情況下具有更小的相對誤差； 對于CFS方法，該方法僅在Gauss子波中Q值估計(jì)的相對誤差小于5%，在其他兩種子波中Q值估計(jì)是不準(zhǔn)確的。綜上所述，無論地震子波的類型如何，所提CVS方法估計(jì)的Q值均與理論Q值最接近。

圖7 不同的地震子波，CVS、CVG、CFS和LSR方法Q值估計(jì)的相對誤差(a)Ricker子波； (b)Gauss子波； (c)GSW子波

2.3 含噪合成地震數(shù)據(jù)算例

為了更進(jìn)一步研究CVS方法的抗噪性，本文估計(jì)了含噪情況下不同方法的Q值。含噪地震記錄如圖8所示。由于加入了高斯白噪聲，實(shí)驗(yàn)結(jié)果統(tǒng)計(jì)了重復(fù)估計(jì)100次Q值結(jié)果的平均值和標(biāo)準(zhǔn)差。模型數(shù)據(jù)理論Q值為40，地震子波采用Ricker子波，其他參數(shù)與上述實(shí)驗(yàn)相同。CVS、CVG、CFS和LSR方法對不同信噪比的含噪地震數(shù)據(jù)所估計(jì)Q值的統(tǒng)計(jì)結(jié)果如表1所示。此外，通過CVS方法分別從信噪比為5、12和20dB的含噪地震數(shù)據(jù)中提取的Q值結(jié)果如圖9所示。可以看出，隨著地震數(shù)據(jù)信噪比的增加，CVS方法計(jì)算結(jié)果均值的精度也在增加，且落在理論Q附近的估計(jì)值越來越多。由表1可得，無論信噪比為5、12還是20dB，本文提出的CVS方法的估計(jì)結(jié)果均是最準(zhǔn)確、穩(wěn)定的，說明CVS能夠通過自適應(yīng)參數(shù)確定合理的積分區(qū)間，從而一定程度上抑制噪聲，提高該方法抗噪能力，最終估計(jì)出穩(wěn)定的Q值。

表1 不同信噪比下CVS、CVG、CFS和LSR的Q值估計(jì)統(tǒng)計(jì)結(jié)果(理論Q為40)

圖8 不同信噪比的含噪地震數(shù)據(jù)

圖9 不同信噪比的含噪地震數(shù)據(jù)CVS方法估計(jì)結(jié)果(a)5dB； (b)12dB； (c)20dB

3 實(shí)際資料應(yīng)用

將CVS、CFS和CVG方法應(yīng)用于實(shí)際地震資料。圖10是來自中國某油田的疊后地震剖面，有2口井位于地震剖面。井1在目標(biāo)層(1.3～1.4s)鉆遇氣砂巖，獲得12.9815m3/d的工業(yè)氣流。井2的儲層不發(fā)育，證實(shí)為干井。

圖10 某油田疊后地震剖面

根據(jù)本文CVS方法，選擇目標(biāo)層附近(1.25～1.45s)的地震數(shù)據(jù)進(jìn)行Q值估計(jì)?？紤]實(shí)際資料的Q一般不是常Q，因此需要利用滑動窗多次計(jì)算Q值。針對實(shí)際數(shù)據(jù)，將非平穩(wěn)的地震數(shù)據(jù)進(jìn)行劃分[20]，然后在每段內(nèi)計(jì)算層間Q。對于局部時窗內(nèi)的地震數(shù)據(jù)，時頻振幅譜的形狀與較大幅值系數(shù)主要由時變子波確定[20]，因此在式(8)中忽略反射系數(shù)序列的時頻譜是成立的。實(shí)際地震資料受多種因素的影響，會出現(xiàn)異常的Q值結(jié)果。剔除異常結(jié)果后，采用樣條擬合方法擬合各地震道的Q值，最終顯示等效吸收剖面(Q-1)，如圖11所示。從圖11可以清楚地看出，CVS、CFS和CVG方法的估計(jì)結(jié)果在總體趨勢上是一致的。在井1位置，CVS方法估計(jì)的Q值最??； 在井2位置，CVS方法估計(jì)的Q

圖11 實(shí)際資料目標(biāo)層附近的等效吸收(Q-1)剖面(a)CVS方法； (b)CFS方法； (c)CVG方法

值最大?？傮w而言，相比于另外兩種方法，CVS方法具有更好的有效性、地震子波適應(yīng)性和穩(wěn)定性。

進(jìn)一步利用CVS計(jì)算的Q值對實(shí)際資料進(jìn)行反Q濾波[6]，補(bǔ)償后的結(jié)果如圖12所示。與圖10相比，藍(lán)色矩形和箭頭指示處的結(jié)果表明原始剖面經(jīng)反Q濾波后，同相軸變地更細(xì)且更加連續(xù)，地震資料的分辨率得到提高。因此，實(shí)際資料驗(yàn)證了CVS方法能夠估計(jì)有效的Q值，并且為后續(xù)反Q濾波提高地震資料的分辨率奠定了堅(jiān)實(shí)的基礎(chǔ)。

圖12 利用CVS方法得到的Q值對實(shí)際資料反Q濾波處理后的結(jié)果

4 結(jié)論

本文提出了一種新的基于S變換和變分法的Q值估計(jì)方法(CVS方法)。模型數(shù)據(jù)和實(shí)際數(shù)據(jù)的計(jì)算結(jié)果表明，應(yīng)用CVS方法能夠從疊后反射地震數(shù)據(jù)中有效提取品質(zhì)因子Q，實(shí)現(xiàn)儲層識別和非平穩(wěn)地震資料的衰減補(bǔ)償。具體結(jié)論如下。

(1)與CVG方法對比，將Gabor變換改為S變換，可以克服CVG方法對窗函數(shù)長度比較敏感的缺點(diǎn)。此外，利用S變換的多尺度分辨率的特點(diǎn)，可以更加精細(xì)地刻畫地震數(shù)據(jù)的非平穩(wěn)性，從而獲得更準(zhǔn)確的Q值。

(2)與CVG方法和LSR方法對比，本文提出的CVS方法能夠根據(jù)時頻譜的能量自適應(yīng)地選擇合適的積分區(qū)間，具有更好的抗噪性能。

(3)與傳統(tǒng)的CFS方法對比，本文提出的方法不需要假設(shè)地震子波的類型，具有更好的子波適應(yīng)性，并且在未知準(zhǔn)確地震子波情況下，提供了準(zhǔn)確估計(jì)Q值的理論基礎(chǔ)。