基于改進的EAL邊界條件的地震波場數值模擬

張平民, 姚剛*

1 中國石油大學(北京)油氣資源與工程全國重點實驗室, 北京 102249 2 中國石油大學(北京)非常規油氣科學與技術研究院, 北京 102249

0 引言

隨著計算機算力的提高和數值方法的發展,波動方程數值求解方法被廣泛用于復雜介質地震波場正演模擬.這也促進了高精度反演和成像方法的實際工業應用,如全波形反演(FWI)(Virieux and Operto, 2009; Warner et al., 2013; da Silva et al., 2016; Yao and Wu, 2017; da Silva and Yao, 2018; Chen et al., 2022a)和最小二乘逆時偏移(LSRTM)(Nemeth et al., 1999; Kühl and Sacchi, 2003; Wu D et al., 2016; Yao et al., 2022; Wu B et al., 2022).波動方程數值求解一直是波動方程正演模擬(Virieux et al., 2011; Tohti et al., 2021)、逆時偏移(Li and Qu, 2022)和全波形反演(Yao et al., 2020; Chen et al., 2022b)的核心問題.波動方程數值解法主要包括有限差分(Alford et al., 1974; Virieux, 1984, 1986; Levander, 1988)、有限元(Cohen and Fauqueux, 2005; Mulder, 2001)、偽譜法(Kosloff and Baysal, 1982)和譜元法(Komatitsch and Vilotte, 1998; Che et al., 2010).其中,有限差分法因其占用內存小,計算效率高和易于實現等優點而被廣泛應用(劉洋, 2014; Liu et al., 2015; Wang et al., 2019; Ren et al., 2022;方修政等, 2023).

在波動方程數值求解的過程中,由于計算區域的截斷,會在計算區域邊界產生強邊界反射,從而很大程度上降低了正演模擬、全波形反演和逆時偏移的精度(Etgen et al., 2009; Wang et al., 2017; Zou et al., 2020; Chen et al., 2022a,b).為了克服這一問題,目前普遍采用吸收邊界條件壓制人工截斷邊界造成的反射.這類方法通常是在計算區域外增加若干層網格對波場進行吸收衰減,從而達到消除人工邊界虛假反射的目的.目前,壓制人工邊界虛假反射的主要方法有:Cerjan吸收邊界條件(Cerjan et al., 1985)、單程波吸收邊界條件(Clayton and Engquist, 1977; Higdon, 1991)、混合吸收邊界條件(Liu and Sen, 2010, 2011, 2012, 2018)和PML(Perfectly Matched Layer)吸收邊界條件(Collino and Tsogka, 2001).

在上述吸收邊界條件中,PML吸收邊界條件效果最佳.該方法最早應用于電磁學領域(Berenger, 1994),由于其具有吸收效果良好和適用范圍廣的特點,故而在地震波場數值模擬中廣泛應用于壓制人工邊界反射(Komatitsch and Tromp, 2003; Hu et al, 2007; 杜啟振等, 2011; Zhao and Shi, 2013;劉延利等, 2021).PML吸收邊界條件本質上是對波動方程做復坐標拉伸變換,即在吸收層將波動方程的實坐標拉伸變換為復數坐標(Chew and Weedon, 1994).經過復拉伸變換后的波動方程存在卷積項,難以在時間域進行數值求解.常規的PML吸收邊界條件采用分裂波場的方法來避免這一問題(Collino and Tsogka, 2001),但同時增加了實現難度、內存消耗和計算量.

Kuzuoglu和Mittra(1996)在PML吸收邊界條件的復拉伸變換中引入復頻移因子和收縮因子形成了C-PML吸收邊界條件.經過前人的研究,C-PML吸收邊界條件有效的克服了地震波進入吸收層時會產生的不均勻波,以及當地震波以大角度入射進吸收層時吸收效果欠佳的問題(Martin et al., 2008; Fang et al., 2022).由于C-PML吸收邊界條件的復坐標拉伸變換形式的復雜性,無法采用原有的分裂波場形式在時間域數值求解.因此,前人做了大量研究分別提出了輔助微分方程法(Martin et al., 2010; Zhang and Shen, 2010)和遞歸積分法(Drossaert and Giannopoulos, 2007; Martin et al., 2008),避免了在時間域進行復雜的卷積計算.然而,這兩種方法都需要引入輔助變量參與計算.此外,Ma等(2018, 2019)結合辛算法提出了一種較為有效且適用于二階波動方程的非分裂C-PML組合算法.

Komatitsch和Martin(2007)利用數值模擬驗證,傳統PML吸收邊界條件和C-PML吸收邊界條件在對某些各向異性介質進行波場模擬中會產生計算不穩定的現象.Zeng等(2011)利用Routh-Horwitz行列式對各向異性介質中PML和C-PML吸收邊界條件的穩定性問題進行深入研究.羅玉欽和劉財(2020)根據方程組波解、特征值靈敏度和模擬結果,分析了C-PML吸收邊界條件的各個因子對計算穩定性和吸收效果的影響.Meza-Fajardo和Papageorgiou(2008)通過嚴格的數學證明,驗證了傳統的PML吸收邊界條件在各向異性介質中都不滿足嚴格的漸近穩定性.Yang等(2003)開創性地提出了n次吸收邊界條件,并且數值實驗結果表明n次吸收邊界條件結合緊致有限差分可以有效地消除各向異性介質波場模擬過程中的邊界數值反射.因此,適用于各向異性介質情況下的吸收邊界條件仍具有很大的研究潛力.

雖然PML類吸收邊界條件因其能夠幾乎完全吸收人工邊界反射波的能量而被廣泛應用,但是實現難度、內存消耗和計算量較大且仍不能穩定地模擬某些各向異性介質中的波場.混合吸收邊界在各向同性介質波場模擬中展現了其良好的吸收效果,但該方法需要根據所采用的波動方程推導相應的單程波動方程,且該方法在各向異性介質波場模擬中的吸收效果和穩定性需進一步研究.

Yao等(2018)通過簡化PML吸收邊界條件中的衰減阻尼項,提出一種更為靈活的吸收邊界條件,即EAL(Effective Absorbing Layer)邊界條件.在PML吸收邊界條件的基礎上,對所有的波場分量采用相同的衰減阻尼項來實現簡化,形成一個新的更簡單的方程系統.采用相同的衰減阻尼項可以將PML吸收邊界條件中的復坐標拉伸變換合并后轉移至波動方程的時間項上,所以該方法不存在復雜的卷積項.該方程系統的解在時間上具有指數衰減的特點,并且適用于一階方程系統和二階方程系統.雖然EAL邊界條件比PML吸收邊界條件的吸收效果略低,但在內存消耗和計算效率方面得到顯著的改善.通過采用特定的反演策略,EAL邊界條件在聲波和彈性波全波形反演的地震波模擬中表現良好.

雖然EAL邊界條件在保證吸收效果的同時兼顧了內存消耗和計算量,但是相對于其他吸收邊界該方法仍存在不足.一方面,該方法的吸收效果略低,特別是模型底部邊界的虛假反射吸收不充分.由于模型底部的速度較高,因此地震波在底部吸收層的傳播時間較短不能得到充分吸收.另一方面,雖然在各向同性介質情況下EAL邊界條件的穩定性已經得到嚴格的數學證明(Yao et al., 2018),但是沒有明確的數值實驗驗證EAL吸收邊界條件在各向異性介質波場模擬中是否穩定.

本文在常規的EAL邊界條件的基礎上,提出了改進方法.通過引入減速阻尼項并應用高斯型空間加權函數構造衰減阻尼項和減速阻尼項形成了MEAL(Modified Effective Absorbing Layer)邊界條件.相比于傳統的指數型空間加權函數,高斯型空間加權函數能夠保證阻尼函數在平滑遞增到外邊界達到相應最大值的同時兼顧了吸收層內外邊界吸收效果.本文通過實驗表明,在相同邊界厚度的情況下,MEAL邊界條件具有更好的吸收效果且適用于復雜構造模型和各向異性模型.特別地,各向異性模型正演模擬結果表明,EAL邊界條件和MEAL邊界條件具有良好的計算穩定性,同時MEAL邊界條件具有更好的吸收效果.

1 方法原理

EAL邊界條件是一種有效的邊界反射吸收衰減方法,與其他吸收類邊界條件類似,其基本思想是在計算區域的四周添加一定厚度的吸收層對進入其中的波以及人工邊界產生的反射波進行吸收.如何利用有限的吸收層厚度達到有效壓制邊界反射是該方法的核心所在.

1.1 EAL邊界條件

EAL技術本質上是將波動方程在EAL層內引入衰減項,并通過控制其衰減阻尼的空間分布,保證其在不影響計算區域模擬結果的同時確保波在吸收層中得到充分的吸收.以二維一階彈性波方程組為例,加載EAL邊界條件的一階彈性波方程組可表示為:

(1)

其中,vi為質點振動速度分量,τij為應力,D為EAL邊界條件的衰減阻尼項,ρ為密度.在各向同性情況下, 介質的彈性常數c11=c33=λ+2μ,c13=λ,c44=μ.其中,λ和μ為拉梅常數,其與縱橫波速度和密度的關系為:

(2)

衰減阻尼D可以由兩個簡單的一維衰減函數組合而成:

(3)

其中,一維衰減函數di(i=x,z)的定義方式為:

(4)

(4)式中:VP是縱波速度,R是邊界理論反射系數,對10層,20層和40層吸收層,R分別取0.001,0.0001,0.00001(Collino and Tsogka, 2001; Zhang and Shen, 2010),l是EAL區域內部計算點到EAL內邊界的距離,L是EAL層厚度.當l>0即位于吸收層時,D>0,衰減項對吸收層的波場起到吸收衰減作用.當l≤0即位于計算區域時,D≡0,衰減項對計算區域的波場不起任何作用.

根據上面的論述,在實施EAL邊界條件時,不需要對波動方程進行分裂或者引入輔助微分方程,并且也不需要改變波動方程的差分格式,整個計算區域可以采用統一的差分格式.因此,相比于PML吸收邊界條件,EAL邊界條件易于實現且計算量和內存消耗更小.這些特征使得EAL邊界條件適用于高性能并行計算(Yao et al., 2018).

1.2 MEAL邊界條件

EAL邊界條件通過引入衰減區實現壓制邊界反射,衰減區與中心研究區域相鄰的邊界稱為內邊界,衰減區的最外層邊界稱為外邊界.衰減型邊界條件的吸收效果取決于內外邊界的殘余反射能量大小.因此,衰減區的衰減函數是EAL邊界條件吸收效果的重要影響因素.由式(4)可知衰減區的衰減阻尼可以看成最大衰減阻尼d0與值域為0～1的空間加權函數的乘積:

di(l)=d0w(l),

(5)

由式(5)可知,空間加權函數w(l)的數值大小決定了波場在衰減區域的衰減程度.由于內邊界更加靠近內部計算區域,波場能量相對較強,如果空間加權函數在靠近內邊界區域變化較大,則會使得波場在內邊界區域衰減過于劇烈,從而產生較強的數值反射,即內邊界反射.但另一方面,如果空間加權函數在靠近內邊界區域變化過于緩慢,則會導致靠近外邊界區域的空間加權函數變化過于劇烈,會在外邊界產生較強的數值反射,即外邊界反射.

EAL原有的空間加權函數屬于指數型函數,需要更多的層數才能保證空間加權函數在內邊界處變化足夠緩慢,使得內邊界區域的波場與中心區域波場更加耦合,過渡更自然,有效壓制內邊界殘余反射.因此,有必要對原有的空間加權函數進行改進,兼顧內外邊界的綜合邊界反射壓制效果.非線性空間加權函數已經被證明能夠有效提升混合吸收邊界條件的吸收效果 (Liu and Sen, 2018; Xie et al., 2020).王紹文等(2021)在此基礎上提出了一種基于高斯型加權函數的混合吸收邊界條件方法,并通過數值實驗驗證其能更好地兼顧內外邊界的吸收效果.本文利用高斯型空間加權函數,對EAL衰減阻尼進行改進,提升EAL邊界條件的吸收效果,其數學表達式為:

(6)

式中,α為高斯型空間加權函數的階數,控制內外邊界處的變化率.

如圖1所示,隨著α的增大,高斯型空間加權函數在內邊界處的變化更加緩慢,而在外邊界處的變化更加劇烈.因此,可以選擇合適的α使得空間加權函數在內邊界處變化比指數型空間加權函數更加緩慢,保證內邊界處波場與中心區域波場的耦合性,同時具有更優的內外邊界反射綜合壓制效果.我們通過大量實驗表明α的最佳取值范圍為2～2.5.

圖1 指數型和高斯型空間加權函數變化曲線

目前,關于衰減型邊界條件的應用存在著一個誤區:在應用吸收邊界時,希望波傳播至外邊界時波場值衰減為0.如果要滿足這一需求,必須增加吸收衰減層的厚度,從而導致計算量的增加.實際上波傳播到衰減層外邊界產生反射后會再次經過衰減層傳播回計算區域,這將會經歷二次衰減.因此,只需要保證波從外邊界反射傳回內邊界時波場能量足夠小即可以保證波場模擬的效果足夠準確.

EAL邊界條件的波場衰減函數D是一個與空間位置相關而與時間不相關的函數.當模型邊界速度較高時,地震波會快速地從吸收層的內邊界傳播至外邊界.由于衰減函數D的值從外邊界至內邊界是遞增的,因此當模型邊界的速度值較高時,吸收層內的波場仍會出現衰減過于劇烈的現象,從而產生殘余邊界反射.所以,EAL的衰減項對波場衰減的程度還與波在吸收層中傳播的時間有關:波在吸收層中的傳播時間越長,波場能量被吸收層吸收得越充分.因此,為了進一步改善EAL邊界條件的吸收效果,本文在原有EAL邊界條件的基礎上增加一個減速阻尼項β,來延長波在高速吸收層傳播的時間.加載減速阻尼項的EAL方程組可表示為:

(7a)

(7b)

其中,β為大于等于1的阻尼參數,當β=1時,它不會對波場起任何特殊作用,當β>1時,會減小波的傳播速度.β可以由兩個簡單的一維函數組合而成:

(8)

參考C-PML(Zhang and Shen, 2010)吸收邊界條件中的收縮因子的定義方式,βi(i=x,z)可被定義為:

(9)

式中,β0為最大減速阻尼,VS為橫波速度,Δh為網格間距,f0是震源子波主頻.

由式(9)可知,空間加權函數w(l)的變化程度決定了波場在衰減區域的傳播速度的變化程度.如果空間加權函數在靠近內邊界區域變化較大,則會增加波在內邊界區域的傳播速度差異,從而產生較強的內邊界反射.另一方面,雖然外邊界區域由于衰減阻尼的作用波場能量相對較小,但是如果外邊界區域的空間加權函數變化過于劇烈,也會產生較強的外邊界反射.

因此,減速阻尼與衰減阻尼均需要一個空間加權函數,使得阻尼函數在平滑遞增到外邊界達到相應最大值的同時兼顧吸收層內外邊界吸收效果.本文仍采用高斯型空間加權函數wG(l)構造減速阻尼函數.

采用上述改進方案形成了本文提出的改進的EAL邊界條件,即MEAL邊界條件.方程(7)即為加載MEAL邊界條件的一階彈性波方程組.

以質點振動速度的水平分量方程為例,加載EAL和MEAL邊界條件的交錯網格差分格式分別為:

(10)

(12)

其中,i、j和n分別為x、z和t軸的樣點序號,Δx、Δz和Δt分別為x、z和t軸的離散步長,Lx和Lz分別為沿x和z方向的空間導數的差分算子,cm為空間導數的差分系數,M為差分階數的1/2.

由式(10)和式(11)可知,EAL和MEAL邊界條件的差分格式基本一致,不同之處是差分方程的系數.這些系數不會在波場差分迭代更新的過程中改變且可以在差分迭代更新波場前計算并保存在計算機內存中.因此,實現MEAL邊界條件的計算量和內存消耗與EAL邊界條件相同.

2 模型試算

為了驗證MEAL邊界條件的吸收效果,下面對均勻模型、Foothills模型、各向異性模型進行了正演模擬,并將正演模擬結果進行對比分析.

2.1 均勻模型

為驗證本文提出的MEAL邊界條件的吸收效果,本文設計VP=6000 m·s-1、VS=3464 m·s-1、ρ=2728 kg·m-3的均勻介質模型.縱橫向空間網格間距為8 m,時間步長為0.5 ms,采用時間二階空間八階精度差分格式.震源位于模型中心,采用主頻為25 Hz的Ricker子波加載于質點振動速度的垂直分量vz上,模擬時長5 s.

本文在常規的EAL邊界條件的基礎上,分別采用指數型與不同階數的高斯型空間加權函數進行模擬.其中,計算20層和30層吸收邊界條件中衰減阻尼的理論反射系數R分別取0.0001和0.000026.經過正演模擬得到了應用20層和30層吸收邊界條件的650 ms時刻的vz分量波場快照,分別如圖2和圖3所示.由圖2和圖3可知,相比于20層EAL邊界條件,30層EAL邊界條件的吸收效果明顯優于前者;不論是20層還是30層吸收邊界,采用指數型空間加權函數的常規EAL邊界條件的模擬結果存在較強的內邊界殘余反射,如圖中黑色箭頭所示;相比于指數型空間加權函數,隨著階數的提高,基于高斯型空間加權函數的EAL邊界條件的模擬結果的內邊界的殘余反射能量逐漸降低;相反地,外邊界的殘余反射能量逐漸增強,如圖中白色箭頭所示.模擬結果表明,通過調節高斯型空間加權函數的階數,可以控制吸收邊界在內外邊界區域的吸收效果;20層和30層吸收邊界分別采用α=2和α=2.5的高斯型空間加權函數更能平衡內外邊界的反射能量.

圖2 應用20層指數型與不同階數高斯型空間加權函數的吸收邊界在 650 ms時刻vz分量波場模擬結果對比

圖3 應用30層指數型與不同階數高斯型空間加權函數的吸收邊界在650 ms時刻vz分量波場模擬結果對比

由于吸收邊界的殘余反射能量遠小于研究區域的有效彈性波場能量,圖2和圖3所示的模擬結果對比并不是很明顯.為了定量地驗證高斯型空間加權函數的吸收效果,本文在模擬過程中計算了應用高斯型空間加權函數與常規指數型空間加權函數EAL邊界條件模擬的每一時刻波場的能量以及兩者間的相對能量差:

(13)

其中,Er為相對能量差,Ee為應用常規的指數型空間加權函數模擬的波場能量,EG為應用高斯型空間加權函數模擬的波場能量.彈性波場能量可被定義為:

(14)

其中,τij為應力,eij為應變,可通過廣義胡克定律換算得到,s為波場模擬的區域.由式(13)可知,波場的相對能量差Er越大,則模擬結果的吸收邊界的殘余反射能量相對越小,即吸收效果相對越好.

圖4為20層和30層應用不同階數的高斯型和傳統的指數型空間加權函數邊界條件模擬的波場能量之間的相對能量差.時間約為750 ms之前,由于震源的激發,波場被注入能量,部分波傳入吸收邊界,能量被吸收層吸收,EAL和MEAL情況下的能量變化基本相同.大約在750 ms之后,研究區域只有吸收層內外邊界的殘余反射波在模型中傳播.值得注意的是,大約在4 s之后,波場相對能量差逐漸減小.這是由于隨著時間的推進其殘余反射能量反復被吸收邊界吸收逐漸趨近于0,使得在4 s之后波場之間的能量差異逐漸減小.由圖4可知,當α為2和2.5時,不論是20層還是30層吸收邊界,吸收層內外邊界的殘余反射能量都低于常規的EAL邊界條件,壓制效果最高提升了約30%～40%;然而,當α=3時,吸收層外邊界區域的空間加權函數變化劇烈,尤其是采用20層吸收邊界情況下產生了較強的外邊界殘余反射,反而降低了吸收效果.雖然,吸收層的剩余反射經過傳播會再次進入吸收層得到吸收,波場能量逐漸接近于α=2.5時的波場能量,但其綜合吸收效果仍然不理想.根據該實驗可知,當α的取值在2～2.5之間時,高斯型EAL邊界條件能夠兼顧對內外邊界的綜合邊界反射的壓制效果.其中,20層和30層吸收邊界分別采用α=2和α=2.5的高斯型空間加權函數時,EAL邊界條件的吸收效果表現最為穩健.

圖4 應用不同階數的高斯型空間加權函數吸收邊界的波場相對能量差對比

為了驗證減速阻尼因子β對邊界吸收的改善作用,本文選取不同階數高斯型空間加權函數的MEAL邊界條件進行模擬.其中,計算20層和30層吸收邊界條件中衰減阻尼的理論反射系數R選取與上述實驗相同的系數,最大減速阻尼β0均取3.464.通過模擬得到應用20層和30層吸收邊界650 ms時刻的波場,分別如圖5和圖6所示.如圖中箭頭所示位置,當衰減函數D和減速阻尼因子β采用相同的空間加權函數時,應用減速阻尼因子后吸收層內外邊界的殘余反射能量明顯降低,且外邊界反射波前面的到達位置明顯滯后.這是由于吸收層引入了減速阻尼因子,使得波傳入吸收層后傳播得更慢,這也會保證波傳播到衰減層外邊界產生反射后會再次經過衰減層傳播回計算區域所需的時間更長,進而得到更好的吸收.為了更直觀的對比兩種吸收邊界條件的殘余邊界反射能量,本文給出了傳播時間相對更長的4500 ms時刻波場快照,分別如圖7和圖8所示.從模擬結果可知,引入減速阻尼因子后計算區域的剩余邊界反射能量明顯減小.因此,減速阻尼因子能夠有效地改善截斷邊界虛假反射波的壓制效果.

圖5 應用20層不同階數高斯型空間加權函數的EAL和MEAL邊界條件在650 ms時刻模擬結果對比

圖6 應用30層不同階數高斯型空間加權函數的EAL和MEAL邊界條件在650 ms時刻模擬結果對比

圖7 應用20層不同階數高斯型空間加權函數的EAL和MEAL邊界條件在4500 ms時刻模擬結果對比

圖8 應用30層不同階數高斯型空間加權函數的EAL和MEAL邊界條件在4500 ms時刻模擬結果對比

同樣地,本文計算了每一時刻應用相同階數的高斯空間加權函數的MEAL邊界條件和EAL邊界條件的波場能量與相對能量差,如圖9所示.可以清晰地看到,應用MEAL邊界條件后,吸收層對殘余反射能量的吸收效果得到顯著提升.引入減速阻尼因子后,吸收效果提高了約40%～70%.與圖4所示實驗結果相同,20層和30層吸收邊界分別采用α=2和α=2.5的高斯型空間加權函數時綜合吸收效果最為均衡.

圖9 應用不同階數的高斯型空間加權函數MEAL邊界條件的波場相對能量差對比

2.2 Foothills模型

Foothills模型描述的是加拿大落基山脈的一個橫截面模型.圖10顯示了其縱波速度模型.橫波速度由泊松體的縱橫波速比計算得出.密度由Gardner關系計算得出.該模型構造十分復雜,可以通過該模型的正演模擬效果來驗證MEAL吸收邊界條件是否具有實際應用價值.本文采用10 m的縱橫向空間網格間距,0.5 ms的時間步長和時間二階空間八階精度差分格式.震源位于模型上邊界中心,采用主頻為15 Hz的Ricker子波加載于質點振動速度的垂直分量vz上.分別選取20層和30層EAL與MEAL邊界條件進行模擬.其中,計算20層和30層吸收邊界條件中衰減阻尼的理論反射系數R分別取0.0001和0.000026,計算MEAL邊界條件中減速阻尼的最大減速阻尼β0均取4.6352.圖11—14為1.5 s時刻的波場快照.

圖10 Foothills模型

圖11 應用20層吸收邊界條件的vx分量波場吸收效果

圖12 應用20層吸收邊界條件的vz分量波場吸收效果

圖13 應用30層吸收邊界條件的vx分量波場吸收效果

圖14 應用30層吸收邊界條件的vz分量波場吸收效果

如圖11—14中箭頭所示,當邊界區域速度場相對較高時,地震波在吸收層的傳播時間較短,且指數型空間加權函數在內邊界區域變化較快,因此當波由研究區域傳播進吸收層時波場能量在空間和時間上都過于劇烈地衰減,產生了較強內外邊界殘余反射能量.

模擬結果表明,20層和30層EAL邊界條件均不能完全吸收模型底部高速區域的邊界反射,其中20層EAL邊界條件的邊界殘余反射能量明顯更強;MEAL邊界條件能夠很好的避免在模型底部高速邊界區域產生內邊界的殘余反射,從而得到清晰可靠的波場模擬結果.特別地,應用30層MEAL邊界條件后,吸收邊界的殘余反射波能量遠遠小于研究區域內的有效波,邊界反射被有效吸收.

為了進一步觀察地震波到達左右邊界時的虛假反射波壓制情況,本文將應用20層與30層的EAL與MEAL邊界條件模擬獲得的地表位置地震記錄進行對比,分別如圖15和圖16所示.模擬獲得的地震記錄表明,EAL邊界條件對左右邊界產生的虛假反射具有一定的壓制效果,但是仍然存在著一些虛假反射沒有得到充分吸收,如圖15a和圖16a中的箭頭所示.可以看出這些殘余虛假反射與其同時被檢波點接收的來自模型深部的反射波能量相近,計算區域的波場模擬精度受到明顯的影響.相反地,MEAL邊界條件能夠有效壓制這些虛假反射,左右邊界的殘余虛假反射能量明顯降低,如圖15b和圖16b所示.特別地,當采用30層MEAL邊界條件時,左右邊界虛假反射得到了有效壓制.

圖15 應用20層吸收邊界條件的vz分量地震記錄

圖16 應用30層吸收邊界條件的vz分量地震記錄

經過對該模型的測試可知,改進后的EAL邊界條件即使在復雜構造區域仍能有效的壓制邊界反射,得到較好的有效波場模擬結果.

2.3 各向異性模型

前人通過嚴格的數學推導和數值實驗已經驗證將PML邊界條件應用于各向異性介質的波場模擬并不能保證數值穩定(Meza-Fajardo and Papageorgiou, 2008).為了驗證MEAL邊界條件在模型各向異性情況下的正確性和有效性,本文選擇Komatitsch和Martin(2007)提出的磷灰石和鋅晶體兩種各向異性介質模型,模型的彈性參數如表1所示.這兩種介質模型已經被證明了PML邊界條件由于存在指數增長解,故而PML邊界條件應用于這兩個模型進行模擬時會造成數值不穩定現象.

表1 各向異性模型彈性參數

本文采用相同的模擬參數分別應用EAL和MEAL邊界條件對磷灰石模型和鋅晶體模型進行模擬.縱橫向空間網格間距為0.0625 cm,時間步長為0.2 μs,采用時間二階空間八階精度差分格式.磷灰石模型和鋅晶體模型分別采用主頻為300 kHz和170 kHz的Ricker子波,震源均位于模型中心且加載于速度垂直分量vz上,模擬時長400 μs,共計40000個時間步長,物理區域四周添加30層吸收層.該組模型的模擬實驗中,計算EAL和MEAL吸收邊界條件中衰減阻尼的理論反射系數R均取0.000026,計算MEAL邊界條件中減速阻尼的最大減速阻尼β0分別取4.853和4.445.圖17和圖18分別為磷灰石模型和鋅晶體模型在不同時刻模擬效果對比.

圖17 磷灰石模型應用30層吸收邊界的vz分量波場吸收效果對比

圖18 鋅晶體模型應用30層吸收邊界的vz分量波場吸收效果對比

Komatitsch和Martin(2007)應用C-PML吸收邊界條件并采用上述模擬參數對磷灰石和鋅晶體模型進行正演模擬.模擬結果表明,磷灰石和鋅晶體模型分別在140 μs和125 μs時刻出現計算不穩定現象,如該文獻中Figure 10所示.

由圖17和圖18可知,EAL邊界條件對各向異性介質具有穩健的適應性,模擬的波場能夠在EAL和MEAL邊界中穩定地傳播,并沒有出現不穩定現象且都得到一定程度的吸收.圖17展示了磷灰石模型加載EAL和MEAL邊界條件的不同時刻模擬結果.不難看出,在25 μs時,彈性波還沒有傳播至模型邊界,加載兩種吸收邊界條件情況下的計算結果幾乎一樣,沒有明顯差別;到了40 μs和50 μs時,EAL邊界條件不能很好地對截斷邊界進行有效的吸收,出現了明顯的強邊界反射,而MEAL邊界條件雖然也存在邊界反射,但是殘余反射能量明顯降低;到了140 μs時,彈性波已經完全離開計算區域,只有吸收邊界的殘余反射波在模型中傳播,但是相比于EAL邊界條件,MEAL條件下的殘余反射波能量明顯降低.同樣地,鋅晶體模型在加載MEAL邊界后,模擬效果明顯提高,來自截斷邊界的反射波能量得到很好的壓制,如圖18所示.

圖19是磷灰石模型和鋅晶體模型的物理區域中每一時刻應用MEAL邊界條件模擬的波場能量與應用EAL邊界條件模擬的波場能量的相對能量差.由圖中可以看出,磷灰石模型和鋅晶體模型分別大約在60 μs和100 μs前,由于震源的激發,波場被注入能量,部分波傳入吸收邊界,能量被吸收層吸收,EAL和MEAL情況下的能量變化基本相同.而在這之后,波已經完全傳播出整個物理區域,兩種邊界條件情況下的波場能量差異逐漸增大,并趨近于一個穩定的區間.相比于EAL邊界條件,MEAL邊界條件對邊界反射吸收效果分別提高了約80%和60%.

圖19 磷灰石和鋅晶體模型應用MEAL邊界條件的波場相對能量差

綜上所述,這兩個各向異性模型的數值實驗結果表明,MEAL和EAL邊界條件在模型各向異性情況下具有一定的穩定性和有效性,且MEAL邊界條件的吸收效果顯著提升.

3 結論

本文以EAL邊界條件為基礎,引入高斯型空間加權函數和減速阻尼因子,對EAL邊界條件進行改進,并采用高階交錯網格有限差分法進行了模型試算,經過正演模擬結果的對比分析得到以下4點結論:

(1) 通過調節高斯型空間加權函數的階數可以保證阻尼項從內邊界平滑遞增到外邊界的同時兼顧內外邊界的綜合邊界反射壓制效果.相比于傳統的指數型空間加權函數,在相同數量吸收層的情況下,能夠更加有效的壓制邊界反射.

(2) 引入減速阻尼因子能夠保證地震波在吸收層中的傳播時間更長,使得地震波在吸收層中更加充分的吸收.同樣地,應用高斯型空間加權函數定義吸收層中減速阻尼因子,能夠確保地震波從吸收層內邊界向外邊界傳播過程中逐步減速的同時,達到平緩過渡的效果,避免由于傳播速度的改變而造成的反射.

(3) 與EAL邊界條件相同,MEAL邊界條件的引入沒有改變控制方程的差分格式,只是改變差分格式的相關的系數.因此,MEAL邊界條件不僅繼承了EAL邊界條件易于實現且計算量和內存消耗小的優點,還提升了對人工邊界反射的吸收效果.

(4) EAL邊界條件在各向異性模型中表現穩定且有效.特別地,相對EAL邊界條件,MEAL邊界條件在各向異性介質中對邊界反射的壓制作用改善明顯,能夠獲得良好的模擬效果.