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

羅玉欽 劉 財

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

0 引言

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

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

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

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

1 理論研究

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

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

(1)

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

采用復坐標伸展變換的形式為

(2)

(3)

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

(4)

(5)

(6)

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

(7)

且有

(8)

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

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

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

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

先分析NPML的平面波解，經過復坐標變換之后方程具有新的平面波解

(9)

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

(10)

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

(11)

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

(12)

(13)

式中：L是PML的層厚；R是理論邊界反射系數；l為PML區計算點與內部邊界的距離；η0的經驗值范圍是0～6；β0的經驗值范圍是1～10； 指數qβ、qη與n相同，通常取1～3。

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

(14)

(15)

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

(16)

=ηx(x)Txx+iωTxx

(17)

轉換到時間域，有

(18)

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

(19)

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

采用交錯網格有限差分數值模擬時，將式(19)離散，得到如下離散格式

(20)

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

(21)

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

圖1 衰減因子PML區域示意圖

將雙衰減剖面應用于NPML，如區域2，有

(22)

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

(23)

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

(24)

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

2 數值模擬

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

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

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

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

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

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

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

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

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

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

圖3 能量衰減曲線

圖4 反射波振動曲線

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

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

3 衰減函數

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

(25)

(26)

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

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

4 結束語

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