基于Lucy-Richardson算法和廣義S變換的Q值提取

李 波,文曉濤,張懿疆,何 健,黃 偉

(1.成都理工大學(xué)地球物理學(xué)院,四川成都610059;2.成都理工大學(xué)油氣藏地質(zhì)及開發(fā)工程重點實驗室,四川成都610059)

地震波在地下傳播時會受到介質(zhì)的吸收而發(fā)生能量衰減。該衰減可分為兩類:一類是由于介質(zhì)非均勻性引起的衰減,包括散射、層狀地層引起的透射等;另一類是由地層吸收引起的衰減特性,即地層本征衰減[1]。地層的吸收特征是一種非常重要的油氣指示標(biāo)志,通常用品質(zhì)因子Q來表征[2]。

利用地震資料估計Q值的方法主要有兩大類。第一類是時間域的方法,如時間域解析信號法、上升時間法、子波模擬法和振幅衰減法等;第二類是頻率域的方法,如譜比法、匹配技術(shù)和頻率域譜模擬法等。TONN[3]利用VSP數(shù)據(jù)比較了7種Q值估計方法,發(fā)現(xiàn)振幅衰減法可靠性最差,譜比法的效果相對較好。對于時間域的方法,最主要的問題是地震資料振幅信息不保真及信噪比較低,因而Q值估算精度不高;頻率域方法的問題是利用傅里葉譜估計一個時窗的地震記錄頻率譜時,對窗函數(shù)和時窗的長度要求較高。盡管兩種方法都各有利弊,但實際應(yīng)用中譜比法仍然是目前最常用的Q值估計方法之一[4]。魏文等[5]利用S變換將地震道由時域變?yōu)闀r頻域,再由時頻域能量分布及中心頻率計算出對應(yīng)時間的中心頻率,進而得到Q值,但該方法的S變換小波基函數(shù)固定,難以適應(yīng)實際地震信號;WANG[4,6]利用時頻函數(shù)將時頻譜轉(zhuǎn)變?yōu)橐痪S數(shù)據(jù),最后進行補償分析得到Q值,此方法雖然穩(wěn)定但繁瑣;王宗俊[7]改進了質(zhì)心頻移法,推導(dǎo)了子波參數(shù)、質(zhì)點頻移與地層Q值之間的關(guān)系式,進而得到地層Q值;袁煥等[8]等提出先利用VSP上、下行波的旅行時間聯(lián)合反演得到地層層速度,再利用VSP上、下波形反演Q值;張繁昌等[9]對地震信號進行子波分解得到匹配子波,將匹配子波中心頻率、中心時間之積與對數(shù)振幅進行投影,利用擬合散點斜率得到Q值,該方法使得自適應(yīng)分解得到的子波能夠應(yīng)用譜比法解決薄層和反射波干涉的問題;范明霏等[10]研究發(fā)現(xiàn)S變換窗函數(shù)具有一定局限性,提出基于廣義S變換模極大值的薄儲層刻畫方法,提高了廣義S變換薄層結(jié)構(gòu)剖面刻畫的可行性與可信度;HAO等[11]推導(dǎo)了衰減地震波的廣義S變換表達式,據(jù)此得到新的譜比法Q值反演公式,消除了窗函數(shù)的頻率響應(yīng),提高了Q值反演精度;金子奇等[12]改進了衰減旅行時層析方法,提高了衰減旅行時的計算精度,解決了上覆地層的影響,使得Q值求取較為準(zhǔn)確。

為了進一步提高Q值求取的精度,本文在傳統(tǒng)廣義S域求取Q值的基礎(chǔ)上,提出了基于Lucy-Richardson算法-廣義S域(以下簡稱LR-GST)的Q值估算方法。文中先介紹了該方法的基本原理,接著利用合成復(fù)雜信號對比了LR-GST算法與廣義S域算法(以下簡稱GST)應(yīng)用效果,分析了兩種算法的頻率分辨率及穩(wěn)定性,再利用粘彈性地層模型進一步試算,將LR-GST算法提取的Q值與GST算法提取的Q值與理論值進行對比,分析了兩種方法的精度差異,最后將LR-GST算法與GST算法應(yīng)用于實際資料的Q值估算,將兩者得到的結(jié)果與實際井資料進行對比,驗證本文方法的精度與可信度。

1 方法原理

陳學(xué)華等[13]提出的兩參數(shù)廣義S變換時域表達式為:

exp(-i2πft)dt(1)

式中:f為頻率;λ和p是高斯窗函數(shù)的尺度調(diào)節(jié)因子;τ為高斯窗函數(shù)在時間軸上滑動時中點對應(yīng)的時間;x(t)表示地震信號。

同時給出了頻率域廣義S變換的表達式:

exp(i2πfaτ)dfa

(2)

式中:X(fa)為地震信號x(t)的傅里葉變換;fa為頻移量。

趙偉等[14]將震源子波的頻域表達式表示成高斯函數(shù)的形式,即S0(f)=exp[-(2πf-2πfd)2/m],其中,fd為震源主頻;m控制震源頻帶寬度,m值越大頻帶越寬,m值趨于無窮大時,震源為脈沖震源;fd和m一般通過對目的層處的子波進行分析得到。根據(jù)Futterman理論[15]和頻率域相移原理,震源子波傳播t*時間后接收的地震子波振幅譜可表示為:

式中:Q表示t*時間內(nèi)地震子波所經(jīng)過地層的Q值;P表示由于幾何擴散、透射、反射等引起的與頻率無關(guān)的能量損失。

將公式(3)代入公式(2),得:

經(jīng)過推導(dǎo)和整理得t*時刻廣義S振幅譜[11]為(推導(dǎo)見附錄A):

高斯窗G(t)的Wigner-Ville分布:

化簡公式(6)得(推導(dǎo)見附錄B):

在此處引入Lucy-Richardson算法,該算法源于貝葉斯的條件概率定理,從最大似然估計的角度出發(fā),假設(shè)圖像服從Possion分布[16-17],通過迭代逼近的方法,使其得到較為精確的值。具體表達式[18-19]為:

式中:k+1表示迭代次數(shù);Sx為x(t)的廣義S振幅譜,令Wx(0)=Sx;Wh為高斯窗的魏格納分布,上標(biāo)*表示相關(guān)算子;?表示褶積算子;Wx(k+1)為基于LR-GST算法經(jīng)過k+1次迭代后的振幅譜。整理公式(8)得:

將公式(5)、公式(7)代入公式(9)得(推導(dǎo)見附錄C):

將不同時刻t1和t2到達的反射子波的振幅譜相比并取對數(shù),有:

式中:Δt=t2-t1;Q為t1和t2時間內(nèi)的品質(zhì)因子;P1和P2分別為兩個反射子波與頻率無關(guān)的能量損失。對地震道進行LR-GST后,取振幅譜局部極大值[20]進行對數(shù)譜比求取Q值時,譜比值與頻率不再符合線性關(guān)系。將公式(11)中的Q和Q2項提出并進行變量替換,得到簡化的Q值反演公式:

依據(jù)公式(12),將對數(shù)譜比S*(f)與γ進行線性擬合,擬合直線斜率值的倒數(shù)即為Q值。

2 方法驗證

2.1 合成復(fù)雜地震信號模型驗證

為對比LR-GST的時頻分辨率與時頻聚集性,本文采用一個合成的復(fù)雜信號進行驗證。該信號包含了低、中、高三個頻段,其中301～510ms包含了20Hz和100Hz兩個頻率分量的信息,詳細頻率信息如表1所示。該信號的表達式為:

圖1為合成信號,圖2和圖3分別為對合成信號進行GST計算和LR-GST計算得到的時頻譜。比較圖2、圖3可看出:①圖2箭頭處端點效應(yīng)明顯,出現(xiàn)了部分假頻;②GST時頻譜上兩種頻率成分交界

表1 合成信號的時頻分布

圖1 合成信號

圖2 GST計算得到的時頻譜

圖3 LR-GST計算得到的時頻譜

處時頻分辨率較差。相對而言,LR-GST時頻譜的整體分辨率較高,不同頻率成分交界處邊界清晰,時頻聚集性相對優(yōu)于GST算法。因此,LR-GST算法對于非平穩(wěn)信號中不同頻率分量的區(qū)分能力更強。

2.2 粘彈性地層模型驗證

圖4a為建立的層狀地層模型,共分5層,其中第3層是目的層(見圖4紅色虛線框),該層具有較小的速度與Q值;每層的厚度分別為100、100、80、70和162m。圖4b和圖4c分別為基于該模型利用褶積加衰減、褶積不加衰減算法得到的疊后地震剖面。圖4d 為圖4b和圖4c中第10道地震信號對比,其中實線為褶積加衰減算法獲得的信號,即考慮了非彈性衰減;虛線為褶積不加衰減算法獲得的信號,即未考慮非彈性衰減。

圖5a和圖5b分別是采用LR-GST算法和GST算法得到的等效Q值剖面。可以看出,兩圖的目的層位置都存在異常,為體現(xiàn)效果,我們提取了單道數(shù)據(jù)進行深入分析。圖6a為提取的第10道地震數(shù)據(jù),圖6b為該道數(shù)據(jù)的LR-GST計算的時頻圖,圖中清晰地顯示了反射層位置的頻譜。圖6c中黑色實線為基于LR-GST算法估算的地層等效Q值,紅色實線為基于GST算法獲得的地層等效Q值,綠色為理論等效Q值。分析該圖發(fā)現(xiàn):①兩種方法都能準(zhǔn)確地反映反射界面位置;②基于LR-GST算法獲得的地層等效Q值與理論值更接近,精度更高。圖6d是S*(f)-γ線性擬合曲線,整體上看,數(shù)據(jù)點的線性關(guān)系比較穩(wěn)定,擬合效果較好,擬合得到的直線斜率即為1/Q。

圖4 層狀模型(a)、利用褶積加衰減合成的疊后地震剖面(b)、利用褶積不加衰減合成的疊后地震剖面(c)以及單道地震信號(第10道)對比(d)

圖5 LR-GST算法得到的等效Q值剖面(a)和GST算法得到的等效Q值剖面(b)

圖6 第10道地震波(a)、第10道的LR-GST時頻圖(b)、LR-GST算法和GST算法得到的單道等效Q值估計值與理論值對比(c)以及S*(f)-γ線性擬合曲線(d)

3 應(yīng)用實例

為了驗證本文提出的基于LR-GST算法Q值估計方法的實用性,我們選取了實際地震數(shù)據(jù)進行試驗,并與基于GST算法的Q值估計方法進行了對比。

研究區(qū)目的層為二疊系中油組,整體為由北西向南東抬高的構(gòu)造背景,區(qū)域上位于鹽邊構(gòu)造帶,主要受北東向雁列狀斷裂帶及鹽體的塑性活動控制,在擠壓應(yīng)力背景下形成北東延伸的大型圈閉群,在此基礎(chǔ)上形成多個獨立的低幅度圈閉、巖性圈閉及復(fù)合圈閉。目的層為砂體,埋深大于4000m,厚度大約為9m,構(gòu)造幅度一般為5～25m,其中A井位于砂體的中間,為高產(chǎn)油氣井。圖7為該地區(qū)的一條過A井疊后地震剖面,其中黑色測井曲線為聲波速度,圖中圓圈中箭頭處為目的層位置(對應(yīng)低速)。對于該剖面,我們分別利用基于LR-GST算法(圖8)和基于GST算法(圖9)求取了Q值。從圖8和圖9可以看出,基于GST提取的低Q值區(qū)域邊界不明顯,而基于LR-GST算法得到的低Q值異常區(qū)邊界清晰,與儲層(黑圈位置)對應(yīng)較好。圖10a為提取的井旁道的地震數(shù)據(jù),圖10b為對應(yīng)的LR-GST時頻譜,時頻譜中時間分辨率和頻率分辨率均較高。圖10c為目的層處Q值的切片,從圖中能看出A井處于明顯的異常低值區(qū),表明地震波在該位置衰減明顯,可能與含油有關(guān),與實際情況相吻合。

圖7 過A井疊后地震剖面

圖8 利用LR-GST算法估算的Q值剖面

圖9 利用GST算法估算的Q值剖面

圖10 井旁道振幅譜(a)、LR-GST時頻譜(b)以及利用LR-GST算法求取的Q值沿層切片(c)

4 結(jié)論

本文在基于GST算法的Q值估計基礎(chǔ)上提出了基于LR-GST算法的Q值估計。與基于GST的Q值估計相比,基于LR-GST算法的Q值估計方法能得到分辨率較高的時頻譜,使Q值估算的精度和穩(wěn)定性更好。在實際地震資料處理中,基于LR-GST算法的Q值估計結(jié)果與鉆井結(jié)果較為吻合,證明了該方法的穩(wěn)定性和實用性。需要注意的是,該方法對參數(shù)λ和p的準(zhǔn)確性要求較苛刻,合理選擇參數(shù)λ和p能夠調(diào)節(jié)LR-GST算法的分辨率,使Q值的估算結(jié)果更準(zhǔn)確。