楊 震,劉俊州,時(shí) 磊,韓 磊,吳高奎

(1.頁巖油氣富集機(jī)理與有效開發(fā)國家重點(diǎn)實(shí)驗(yàn)室,北京100083;2.中國石化彈性波理論與探測技術(shù)重點(diǎn)實(shí)驗(yàn)室,北京100083;3.中國石油化工股份有限公司石油勘探開發(fā)研究院,北京102206)

我國致密砂巖氣地質(zhì)資源量為17.4×1012～25.1×1012m3,可采儲量為8.8×1012～12.1×1012m3[1],目前已形成四川盆地須家河組和鄂爾多斯盆地上古生界等2個(gè)致密砂巖氣大氣區(qū)。自2009年以來,四川盆地須家河組天然氣勘探潛力逐漸引起關(guān)注[2-5],該致密砂巖儲層成層性明顯且具有儲層厚度薄的特點(diǎn),優(yōu)質(zhì)儲層的單層厚度在4～9m。使得基于精確Zoeppritz方程和相應(yīng)近似正演模擬不再適用于致密砂巖薄儲層。這是因?yàn)閆oeppritz方程只適用于單一反射界面,對于層狀模型,尤其是含有薄互層的儲層,地震記錄中的多次波、轉(zhuǎn)換波和透射損失等更為復(fù)雜,Zoeppritz方程難以適用。

為解決上述問題,需要能模擬薄儲層波場響應(yīng)的正演模擬算法。反射率法使用了波動方程的解析解,可模擬層內(nèi)的復(fù)雜波場響應(yīng),包括反射波、透射波、轉(zhuǎn)換波和多次波。MALLICK等[6]和KENNETT[7-8]指出,雖然反射率法基于一維模型假設(shè),但由于采用了波動方程的解析解,因此該類方法比有限差分法等正演方法有更高的模擬精度且計(jì)算速度更快。FUCHS[9]首次提出將反射率法用于合成地震記錄,該方法可模擬層內(nèi)的復(fù)雜地震波場響應(yīng),包括反射波、透射波、轉(zhuǎn)換波和多次波。在此基礎(chǔ)上,KENNETT等[7-8,10-11]首次提出了用于計(jì)算層狀地層的總反射和透射系數(shù)的遞歸算法,該方法的特點(diǎn)是對所有頻率和慢度無條件穩(wěn)定。為提高計(jì)算速率,PHINNEY等[12]在KENNETT研究成果的基礎(chǔ)上提出了快速反射率法,即將原有的反射率法進(jìn)行重組,實(shí)現(xiàn)矢量化,使計(jì)算速度提高了約20倍。YANG等[13]基于反射率法推導(dǎo)了針對單薄層的二階近似公式,并利用遞歸算法將其推廣到薄互層。

在地震疊前反演方面,國內(nèi)外學(xué)者已經(jīng)嘗試將反射率法用于薄儲層AVA反演。SEN等[14]比較了KENNETT的反射率法與快速反射率法,并建立基于Gauss-Newton的反演方法用于疊前波形反演。王建花[15]在層參數(shù)等效的情況下,推導(dǎo)了P波入射到單薄層時(shí)的波場能量再分配的簡化公式,并進(jìn)行了初步的P波反演試驗(yàn)。在這些研究中,反演均是在截距時(shí)間和射線參數(shù)域(即τ-p域)進(jìn)行的。然而,當(dāng)可用偏移距數(shù)據(jù)較少時(shí),τ-p域變換會受到混疊現(xiàn)象的影響。為了避免混疊現(xiàn)象,MALLICK等[16]建議可生成角道集與實(shí)際道集數(shù)據(jù)相匹配。LIU等[17]針對P波AVA反演,提出基于快速反射率法的Bayesian反演方法。YANG等[18]基于快速反射率法對PP波進(jìn)行反演,并在反演中考慮到稀疏約束。LIU等[19]提出基于快速反射率法的遺傳算法,用于AVA疊前聯(lián)合反演。盡管已有多位學(xué)者發(fā)表相關(guān)文章,但受到公式本身復(fù)雜性的制約,基于反射率法的疊前反演并沒有被廣泛應(yīng)用,仍需要進(jìn)一步的研究。

本文將基于快速反射率法的AVA反演應(yīng)用于致密砂巖儲層。在快速反射率法理論的基礎(chǔ)上,利用Gauss-Newton法建立目標(biāo)函數(shù)。為了估計(jì)彈性參數(shù)(縱波和橫波速度以及密度),目標(biāo)函數(shù)是通過角度域中模擬數(shù)據(jù)和觀測數(shù)據(jù)之間的最小差異實(shí)現(xiàn)。最后,分別使用薄互層模型和實(shí)際數(shù)據(jù)進(jìn)行基于反射率法的反演,并與基于精確Zoeppritz方程的反演方法進(jìn)行比較,結(jié)果表明,基于快速反射率法的反演結(jié)果與真實(shí)值更加接近。

1 技術(shù)方法

1.1 反射系數(shù)的計(jì)算

KENNETT[7-8]推導(dǎo)出了反射率法的迭代公式,可以計(jì)算目標(biāo)區(qū)域的總反射系數(shù)及總透射系數(shù),包括多次波和轉(zhuǎn)換波等。為提高運(yùn)算速度,PHINNEY等[12]對KENNETT[7-8]反射率法進(jìn)行了改進(jìn),將最內(nèi)層的循環(huán)改為計(jì)算頻率,并拓展到整個(gè)頻率域,實(shí)現(xiàn)矢量化,不僅加快了計(jì)算速度,還能實(shí)現(xiàn)對彈性參數(shù)偏導(dǎo)數(shù)的解析計(jì)算,提高了后續(xù)反演工作的精度和效率。據(jù)PHINNEY等[12]的研究成果,頻率域PP波的總反射系數(shù)公式可以寫為:

(1)

式中:p為慢度;ω為角頻率;v01、v04分別為數(shù)組v0的第1個(gè)和第4個(gè)元素;v0為PHINNEY基于Haskell-Dunkin minor matrix方法在慢度-頻率域(p,ω)定義的含有6個(gè)元素的數(shù)組,則有:

v0=

(2)

式中:Δ為比例系數(shù);Rij(i=P或S,j=P或S)為對應(yīng)波類型的總反射系數(shù);detR為反射系數(shù)行列式。

v0可以通過遞歸法計(jì)算得到:

v0=Q0Q1…Qn…QN-1vN

(3)

其中,vN為初始項(xiàng)。

(4)

Qn為地層n的傳播矩陣。

Qn=EnFn

(5)

式中:En為層穿越矩陣,描述了平面波在第n層地層中傳播時(shí)相位的移動;Fn為界面矩陣。En和Fn可以具體寫為:

(6)

(7)

式中:qP,qS分別為縱波、橫波的垂直慢度;Tn源自Dunkin矩陣,是第n層大小為6×6的矩陣,與頻率無關(guān),依賴于第n層的水平慢度p和彈性參數(shù),具體形式可以參看PHINNEY等[12]的文獻(xiàn)。

1.2 反演原理

本文的反演基于Gauss-Newton法,其核心思想是利用快速反射率法計(jì)算層狀地層PP波反射系數(shù),最大程度地模擬目標(biāo)地震記錄,以此來估算目標(biāo)區(qū)域的彈性參數(shù)(縱波速度α,橫波速度β和密度ρ)。為得到合適的模型參數(shù),使得模型響應(yīng)和觀測數(shù)據(jù)的差異最小,反演的目標(biāo)函數(shù)定義為:

Q(M)=‖Φ-Φobs‖2

(8)

式中:M=(α,β,ρ)為模型彈性參數(shù)向量;Φobs為實(shí)際地震記錄;Φ為模擬地震記錄。

需要注意的是,這里使用的Φobs和Φ均為時(shí)間域角道集,而本文所述的快速反射率法在慢度-頻率域計(jì)算反射系數(shù),因此需要進(jìn)行域的轉(zhuǎn)換。由慢度-頻率域轉(zhuǎn)換到時(shí)間-角度域的過程可以分成以下兩個(gè)步驟:

1) 慢度轉(zhuǎn)換為角度。在反射系數(shù)R(p,ω)的計(jì)算中有3個(gè)循環(huán),即最外層循環(huán)是慢度p,中間循環(huán)為層循環(huán)(n)和最內(nèi)層循環(huán)為頻率循環(huán)。一般來說,慢度在每一個(gè)最外層循環(huán)上是固定的,入射角在每一層上都會發(fā)生變化,即θn=arcsin(αnp)。但將入射角設(shè)置在最外層循環(huán),每層的慢度都會發(fā)生變化,因此,通過公式(9),R(p,ω)被重采樣為R(θ,ω),最后得到角度道集。

(9)

式中:下角標(biāo)n表示第n層。

2) 頻率域轉(zhuǎn)到時(shí)間域。類比FRYER[20],KENNETT[7-8]的處理方法,對反射系數(shù)做一次逆Fourier變換,最終得到截距時(shí)間-角度域(τ-θ)的反射系數(shù)為:

(10)

得到時(shí)間-角度域的反射系數(shù)后,地震記錄Φ可以通過反射系數(shù)R與子波W的褶積得到。

Φ=W*R

(11)

給定符合假設(shè)條件的初始模型M0(能反映出地層的大致形態(tài)),利用Gauss-Newton法[21-24],模型的更新參數(shù)公式可以寫為:

ΔM=[JTJ]-1JT(Φ-Φobs)

(12)

式中:J為雅可比矩陣。

(13)

通過迭代ΔM更新初始彈性參數(shù)向量M0直到目標(biāo)函數(shù)Q足夠小,得到反演結(jié)果。公式(12)的難點(diǎn)在于其不穩(wěn)定性,這里應(yīng)用Levenberg-Marquardt方法修改公式(12),得到:

ΔM=[H+λI]-1JT(Φ-Φobs)

(14)

式中:H=JTJ,為海森矩陣;λ為比例系數(shù);I為單位矩陣。

根據(jù)公式(10)和公式(13),雅可比矩陣J具體可以寫為:

(15)

其中,

(16)

其中,v0對Mn的偏導(dǎo)數(shù)為:

(17)

其中,

(18)

(19)

式中:Mn代表第n層的α,β,ρ。

對于Fn的偏導(dǎo)數(shù)具體可以寫為:

(20)

其中,

(21)

En對M的偏導(dǎo)數(shù)可以求出解析解,具體可以寫成以下形式:

(22)

式中:A=qP-p2/qS,B=qP+p2/qS。

(23)

(24)

在E拓展后,對應(yīng)的v0也需要做同樣的拓展,拓展方法可參考PHINNEY[12]。

2 模型測試

2.1 反射系數(shù)對比

鑒于實(shí)際地層中薄層不能單獨(dú)存在的原因,為驗(yàn)證快速反射率法在AVA建模和反演中的有效性,采用多層模型測試其中薄層的AVA曲線特點(diǎn),模型參數(shù)如表1所示。PP波合成角道集由快速反射率法和精確Zoeppritz方程分別進(jìn)行模擬,子波采用40Hz的Ricker子波。與基于精確Zoeppritz方程的AVA道集相比(圖1),基于快速反射率法的AVA道集包含更完整的信息,例如圖1c中所示的多次波。這是因?yàn)槭褂镁_Zoeppritz方程的合成地震記錄沒考慮多次波,可能會導(dǎo)致解釋和反演中出現(xiàn)錯(cuò)誤。反射率法考慮了波在地層中傳播時(shí)的傳播效應(yīng),可以計(jì)算傳播過程中所有的反射系數(shù),由于現(xiàn)有的地震資料處理技術(shù)難以在不影響初級反射振幅的情況下去除波在地層中的傳播效應(yīng),如透射損失、層間多次波和幾何擴(kuò)散等,因此,采用快速反射率法可以更逼真地模擬地震記錄。

為了進(jìn)一步顯示快速反射率法和精確Zoeppritz方程方法之間的差異,沿圖1中黃線位置提取瞬時(shí)振幅隨入射角的變化數(shù)值,如圖2所示。受到傳播效應(yīng)的影響,使用快速反射率法計(jì)算兩個(gè)界面的振幅隨深度的增加而逐漸減小,比使用Zoeppritz方程模擬的振幅小,主要是由于該方法考慮了波在傳播時(shí)的透射損失。而基于Zoeppritz方程計(jì)算的AVA曲線不會隨深度的變化而變化。如果在地震數(shù)據(jù)處理程序中未正確校正傳輸損耗和層間多次波,則不建議基于Zoeppritz進(jìn)行正演模擬,利用快速反射率法模擬的道集將與現(xiàn)場數(shù)據(jù)更匹配。

表1 等厚模型參數(shù)

圖1 PP波AVA道集對比(紅色曲線代表P波速度)a 基于快速反射率法的AVA道集; b 基于精確Zoeppritz方程的AVA道集; c 兩種方法模擬的AVA道集之間的差異

圖2 AVA曲線對比(Phinney 1和Phinney 2是圖1a中的兩個(gè)界面,Zoeppritz 1和Zoeppritz 2是圖1b中的兩個(gè)界面)

2.2 薄互層模型反演

將薄互層模型的單個(gè)薄層的厚度設(shè)置為8m,具體參數(shù)如表2所示。利用快速反射率法計(jì)算反射系數(shù)時(shí),頻率的計(jì)算范圍為0～125Hz,合成地震記錄所用子波為40Hz Ricker子波,采樣率為1ms。考慮到實(shí)際數(shù)據(jù)的角道集數(shù)據(jù)入射角范圍通常較窄,該薄互層模型角道集的角度范圍是5°～30°,如圖3a所示。圖3b為對地震記錄添加15%隨機(jī)噪聲后的地震記錄。此外,該模型為深度域模型,后續(xù)反演得到的模型更新量是時(shí)間域模型,需將深度域模型轉(zhuǎn)換到時(shí)間域(PP波時(shí)間域),各個(gè)界面對應(yīng)于表2中的雙程旅行時(shí)。

表2 薄互層模型參數(shù)

圖3 薄互層模型PP波角道集記錄(入射角為5°～30°,紅色曲線為縱波速度)a 不加噪聲的地震記錄; b 添加15%噪聲的地震記錄

以平滑后的真實(shí)模型為初始模型,反演的目標(biāo)數(shù)據(jù)為圖3所示地震記錄。分別采用基于快速反射率法和精確Zoeppritz的反演方法進(jìn)行反演,反演結(jié)果如圖4所示。由反演結(jié)果可以看出,基于快速反射率法的反演結(jié)果基本與真實(shí)值重合。而基于Zoeppritz方程的反演結(jié)果會有一些波動,與地震記錄中的多次波有關(guān)。圖5為反演迭代次數(shù)收斂曲線,紅色為基于快速反射率法的收斂曲線;綠色為基于Zoeppritz方程的收斂曲線。若是將收斂精度定義在剩余誤差絕對值小于1時(shí)反演收斂,那么基于快速反射率法的反演在迭代3次后收斂,基于Zoeppritz方程的反演在迭代4次后收斂。

圖4 基于快速反射率法和精確Zoeppritz方程的反演結(jié)果

與之前的反演類似,均以平滑模型為初始模型,圖6為兩種反演方法的反演結(jié)果。受地震記錄中隨機(jī)噪聲的影響,兩種反演方法的反演結(jié)果中均含有隨機(jī)噪聲。但是在目標(biāo)層位置處,特別是單獨(dú)的PP波反演,反演結(jié)果仍然能較好地反演地層的真實(shí)參數(shù)。同時(shí),基于快速反射系數(shù)法的反演結(jié)果更逼近于真實(shí)值,說明該方法具有更強(qiáng)的抗噪性。

圖5 不添加噪聲的薄互層模型AVA反演迭代次數(shù)收斂曲線

圖6 基于快速反射率法和精確Zoeppritz方程的反演結(jié)果(對模型地震記錄添加15%隨機(jī)噪聲后)

需要強(qiáng)調(diào)的是,這里的模擬地震記錄包含層間多次波,并考慮了透射損失,因此將該方法應(yīng)用到實(shí)際數(shù)據(jù)的時(shí)候,實(shí)際數(shù)據(jù)含透射損失和多次波。此外,本反演方法沒有反演地層厚度參數(shù),未討論厚度對反演結(jié)果的影響,這要求使用的初始模型能大體反應(yīng)地層厚度的變化。

3 應(yīng)用實(shí)例

實(shí)際數(shù)據(jù)選自四川某工區(qū),其目的層段為須家河組的須三段。由于龍門山造山帶北段和西部山系的持續(xù)隆升,造成須三段沉積時(shí)期盆地邊緣大量沉積物快速堆積,物源以近南北向?yàn)橹?在盆地中央深處發(fā)育湖相沉積,目的層巖性組合為致密砂巖、泥巖互層,埋深4000～5000m。

在應(yīng)用本文反演方法之前,首先,要確保數(shù)據(jù)已經(jīng)被合理地處理過。由于反射率法在計(jì)算反射系數(shù)時(shí)會考慮透射損失以及多次波,因此地震數(shù)據(jù)處理時(shí),可以不用考慮透射損失以及多次波的處理。其次,利用井?dāng)?shù)據(jù)做正演模擬,然后對地震數(shù)據(jù)進(jìn)行標(biāo)定,這個(gè)過程是將井?dāng)?shù)據(jù)(深度域)與地震數(shù)據(jù)(時(shí)間域)進(jìn)行對應(yīng),對于該數(shù)據(jù),標(biāo)定是基于商業(yè)軟件完成,所以標(biāo)定的過程是在Zoeppritz方程的基礎(chǔ)上完成的。標(biāo)定完成后,模擬的地震記錄與實(shí)際數(shù)據(jù)之間的相似度約0.78?；赯oeppritz和基于反射率法的合成地震記錄與實(shí)際數(shù)據(jù)的對比如圖7所示,為了更細(xì)致地觀察二者的區(qū)別,在圖7中選3個(gè)時(shí)間點(diǎn)(T1,T2,T3)提取AVO曲線,AVO曲線對比如圖8所示。由圖8可以看出,無論是數(shù)據(jù)振幅的幅值,還是AVO曲線的走向趨勢,均是基于反射率法的地震記錄與實(shí)際地震記錄更加接近。最后是初始模型的選擇,將實(shí)際測井?dāng)?shù)據(jù)限頻后的曲線作為初始模型,頻率的選取范圍為0～60Hz。需要注意的是,反射率法的基本假設(shè)條件是數(shù)據(jù)為局部一維數(shù)據(jù)或是水平層狀。結(jié)合疊前時(shí)間偏移,該方法可以適用于地質(zhì)條件簡單的地區(qū)。在地質(zhì)條件復(fù)雜的地區(qū),需要進(jìn)行疊前時(shí)間偏移或逆時(shí)偏移。本文數(shù)據(jù)因經(jīng)過疊前偏移處理,所以是水平地層的剖面。

圖7 井旁地震道及不同方法合成地震記錄(紅色曲線為縱波速度測井曲線)a 實(shí)際地震數(shù)據(jù); b 基于快速反射率法的合成地震記錄; c 基于Zoeppritz方程的合成地震記錄

圖8 圖7中T1,T2,T3位置處的AVO曲線對比(散點(diǎn)代表振幅數(shù)據(jù)點(diǎn);曲線代表回歸曲線)a T1時(shí)間; b T2時(shí)間; c T3時(shí)間

為了驗(yàn)證基于反射率法反演的優(yōu)越性,分別基于精確Zoeppritz方程和反射率法對實(shí)際資料進(jìn)行AVA反演。反演的目標(biāo)區(qū)是2150～2350ms,井旁角道集數(shù)據(jù)如圖9所示。共有6個(gè)角道集,范圍為3°～28°。該區(qū)域主要的地層特征為砂泥巖互層,圖中綠色方框所圈出區(qū)域?yàn)槟繕?biāo)層段。圖10為井旁道地震數(shù)據(jù)反演結(jié)果,由反演結(jié)果可以看出,基于反射率法的反演結(jié)果與測井曲線更加貼合,能顯示出更多的細(xì)節(jié)。值得注意的是,相對于縱波速度反演結(jié)果,橫波速度和密度的反演結(jié)果相對較差,這是因?yàn)榉瓷湎禂?shù)公式本身對縱波速度更加敏感[25]。圖11和圖12 分別為過井的基于快速反射率法和基于精確Zoeppritz方程的反演剖面。反演中所使用的子波需要對每個(gè)入射角分別提取。盡管受到地震分辨率的限制,反演結(jié)果與測井相比頻率較低,但從反演結(jié)果可以看出,基于反射率法的反演結(jié)果要更接近于測井曲線,能更準(zhǔn)確地反映地層變化的趨勢,對識別有效的致密砂巖儲層非常有幫助。

圖9 目標(biāo)區(qū)域井旁角道集(紅色曲線為伽馬測井曲線;綠色的方框?yàn)槟繕?biāo)層)

圖10 井旁道反演結(jié)果對比

圖11 基于反射率法的反演結(jié)果a 縱波速度反演結(jié)果; b 橫波速度反演結(jié)果; c 密度反演結(jié)果

圖12 基于Zoeppritz方程的反演結(jié)果a 縱波速度反演結(jié)果; b 橫波速度反演結(jié)果; c 密度反演結(jié)果

4 結(jié)論

討論了基于快速反射率法的疊前反演理論,展示了其在致密砂巖儲層中的應(yīng)用效果。該反演是在時(shí)間-角度域而不是τ-p域中實(shí)現(xiàn),從而避免了積分變換中的混疊。

實(shí)際地震數(shù)據(jù)與基于快速反射率法和精確Zoeppritz方程合成地震記錄的比較,證明在地震數(shù)據(jù)的幅值及AVO曲線趨勢方面,基于快速反射率法的模擬結(jié)果與實(shí)際地震記錄更加接近。

模型數(shù)據(jù)的反演結(jié)果表明,基于快速反射率法的反演方法收斂速度更快且具有一定的抗噪能力。此外,基于快速反射率法的反演結(jié)果更加準(zhǔn)確。實(shí)際數(shù)據(jù)的反演結(jié)果顯示,基于快速反射率法的反演結(jié)果在精度和連續(xù)性方面優(yōu)于基于精確Zoeppritz方程的反演結(jié)果。

雖然展示了PP波反演的應(yīng)用結(jié)果,但該反演方法也可應(yīng)用于PP和PS波的聯(lián)合反演,這需對正演模擬算子進(jìn)行修改,并考慮如何匹配PP波和PS波。需要強(qiáng)調(diào)的是,反演前的處理過程應(yīng)包括與地表相關(guān)的多次波衰減、幾何擴(kuò)展補(bǔ)償?shù)?但應(yīng)排除透射損失補(bǔ)償和層間多次波。