徐振旺

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

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

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

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

衰減吸收層技術(shù)是在人工邊界附近設(shè)置一定厚度的吸收層,使波在傳播到吸收層時(shí),按指數(shù)規(guī)律逐步衰減。目前最為成功的衰減吸收層技術(shù)是BéRENGER[6]在求解Maxwell方程時(shí)提出的完全匹配層(perfectly matched Layer,PML)技術(shù)。理論上PML技術(shù)的吸收層內(nèi)無波阻抗差異,因此反射系數(shù)為零。相較于透射邊界條件,PML技術(shù)適用于較大范圍入射角的地震波吸收,同時(shí)不影響時(shí)間遞推穩(wěn)定性條件。CHEW等[7]率先驗(yàn)證了PML技術(shù)在地震波數(shù)值模擬中的有效性;王守東[8]研究了PML技術(shù)在聲波方程中的應(yīng)用;陳可洋[9-10]提出了基于聲波方程數(shù)值模擬的PML改進(jìn)算法;高剛等[11]對(duì)PML的衰減因子進(jìn)行了詳細(xì)討論,并對(duì)PML衰減因子的設(shè)定進(jìn)行了理論研究;CHEN等[12]將PML技術(shù)應(yīng)用于任意角度標(biāo)量波動(dòng)方程的數(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)行波場(chǎng)分離,可節(jié)省計(jì)算機(jī)內(nèi)存。目前,非分裂形式的CFS-PML實(shí)現(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方程的時(shí)間域有限差分法求解;HE等[27]對(duì)二階Biot方程的有限元法求解提出了一種新的非分裂形式PML,并借助COMSOL軟件平臺(tái)討論了其吸收效果。值得關(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能更直接地實(shí)現(xiàn)復(fù)頻移拉伸算子,因此MZT-PML技術(shù)可被進(jìn)一步應(yīng)用于粘滯聲波方程數(shù)值模擬[29]。

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

1 方法原理

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

二維孔隙介質(zhì)彈性波波動(dòng)方程表示如下[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]對(duì)(4)式進(jìn)行改進(jìn)并提出了如下的復(fù)頻移拉伸算子(CFS-PML邊界條件下),改善了大角度入射波的吸收效果:

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

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

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

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

式中:Δt表示時(shí)間采樣間隔?？紤]到(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表示離散時(shí)間域一個(gè)步長的延遲,那么(14)式、(15)式、(16)式和(17)式可以表示為如下的有限差分時(shí)間遞推格式:

(21)

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

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

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

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

均勻孔隙介質(zhì)彈性波數(shù)值模擬包括兩部分,圖1中的陰影部分為MZT-PML吸收層區(qū)域,中間的空白部分為求解區(qū)域,模型為長條狀,有利于檢驗(yàn)MZT-PML吸收層對(duì)近似掠射入射波的吸收效果。圓圈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,空間離散時(shí),x方向和y方向的相鄰網(wǎng)格點(diǎn)間距均為0.5m,全區(qū)域離散網(wǎng)格點(diǎn)為622×142個(gè)。依據(jù)時(shí)間遞推穩(wěn)定性條件[26],時(shí)間步長應(yīng)小于0.11ms,故此處設(shè)定時(shí)間步長為0.1ms,采樣時(shí)間長度為600ms,時(shí)間采樣點(diǎn)為6000個(gè)。MZT-PML吸收層的厚度為5m,包含10個(gè)有限差分單元。將震源置于頂部吸收邊界附近,接收點(diǎn)R1和R2置于底部和頂部吸收邊界附近。MZT-PML吸收層的吸收效果受衰減參數(shù)dη、κη和αη的控制。這3個(gè)衰減參數(shù)通常由以下3個(gè)公式確定[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時(shí)刻的波場(chǎng)快照。圖2a和圖2b分別為40ms,100ms和200ms時(shí)刻固相和流相的x分量,圖2c和圖2d分別為40ms,100ms和200ms時(shí)刻固相和流相的y分量?？梢钥闯?MZT-PML技術(shù)對(duì)各個(gè)角度的入射波場(chǎng)都有良好的吸收效果。對(duì)比圖2中固相和流相的波場(chǎng)快照,可發(fā)現(xiàn)慢縱波振幅比快縱波振幅大(波前超前的為快縱波,滯后的為慢縱波),流相和固相中的慢縱波相位相反,流相中的快縱波振幅相對(duì)更小,這種振幅上的差異導(dǎo)致圖2b和圖2d中并未顯示出快縱波。

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

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

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

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

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

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

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

3 結(jié)論

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