楊茜娜 張 偉

(①中國科學技術大學地球和空間科學學院，安徽合肥 230026； ②南方科技大學地球與空間科學系，廣東深圳 518055； ③南方科技大學深圳市深遠海油氣勘探技術重點實驗室，廣東深圳 518055)

0 引言

完全匹配層(Perfectly matched layer，PML)被廣泛應用于波動傳播的數值模擬，以解決數值模型的邊界吸收問題。其初期應用是Bérenger[1]在計算電磁學領域提出的一種有效吸收的邊界條件。由于PML具有吸收效果好、適用范圍廣的特點，Chew等[2]將其引入地震波數值模擬。此后，PML在地震波模擬中被廣泛應用[3-7]。

后續(xù)研究發(fā)現PML對掠射到吸收層的虛假反射的吸收效果欠佳。在計算電磁學領域，為了增強對掠射波和非均勻波的吸收效果，Kuzuoglu等[8]提出將PML在復數域的坐標拉伸方程的極點移到非零區(qū)域。在計算地震學領域，Festa等[5]將基于復頻移的PML應用于地震波數值模擬，證實這樣可提高掠射波的吸收。Roden等[9]將這種PML的改進形式命名為復頻移完全匹配層(Complex frequency-shifted PML，CFS-PML)。

在此基礎上，近年也有不少學者對已提出的PML方法進行對比和改進[10-12]。李佩笑等[10]對比分析了非分裂與分裂兩種實現方法； 廉西猛等[11]研討了多種吸收方式的異同； Liu等[12]剖析了影響衰減的各因子。為了在地震波數值模擬程序中應用CFS-PML，人們研發(fā)了多種具體實現形式，包括卷積法PML(C-PML)[9，13]、積分法[14]、Z變化法[15]和輔助微分方程法(ADE CFS-PML)[16-17]等。其中ADE CFS-PML將CFS-PML轉換成偏微分方程，可使用與內部相同的數值格式求解，具有實現簡單、非分裂形式、不依賴于時間積分格式等優(yōu)點，獲得廣泛應用[16-17]。

實際應用CFS-PML需合理設置三個參數：衰減參數d、頻移參數α和坐標實拉伸參數β。d是地震波在PML傳播中最主要的衰減因子，表征地震波在PML內傳播時的振幅和能量衰減。α使PML的衰減變成與頻率相關，即低頻(頻率小于α)是衰減弱→不衰減，高頻(頻率大于α)是正常衰減。該參數使CFS-PML對地震波傳播的作用類似于低通濾波器。α的引入盡管破壞了標準PML對所有頻段都相同吸收的優(yōu)異性能，但實際模擬中發(fā)現該參數可降低掠射波入射到PML界面時產生的非均勻波[9，13]。補償該參數對低頻吸收的破壞的一種方案是采用雙極點CFS-PML[18-19]。β使垂向空間長度拉伸β倍，等效于長度不拉伸、但介質變成垂向VTI，使掠射波在PML中波前面向界面法線方向彎曲，增大方程中的角度項，增強掠射波吸收效果[19]。

CFS-PML的參數設置包括兩個方面： 一是參數沿PML的變化形式，即空間分布； 二是參數在最內或最外層的最大值。前人對d提出過多種空間分布(函數)形式，如冪函數[3，13，16，20-21]、對稱函數[22]、正余弦函數[23-24]、復合函數[25-26]等形式。α、β的空間分布以冪函數形式為主[13，16]。

衰減參數最大值d0的取值本質上取決于預期最小反射系數R，通常采用Collino等[3]給出的經驗設置。Festa等[5]給出頻移參數最大值取值公式α0=πfc，其中fc是震源子波中心頻率。Zhang等[16]指出β0可允許坐標拉伸后中心波長對應的每波長網格數降為色散誤差允許的每波長網格數的一半。

以上各參數的最大值都是基于均勻速度模型設置的。復雜速度模型存在橫向和縱向速度變化，對衰減起關鍵作用的參數d0與吸收層內速度值相關； 在復雜速度模型情況下，是按每個點速度值估算參數最大值，還是在整個吸收層按統(tǒng)一參考速度設置衰減參數最大值，尚難覓相關文獻報道及對此的系統(tǒng)評估。本文充分調研了復雜速度模型CFS-PML參數設置方法，通過對多種復雜速度模型做數值測試，總結出復雜速度模型CFS-PML優(yōu)化取值方案，為CFS-PML的實際應用提供參考。

1 CFS-PML公式形式和參數取值

1.1 ADE CFS-PML實現形式

從CFS-PML的各種實現方式看，ADE CFS-PML具有適用于各種時間積分法的優(yōu)勢，本文以ADE CFS-PML為例介紹CFS-PML的實現[16]。

(1)

式中：sx(η)表示拉伸方程；η表示x方向積分變量。

不同的拉伸方程對應不同類型的PML。對式(1)等號兩邊分別求導可得復坐標下的空間導數與實坐標下的空間導數的關系

(2)

以二維P-SV波波動方程為例，方程如下

ρvx,t=τxx,x+τxz,z

(3)

ρvz,t=τzx,x+τzz,z

(4)

τxx=(λ+2μ)vx,x+λvz,z

(5)

τzz=λvx,x+(λ+2μ)vz,z

(6)

τxz=μ(vx,z+vz,x)

(7)

式中：ρ是密度；λ、μ是拉梅常數；vx、vz為時間域速度分量；τxx、τzz、τxz為時間域應力分量，下標中逗號后的t、x、z表示對t、x、z的導數。以式(3)為例，將拉伸方程代入得到

(8)

(9)

(10)

式中n表示x、y和z中的一個方向，則dn、αn、βn分別是對應方向上的衰減參數、 頻移參數和坐標實拉伸參數。以下簡述ADE方法實現方式。

式(8)等號右側代入拉伸方程后的微分項為

(11)

(12)

將式(12)代入式(8)，可得

(13)

將式(11)和式(12)從頻率域轉換到時間域，就得到了輔助方程CFS-PML方程

(14)

(15)

將以上方法用于整個一階速度—應力方程組就是ADE CFS-PML的實現方式。

1.2 CFS-PML性質分析

以沿+x方向傳播的行波(traveling wave)、非均勻波(evanescent wave)兩類平面地震波為例，通過二維平面(x，z)波動方程的傳播具體分析地震波在PML中的傳播和衰減特性。行波可表示為

(16)

平行于介質—PML下界面?zhèn)鞑?，沿界面法線方向衰減的非均勻波可表示為

(17)

(18)

將式(18)代入式(16)，則可轉化為

(19)

再將式(18)代入式(17)，可轉化為

(20)

CFS-PML復拉伸方程對應的實部和虛部分別是

(21)

(22)

行波和非均勻波在CFS-PML中的方程分別是

(23)

(24)

從式(23)、式(24)易知，CFS-PML存在相應的虛部對行波衰減、相應的實部對非均勻波衰減，吸收效果大大提升。

1.3 均勻速度模型CFS-PML參數取值

由于參數的空間分布形式對結果影響不大，這里采用最常用的冪函數形式[16，20]

(25)

(26)

(27)

式中：指數pd通常取2、3、4，現有研究尚未發(fā)現這三個不同取值對結果具有本質影響； 指數pα通常取1、pβ通常取2；d0、α0、β0分別是衰減、頻移、坐標實拉伸參數最大值；D是波傳播到吸收層內的距離；L是吸收層的厚度。

上述三種參數最大值中，d0是控制地震波在PML中衰減程度的參數，即是決定吸收效果的一個重要因子。d0值通常依據對于垂直入射的平面波經過PML后的預期反射系數R確定：理論上，越小的預期反射系數越符合實際需求； 但預期反射系數太小，會導致d在PML內的空間分布變化過快，有可能由于過快的空間分布變化導致反射波。因此對不同的PML層數通常設置不同的預期反射系數。不同層數對應的最優(yōu)R值大多是按經驗設定的，并往往通過數值測試進行測定。Collino等[3]經驗性地給出5層R=10-2，10層R=10-3，20層R=10-4。為了能對任意PML層數設置R值，Marcinkovich等[21]對Collino等[3]的離散值進行了擬合。Zhang等[16]通過分析離散值之間的關系，得到如下更具一般性的層數與R值的經驗關系式

(28)

式中N為設定的吸收層層數。給定預期反射系數R后，根據垂直入射的平面波，按d的空間函數分布穿過PML再反射回物理區(qū)間的衰減預取反射系數R，給出最大值d0。沿x傳播的平面波經坐標拉伸后表示為

(29)

式中u0為初始時刻振幅。

結合式(1)、式(21)、式(22)可得地震波在垂直入射時標準PML傳播方程

(30)

該式第二個指數項是衰減項?？紤]雙程衰減，R可用衰減項表示為

(31)

針對式(25)中不同的pd形式，利用垂直入射測試模型進行數值測試，發(fā)現Collino等[3]給出的預期反射系數不是最優(yōu)結果，相同層數下的預期反射系數可更低。按照原始公式進行設置：N=5時R=10-2，N=10時R=10-3，N=20時R=10-4； 不同吸收層數對應的誤差量級分別為：N=5時為10-2，N=10時為10-3，N=20時為10-4。而實際上按照N=5時R=10-3、N=10時R=10-4、N=20時R=10-5設置，誤差在原有基礎上減小一個量級，即N=5時為10-3、N=10時為10-4、N=20時為10-5，顯然其吸收效果更優(yōu)。

通過對不同層數對應的最優(yōu)R值的測試結果進行擬合，得到如下的層數與預期反射系數關系

(32)

將相應的衰減參數d的取值公式代入式(31)，通過積分變換就可得到用R表示的d0取值公式[3]

(33)

式中VPr是吸收層內設置的參考P波速度。

α0是CFS-PML作為低通濾波器時的拐角角頻率[4]，取值太低時CFS-PML對地震波傳播等效于常規(guī)PML，不能發(fā)揮抑制非均勻波的作用； 取值太高則導致大部分頻率不吸收，沒有起到PML的作用。通常設置為震源中心角頻率fc的一半[5]，即

(34)

β0是垂直坐標的最大拉伸程度，拉伸程度越大就越有利于將波前向界面法線方向彎曲； 但過大的β0會導致等效網格的尺度過大，超過數值格式分辨能力，進而導致強反射。可以預期，網格的空間步長越小(相比于色散誤差要求的網格大小)，β0越大。通過數值實驗，Zhang等[16]指出最大β0可允許坐標拉伸后中心波長對應的每波長網格數降為色散誤差允許的每波長網格數的一半。馮海珂[19]進一步將該準則明確表示為如下形式

(35)

式(25)～式(27)和式(32)～式(35)構成CFS-PML參數的最優(yōu)取值方案。從上述公式易知，這些取值都只是針對均勻速度模型。

1.4 復雜速度模型CFS-PML參數取值

從上述針對CFS-PML各參數的取值公式及最大值設置的討論得知： 預期反射系數R和頻移參數α0不依賴于速度模型，反射系數取值(式(32))僅與吸收層數相關，頻移參數取值(式(34))僅與震源中心頻率相關； 但衰減參數d0和坐標實拉伸參數β0除了與固定的震源性質和網格性質相關外，還同吸收層內速度值和波長相關，即依賴于速度模型，其中d0與P波速度相關而β0與S波速度相關，d0的設置對吸收效果最為重要。從式(33)可看出，復雜速度模型下關鍵是如何設置參考速度VPr。所謂PML層內參考速度，即按該速度垂直傳播進入邊界的波，正好可實現預計的反射系數R的吸收。理論上若實際傳播速度小于該速度，則在吸收層內經歷過多衰減； 若實際傳播速度大于該速度，則在吸收層內吸收不足，在邊界處產生較大反射。

根據實際情況，本文歸納出三種實用的復雜模型參考速度取值方法： ①PML層內逐點設置不同的值，每點取所處物理區(qū)間P波速度； ②針對復雜速度模型取一個等效速度值作為參考速度值，可以是最大值(即相應PML層內所有速度的最大值)或中間值(即PML層內所有速度的平均值)； ③對物理區(qū)間真實垂向速度值分布進行光滑，每點取光滑后逐點對應的不同的速度值。

2 模型實測結果分析與對比

采用上述各方法針對多種速度模型進行實測和分析，以對比和遴選針對復雜模型的參考速度取值方案。

2.1 雙層速度模型

首先考察存在速度劇烈變化的雙層模型，本次測試所用雙層模型的速度設定為： 上層VP1=0.5km/s、VS1=0.3km/s，下層VP2=10.0km/s、VS2=6.0km/s，速度變化幅度達20倍(圖1)。各衰減參數空間分布參照式(25)～式(27)進行選取，與VP相關的d0和β0通過選取VP而確定； 反射系數R按式(32)選取為10-5； 頻移參數按式(34)選取為α0=πfc≈4.71。

采取以下四種方式(圖2)設置參考速度： ①逐點取邊界介質速度對應值； ②統(tǒng)一取邊界介質速度中間值(5km/s)； ③統(tǒng)一取邊界介質速度最大值； ④計算一維垂向分布的每層平均速度，光滑后作為該層參考速度。

圖1 高波速比雙層模型

模型尺寸為5.7km×3.0km，網格單元為10m×10m，吸收層數N=20； 震源是中心頻率為1.5Hz雷克子波，位于(0.4km，-1.0km)，與邊界距離為0.2km； 三個檢波點由上至下編號為1～3，所處位置橫坐標均為0.3km，縱坐標為-0.5～-2.5km，深度間隔為1.0km測試的參考解通過在原模型上保持震源、檢波點相對位置不變，擴大計算區(qū)間(7km×7km)設計而成。擴大后形成的參考模型保證在觀測時窗內、最外層邊界產生的誤差不會被檢波器接收到，從而使檢波器在固定時窗內接收到的波形等效為在該時間范圍內“無限大理想模型”接收到的地震波。

圖2 雙層模型不同參考速度設置方式下垂直速度分布

在t=2.4s時通過各取值方式與參考模型(圖3e)波場快照的對比可明顯看出，設置為對應值(圖3a)或垂向速度值(圖3d)在左側下半部邊界處吸收不干凈并有一定程度的反射，取統(tǒng)一值(圖3b、圖3c)吸收效果更好。而t=6.0s時，對于上部低速層，取對應值(圖3f)時反射較弱，吸收效果較好，其余方式(圖3g～圖3j)均有不同程度的反射。因上、下層速度差異很大，界面反射程度高，上層的波場振幅和能量更大。為了更清楚地顯示下層波動，所以波場快照范圍較小，顯示的反射波對主波的影響無法定量估量。為定量地考察上、下層吸收效果，通過波形圖和誤差圖進行進一步分析、討論。

圖3 吸收層參考速度不同取值方式下雙層模型的波場快照

從位于不同深度檢波器接收的波形對比圖(圖4)不難看出，參考速度設置為對應值和垂向平滑值時在下部高速層(圖4c)呈現與參考波形的大幅度差異，但在兩層界面處(圖4b)參考速度設置為對應值和垂向平滑值在主波幅值上與其余兩種方法差異較小，只是吸收效果仍不理想，且設置為垂向平滑值在上部低速層(圖4a)主波后也出現較明顯差異。統(tǒng)一設置為中間值或最大值時則與參考波形更契合，但在上部低速層(圖4a)和兩層界面處(圖4b)主波后也出現小幅度反射現象。

圖4 雙層模型四種參數設置方法記錄波形與參考波形對比

圖5 雙層模型四種參數設置方法的測試波形與參考波形的相對誤差對比

圖5用具象化誤差對比展示四種方法吸收效果的差異。以參考模型為標準，參考波形最大值為歸一化值對參考模型歸一化。測試模型波形除以該值與參考模型歸一化的差，即兩者間的相對誤差。不同參數設置方法的誤差不同，從誤差值的差異評估不同方法的吸收效果。統(tǒng)一設置為中間值或最大值在波形上都較為契合，但從誤差的數值分布可見設置為中間值或最大值的誤差遠小于另兩種方法(圖5c)，那兩種方法的誤差遠超出圖示范圍。而相較于統(tǒng)一設置為最大值，統(tǒng)一設置為中間值時在各深度處(檢波器)的吸收效果都更好(圖5a～圖5c)，吸收誤差(范圍約為10-3～10-2)相對更小、更平穩(wěn)。

2.2 多層模型

雙層模型中設定了一個極端的上下速度比為20的例子說明統(tǒng)一設置參考速度吸收效果的優(yōu)越性，在實際應用中很少出現類似情況，因此又設計了多層且速度逐漸變化的模型，以更契合地模擬實際復雜地形情況下的吸收效果。

圖6 多層速度模型

各層速度和對應界面深度由淺至深分別為1.0、2.5、4.0、5.5、7.0km/s，-0.7、-1.2、-1.8、-2.3km； 震源是中心頻率為5Hz的雷克子波，其位置及距邊界的距離、三個檢波點所處位置、深度間隔及編號等同圖1

參考速度設置(圖7)采取與雙層模型相同的四種方式，只是第二種取邊界介質速度中間值時改用4.0km/s作為該層參考速度值。

圖7 多層模型參考速度不同設置方式下垂直速度分布

t=0.84s時刻不同設置方式波場快照與參考模型(圖8e)對比可顯示四種方法差異，設置為對應值(圖8a)和垂向速度值(圖8d)在紅圈區(qū)域產生較明顯反射波，而統(tǒng)一設置(圖8b、圖8c)則吸收得較干凈。

圖8 吸收層參考速度不同取值方式下多層模型t=0.84s時刻波場快照

類似于雙層模型，參考速度設置為對應值和垂向速度值在上部低速層(圖9a)與參考波形較吻合，而隨著深度(速度)增加波形的不匹配性逐漸顯著(圖9b)，到下部高速層(圖9c)出現遠高于主波的反射。統(tǒng)一設置為中間值和最大值，在3個檢波點與參考波形的對比中都展現出高度匹配性，兩者幾乎重合。統(tǒng)一設置為中間值和最大值的相對誤差低至10-4～10-3，而其余兩種方法的誤差為統(tǒng)一設置速度值的102～103倍，因此在誤差對比圖中遠超出圖幅范圍。統(tǒng)一設置為中間值和最大值兩者吸收差異在中部較高速層(圖10b、圖10c)相對較小，但在上部低速層(圖10a)統(tǒng)一設置為最大值時則呈現較大誤差值。雖然統(tǒng)一設置為中間值在高速層(圖10c)的誤差略高于設定為最大值，但從整體吸收效果看，參考速度統(tǒng)一設置為中間值吸收效果更優(yōu)異。

圖9 多層模型四種參數設置方法記錄波形與參考波形對比

為了驗證上述結果的普適性，針對不同方向(z)、不同吸收層數(N=10、15)也進行了相同測試，所得結果基本一致。通過對雙層、多層模型不同設置方法的波場快照、波形和誤差的對比及綜合分析可知，對于復雜速度模型而言，參考速度統(tǒng)一設置為中間速度值對衰減參數進行取值的吸收效果更優(yōu)越。

2.3 Marmousi模型

與上述模型相比，Marmousi模型除了存在多層變化速度外，還存在局部速度陡增(如斷層)等情況。震源參數除子波中心頻率為10Hz外，其余與上述相同。

圖10 多層模型四種參數設置方法的測試波形與參考波形的相對誤差對比

由于Marmousi模型無完全參考解，所以測試模型僅選取圖11的中間部分，吸收層數為20。而參考解則是略微擴大模型將吸收層數增至100。

參考速度設置(圖12)仍采用與上述模型相同的四種方式，只是第二種取邊界介質速度中間值時改用2.0km/s，且第三種取邊界介質速度最大值是4.0km/s，作為對應層位的參考速度值。

圖11 Marmousi模型示意圖

對比參考速度四種選取方法波場快照(圖13a)，發(fā)現這四者基本一致；再與參考模型(圖13b)相比，波場快照也基本一致。說明這四種方法在吸收過程中波場殘留較弱，均未出現明顯的非均勻轉換波和強不穩(wěn)定雜波而大范圍、大幅度地影響正常的波傳播。

圖12 Marmousi模型不同參考速度設置方式下垂直速度分布

圖13a Marmousi測試模型取邊界介質速度中間值在不同時刻(t=0.30、0.66、1.14、1.50s)的波場快照

圖13b Marmousi參考模型(速度場Vx)在不同時刻(t=0.30、0.66、1.14、1.50s)的波場快照

圖14 Marmousi模型四種參數設置方法記錄波形與參考波形對比

圖15 Marmousi模型四種參數設置方法的測試波形與參考波形的相對誤差對比圖

對于Marmousi模型，無論是波場快照(圖13)，還是波形圖(圖14)，四種方法的吸收效果基本一致，誤差也基本在一個量級(10-2)。因為Marmousi模型雖較復雜，但其層間速度變化和整體速度變化相對較小，因此四種方法的吸收效果也很接近。盡管四種方法在1號(圖15a)、2號(圖15b)檢波點的誤差基本一致，但在3號檢波點(圖15c)中t=0.8～1.5s時段統(tǒng)一設置為中間值的誤差比其他方法略小(誤差范圍約為±3×10-2)。綜合考慮速度差異大和速度差異小的情況，建議統(tǒng)一設置為中間取值方式。

3 結論

在數值模擬中，CFS-PML對虛假反射的吸收效果尤為突出。本文對各參數設置和作用進行了深入討論，歸納出各參數的取值方式并進行了優(yōu)化，對反射系數R提出了新的取值公式； 針對復雜速度模型，其參數設置方式(取值公式)的關鍵在參考速度的選取，與參考速度無關的衰減參數(反射系數、頻移參數)可按模型相關參數代入取值公式進行設置，與參考速度相關的衰減參數(衰減參數d0、坐標實拉伸參數)則采取先選定參考速度值再做設置。

通過雙層高速度比模型、多層漸變速模型、Marmousi模型，對參考速度的幾種設置方法進行測試，對比分析了各方法所得波場快照、波形圖和誤差圖，發(fā)現“統(tǒng)一設置參考速度值為邊界介質速度中間值”的方法簡捷、吸收效果穩(wěn)定，可為實際復雜速度模型的地震波數值模擬提供參考和借鑒。