徐振旺

(中國石油天然氣股份有限公司遼河油田分公司勘探開發(fā)研究院,遼寧盤錦124010)

真實的油氣儲層被認(rèn)為最接近于含流體的孔隙介質(zhì),即主要由固體骨架顆粒和孔隙流體(如油、氣和水等)組成。研究地震波在此類介質(zhì)中的傳播規(guī)律,對油氣勘探開發(fā)具有重要的意義。近年來,孔隙介質(zhì)彈性波傳播的數(shù)值模擬研究屢有發(fā)表[1-3]。

地震波數(shù)值模擬是描述和認(rèn)識地震波傳播規(guī)律的主要手段之一。考慮數(shù)值模擬的計算可行性,將無限的介質(zhì)區(qū)域人為地截斷為有限區(qū)域,截斷的邊界稱為人工邊界。然而,人工邊界的引入會導(dǎo)致強(qiáng)烈的邊界虛假反射,影響數(shù)值模擬的真實性和準(zhǔn)確性。為解決這一問題,人們研究了一系列的方法技術(shù),主要包括兩類:透射邊界條件和衰減吸收層技術(shù)。

透射邊界條件的本質(zhì)是在人工邊界處改用與原定解問題中不同的波動方程來求解全區(qū)域波場。最經(jīng)典的透射邊界條件是由CLAYTON等[4]提出的基于單程波動方程的旁軸近似吸收邊界條件。由于該邊界條件僅利用人工邊界處的單層介質(zhì)來吸收邊界反射,故只需引入額外一層的計算量;而衰減吸收層技術(shù)則需要開辟與層數(shù)相關(guān)的內(nèi)存進(jìn)行計算,因此透射邊界條件的計算效率相對更高。CLAYTON等提出的吸收邊界條件所采用的低階旁軸近似,能有效吸收近垂直入射的波場,但對大角度入射波場的吸收效果不佳。為提高對大角度入射波場的吸收能力,此后又發(fā)展了高階旁軸近似吸收邊界條件[5],然而高階旁軸近似吸收邊界條件要求時間步長必須滿足某一條件的約束,這意味著時間遞推穩(wěn)定性條件變得更為苛刻,在時間域顯式遞推求解時需引入小時間步長,從而導(dǎo)致計算效率的降低。

衰減吸收層技術(shù)是在人工邊界附近設(shè)置一定厚度的吸收層,使波在傳播到吸收層時,按指數(shù)規(guī)律逐步衰減。目前最為成功的衰減吸收層技術(shù)是BéRENGER[6]在求解Maxwell方程時提出的完全匹配層(perfectly matched Layer,PML)技術(shù)。理論上PML技術(shù)的吸收層內(nèi)無波阻抗差異,因此反射系數(shù)為零。相較于透射邊界條件,PML技術(shù)適用于較大范圍入射角的地震波吸收,同時不影響時間遞推穩(wěn)定性條件。CHEW等[7]率先驗證了PML技術(shù)在地震波數(shù)值模擬中的有效性;王守東[8]研究了PML技術(shù)在聲波方程中的應(yīng)用;陳可洋[9-10]提出了基于聲波方程數(shù)值模擬的PML改進(jìn)算法;高剛等[11]對PML的衰減因子進(jìn)行了詳細(xì)討論,并對PML衰減因子的設(shè)定進(jìn)行了理論研究;CHEN等[12]將PML技術(shù)應(yīng)用于任意角度標(biāo)量波動方程的數(shù)值模擬。上述PML研究均基于有限差分法,劉有山等[13]將PML技術(shù)應(yīng)用于三角網(wǎng)格剖分的顯式有限元地震波數(shù)值模擬。目前傳統(tǒng)的分裂式PML技術(shù)研究較之前有所減少[14],張衡等[15]完成了分裂式的PML在TTI介質(zhì)聲波方程模擬中的加載;馬銳等[16]提出在偽譜法彈性波模擬中使用SPML和海綿邊界的混合邊界條件,并取得了良好的吸收效果。

近20年來,PML技術(shù)不斷改進(jìn),其中較為成功的改進(jìn)技術(shù)是將經(jīng)典PML中的復(fù)數(shù)坐標(biāo)變換為復(fù)頻移拉伸算子[17],稱之為復(fù)頻移拉伸(complex frequency shifted,CFS)PML(CFS-PML)技術(shù)。相較于傳統(tǒng)的PML,CFS-PML能更有效地吸收近似掠射的入射波;非分裂形式的CFS-PML無需進(jìn)行波場分離,可節(jié)省計算機(jī)內(nèi)存。目前,非分裂形式的CFS-PML實現(xiàn)方法主要有卷積法(convolutional PML,CPML)[18]、積分法[19]和輔助微分方程法[20]。其中CPML應(yīng)用廣泛,最初應(yīng)用于求解Maxwell方程組,而后逐步應(yīng)用于地震波數(shù)值模擬[21-25]。MARTIN等[26]率先將CPML應(yīng)用于一階Biot方程的時間域有限差分法求解;HE等[27]對二階Biot方程的有限元法求解提出了一種新的非分裂形式PML,并借助COMSOL軟件平臺討論了其吸收效果。值得關(guān)注的是,SHI等[28]基于CFS-PML提出了一種匹配Z變換(matched Z-transform)PML(MZT-PML)技術(shù),并將其應(yīng)用于彈性波數(shù)值模擬。MZT-PML繼承了CFS-PML所有優(yōu)點(diǎn),而采用Z變換技術(shù)使MZT-PML相比CPML能更直接地實現(xiàn)復(fù)頻移拉伸算子,因此MZT-PML技術(shù)可被進(jìn)一步應(yīng)用于粘滯聲波方程數(shù)值模擬[29]。

本文采用時域交錯網(wǎng)格有限差分法,將MZT-PML技術(shù)擴(kuò)展應(yīng)用于孔隙介質(zhì)彈性波傳播的數(shù)值模擬,然后給出了在Biot方程加載MZT-PML的一般推導(dǎo)過程,最后通過數(shù)值模擬算例證明了MZT-PML技術(shù)的有效性。

1 方法原理

1.1 孔隙介質(zhì)彈性波動方程

二維孔隙介質(zhì)彈性波波動方程表示如下[30]:

式中:自由指標(biāo)i和j在二維情形下分別可取x或y;啞指標(biāo)k遵守愛因斯坦求和法則;δij表示Kronecker函數(shù);μ表示剪切模量;λs=λ-α2M表示固體骨架的拉梅系數(shù),其中λ表示飽含流體骨架的拉梅系數(shù)。

1.2 匹配Z變換完全匹配層

將(1)式和(2)式寫成二階位移格式和一階速度-應(yīng)力格式,并變換到頻率域,有:

KUZUOGLU等[17]對(4)式進(jìn)行改進(jìn)并提出了如下的復(fù)頻移拉伸算子(CFS-PML邊界條件下),改善了大角度入射波的吸收效果:

式中:衰減參數(shù)dη,κη和αη均表示沿η坐標(biāo)軸的實函數(shù),具體公式見后文,dη≥0,κη≥1,αη≥0,其中η為自由指標(biāo),在二維情形下分別可取x或y。當(dāng)κη=1和αη=0時,CFS-PML退化為傳統(tǒng)的PML。將PML的復(fù)拉伸坐標(biāo)代入(3)式中的第1個方程(第2個方程的處理跟第1個方程的處理過程相同,不再贅述),可得:

(5)式為關(guān)于頻率的函數(shù),如果將(6)式和(7)式進(jìn)行傅里葉反變換到時間域,得到的波動方程將會產(chǎn)生卷積項。CPML的基本原理即為引入記憶變量,并采用遞推卷積技術(shù)來計算卷積項,從而實現(xiàn)復(fù)頻移拉伸算子的加載。我們注意到無論是時間域中的卷積或者頻率域的乘法運(yùn)算在Z域均簡化為乘法運(yùn)算,于是將(6)式和(7)式變換到Z域能相對容易地應(yīng)用CFS-PML技術(shù),這便是本文所采用的MZT-PML技術(shù)的基本思想。

MZT-PML在Biot方程中實現(xiàn)復(fù)頻移拉伸算子的基本推導(dǎo)過程如下。將(5)式的復(fù)頻移拉伸算子變換到Z域,同時將式中所有速度-應(yīng)力方程變換到Z域。首先考慮(5)式在S域的倒數(shù)形式為:

式中:s表示復(fù)頻率參數(shù)。對于任意變量P,從S變換到Z變換有如下對應(yīng)關(guān)系:(s-P)?(1-ePΔtz-1),故(8)式對應(yīng)的Z變換結(jié)果為:

式中:Δt表示時間采樣間隔。考慮到(1-z-1)/Δt為iω的Z變換,那么(6)式的Z域形式可表示為:

(11)式,(12)式和(13)式可進(jìn)一步改寫為:

其中,中間變量bη=e-(αη+dη/κη)Δt(η=x,y)。將(11)式、(12)式和(13)式代入(10)式,可得:

Ψσxy,y-ayz-1Ψσxy,y)+ρf(ΨPf,x-

axz-1ΨPf,x)+ρfKvfx

(17)

式中:aη=e-αηΔt(η=x,y);考慮到z-1表示離散時間域一個步長的延遲,那么(14)式、(15)式、(16)式和(17)式可以表示為如下的有限差分時間遞推格式:

(21)

(18)式、(19)式、(20)式和(21)式即為包含了MZT-PML的時間域有限差分遞推式。對于(3)式中的其它方程采用同樣的方法處理即可以完成MZT-PML的加載,然后對這些方程采用空間四階時間二階交錯網(wǎng)格有限差分法離散求解,則可算出全部分量的波場。

2 數(shù)值模擬及分析

為驗證MZT-PML技術(shù)在孔隙介質(zhì)彈性波數(shù)值模擬中的有效性及長時間模擬的穩(wěn)定性,設(shè)計如圖1 所示的均勻孔隙介質(zhì)模型進(jìn)行數(shù)值模擬實驗。將得到的數(shù)值結(jié)果與參考解進(jìn)行對比分析,并與采用CPML技術(shù)得到的數(shù)值解進(jìn)行比較分析,以驗證MZT-PML技術(shù)的可靠性。

圖1 二維均勻孔隙介質(zhì)模型

均勻孔隙介質(zhì)彈性波數(shù)值模擬包括兩部分,圖1中的陰影部分為MZT-PML吸收層區(qū)域,中間的空白部分為求解區(qū)域,模型為長條狀,有利于檢驗MZT-PML吸收層對近似掠射入射波的吸收效果。圓圈S為震源激發(fā)位置,坐標(biāo)為(30m,5.5m),震源為固相縱波源,主頻為80Hz,三角形R1和R2為接收點(diǎn),坐標(biāo)分別為(40m,60m)和(300m,5.5m),該均勻孔隙介質(zhì)模型的物性參數(shù)如表1所示。模型大小為310.5m×70.5m,空間離散時,x方向和y方向的相鄰網(wǎng)格點(diǎn)間距均為0.5m,全區(qū)域離散網(wǎng)格點(diǎn)為622×142個。依據(jù)時間遞推穩(wěn)定性條件[26],時間步長應(yīng)小于0.11ms,故此處設(shè)定時間步長為0.1ms,采樣時間長度為600ms,時間采樣點(diǎn)為6000個。MZT-PML吸收層的厚度為5m,包含10個有限差分單元。將震源置于頂部吸收邊界附近,接收點(diǎn)R1和R2置于底部和頂部吸收邊界附近。MZT-PML吸收層的吸收效果受衰減參數(shù)dη、κη和αη的控制。這3個衰減參數(shù)通常由以下3個公式確定[31]:

式中:η表示PML吸收層內(nèi)一點(diǎn)到交界面的距離;L表示PML吸收層的厚度。本算例中,各衰減參數(shù)取值如下:dmax=-(pd+1)cplg(Rc)/2L,amax=πf0,κmax=-(pκ+1)blg(Rc)/2L,其中,Rc=10-5,pd=pκ=2,pα=1,b=10Δh,Δh為差分網(wǎng)格點(diǎn)最大間距,cp表示最大縱波速度,本算例中為快縱波速度;f0為子波主頻。

表1 均勻孔隙介質(zhì)模型的物性參數(shù)

圖2為均勻孔隙介質(zhì)模型加載MZT-PML后數(shù)值模擬得到的40ms,100ms和200ms時刻的波場快照。圖2a和圖2b分別為40ms,100ms和200ms時刻固相和流相的x分量,圖2c和圖2d分別為40ms,100ms和200ms時刻固相和流相的y分量。可以看出,MZT-PML技術(shù)對各個角度的入射波場都有良好的吸收效果。對比圖2中固相和流相的波場快照,可發(fā)現(xiàn)慢縱波振幅比快縱波振幅大(波前超前的為快縱波,滯后的為慢縱波),流相和固相中的慢縱波相位相反,流相中的快縱波振幅相對更小,這種振幅上的差異導(dǎo)致圖2b和圖2d中并未顯示出快縱波。

圖2 固相x分量(a)、流相x分量(b)、固相y分量(c)和流相y分量(d)在40ms,100ms和200ms時刻的波場快照(紅色直線表示MZT-PML吸收層的內(nèi)邊界)

圖3和圖4分別為CPML和MZT-PML邊界條件下接收點(diǎn)R1和R2處的固相x分量記錄,其中的參考解利用擴(kuò)邊方法獲得。從圖3a可看出,CPML和MZT-PML邊界條件下得到的數(shù)值解與參考解高度吻合,且擬合誤差在10-2量級,說明CPML和MZT-PML邊界條件下近垂直方向入射波吸收效果好。圖3b 為CPML和MZT-PML邊界條件下得到的數(shù)值解和參考解的差值對比,可以看出二者僅存在微小的差異,這是CPML和MZT-PML邊界條件的離散格式差異造成的。圖4進(jìn)一步展現(xiàn)了CPML和MZT-PML邊界條件下近似掠射的入射波吸收效果,CPML和MZT-PML邊界條件下得到的數(shù)值解與參考解幾乎完全吻合,說明CPML和MZT-PML邊界條件下近似掠射的入射波均吸收效果好。

長時間數(shù)值模擬的穩(wěn)定性也是評價各種完全匹配層的重要指標(biāo)之一。圖5顯示了CPML和MZT-PML邊界條件下波場總能量隨時間的衰減情況。波場總能量的計算公式如下[26]:

式中:點(diǎn)號表示對時間的一階偏導(dǎo)。本算例中將傳播時間延長至10s,即10000倍時間步長,傳播時長分別為0.6s和10.0s時波場總能量衰減情況如圖5所示。由圖5a可見約0.4s之后能量陡降,這是因為此時直達(dá)波完全離開了求解區(qū)域;0.4s之后能量逐漸趨于穩(wěn)定,理論上此時殘留的能量全部為虛假反射能量。由圖5b可知即使是在6.0s之后,能量仍呈微弱的下降趨勢,說明了MZT-PML和CPML邊界條件均具有較好的長時間數(shù)值模擬穩(wěn)定性,經(jīng)MZT-PML和CPML邊界條件吸收后殘留的總能量基本處于同一數(shù)量級。

圖3 CPML和MZT-PML邊界條件下接收點(diǎn)R1處的固相x分量記錄a 數(shù)值解與參考解對比; b 數(shù)值解與參考解的差值對比

圖4 CPML和MZT-PML邊界條件下接收點(diǎn)R2處的固相x分量記錄a 數(shù)值解與參考解對比; b 數(shù)值解與參考解的差值對比

圖5 CPML和MZT-PML邊界條件下不同傳播時長的波場總能量衰減情況a 傳播時長0.6s; b 傳播時長10.0s

3 結(jié)論

本文詳細(xì)推導(dǎo)了非分裂MZT-PML在Biot方程中實現(xiàn)復(fù)頻移拉伸算子的一般過程,并采用時域交錯網(wǎng)格有限差分法對加載了MZT-PML的Biot方程進(jìn)行離散求解。與基于傅里葉變換和卷積算子的CPML不同的是,采用匹配Z變換技術(shù)的MZT-PML可更直接地實現(xiàn)復(fù)頻移拉伸算子。數(shù)值模擬結(jié)果表明,MZT-PML對孔隙介質(zhì)中固相和流相的各分量均具有良好的吸收效果。與CPML類似,MZT-PML也能有效吸收各個角度入射的地震波,尤其對于近似掠射的入射波的吸收效果顯著。由長時間能量衰減計算結(jié)果可知,MZT-PML在Biot方程求解中具有長時間數(shù)值模擬穩(wěn)定性。