覃發(fā)兵,高志偉,解皓楠,徐振旺

(1.長(zhǎng)江大學(xué)地球物理與石油資源學(xué)院,湖北武漢430100;2.長(zhǎng)江大學(xué)管理學(xué)院,湖北荊州434023;3.中國(guó)石油天然氣股份有限公司新疆油田分公司采油一廠,新疆克拉瑪依834000;4.中國(guó)石油天然氣股份有限公司遼河油田分公司勘探開(kāi)發(fā)研究院,遼寧盤(pán)錦124010)

高精度彈性波數(shù)值模擬是彈性波動(dòng)方程成像和反演的核心步驟之一。彈性波數(shù)值模擬時(shí),從計(jì)算可行性考慮,無(wú)界的地下介質(zhì)需加人工邊界截?cái)酁橛薪鐓^(qū)域,但在人工邊界處會(huì)出現(xiàn)虛假反射,從而嚴(yán)重影響數(shù)值模擬的精度。因此,需要在人工邊界處進(jìn)行特殊處理以消除虛假反射對(duì)數(shù)值模擬結(jié)果的影響。

根據(jù)地震波場(chǎng)數(shù)值模擬時(shí)是否在求解區(qū)域最外層附加額外的幾何空間,可將近四十年發(fā)展起來(lái)的人工邊界處理技術(shù)分為兩類:吸收邊界方法和吸收衰減層方法。吸收邊界方法是在計(jì)算區(qū)域最外層網(wǎng)格節(jié)點(diǎn)上引入額外變量壓制虛假反射。由于沒(méi)有增加計(jì)算區(qū)域的幾何空間,吸收邊界方法產(chǎn)生的額外計(jì)算量較少。經(jīng)典的吸收邊界方法為CLAYTON等[1]于1977年提出的旁軸近似法,該方法能有效吸收垂直入射的波場(chǎng)能量,但對(duì)于大角度入射能量吸收效果較差,難以滿足精度要求。針對(duì)該問(wèn)題,HIGDON[2]將旁軸近似法推廣為高階形式,使其能吸收更大范圍角度入射的波場(chǎng)能量,但高階形式的旁軸近似實(shí)現(xiàn)過(guò)程較復(fù)雜,且引入的額外計(jì)算量較大。不同于吸收邊界方法,吸收衰減層方法是在實(shí)際物理求解區(qū)域外額外設(shè)置多個(gè)衰減層以逐漸吸收傳入到衰減層中的波場(chǎng)能量。該方法顯然會(huì)引入額外的計(jì)算量,但往往能吸收較大范圍入射角的波場(chǎng)能量。1994年,BERENGER[3]在求解Maxwell方程組時(shí)提出了目前最經(jīng)典的吸收衰減層方法——完全匹配層(PML)。該方法理論上能使任意角度任意頻率成分的邊界入射波的反射系數(shù)為零,并因此而得名。1996年,CHEW等[4]首次將PML推廣到彈性波方程的求解。早期PML需要分裂波場(chǎng)來(lái)求解,這對(duì)于大規(guī)模數(shù)值模擬是不容忽視的計(jì)算量。2007年,KOMATITSCH和MARTIN在地震波場(chǎng)模擬中發(fā)展了不需要分裂波場(chǎng)的卷積完全匹配層(CPML)[5],顯著提高了PML的計(jì)算效率及對(duì)近掠射的邊界入射波的吸收效果。隨后,CPML技術(shù)被迅速推廣到聲波方程模擬及更為復(fù)雜介質(zhì)的數(shù)值模擬中[6-9]。近期,高正輝等[10]在2.5維聲波近似方程數(shù)值模擬中提出了一種改進(jìn)的PML方法,該方法可使邊界入射波在PML吸收層中進(jìn)行兩次能量衰減,從而進(jìn)一步提高了PML的吸收效果。

然而,關(guān)于PML的應(yīng)用研究大多在有限差分法數(shù)值模擬方面[10-15],這是因?yàn)镻ML在一階形式的波動(dòng)方程求解中較容易加載。有限元法數(shù)值模擬主要基于二階波動(dòng)方程,對(duì)于PML的應(yīng)用存在一定的障礙。2003年,KOMATITSCH等[16]提出了一組適用于二階波動(dòng)方程的PML方程組,并將其成功應(yīng)用于譜元法地震波場(chǎng)數(shù)值模擬。盡管如此,關(guān)于PML在有限元法數(shù)值模擬中的應(yīng)用研究仍然較少[17-18]。

本文基于LIU等[19]提出的PML在譜元法中的加載方法,推導(dǎo)了分裂式PML在二維彈性波時(shí)間域有限元數(shù)值模擬中的加載過(guò)程,說(shuō)明了理論反射系數(shù)選取對(duì)PML吸收效果的影響,討論了狄利克雷邊界在PML外邊界層和自由邊界處的正確使用方式及其對(duì)PML吸收效果的影響。

1 PML的基本原理

頻率域各向異性彈性波動(dòng)方程可表示為:

(1)

圖1 目標(biāo)區(qū)域和PML邊界層區(qū)域圖示

(2)

依據(jù)鏈?zhǔn)椒ㄓ?

(3)

?x表示對(duì)x坐標(biāo)的偏導(dǎo)(下文同)。則:

(4)

其中

(5)

PML從純數(shù)學(xué)上講,是將對(duì)實(shí)坐標(biāo)的偏導(dǎo)做如(4)式所示對(duì)復(fù)數(shù)坐標(biāo)偏導(dǎo)的復(fù)數(shù)變換。對(duì)應(yīng)地,PML的物理意義可借助平面波解來(lái)很好地理解。將(2)式代入平面波解:

(6)

2 含PML的有限元控制方程

2.1 含PML的彈性波方程

以二維各向同性完全彈性介質(zhì)為例,其本構(gòu)方程為:

(7)

(8)

將(7)式代入(8)式,有:

(9)

(10)

又知:

(11)

將(11)式代入(10)式得:

將位移場(chǎng)x分量Ux分裂為Ux=Ux1+Ux2+Ux3,并將其代入(12)式中的第一個(gè)方程,整理后有:

(13)

同理,將Uz分裂為Uz=Uz1+Uz2+Uz3,并將其代入(12)式中第二個(gè)方程,整理后有:

(14)

(13)式和(14)式中均含有高次項(xiàng),故引入以下4個(gè)中間變量:

(15)

將(15)式代入(13)式和(14)式中,有:

(16)

和

(17)

對(duì)(16)式和(17)式進(jìn)行時(shí)間傅里葉反變換,可得包含了PML的彈性波動(dòng)方程:

(18)

和

(19)

式中:uj(j=x,z)表示對(duì)應(yīng)Uj(j=x,z)的時(shí)間域位移場(chǎng)量,pxx、pxz、pzx、pzz為時(shí)間域中間變量。

2.2 彈性波方程的等效積分弱形式

設(shè)Ω為全區(qū)域,Γ為PML外邊界,nx為單位外法向量的x分量,nz為單位外法向量的z分量。采用Galerkin法,取試探函數(shù)φx和φz分別為位移x分量和z分量的變分,將(18)式和(19)式分別乘以φx和φz并求其加權(quán)積分。根據(jù)分部積分性質(zhì)有:

(20)

和

φz+C12?xuxnzφz)dΓ-

dzuz2φzdΩ

(21)

(20)式和(21)式即為包含PML的彈性波方程等效積分弱形式。

2.3 空間離散及有限元控制方程

有限元的主要思想是將全區(qū)域Ω的求解轉(zhuǎn)變?yōu)橐幌盗凶訁^(qū)域Ωe(或單元區(qū)域)上的局部求解,將這些局部解組合即得到全局解。這里將全區(qū)域分解為若干線性三角形單元,即每個(gè)單元內(nèi)有3個(gè)節(jié)點(diǎn),進(jìn)行雙線性插值,則單元內(nèi)部的位移場(chǎng)可近似為:

(22)

(23)

其中

(24)

(25)

其中

(26)

符號(hào)“?”表示兩個(gè)向量之間對(duì)應(yīng)元素相乘。綜上所述,可將(20)式和(21)式寫(xiě)成有限元控制方程形式:

(27)

和

(28)

式中:M為全局質(zhì)量矩陣;Kx1、Kx2、Kx3、Kx4、Kx5、Kz1、Kz2、Kz3、Kz4和Kz5均為全局剛度矩陣;Cx、Cz、Cxz、Cxx和Czz為全局阻尼矩陣。

3 有限元控制方程的時(shí)間離散

觀察(27)式和(28)式可見(jiàn),加載PML的彈性波有限元控制方程有如下一般形式:

采用Euler方法進(jìn)行時(shí)間離散求解類似(30)式的方程,可得到四個(gè)中間變量pxx、pxz、pzx和pzz的時(shí)間遞推格式:

(31)

其中,Δt為時(shí)間步長(zhǎng)。對(duì)于類似(29)式的方程,可用中心差分法得到:

(32)

和

(33a)

(33b)

(33c)

(34)

4 模型算例

4.1 全空間均勻介質(zhì)模型

為了說(shuō)明理論反射系數(shù)的設(shè)定對(duì)PML吸收效果的影響,進(jìn)行了二維彈性波均勻全空間數(shù)值模擬。圖2a為均勻各向同性彈性介質(zhì)模型,密度為2000kg/m3,縱波速度為1800m/s,橫波速度為1100m/s,泊松比約為0.202。震源為正z方向激發(fā)的方向力源,子波函數(shù)為雷克子波,主頻為5Hz,位于(2000m,1000m)點(diǎn)處。采用對(duì)稱結(jié)構(gòu)化線性三角形網(wǎng)格(圖2b),先將全區(qū)域用邊長(zhǎng)為10m的正方形單元剖分,然后將每個(gè)正方形單元?jiǎng)澐殖?個(gè)等腰直角三角形。在全空間條件下(即無(wú)面波情形),最短波長(zhǎng)為1100/5/2.5=88.0m,三角形最大邊長(zhǎng)為10m,則一個(gè)波長(zhǎng)內(nèi)的采樣點(diǎn)數(shù)至少為8.8。此外,三角形單元的高最小為5m,線性三角形單元CFL(courant-friedrichs-lewy)數(shù)約為0.67,時(shí)間步長(zhǎng)應(yīng)滿足穩(wěn)定性條件Δt≤0.67×5/1800=1.86ms,本算例中取時(shí)間步長(zhǎng)為1ms。

根據(jù)前人研究[5],阻尼函數(shù)dx可表示為:

(35)

圖2 均勻各向同性彈性介質(zhì)模型(a)與對(duì)稱結(jié)構(gòu)化線性三角形網(wǎng)格(b)

4.2 半空間均勻介質(zhì)模型

狄利克雷邊界條件不同于自由邊界條件,是一種剛性邊界條件,即波傳播至此會(huì)被完全反射回去。假設(shè)PML吸收層最外層區(qū)域?yàn)棣竨,那么狄利克雷邊界條件可表示為:

(36)

下面討論半空間均勻介質(zhì)情況下,在PML吸收層最外層施加狄利克雷邊界條件對(duì)PML吸收效果的影響。震源位置設(shè)為(2000m,50m),z=0處為自由表面。圖5和圖6分別為未加載狄利克雷邊界條件和不包括地表自由表面部分的PML最外吸收層上加載了狄利克雷邊界條件的波場(chǎng)快照。從圖5a和圖6a 均可看到自由表面處有明顯的面波產(chǎn)生。對(duì)比圖5b和圖6b可見(jiàn),PML最外吸收層上不加載狄利克雷邊界條件會(huì)降低PML邊界的穩(wěn)定性,產(chǎn)生可見(jiàn)的虛假波場(chǎng)(圖5b黃色虛線圈內(nèi))。因此,在半空間均勻介質(zhì)情形下,PML吸收層最外層加載狄利克雷邊界條件能增強(qiáng)PML的數(shù)值穩(wěn)定性。

圖3 不同理論反射系數(shù)1.5s時(shí)刻的總速度場(chǎng)波場(chǎng)快照對(duì)比a 理論反射系數(shù)R=10-2; b 理論反射系數(shù)R=10-3; c 理論反射系數(shù)R=10-4; d 理論反射系數(shù)R=10-5

圖4 不同理論反射系數(shù)條件下(1700m,1300m)點(diǎn)處地震記錄a 速度場(chǎng)x分量; b 速度場(chǎng)z分量

圖5 未加狄利克雷邊界條件的總速度場(chǎng)波場(chǎng)快照a 2.0s時(shí)刻; b 2.5s時(shí)刻

需要指出的是,在PML吸收層最外層加載狄利克雷邊界條件存在一個(gè)嚴(yán)重的數(shù)值陷阱。圖7為在包括地表自由表面部分的所有PML最外吸收層上加載了狄利克雷邊界條件模擬的波場(chǎng)快照?？梢钥闯?在自由表面左右兩端產(chǎn)生了非常嚴(yán)重的面波虛假反射,嚴(yán)重干擾了目標(biāo)區(qū)域內(nèi)的波場(chǎng)求解,說(shuō)明處于地表自由表面的PML吸收層最外層部分不可加載狄利克雷邊界條件。

圖6 加載了狄利克雷邊界條件的總速度場(chǎng)波場(chǎng)快照a 2.0s時(shí)刻; b 2.5s時(shí)刻

圖7 在處于自由表面上的PML吸收層最外層上加載狄利克雷邊界條件的總速度場(chǎng)波場(chǎng)快照a 2.0s時(shí)刻; b 2.5s時(shí)刻

4.3 起伏地表非均勻介質(zhì)模型

為了檢驗(yàn)PML吸收層中包含狄利克雷邊界在復(fù)雜介質(zhì)條件下的吸收效果,采用如圖8所示的起伏地表非均勻介質(zhì)模型進(jìn)行了有限元數(shù)值模擬測(cè)試。模型x方向3000m,最大深度為1000m,地表最大高程為210 m。上層介質(zhì)密度為2100kg/m3,縱波速度2800m/s,橫波速度1500m/s;下層介質(zhì)密度為2250kg/m3,縱波速度3400m/s,橫波速度1800m/s。震源函數(shù)為雷克子波,主頻20Hz,位于(1500m,100m)點(diǎn)處。模型全區(qū)域被剖分為686025個(gè)非結(jié)構(gòu)化線性三角形單元,其最大邊長(zhǎng)約為5.4m,高最小為3.1m。由于穩(wěn)定性要求,時(shí)間采樣間隔取為0.5ms。PML吸收層厚度設(shè)置為170m,對(duì)應(yīng)于一個(gè)最大主波長(zhǎng)(約20個(gè)網(wǎng)格點(diǎn))。理論反射系數(shù)取10-6,以保證復(fù)雜模型中PML吸收層能充分吸收不同方向不同波長(zhǎng)的能量。

圖8 起伏地表非均勻介質(zhì)模型

圖9為0.4s、1.0s和2.0s時(shí)刻的總速度場(chǎng)波場(chǎng)快照。其中,圖9a、圖9c和圖9e為PML吸收層最外層不加狄利克雷邊界時(shí)的模擬結(jié)果;圖9b、圖9d和圖9f為PML吸收層最外層加載了狄利克雷邊界條件時(shí)的模擬結(jié)果。從圖9可以看出,在0.4s和1.0s時(shí)刻,PML吸收層最外層加載和不加載狄利克雷邊界條件得到的波場(chǎng)信息一致;但在2.0s時(shí)刻,PML吸收層最外層不加狄利克雷邊界條件時(shí),模擬結(jié)果出現(xiàn)了不穩(wěn)定現(xiàn)象(見(jiàn)圖9e黃色虛線框部分),這進(jìn)一步說(shuō)明了在PML吸收層最外層加載狄利克雷邊界條件能增強(qiáng)PML的數(shù)值穩(wěn)定性。

圖9 不同時(shí)刻總速度場(chǎng)波場(chǎng)快照(從上到下波場(chǎng)快照時(shí)刻依次為0.4s、1.0s和2.0s)a,c,e PML吸收層最外層未加狄利克雷邊界條件; b,d,f PML 吸收層最外層加載狄利克雷邊界條件

4.4 討論

在全空間均勻介質(zhì)算例中,PML吸收層厚度一定時(shí),理論反射系數(shù)越小,PML吸收效果越好。由于PML是一種衰減性質(zhì)的邊界條件,這類加衰減層的邊界條件經(jīng)驗(yàn)上遵循衰減層越厚吸收效果越好的規(guī)律,因此推斷理論反射系數(shù)一定時(shí),PML厚度越大,PML吸收效果越好。總之,理論反射系數(shù)越小,PML吸收層越厚,PML吸收效果越好。針對(duì)不同復(fù)雜程度的速度模型,理論反射系數(shù)和PML吸收層厚度的選擇可根據(jù)數(shù)值試驗(yàn)確定。觀察分析PML吸收層最外層加載狄利克雷邊界條件對(duì)PML吸收效果的影響,我們認(rèn)為其與狄利克雷邊界本身性質(zhì)有關(guān)。狄利克雷邊界條件是一種剛性邊界條件,在有限區(qū)域內(nèi)向外傳播的波場(chǎng)經(jīng)過(guò)PML衰減層逐步衰減,傳播到PML最外層時(shí)遇到狄利克雷邊界,未被PML層完全吸收的波場(chǎng)會(huì)被完全反射回PML衰減層,進(jìn)行第二次指數(shù)衰減。因此,理論分析認(rèn)為,在PML最外層添加狄利克雷邊界條件可使傳播到人工邊界處的波場(chǎng)經(jīng)歷兩倍厚度的PML吸收層的指數(shù)衰減,減小了時(shí)間遞推過(guò)程中計(jì)算誤差的傳遞,從而增強(qiáng)了PML的數(shù)值穩(wěn)定性。但在自由表面PML吸收層區(qū)域加載狄利克雷邊界條件時(shí),產(chǎn)生了嚴(yán)重的面波虛假反射,這是因?yàn)樽杂杀砻娴腞ayleigh面波沿自由表面?zhèn)鞑ブ罰ML吸收層時(shí),遇到狄利克雷邊界產(chǎn)生了全反射,繼續(xù)傳播至PML外邊界自由表面角落時(shí),如同產(chǎn)生了二次震源,最終造成強(qiáng)虛假反射。

5 結(jié)束語(yǔ)

本文推導(dǎo)和實(shí)現(xiàn)了PML吸收邊界條件在二階彈性波方程有限元數(shù)值解法中的變分形式,通過(guò)全空間均勻介質(zhì)數(shù)值實(shí)驗(yàn),測(cè)試了理論反射系數(shù)對(duì)PML吸收效果的影響,發(fā)現(xiàn)當(dāng)PML吸收層厚度約為半個(gè)最大主波長(zhǎng)、理論反射系數(shù)設(shè)置小于(或等于)10-5時(shí),PML吸收效果達(dá)到最優(yōu)。

半空間均勻介質(zhì)數(shù)值實(shí)驗(yàn)結(jié)果表明,在PML吸收層最外層加載狄利克雷邊界條件可提高PML的數(shù)值穩(wěn)定性,但是處于自由表面的PML吸收層最外層部分不可加載狄利克雷邊界條件,否則會(huì)產(chǎn)生嚴(yán)重的虛假反射。

起伏地表復(fù)雜介質(zhì)數(shù)值模擬結(jié)果進(jìn)一步表明,在PML吸收層最外層加載狄利克雷邊界條件對(duì)PML數(shù)值穩(wěn)定性改善效果顯著,說(shuō)明狄利克雷邊界條件的引入可作為提高PML數(shù)值穩(wěn)定性的有效手段之一。