頻譜一致性Q因子反演方法

趙東鳳,梁伽福,馬 雄,李國發(fā),范 旭

(1.中國石油大學(xué)(北京)油氣資源與探測國家重點實驗室,北京102249;2.中國石油新疆油田分公司勘探開發(fā)研究院,新疆克拉瑪依834000)

品質(zhì)因子Q是地球物理勘探中的重要地層參數(shù),其不僅可用于基于吸收補償?shù)母叻直媛实卣鹳Y料處理[1],還可用于指示地下介質(zhì)的流體特征、孔隙度、滲透率等信息[2]。自FUTTERMAN[3]提出地層吸收衰減理論和相關(guān)計算方法以來,許多學(xué)者在時間域和頻率域發(fā)展了一系列的Q因子估計方法。時間域Q因子估計方法根據(jù)地震子波形態(tài)變化與衰減的關(guān)系進行Q值估計,如振幅衰減法、上升時間法等[4]。由于幾何擴散、透射損失等非地層吸收因素會引起子波振幅的改變,降低了時間域Q因子估算精度,因此,目前很少采用時間域方法進行Q值估算。相對而言,頻率域Q因子估計方法能夠很方便地分離出幾何擴散等與頻率無關(guān)的影響因素,逐漸成為Q值估算的主流方法[5]。譜比法[6]和質(zhì)心頻移法[7]是目前頻率域Q因子估算最為常用的方法。譜比法基于地震信號衰減函數(shù)的斜率對Q進行估算,該方法的物理意義明確,計算方法直觀,但其估算結(jié)果易受噪聲干擾和層間干涉等因素的影響。質(zhì)心頻移法基于地震信號質(zhì)心頻移的變化對Q進行估算。質(zhì)心頻移的估算具有一定的統(tǒng)計效應(yīng),其抗噪性要優(yōu)于譜比法。但是,質(zhì)心頻移法假設(shè)地震子波的頻譜為高斯函數(shù),當(dāng)?shù)卣鹱硬ǖ念l譜偏離高斯函數(shù)時,估算結(jié)果也隨之產(chǎn)生誤差[8]。為提高這兩種方法的穩(wěn)定性和適應(yīng)性,許多學(xué)者從不同角度對以上方法進行了修改和完善。ZHANG等[9]將質(zhì)心頻率修改為峰值頻率,提出了基于峰值頻率的Q估算方法。WANG等[10]提出了基于對數(shù)譜面積差值法的Q值估計方法,提高了Q估計的穩(wěn)定性和抗噪性。李擁軍等[11]將非零相位雷克子波的復(fù)數(shù)域快速匹配追蹤分解法與對數(shù)譜比法相結(jié)合,進行了微測井?dāng)?shù)據(jù)近地表Q值估計。劉國昌等[12]以地震波射線傳播為基礎(chǔ),根據(jù)同相軸局部斜率和射線參數(shù)的映射關(guān)系,將多射線頻譜同時代入譜比法進行聯(lián)合反演,有效避免了速度誤差對Q值估計的影響。李波等[13]將Lucy-Richardson算法與廣義S變換相結(jié)合,提高了Q因子估算精度。LI等[14]提出了一種雙參數(shù)回歸方法,從非線性地震衰減中估計與頻率相關(guān)的Q因子。另外,LEI等[15]借助于小波變換和時頻分析技術(shù),提高了譜比法和質(zhì)心頻移法進行Q值估計的精度。趙秋芳等[16]對近地表品質(zhì)因子估算方法進行了總結(jié)和梳理,指出井地聯(lián)合及多源信息融合是未來近地表Q反演的主要研究方向。云美厚等[17-18]就品質(zhì)因子地震反演問題進行了系統(tǒng)分析,并指出了基于品質(zhì)因子進行儲層預(yù)測的風(fēng)險。

常用Q值估計方法利用兩個不同位置的地震信號或者相同位置兩個不同時間的地震信號進行Q值估算,當(dāng)利用多個地震信號同時對地下吸收參數(shù)進行估算時,這種方法不能充分利用所有地震信號的頻率變化及其相互關(guān)系進行Q參數(shù)反演。隨著速度層析反演技術(shù)的發(fā)展[19],Q值估算逐漸與各種層析反演技術(shù)結(jié)合,通過衰減信號與Q值的關(guān)系建立矩陣方程進行Q值估計。WONG等[20]首先將VSP數(shù)據(jù)與層析反演相結(jié)合進行Q值估算。BREGMAN等[21]利用初至波信息,結(jié)合基于射線理論的層析反演與譜比法來進行Q值估計,并分析了大量實例,指出Q值層析反演中數(shù)據(jù)預(yù)處理的重要性。JAMES[22]利用不同深度地震子波的對數(shù)譜構(gòu)建層析反演方程來估計層間衰減,并將數(shù)據(jù)的高信噪比部分和模型約束作為先驗信息加入到反演方程,提高了反演的抗噪能力和穩(wěn)定性。HU等[23]基于質(zhì)心頻移法,加入一個源振幅譜函數(shù)使Q層析成像很好地適應(yīng)實際地震數(shù)據(jù),提高了衰減層析反演的精度。基于全波形反演的吸收參數(shù)建模是近年來的熱點研究內(nèi)容[24],這類方法具有明顯的理論優(yōu)勢,但在實際應(yīng)用中也遇到了諸多困難。

就工業(yè)界應(yīng)用現(xiàn)狀而言,譜比法與質(zhì)心頻移法依然是目前應(yīng)用最為廣泛的Q值估算方法。但是,這兩種方法只是利用了地震頻譜的某個單一屬性(衰減斜率或者質(zhì)心頻移),未能充分考慮地震頻譜及其空間變化的完整信息,估算結(jié)果很容易受干擾因素的影響,另外,質(zhì)心頻移法還需要高斯譜等基本假設(shè)。為此,本文提出了一種頻譜一致性Q因子估算方法,其基本思想是:基于粘彈性介質(zhì)地震波傳播理論,將不同位置或不同時間的多個地震信號反向傳播到同一位置,以反向傳播之后地震頻譜的一致性作為目標函數(shù),采用全局優(yōu)化算法對Q模型進行非線性反演。該方法既不需要已知震源子波函數(shù),也不需要對震源子波的性質(zhì)進行任何假設(shè),其估算結(jié)果不依賴于地震頻譜的某個單一屬性,而是基于地震頻譜的整體變化,因而增強了Q估算的穩(wěn)定性、抗噪性和適應(yīng)性。本文對模型數(shù)據(jù)采用頻譜一致性Q因子反演方法進行實驗分析和效果對比,并利用實際數(shù)據(jù)進行了Q因子反演和吸收補償?shù)脑囼炋幚怼?/p>

1 方法原理

根據(jù)粘彈性介質(zhì)波動理論,地震波的傳播過程在頻率域表示[3]為:

R(f)=G·H(f)·S(f)

(1)

式中:f為頻率;S(f)為震源信號振幅譜;R(f)是接收點信號的振幅譜;G表示不依賴頻率的衰減,包括幾何擴散、反射、透射損失等;H(f)為地層吸收衰減函數(shù)。H(f)表示為[25]:

(2)

式中:N為地層數(shù)目;Qj是第j個地層的品質(zhì)因子;Δtj是地震波在第j個地層的旅行時間。將公式(2)代入公式(1)中,忽略不依賴于頻率的衰減因子G的影響,有:

(3)

(3)式完整地表達了地層吸收對地震信號頻譜的改造效應(yīng)。

為了更加直觀地展示地層吸收對地震頻譜的影響,進而引出本文方法的基本原理,在一維層狀介質(zhì)假設(shè)條件下,建立了一個4層吸收介質(zhì)模型,模型參數(shù)見表1。假設(shè)震源子波為50Hz的雷克子波,利用(3)式計算地層吸收后每個界面反射信號的振幅譜,圖1a給出了這4個反射信號的振幅譜。再利用傅里葉反變換得到圖1b所示的合成地震記錄?？梢钥闯?從淺至深,來自4個不同界面(R1,R2,R3,R4)地震信號的主頻逐漸降低,有效頻寬逐漸變窄,振幅譜的形態(tài)發(fā)生了較大變化,相應(yīng)時間反射記錄振幅逐漸減弱,主波瓣寬度不斷增大。

可以設(shè)想,若將不同深度的地震信號反向傳播到同一位置,在反向傳播的過程中,依據(jù)地層Q模型對地震信號進行吸收補償,此時,所采用的Q模型與實際模型越接近,則這些信號反向傳播之后所重構(gòu)的振幅譜一致性越好,否則,不同信號重構(gòu)振幅譜的差異性越大。此即本文方法的基本思想,稱之為頻譜一致性Q因子反演方法。通常,可以采用不同的數(shù)學(xué)準則對重構(gòu)頻譜的一致性進行評價和分析,本文采用方差函數(shù)作為評價準則,其目標函數(shù)為:

F(Q)=

(4)

式中:N為不同深度或時間地震信號數(shù)目;Q=(Q1,Q2,…,QN)T為待反演的吸收模型;Si(f,Q)表示第i個地震信號的重構(gòu)振幅譜;Sj(f,Q)表示第j個地震信號的重構(gòu)振幅譜;fmin和fmax分別為地震信號的最低和最高有效頻率。

(4)式給出了頻譜一致性Q因子反演的一般性目標函數(shù),適用于地面地震和VSP地震等各類野外觀測數(shù)據(jù)的Q模型反演。如將不同位置接收的地震信號反向傳播到同一位置,需要在已知速度模型的情況下進行射線追蹤以確定地震信號衰減的軌跡和大小。為簡化問題,在后面的模型試驗和實際資料試算中均采用了一維層狀介質(zhì)模型假設(shè)進行反向逆推計算。

表1 合成地震記錄的模型參數(shù)(反射時間、反射系數(shù)和Q值)

圖1 不同反射界面地震信號的振幅譜(a)及其合成地震記錄(b)

基于(4)式所示的目標函數(shù)對Q模型進行估算是典型的非線性反演問題,該目標函數(shù)呈現(xiàn)明顯的非線性特征,模擬退火法、遺傳算法、蟻群算法和粒子群算法等是目前非線性反演問題的常用算法,相對而言,粒子群算法的基本思想更加簡潔,其結(jié)果也較穩(wěn)定[26]。

粒子群算法(particle swarm optimization,PSO)是KENNEDY等[27]受大自然中鳥類覓食的啟示而提出的一種全局尋優(yōu)智能算法。該算法的核心在于,鳥群在覓食時每只鳥作為一個單獨的信息源,它將自己所獲得的信息反饋到鳥群中去,通過每只鳥不停移動并反饋自身認知、更新群體認知,從而加快了整個群體獲得最佳食物的進程。PSO對鳥類覓食過程進行了高度概括和總結(jié),數(shù)學(xué)模型簡單、適用范圍大、收斂較快,被廣泛應(yīng)用于解決全局優(yōu)化問題[28]。

(5)

(6)

式中:i=1,2,…,M;j=1,2,…,N;ω為慣性因子,用于控制算法的收斂;pi,pt分別為局部和全局最優(yōu)位置;c1,c2為學(xué)習(xí)因子;r1,r2為0～1內(nèi)的隨機數(shù);α為控制速度權(quán)重的約束因子。

本文頻譜一致性Q反演的基本流程表述如下。

1) 對地震記錄進行頻譜分析,獲得不同反射時間tk及其對應(yīng)的振幅譜Sk(f)。

2) 給定ω,c1,c2,r1和r2等粒子群算法的控制參數(shù)。

3) 粒子群初始化。給定Q的最大值和最小值,隨機生成M個Q模型。

4) 計算每個粒子的適應(yīng)度值。亦即將反射時間tk、振幅譜Sk(f)和Q模型代入(4)式,計算頻譜一致性函數(shù)。

5) 計算每個粒子在當(dāng)前時刻的最優(yōu)值,也就是說,每個初始Q模型迭代之后的當(dāng)前最優(yōu)值。然后,由每個粒子的最優(yōu)值再篩選出整個粒子群在當(dāng)前時刻的最優(yōu)值。

6) 利用(5)式和(6)式更新粒子的速度與位置,即計算下一個時刻的Q模型及其更新速度。

7) 更新每個粒子的局部最優(yōu)值與整個粒子群的全局最優(yōu)值。

8) 重復(fù)步驟5)和步驟6),直至滿足搜索精度或到達最大迭代次數(shù)。

圖2給出了利用PSO反演Q值的流程。

圖2 利用PSO反演Q值的流程

2 模型試驗

基于一維層狀介質(zhì)模型對頻譜一致性Q因子反演方法的結(jié)果與譜比法和質(zhì)心頻移法Q因子估算結(jié)果進行對比分析。

2.1 非高斯譜源子波試驗

質(zhì)心頻移法假設(shè)地震子波的頻譜為高斯函數(shù),頻譜一致性Q因子反演方法(本文方法)基于重構(gòu)頻譜的一致性對Q因子進行反演估算,無需對震源子波及其頻譜進行任何假設(shè)。為此采用表1所示的模型基于非高斯譜對兩種估算方法的影響進行測試。俞氏子波的振幅譜是非高斯的[29],因此,選擇圖3所示的主頻為50Hz的俞氏子波作為震源函數(shù)。圖4顯示了采用本文方法和質(zhì)心頻移法的Q值估算結(jié)果,其中,黑虛線代表真實值,藍實線代表質(zhì)心頻移法估算結(jié)果,紅實線代表本文方法的估算結(jié)果。采用質(zhì)心頻移法在第1,2,3層的誤差分別為8.50%,4.10%,3.70%;采用本文方法在第1,2,3層的誤差分別為1.20%,1.70%,0.92%,由此可知,當(dāng)震源子波不滿足高斯譜假設(shè)時,采用質(zhì)心頻移法進行Q估計的準確性會下降,而采用本文方法則準確性會高很多。

圖3 俞氏子波(a)及其歸一化振幅譜(b)

圖4 采用本文方法與質(zhì)心頻移法的Q值估算結(jié)果

2.2 抗噪性試驗

對圖1b所示的雷克子波合成地震記錄分別加入振幅強度為5%,10%,15%的高斯噪聲,獲得如圖5所示的含噪地震記錄。在這3種強度的噪聲下分別計算100次,圖6和圖7分別給出了10%和15%噪聲環(huán)境下采用譜比法、質(zhì)心頻移法和本文方法得到的Q值估計結(jié)果的概率分布,表2給出了3種噪聲環(huán)境下估算結(jié)果的均值和方差。

圖5 加入不同振幅強度隨機噪聲后的合成地震記錄

分析圖6、圖7和表2可以發(fā)現(xiàn):在5%噪聲環(huán)境下,采用這3種方法均能得到可接受的反演結(jié)果;在10%噪聲環(huán)境下,采用譜比法、質(zhì)心頻移法和本文方法在第3層的方差分別為257.12,142.35,63.37;在15%噪聲環(huán)境下,采用譜比法、質(zhì)心頻移法和本文方法在第3層的方差分別為361.15,183.20,97.96。可見,隨著噪聲強度的增大,譜比法和質(zhì)心頻移法的反演誤差迅速增大,本文方法的優(yōu)勢逐漸凸顯出來。

圖6 加入10%隨機噪聲時不同方法得到的Q因子概率分布結(jié)果a～c,d～f,g～i分別為3種方法對第1層、第2層、第3層估算的Q值(這3層的Q理論值分別為40,70,100)

圖7 加入15%隨機噪聲時不同方法得到的Q因子概率分布結(jié)果a～c,d～f,g～i分別為3種方法對第1層、第2層、第3層估算的Q值(這3層的Q理論值分別為40,70,100)

表2 不同噪聲環(huán)境下3種方法Q值估計結(jié)果的均值與方差

3 實際數(shù)據(jù)測試

圖8是疊前時間偏移后的地震剖面,該數(shù)據(jù)在疊前經(jīng)過了地表一致性反褶積處理,但沒有進行預(yù)測反褶積等提高分辨率的處理,因此,地震數(shù)據(jù)一定程度上保持了實際地層的吸收特征。

不同于表1所示的模型數(shù)據(jù),實際地震數(shù)據(jù)中不同時間的地震反射相互干涉,很難得到不同時間的孤

圖8 疊前時間偏移地震剖面

立子波及其頻譜,需要采用Gabor變換等時頻分析方法盡量降低子波干涉對頻譜分析的影響。另外,吸收估算所依據(jù)的是不同時間地震頻譜的差異,只有在時間間隔較大時,其頻譜差異才能保證Q反演的穩(wěn)定性。為此,本文試驗選擇了圖8中黃線所示的6個反射層(H1～H6),其時間間隔大致在500ms。以CDP3358的地震道為例,說明Q反演的具體過程。首先對地震道進行時頻分析,然后再以反射層為中心,在200ms的時窗內(nèi),沿時間方向進行三角加權(quán)平均,得到反射時間分別為0.42,0.82,1.41,1.93,2.38,3.02s這6個反射層的振幅譜(圖9),可見,主頻隨反射時間增大而減小的趨勢十分明顯。然后,進行頻譜一致性Q因子反演,反演結(jié)果如圖10所示。進一步利用反演的Q因子對該地震道進行吸收補償,圖11顯示了補償之后6個反射時間的振幅譜,可以看出,不同深度地震反射在吸收補償之后,其振幅譜在整體趨勢上趨于一致,表明反演的Q模型較好地刻畫了實際地層的吸收特征。圖12是采用本文方法反演得到的Q模型。圖13是基于圖12的Q模型進行吸收補償前、后的地震剖面,可以看出,補償之后地震剖面的分辨率得到了明顯改善,更好地揭示了大套地層之間的層序結(jié)構(gòu)及其接觸關(guān)系。圖14是吸收補償前、后的振幅譜對比結(jié)果,振幅譜是圖13所有地震道的平均值,可以看出,地震信號的主頻和頻寬得到了明顯提升。為了更好地展示吸收補償對地震反射特征的恢復(fù)能力,圖15展示了吸收補償前、后地震數(shù)據(jù)與合成地震記錄的標定結(jié)果,其中,補償前、后合成地震記錄的子波分別從其地震剖面統(tǒng)計估算,可以看出,補償之后的地震記錄與合成地震記錄在反射特征上具有較好的一致性,表明采用頻譜一致性反演方法估算的Q模型較好地描述了實際地層的吸收特征。

圖9 CDP3358位置處的6個不同反射時間的振幅譜

圖10 CDP3358位置處反演的Q值曲線

圖11 CDP3358位置處的6個不同反射時間的振幅譜(補償之后)

圖12 采用本文方法反演得到的Q模型

圖13 吸收補償前(a)、后(b)的地震剖面

圖14 吸收補償前、后的振幅譜

圖15 實際數(shù)據(jù)吸收補償前(a)、后(b)與合成地震記錄對比結(jié)果

4 結(jié)論與認識

譜比法和質(zhì)心頻移法基于地震頻譜的某個單一屬性對Q參數(shù)進行反演和估算,反演結(jié)果易受噪聲干擾和地震干涉等因素的影響。不同于以上兩種方法,頻譜一致性反演方法基于吸收補償之后地震頻譜在整體形態(tài)上的一致性對Q因子進行評價和反演,反演結(jié)果的穩(wěn)定性和抗噪性更好。另外,本文采用的粒子群算法能較好地適應(yīng)Q因子的全局尋優(yōu)過程。

需要指出的是,盡管本文方法的模型試驗和實際資料分析都采用了疊后地震數(shù)據(jù),但是,無論是理論方法還是實現(xiàn)過程,頻譜一致性反演方法很容易擴展到疊前地震數(shù)據(jù)、VSP數(shù)據(jù)甚至微測井地震數(shù)據(jù)。