楊茜娜 張 偉

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

0 引言

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

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

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

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

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

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

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

1 CFS-PML公式形式和參數(shù)取值

1.1 ADE CFS-PML實現(xiàn)形式

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

(1)

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

不同的拉伸方程對應(yīng)不同類型的PML。對式(1)等號兩邊分別求導(dǎo)可得復(fù)坐標(biāo)下的空間導(dǎo)數(shù)與實坐標(biāo)下的空間導(dǎo)數(shù)的關(guān)系

(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)

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

(8)

(9)

(10)

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

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

(11)

(12)

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

(13)

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

(14)

(15)

將以上方法用于整個一階速度—應(yīng)力方程組就是ADE CFS-PML的實現(xiàn)方式。

1.2 CFS-PML性質(zhì)分析

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

(16)

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

(17)

(18)

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

(19)

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

(20)

CFS-PML復(fù)拉伸方程對應(yīng)的實部和虛部分別是

(21)

(22)

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

(23)

(24)

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

1.3 均勻速度模型CFS-PML參數(shù)取值

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

(25)

(26)

(27)

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

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

(28)

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

(29)

式中u0為初始時刻振幅。

結(jié)合式(1)、式(21)、式(22)可得地震波在垂直入射時標(biāo)準(zhǔn)PML傳播方程

(30)

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

(31)

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

通過對不同層數(shù)對應(yīng)的最優(yōu)R值的測試結(jié)果進(jìn)行擬合，得到如下的層數(shù)與預(yù)期反射系數(shù)關(guān)系

(32)

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

(33)

式中VPr是吸收層內(nèi)設(shè)置的參考P波速度。

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

(34)

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

(35)

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

1.4 復(fù)雜速度模型CFS-PML參數(shù)取值

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

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

2 模型實測結(jié)果分析與對比

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

2.1 雙層速度模型

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

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

圖1 高波速比雙層模型

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

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

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

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

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

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

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

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

2.2 多層模型

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

圖6 多層速度模型

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

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

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

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

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

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

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

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

2.3 Marmousi模型

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

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

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

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

圖11 Marmousi模型示意圖

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

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

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

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

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

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

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

3 結(jié)論

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

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