多軸復(fù)頻移近似完全匹配層彈性波模擬

羅玉欽 劉 財

(吉林大學(xué)地球探測科學(xué)與技術(shù)學(xué)院，吉林長春 130026)

0 引言

為滿足地震波正、反演的需要，邊界條件的研究引起人們的極大關(guān)注。目前應(yīng)用最廣泛的邊界條件是由Berenger[1-2]在研究麥克斯韋方程時提出的完全匹配層(Perfectly matched layers，PML)。Chew等[3]和Collino等[4]將其解釋為復(fù)坐標(biāo)伸展變換的結(jié)果。當(dāng)入射波角度很大時，衰減系數(shù)會變得很小，因此常規(guī)PML無法吸收大角度入射波，且在極低頻時產(chǎn)生奇異值。復(fù)頻移(Complex frequency shifted，CFS)技術(shù)對復(fù)坐標(biāo)伸展變換(CCS)做進(jìn)一步處理，原來為1的項被尺度因子替換并在大角度入射波的吸收中起作用，在分母中加上頻移因子以消除低頻奇異值[5]。Roden等[6]提出基于CFS卷積的PML(即C-PML)； Komatitsch等[7]進(jìn)一步將其推廣到彈性介質(zhì)中。

此外，近似PML(Nearly PML，NPML)[8]不同于上述常規(guī)PML，它是直接對波場進(jìn)行變換而得到，不改變波場形式，也不需做卷積處理。該方法[8]自提出后很快被推廣試用于聲波介質(zhì)、彈性介質(zhì)以及雙相介質(zhì)中[9-11]。但是，由于該方法在大入射角時吸收能力較弱且在復(fù)雜介質(zhì)的模擬中存在不穩(wěn)定現(xiàn)象，因此，在地震正、反演中現(xiàn)今主流方法仍舊是C-PML。

對于邊界條件的穩(wěn)定性，F(xiàn)esta等[12]指出PML邊界與面波之間會降低吸收效果并產(chǎn)生不穩(wěn)定性。Mezafajardo等[13]認(rèn)為常規(guī)PML在各向同性和各向異性中都不滿足嚴(yán)格的漸進(jìn)穩(wěn)定性，并引入多軸技術(shù)，即在多個正交方向上同時引入衰減因子。同時，衰減因子對邊界吸收效果有很大影響。Collino等[14]提出內(nèi)截斷邊界與PML層的距離呈指數(shù)關(guān)系的衰減函數(shù)； 此后Groby等[15]對其進(jìn)行了改進(jìn)和整理。近年來，中國學(xué)者也對PML做了許多研究[16-28]。田坤等[29]將多軸技術(shù)應(yīng)用于C-PML中； 陳可洋[30]提出了正弦型和余弦型衰減函數(shù)； 羅玉欽等[31]對衰減因子進(jìn)行了系統(tǒng)研究，提出的衰減因子將吸收效果提升了20%～60%。

雖然NPML擁有其自身優(yōu)勢，但它仍存在穩(wěn)定性與吸收效果上的問題。本文基于NPML并引入頻移因子、尺度因子及穩(wěn)定性因子，提出了多軸復(fù)頻移NPML(MCFS-NPML)吸收邊界，旨在提升NPML對大角度入射波的吸收且增強(qiáng)其穩(wěn)定性； 在此基礎(chǔ)上采用文獻(xiàn)[31]提出的衰減函數(shù)，進(jìn)一步提升MCFS-NPML的吸收能力。本文詳細(xì)分析MCFS-NPML的穩(wěn)定性和吸收效果，探討各因子在改善吸收效果及穩(wěn)定性方面所起的作用，證實(shí)了MCFS-NPML擁有更好的吸收效果和穩(wěn)定性； 同時，將該衰減函數(shù)與廣泛運(yùn)用的指數(shù)型衰減函數(shù)進(jìn)行比對，結(jié)果表明在該衰減函數(shù)下吸收效果得到進(jìn)一步改善。

1 理論研究

1.1 近似完全匹配層(NPML)吸收邊界

本文采用交錯網(wǎng)格有限差分法，因此采用一階速度—應(yīng)力形式方程組。將方程組變換到頻率域

(1)

式中：ω是角頻率；λ、μ是拉梅常數(shù)；ρ是密度；Vx、Vz為頻率域彈性波場速度分量；Txx、Tzz、Txz為頻率域彈性波場應(yīng)力分量。

采用復(fù)坐標(biāo)伸展變換的形式為

(2)

(3)

將式(3)代入式(1)，并直接作用在速度分量和應(yīng)力分量上，得到復(fù)坐標(biāo)伸展變換后的方程組

(4)

(5)

(6)

同理，采用相同方法也將其他變量變換到時間域，有

(7)

且有

(8)

ξ=τxx,τzx,τzz,vx,vzm=x,z

式中：τxx、τzx、τzz是時間域應(yīng)力；vx、vz為時間域速度參數(shù)。分別與頻率域中的Txx、Tzx、Tzz和Vx、Vz相對應(yīng)。

從式(7)可見主控方程的形式并未發(fā)生變化。這里需引入8個輔助變量替換原主控方程中的變量。在程序編寫時不需改動主代碼，易于實(shí)現(xiàn)。式(8)中的原始控制方程中的輔助變量是原變量經(jīng)伸展變換得到，是將原波場直接變換后再代入相同形式的控制方程求解，因此易于將NPML推廣到其他介質(zhì)。在計算區(qū)域中衰減因子值為0，故式(8)中原變量與新變量是相同的，因此在編程時可只在PML區(qū)域中進(jìn)行波場變換以減少存儲空間。

1.2 復(fù)頻移近似完全匹配層(CFS-NPML)吸收邊界

先分析NPML的平面波解，經(jīng)過復(fù)坐標(biāo)變換之后方程具有新的平面波解

(9)

修改式(3)，并引入尺度因子βx(x)和頻移因子ηx(x)，可得

(10)

得到新的平面波解的衰減項為

(11)

尺度因子βx(x)使大角度入射波的傳播方向向法線方向彎曲，導(dǎo)致衰減系數(shù)增加； 頻移因子ηx(x)對入射角的影響不大，但可有效地避免奇異值的出現(xiàn)。當(dāng)入射角較大時，在常規(guī)PML中波只能在很淺的介質(zhì)中傳播，無法達(dá)到較好的吸收效果，尺度因子的引入增強(qiáng)了邊界對這一類波的吸收[7]。衰減函數(shù)、尺度因子與頻移因子的表達(dá)式為

(12)

(13)

式中：L是PML的層厚；R是理論邊界反射系數(shù)；l為PML區(qū)計算點(diǎn)與內(nèi)部邊界的距離；η0的經(jīng)驗(yàn)值范圍是0～6；β0的經(jīng)驗(yàn)值范圍是1～10； 指數(shù)qβ、qη與n相同，通常取1～3。

以式(10)替換式(3)，如Txx

(14)

(15)

=[ηx(x)+iω]Txx

(16)

=ηx(x)Txx+iωTxx

(17)

轉(zhuǎn)換到時間域，有

(18)

同理，求得其余變量的微分方程且替換式(8)，則可簡寫為

(19)

ξ=τxx,τzx,τzz,vx,vzm=x,z

采用交錯網(wǎng)格有限差分?jǐn)?shù)值模擬時，將式(19)離散，得到如下離散格式

(20)

1.3 多軸復(fù)頻移NPML(MCFS-NPML)吸收邊界

(21)

式中αx的取值不再是0，而與原衰減函數(shù)成比例，比例系數(shù)P(x/z)即是穩(wěn)定性因子 。

圖1 衰減因子PML區(qū)域示意圖

將雙衰減剖面應(yīng)用于NPML，如區(qū)域2，有

(22)

對于 CFS-NPML，在區(qū)域1和區(qū)域2中的處理方法與上面不一樣。同樣以區(qū)域1為例，變換函數(shù)為

(23)

將式(23)代入式(14)中，得到

(24)

這樣就避免了單剖面衰減引起的平面波解呈指數(shù)增長。采用多軸衰減剖面引入了穩(wěn)定性因子，使系數(shù)矩陣特征值向負(fù)半軸移動，消除了不穩(wěn)定性。穩(wěn)定性因子應(yīng)小于1，且隨介質(zhì)的復(fù)雜程度而改變，不宜太大，否則會導(dǎo)致虛假反射增強(qiáng)。

2 數(shù)值模擬

在地震波模擬中分別實(shí)現(xiàn)NPML、M-NPML、CFS-NPML及MCFS-NPML。選取的時間步長為1ms，加載震源是主頻為25Hz的雷克子波。構(gòu)建模型的尺寸是5200m×1200m(包含邊界厚度)，x方向和z方向網(wǎng)格步長都是10m?？v波速度為3000m/s，橫波速度為1400m/s，密度為2000kg/m3，PML層數(shù)L設(shè)為10，R取10～6。將震源加載于x=2610m、z=110m處，貼近于PML邊界，往兩側(cè)傳播的波近似平行于上邊界傳播，截取上述幾種近似完全匹配層下的地震波波場快照。

圖2a和圖2b分別顯示上述四種PML在0.90s和0.95s時刻地震波傳播情況。從上到下依次呈現(xiàn)穩(wěn)定性因子、尺度因子、頻移因子對匹配層的影響及三個因子共同作用下的波場； 圖2a第一排標(biāo)示的①、②和③分別對應(yīng)直達(dá)波及在下、上邊界形成的反射波； 圖2b第一、第四排的第一列中黑色橢圓即匹配層中殘余能量團(tuán)，該能量團(tuán)是因波在邊界處衰減不完全產(chǎn)生的損耗波，且入射角度越大，能量團(tuán)越明顯。該能量團(tuán)的存在會導(dǎo)致虛假反射回傳到計算區(qū)域，且能量團(tuán)的積累使計算不穩(wěn)定。

通過圖2a第一排可見增加穩(wěn)定性因子時反射波能量被增強(qiáng)； 從圖2b第一排中可看出穩(wěn)定性因子值越大，邊界區(qū)域越干凈，消除了大部分的殘余能量團(tuán)。對于CFS-NPML，選取多組參數(shù)分析頻移因子和尺度因子對近似完全匹配層的影響。通過圖2a對比NPML與η0=0、β0=2 時的CFS-NPML，可見殘余能量團(tuán)稍有減弱，且下邊界處的虛假反射也得到削減，但幅度較小。

對比η0、β0分別取1、2和5這幾組波場快照可知：改變η0值會顯著影響吸收效果，且η0取值越大吸收效果越好，但在邊界中的波“拖尾”現(xiàn)象更嚴(yán)重(圖2b第三排)； 但β0取值過大，反而會影響吸收效果。對于MCFS-NPML、CFS-NPML中邊界處未能完全吸收的殘余能量被再次削弱，MCFS-NPML能進(jìn)一步清除PML層中的能量以保持匹配層的穩(wěn)定且再次壓制了因大角度入射所帶來的虛假反射(圖2b第四排)。

圖3中0時刻介質(zhì)中的能量迅速增加，在約0.30s時因波的下側(cè)已到達(dá)邊界，波的能量逐漸減弱。圖中實(shí)線和虛線分別對應(yīng)內(nèi)區(qū)間(不包括PML邊界)和整體區(qū)域的能量曲線。在0.83s時刻波抵達(dá)左右兩側(cè)，經(jīng)過PML層作用能量迅速降低，在約2.50s時殘余的反射波再次抵達(dá)左、右兩邊界，能量再次減少。圖3a中NPML下虛線與實(shí)線的差異較大，表明有大量能量殘留在邊界中。該殘留能量會使PML層變得不穩(wěn)定，同時過多的能量積累會影響邊界的吸收性能。M-NPML曲線差隨穩(wěn)定性因子取值增大而逐漸減小。當(dāng)穩(wěn)定性因子取0.1時，虛線與實(shí)線幾乎重合，證明邊界中的大部分能量均被消除。對于CFS-NPML，從β0取2、η0取0、1、2、3、5這幾組能量曲線(圖3c)看出，η0值越大，在約0.83s時能量曲線最低。

正演模擬中，期望地震波在第一次達(dá)到邊界處時就被盡可能地吸收，防止地震波反射回計算區(qū)域。頻移因子改善了PML吸收性能，但隨著η0增大，左右兩側(cè)初次未被吸收的微弱反射波再次到達(dá)邊界時，整體區(qū)間的能量曲線擾動劇烈，是不穩(wěn)定的。β0取值過大也會有這樣的擾動，且CFS-NPML中曲線差依舊很大； 而在引入穩(wěn)定性因子之后(MCFS-NPML)，能量擾動和曲線差過大將被壓制。

圖2 三種因子不同組合下的波場快照

圖4是模型中x=110m、z=110m點(diǎn)處0～3.20s時段的x方向速度分量的變化曲線，可清晰地看到上述三個因子對常規(guī)入射波吸收效果的影響。穩(wěn)定性因子取0.1時，來自底界的虛假反射波最強(qiáng)。從反射波振幅看，CFS-NPML比NPML和M-NPML低很多，來自另一側(cè)的反射波隨η0增大而減?。?對于MCFS-NPML，不僅反射波能量微弱，來自另一側(cè)的反射波衰減也特別明顯。證實(shí)了M-NPML會增強(qiáng)虛假反射，CFS-NPML和MCFS-NPML卻有利于入射波的吸收。

為進(jìn)一步測試邊界的穩(wěn)定性與吸收性能，采用了復(fù)雜Marmousi模型。在該模型中不穩(wěn)定性被放大，震源靠近邊界，大角度的入射波進(jìn)入上邊界并截取0.5s和4.0s時的波場快照(圖5)。

在NPML中，上邊界存在強(qiáng)烈的地震波殘留(圖5a，紅色橢圓所示)，CFS-NPML也是如此。從4.0s波場快照可見NPML下整個計算區(qū)域已因誤差積累而發(fā)散。至于M-NPML(圖5b、圖5c)，隨著穩(wěn)定性因子的增加，波場殘留減弱。β0取2、η0取0時的CFS-NPML(圖5d)與NPML對比，雖然產(chǎn)生累積誤差，但污染區(qū)域相對較小，表明β0有微弱壓制作用； 當(dāng)η0取2時(圖5e)也未見波場被污染，表明頻移因子對穩(wěn)定性有一定作用。而MCFS-NPML(圖5f)因?yàn)榉€(wěn)定性因子的引入，上邊界比CFS-NPML更干凈。

圖3 能量衰減曲線

圖4 反射波振動曲線

圖5 0.5s (左)和4.0s(右)時刻Marmousi模型下的波場快照

簡而言之，在提高穩(wěn)定性方面，三個因子都有貢獻(xiàn)，但機(jī)理不同。穩(wěn)定性因子沿平行匹配層方向增加一個衰減剖面，吸收沿該方向傳播的耗散波； 頻移因子避免了低頻奇異值的產(chǎn)生； 尺度因子使波向法線方向彎曲，增強(qiáng)對匹配層中波的吸收能力。Marmousi模型模擬中的不穩(wěn)定隨頻移因子的引入而消除，故該不穩(wěn)定是低頻奇異值所致。雖然CFS-NPML下計算保持穩(wěn)定，但上界殘余能量團(tuán)十分明顯，因此需同時引入三種因子使NPML具有更好的穩(wěn)定性和吸收效能。

3 衰減函數(shù)

嘗試通過修改衰減函數(shù)，進(jìn)一步提升邊界吸收能力。用網(wǎng)格算法模擬地震波傳播，衰減函數(shù)不再是一條連續(xù)的曲線，匹配層之間的衰減函數(shù)值差異較大，引起虛假反射。通過增加層數(shù)縮小該差異，但會增加存儲量，降低計算效率。探討在不增加層數(shù)的基礎(chǔ)上進(jìn)一步削弱離散差異帶來的虛假反射。設(shè)計下式

(25)

(26)

波場模擬采用MCFS-NPML，且β0=2、η0=2及P=0.005，模型尺寸為2000m×2000m。當(dāng)層數(shù)取8時，γ=0.065、δ=0.094這一組系數(shù)下波反射最弱。圖6b中藍(lán)線和黑線對應(yīng)傳統(tǒng)指數(shù)型衰減函數(shù)和式(26)給出的函數(shù)下的反射情況，可看出新函數(shù)下的虛假反射更弱。新函數(shù)取8層達(dá)到的吸收效果與原函數(shù)取10層相當(dāng)。

圖6 吸收效果圖(a)與波形曲線(b)

4 結(jié)束語

為將NPML更好地應(yīng)用于正反演，本文推導(dǎo)出MCFS-NPML吸收邊界條件。該過程中引入了頻移因子、尺度因子和穩(wěn)定性因子，這三種因子保證了波場模擬穩(wěn)定有效地進(jìn)行。將它們同時引入彌補(bǔ)了僅引入某種單一因子的不足，相互彌補(bǔ)了各自固有的不利影響。通過分析衰減函數(shù)的引入過程以及分布方式，改進(jìn)了衰減函數(shù)，最終達(dá)到進(jìn)一步提升吸收效果的目的。