楊 震,劉俊州,時 磊,韓 磊,吳高奎

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

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

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

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

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

1 技術方法

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

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

(1)

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

v0=

(2)

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

v0可以通過遞歸法計算得到:

v0=Q0Q1…Qn…QN-1vN

(3)

其中,vN為初始項。

(4)

Qn為地層n的傳播矩陣。

Qn=EnFn

(5)

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

(6)

(7)

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

1.2 反演原理

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

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

(8)

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

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

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

(9)

式中:下角標n表示第n層。

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

(10)

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

Φ=W*R

(11)

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

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

(12)

式中:J為雅可比矩陣。

(13)

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

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

(14)

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

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

(15)

其中,

(16)

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

(17)

其中,

(18)

(19)

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

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

(20)

其中,

(21)

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

(22)

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

(23)

(24)

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

2 模型測試

2.1 反射系數(shù)對比

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

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

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

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

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

2.2 薄互層模型反演

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

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

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

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

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

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

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

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

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

3 應用實例

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

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

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

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

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

圖9 目標區(qū)域井旁角道集(紅色曲線為伽馬測井曲線;綠色的方框為目標層)

圖10 井旁道反演結果對比

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

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

4 結論

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

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

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

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