羅玉欽 劉 財(cái)

(吉林大學(xué)地球探測(cè)科學(xué)與技術(shù)學(xué)院，吉林長(zhǎng)春 130021)

1 引言

為了消除正演模擬中邊界反射，需要引入吸收邊界條件。Berenger[1]在研究麥克斯韋方程時(shí)提出了完全匹配層(PML)邊界條件，是目前吸收效果最好的邊界條件。在地震波場(chǎng)模擬中，最早采用的是分裂PML邊界條件。Chew等[2]和Collino等[3]對(duì)PML邊界條件進(jìn)行了一般性地推導(dǎo)，將其解釋為復(fù)坐標(biāo)伸展變換的結(jié)果；Chew等[4]和Hastings等[5]首先將PML邊界條件引入地震波模擬。Komatitsch等[6]將該方法推廣到二階彈性波方程中，但需更大的存儲(chǔ)量并且降低了計(jì)算效率，因此發(fā)展了非分裂PML邊界條件，但如何計(jì)算或者避開(kāi)卷積項(xiàng)一直是研究重點(diǎn)。Drossaert等[7]采用遞歸積分的方法避開(kāi)卷積項(xiàng)；Cummer[8]提出了近似PML(Nearly PML，NPML)邊界條件；NPML邊界條件應(yīng)用了波場(chǎng)變換，既不改變波場(chǎng)形式，也不需要進(jìn)行卷積處理，很快就被推廣到聲波介質(zhì)、彈性介質(zhì)和雙相介質(zhì)的地震波模擬[9-11]。但當(dāng)波近平行入射PML邊界時(shí)，衰減系數(shù)會(huì)變得很小，無(wú)法吸收大角度入射波和瞬逝波。Kuzouglu等[12]提出了復(fù)頻移技術(shù)(CFS)，將進(jìn)入PML區(qū)域的地震波傳播方向向法線方向彎曲，增強(qiáng)了對(duì)大角度入射波的吸收。Roden等[13]提出了基于CFS的卷積完全匹配層(C-PML)邊界條件；Komatitsch等[14]將C-PML邊界條件推廣到彈性介質(zhì)。近年來(lái)，中國(guó)學(xué)者也對(duì)完全匹配層做了許多研究[15-22]，王守東[23]給出了聲波模擬中PML邊界條件的基本原理；秦臻等[24]使用輔助變量推導(dǎo)了微分形式下的PML邊界條件，在此基礎(chǔ)上，熊章強(qiáng)等[25]采用CFS技術(shù)提出了輔助微分方程完全匹配層邊界條件。

在完全匹配層中，衰減因子占有十分重要的地位，決定著匹配層的吸收能力。Collino等[26]提出了與內(nèi)邊界距離呈指數(shù)關(guān)系的衰減函數(shù)，后經(jīng)Groby等[27]改進(jìn)整理，一直沿用至今。陳可洋[28]指出指數(shù)衰減函數(shù)隨著階數(shù)的變化，在內(nèi)邊界處過(guò)于平緩或者劇烈，不利于波的傳入，并提出了正弦型和余弦型衰減函數(shù)，但效果有限。理論上指數(shù)型和正余弦型衰減函數(shù)都能完全吸收入射波，但離散之后，反射是由每一層的變化差以及整體邊界層的反射系數(shù)共同決定，因此要進(jìn)一步提升匹配層的吸收能力則十分困難。為此，本文設(shè)計(jì)一種更為靈活且能夠控制衰減因子在各層增速的衰減函數(shù)，在不增加PML層數(shù)且不減小衰減因子最大值的情況下，壓制離散差異帶來(lái)的虛假反射，提升邊界的吸收能力。

2 理論

2.1 近似完全匹配層邊界條件

本文基于速度—應(yīng)力方程組采用時(shí)間二階、空間十二階的交錯(cuò)網(wǎng)格有限差分法進(jìn)行正演模擬。二維頻率域一階彈性波動(dòng)方程為

(1)

式中：ω是圓頻率；λ、μ是拉梅常數(shù)；ρ是密度；Vx、Vz為頻率域彈性波場(chǎng)速度分量；Txx、Tzz、Txz為頻率域彈性波場(chǎng)應(yīng)力分量。復(fù)坐標(biāo)伸展變換(CCS)的形式為

(2)

有

(3)

式中α是衰減函數(shù)。將式(3)代入式(1)，可得復(fù)坐標(biāo)伸展變換后的方程組

(4)

然后將上式進(jìn)行分裂，將應(yīng)力與速度分裂為x方向和z方向分量，并轉(zhuǎn)換到時(shí)間域，得到在PML計(jì)算區(qū)域波場(chǎng)迅速衰減的方程組

(5a)

(5b)

(5c)

(5d)

(5e)

式中：上標(biāo)“x”、“z”分別為波場(chǎng)沿x和z方向的分解量；vx、vz為時(shí)間域彈性波場(chǎng)速度分量；τxx、τzz、τxz為時(shí)間域彈性波場(chǎng)應(yīng)力分量。式(5)即為最早經(jīng)典的分裂PML邊界條件。該方法將五個(gè)變量分解為十個(gè)，方程的數(shù)量也為原來(lái)的兩倍，計(jì)算效率低，占用存儲(chǔ)空間大。而近似完全匹配層是對(duì)式(4)進(jìn)行變換，有

(6a)

(6b)

(7)

(8)

對(duì)其他變量采取相同的方法并變換到時(shí)間域中，有

(9)

式中

(10)

式(9)即為時(shí)間域中基于近似完全匹配層的波動(dòng)方程。對(duì)比式(9)與式(5)，二者形式一致，但前者不需采取分裂的方式實(shí)現(xiàn)完全匹配層，只需要引入8個(gè)輔助變量來(lái)替換原主控方程中的變量，而且這8個(gè)輔助變量可通過(guò)式(10)的微分方程求取。在程序編寫(xiě)時(shí)不必改動(dòng)原代碼，實(shí)現(xiàn)簡(jiǎn)單。對(duì)于這種不需要改變方程形式的完全匹配層也就十分容易推廣到其他復(fù)雜介質(zhì)。而分裂完全匹配層需將變量分裂為兩個(gè)或者多個(gè)分量并改變方程形式，實(shí)現(xiàn)起來(lái)十分復(fù)雜且計(jì)算效率低的缺點(diǎn)被放大。

2.2 吸收機(jī)理及衰減因子研究

采用復(fù)坐標(biāo)伸展變換之后轉(zhuǎn)換到時(shí)間域中的波動(dòng)方程的解(以x方向?yàn)槔?為

(11)

(12)

式中：L為完全匹配層的厚度；R為理論邊界反射系數(shù)；vP為縱波速度；l為PML區(qū)計(jì)算點(diǎn)與內(nèi)部邊界的距離。Groby等[27]改進(jìn)的衰減函數(shù)為

(13)

式(13)是式(12)的拓展，當(dāng)n=2時(shí)二者等價(jià)。當(dāng)取n=1時(shí)式(13)是線性函數(shù)，在內(nèi)邊界處的衰減梯度較大，易產(chǎn)生邊界反射。當(dāng)n值增加時(shí)內(nèi)邊界衰減因子變化緩慢、在外邊界處迅速增加，也不利于反射的吸收。陳可洋[28]提出的正弦型和余弦型衰減函數(shù)提升了內(nèi)邊界處的增速，抑制了外邊界處的突變。余弦型衰減函數(shù)為

(14)

正弦型衰減函數(shù)為

(15)

式中B為衰減幅度因子。B值不宜過(guò)大，否則阻礙了波從內(nèi)計(jì)算區(qū)域傳播到PML區(qū)域。

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

圖2為指數(shù)型衰減函數(shù)和余弦型衰減函數(shù)曲線及其梯度曲線的對(duì)比。由圖可見(jiàn)：①在相同的匹配層層數(shù)情況下，隨著n的增加，指數(shù)型衰減因子在內(nèi)邊界處變小，在外邊界附近增長(zhǎng)劇烈。②余弦型衰減函數(shù)在內(nèi)邊界處的增加速度介于n=1與n=2指數(shù)型衰減因子之間。在內(nèi)邊界處，余弦型函數(shù)的增速較快，在外邊界處增速較慢。③對(duì)于指數(shù)型函數(shù)通過(guò)增加n的值使內(nèi)邊界與計(jì)算區(qū)域接合更好，但在末端突變較大。④當(dāng)匹配層層數(shù)縮減時(shí)，如果保持衰減函數(shù)的系數(shù)不變，余弦型衰減函數(shù)最大值是不變的，這樣余弦型函數(shù)的梯度整體增加，層與層之間的離散差異增大。⑤對(duì)于指數(shù)型衰減函數(shù)，隨著層數(shù)的減少，由于分母L存在，α0的最大值是增大的，使指數(shù)型函數(shù)離散差異變得更為劇烈。

圖2 L=30(a)和L=20(b)的衰減因子曲線(上)、衰減梯度曲線(下)對(duì)比

為了更為靈活地控制衰減因子在各層的增速，需要找到一個(gè)函數(shù)作為修正項(xiàng)，要求該函數(shù)的梯度值先增加再減小并且能夠根據(jù)系數(shù)的改變控制其最大值的分布，以便控制衰減函數(shù)在各部分的梯度值。再將這個(gè)修正函數(shù)與原指數(shù)衰減函數(shù)相結(jié)合，可得修正后的指數(shù)型衰減函數(shù)

α(x)=K[β(l/L)n+γexp(-Lδ/l)]

(16)

上式第一部分β(l/L)n使該函數(shù)滿足完全匹配層中衰減函數(shù)的要求，第二部分γexp(-Lδ/l)作為修正函數(shù)控制衰減函數(shù)在各完全匹配層之間的變化調(diào)整，其中K、γ、β和δ是可調(diào)整的系數(shù)。如圖3所示，該修正函數(shù)通過(guò)改變?chǔ)牡闹悼刂浦荻茸畲簏c(diǎn)的位置，δ值越小，最大值點(diǎn)越接近內(nèi)邊界，且最大值越大。這樣的特性能夠很好地提升衰減梯度過(guò)小處的增速，特別是對(duì)衰減函數(shù)接近內(nèi)邊界區(qū)域的控制。當(dāng)δ的值不變，改變?chǔ)弥担弥翟酱?，梯度最大值越大，以此控制梯度的最大值。該修正函?shù)不用擔(dān)心外邊界處的增速過(guò)大，因?yàn)樵谕膺吔缣幤涮荻榷际菧p小的，正好起到了壓制衰減函數(shù)在外邊界處增速過(guò)快的效果?；谶@樣的特性，將修正函數(shù)與傳統(tǒng)的指數(shù)型衰減函數(shù)或者余弦衰減函數(shù)相加，適當(dāng)?shù)貙?duì)衰減因子做出細(xì)微的改變，以達(dá)到更好的吸收效果。

修正后衰減函數(shù)依舊保留了指數(shù)衰減函數(shù)的特點(diǎn)(圖4)，n越大，內(nèi)邊界因子值增速越小，而外邊界處變化劇烈。在n=2時(shí)，其梯度曲線不再是一條直線。當(dāng)增大β值時(shí)，衰減函數(shù)值整體隨之增加(圖5)，與增大衰減函數(shù)α0的值效果相同。當(dāng)增大γ值時(shí)(圖6)，衰減因子梯度值在匹配層中間部分得到提升，其位置還與δ的值有關(guān)。在外邊界附近，無(wú)論β值與γ值如何取值，梯度在外邊界處都有收縮的趨勢(shì)(圖7)，與傳統(tǒng)衰減函數(shù)相比，一定程度地壓制了衰減因子在邊界處的增速。

圖3 不同參數(shù)時(shí)的修正函數(shù)曲線(上)及其梯度曲線(下)對(duì)比

圖4 n變化、其他參數(shù)不變時(shí)修正后的衰減函數(shù)曲線(左)及其梯度曲線(右)

圖5 β變化、其他參數(shù)不變時(shí)修正后的衰減函數(shù)曲線(左)及其梯度曲線(右)

圖6 γ變化、其他參數(shù)不變時(shí)修正后的衰減函數(shù)曲線(左)及其梯度曲線(右)

圖7 δ變化、其他參數(shù)不變時(shí)修正后的衰減函數(shù)曲線(左)及其梯度曲線(右)

3 數(shù)值模擬分析

本文采用交錯(cuò)網(wǎng)格有限差分進(jìn)行正演模擬，時(shí)間步長(zhǎng)為1ms，震源是主頻為25Hz的雷克子波。模型尺寸為2000m×2000m(不包括邊界厚度)，兩個(gè)方向網(wǎng)格步長(zhǎng)均為10m，縱波速度為3000m/s，橫波速度為1400m/s，密度為2000kg/m3，PML層數(shù)L分別設(shè)為5、10、15和20，純縱波震源加載于(1000m，1000m)處。R取為10-6，余弦衰減函數(shù)中的B值取220。為了方便與指數(shù)型衰減函數(shù)做比較，將修正后的函數(shù)與指數(shù)型衰減函數(shù)取相同的系數(shù)。令

(17)

依次改變?chǔ)煤挺牡闹狄约巴耆ヅ鋵訉訑?shù)，分析其對(duì)邊界吸收效果的影響。同時(shí)將修正函數(shù)加入余弦衰減函數(shù)

(18)

當(dāng)層數(shù)取5、10、15和20的時(shí)候，記錄不同系數(shù)情況下的邊界反射的振幅，以衡量PML邊界的吸收效果。

圖8是時(shí)間采樣間隔為1ms、γ的取值范圍為0～0.08、δ的取值范圍為0～0.16時(shí)應(yīng)力τxx的邊界反射振幅，可見(jiàn)： ①當(dāng)γ為零時(shí)，顯示的是沒(méi)有經(jīng)過(guò)修正的傳統(tǒng)衰減函數(shù)下的邊界反射波振幅，此時(shí)的反射波能量很大；邊界反射振幅隨著γ值的增加先減小后增大，在某一個(gè)范圍內(nèi)吸收效果最好，而該范圍會(huì)隨著PML層數(shù)的改變而移動(dòng)。②無(wú)論指數(shù)型衰減函數(shù)還是余弦型衰減函數(shù)，隨著匹配層層數(shù)的增加，吸收效果最好區(qū)間的δ值先增加再減小，而γ值在余弦型函數(shù)中逐漸減小，在指數(shù)型函數(shù)中先增加再降低。③如圖2所示，因?yàn)樗p函數(shù)最大值是不變的，層數(shù)越少，衰減因子在內(nèi)邊界以外的區(qū)域增長(zhǎng)越快，兩點(diǎn)之間的離散差異變大。而修正函數(shù)能夠提升衰減因子在內(nèi)邊界處的增速，縮小其他區(qū)域的離散差異。因此當(dāng)層數(shù)為5時(shí)，吸收效果最好的點(diǎn)都在δ=0的位置。④當(dāng)層數(shù)增加時(shí)，新的衰減函數(shù)通過(guò)增加δ值和減小γ值降低內(nèi)邊界處的提升速度，并使修正函數(shù)梯度最大值遠(yuǎn)離內(nèi)邊界區(qū)域，提高匹配層中間位置梯度值。使衰減梯度值合理分布，以此達(dá)到更好的吸收效果。⑤當(dāng)取20層時(shí)，由于離散差異引起的誤差反射十分微弱，需要修正的強(qiáng)度也隨之減弱，因此γ取值在20層時(shí)是最小的。

圖9給出了在不同層數(shù)、不同衰減函數(shù)下地震波在0.5s時(shí)的vx分量波場(chǎng)快照。由圖可見(jiàn)余弦型衰減函數(shù)比指數(shù)型衰減函數(shù)的吸收效果要好一些，修正余弦衰減函數(shù)和修正指數(shù)型衰減函數(shù)吸收效果更好。

為了更清晰地分析邊界的衰減效果，給出了模型整體區(qū)間的能量隨時(shí)間變化曲線(圖10)。由圖10可見(jiàn)，地震波在0.34s左右到達(dá)邊界位置，然后迅速衰減，因此能量曲線在這一時(shí)刻驟降。殘余的反射波返回計(jì)算區(qū)域，能量已經(jīng)得到很大程度的削弱，當(dāng)時(shí)間為1s左右時(shí)，反射波再次到達(dá)邊界處，能量再次被吸收，幾乎將能量吸收完全。對(duì)于余弦衰減函數(shù)，當(dāng)層數(shù)較少時(shí)出現(xiàn)了吸收不穩(wěn)定現(xiàn)象，吸收能力較弱，而指數(shù)型衰減函數(shù)沒(méi)有此現(xiàn)象。就能量衰減情況來(lái)看，當(dāng)PML層數(shù)為10和15時(shí)，余弦型衰減函數(shù)的吸收效果優(yōu)于指數(shù)型衰減函數(shù)。當(dāng)PML層數(shù)為20時(shí)兩種邊界條件吸收效果相當(dāng)。由圖中虛線所示的衰減函數(shù)修正后的能量衰減曲線可見(jiàn)，無(wú)論P(yáng)ML層數(shù)為多少，吸收效果都得到了極大的改進(jìn)，其對(duì)應(yīng)的δ和γ的取值為圖8中吸收效果最佳點(diǎn)。

圖8 余弦型(a)和指數(shù)型(b)衰減因子的吸收效果分析

圖9 應(yīng)用不同衰減因子PML邊界條件的0.5s時(shí)刻波場(chǎng)vx分量模擬快照

對(duì)比(1000m，300m)處使用不同衰減函數(shù)下的應(yīng)力τxx分量地震模擬記錄，并求取使用傳統(tǒng)衰減函數(shù)與修正衰減函數(shù)模擬記錄τxx分量的差值曲線(圖11)。對(duì)于傳統(tǒng)的衰減函數(shù)，余弦型的吸收效果較好，只是當(dāng)PML層數(shù)較少時(shí)在0.9s時(shí)接收到了十分強(qiáng)烈的虛假反射。差值曲線多分布于正半軸，且值較大，表明修正衰減函數(shù)壓制虛假反射的效果十分明顯。采用修正衰減函數(shù)PML邊界條件的邊界反射振幅最弱，PML層數(shù)越少，壓制效果越明顯，這是因?yàn)閷訑?shù)的增加減弱了離散差異帶來(lái)的誤差。為對(duì)修正之后壓制效果進(jìn)行定量描述，計(jì)算各系數(shù)組下修正衰減函數(shù)與傳統(tǒng)衰減函數(shù)的邊界反射振幅比，使其歸一化，結(jié)果如圖12所示。當(dāng)γ=0時(shí)為傳統(tǒng)衰減函數(shù)下的反射強(qiáng)度且值為100%，隨著γ的增大，在20層時(shí)比值減小了20%，隨著PML層數(shù)的減小比值甚至能夠減小60%，不到原來(lái)的一半。

圖10 應(yīng)用不同衰減因子PML邊界條件的能量衰減曲線

圖11 使用不同衰減函數(shù)PML邊界條件的(1000m，300m)處振動(dòng)曲線以及差值曲線對(duì)比

圖12 不同類型衰減函數(shù)的修正前、后邊界反射振幅比曲線

4 結(jié)束語(yǔ)

本文在理論吸收系數(shù)不變以及固定的PML層數(shù)情況下，針對(duì)離散差異所導(dǎo)致的誤差，設(shè)計(jì)了合理修正函數(shù)，使修正后的衰減函數(shù)在層數(shù)少時(shí)增加內(nèi)邊界處衰減函數(shù)的梯度值，減小在外邊界處的負(fù)擔(dān)，在層數(shù)較多時(shí)又能減小在內(nèi)邊界處的梯度值并將梯度最大值外移，使梯度值在各部分合理分布。正演模擬結(jié)果顯示，修正后的衰減函數(shù)大大提升了PML邊界的吸收能力。