匹配Z變換完全匹配層在孔隙介質彈性波數值模擬中的應用

徐振旺

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

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

地震波數值模擬是描述和認識地震波傳播規律的主要手段之一?？紤]數值模擬的計算可行性,將無限的介質區域人為地截斷為有限區域,截斷的邊界稱為人工邊界。然而,人工邊界的引入會導致強烈的邊界虛假反射,影響數值模擬的真實性和準確性。為解決這一問題,人們研究了一系列的方法技術,主要包括兩類:透射邊界條件和衰減吸收層技術。

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

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

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

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

1 方法原理

1.1 孔隙介質彈性波動方程

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

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

1.2 匹配Z變換完全匹配層

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

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

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

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

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

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

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

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

其中,中間變量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的加載,然后對這些方程采用空間四階時間二階交錯網格有限差分法離散求解,則可算出全部分量的波場。

2 數值模擬及分析

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

圖1 二維均勻孔隙介質模型

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

式中:η表示PML吸收層內一點到交界面的距離;L表示PML吸收層的厚度。本算例中,各衰減參數取值如下: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為差分網格點最大間距,cp表示最大縱波速度,本算例中為快縱波速度;f0為子波主頻。

表1 均勻孔隙介質模型的物性參數

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

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

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

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

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

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

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

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

3 結論

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