復(fù)頻移完全匹配層在復(fù)雜速度模型中的應(yīng)用

楊茜娜 張 偉

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

0 引言

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

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

在此基礎(chǔ)上，近年也有不少學(xué)者對(duì)已提出的PML方法進(jìn)行對(duì)比和改進(jìn)[10-12]。李佩笑等[10]對(duì)比分析了非分裂與分裂兩種實(shí)現(xiàn)方法； 廉西猛等[11]研討了多種吸收方式的異同； Liu等[12]剖析了影響衰減的各因子。為了在地震波數(shù)值模擬程序中應(yīng)用CFS-PML，人們研發(fā)了多種具體實(shí)現(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ù)值格式求解，具有實(shí)現(xiàn)簡(jiǎn)單、非分裂形式、不依賴于時(shí)間積分格式等優(yōu)點(diǎn)，獲得廣泛應(yīng)用[16-17]。

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

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

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

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

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

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

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

(1)

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

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

(2)

以二維P-SV波波動(dòng)方程為例，方程如下

ρ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為時(shí)間域速度分量；τxx、τzz、τxz為時(shí)間域應(yīng)力分量，下標(biāo)中逗號(hào)后的t、x、z表示對(duì)t、x、z的導(dǎo)數(shù)。以式(3)為例，將拉伸方程代入得到

(8)

(9)

(10)

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

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

(11)

(12)

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

(13)

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

(14)

(15)

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

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

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

(16)

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

(17)

(18)

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

(19)

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

(20)

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

(21)

(22)

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

(23)

(24)

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

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

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

(25)

(26)

(27)

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

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

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

(30)

該式第二個(gè)指數(shù)項(xiàng)是衰減項(xiàng)。考慮雙程衰減，R可用衰減項(xiàng)表示為

(31)

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

通過對(duì)不同層數(shù)對(duì)應(yīng)的最優(yōu)R值的測(cè)試結(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作為低通濾波器時(shí)的拐角角頻率[4]，取值太低時(shí)CFS-PML對(duì)地震波傳播等效于常規(guī)PML，不能發(fā)揮抑制非均勻波的作用； 取值太高則導(dǎo)致大部分頻率不吸收，沒有起到PML的作用。通常設(shè)置為震源中心角頻率fc的一半[5]，即

(34)

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

(35)

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

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

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

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

2 模型實(shí)測(cè)結(jié)果分析與對(duì)比

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

2.1 雙層速度模型

首先考察存在速度劇烈變化的雙層模型，本次測(cè)試所用雙層模型的速度設(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è)置參考速度： ①逐點(diǎn)取邊界介質(zhì)速度對(duì)應(yīng)值； ②統(tǒng)一取邊界介質(zhì)速度中間值(5km/s)； ③統(tǒng)一取邊界介質(zhì)速度最大值； ④計(jì)算一維垂向分布的每層平均速度，光滑后作為該層參考速度。

圖1 高波速比雙層模型

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

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

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

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

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

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

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

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

2.2 多層模型

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

圖6 多層速度模型

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

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

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

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

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

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

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

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

2.3 Marmousi模型

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

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

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

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

圖11 Marmousi模型示意圖

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

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

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

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

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

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

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

3 結(jié)論

在數(shù)值模擬中，CFS-PML對(duì)虛假反射的吸收效果尤為突出。本文對(duì)各參數(shù)設(shè)置和作用進(jìn)行了深入討論，歸納出各參數(shù)的取值方式并進(jìn)行了優(yōu)化，對(duì)反射系數(shù)R提出了新的取值公式； 針對(duì)復(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ù))則采取先選定參考速度值再做設(shè)置。

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