近似完全匹配層邊界條件吸收效果分析及衰減函數的改進

羅玉欽 劉 財

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

1 引言

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

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

2 理論

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

本文基于速度—應力方程組采用時間二階、空間十二階的交錯網格有限差分法進行正演模擬。二維頻率域一階彈性波動方程為

(1)

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

(2)

有

(3)

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

(4)

然后將上式進行分裂，將應力與速度分裂為x方向和z方向分量，并轉換到時間域，得到在PML計算區域波場迅速衰減的方程組

(5a)

(5b)

(5c)

(5d)

(5e)

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

(6a)

(6b)

(7)

(8)

對其他變量采取相同的方法并變換到時間域中，有

(9)

式中

(10)

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

2.2 吸收機理及衰減因子研究

采用復坐標伸展變換之后轉換到時間域中的波動方程的解(以x方向為例)為

(11)

(12)

式中：L為完全匹配層的厚度；R為理論邊界反射系數；vP為縱波速度；l為PML區計算點與內部邊界的距離。Groby等[27]改進的衰減函數為

(13)

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

(14)

正弦型衰減函數為

(15)

式中B為衰減幅度因子。B值不宜過大，否則阻礙了波從內計算區域傳播到PML區域。

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

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

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

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

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

(16)

上式第一部分β(l/L)n使該函數滿足完全匹配層中衰減函數的要求，第二部分γexp(-Lδ/l)作為修正函數控制衰減函數在各完全匹配層之間的變化調整，其中K、γ、β和δ是可調整的系數。如圖3所示，該修正函數通過改變δ的值控制著梯度最大點的位置，δ值越小，最大值點越接近內邊界，且最大值越大。這樣的特性能夠很好地提升衰減梯度過小處的增速，特別是對衰減函數接近內邊界區域的控制。當δ的值不變，改變γ值，γ值越大，梯度最大值越大，以此控制梯度的最大值。該修正函數不用擔心外邊界處的增速過大，因為在外邊界處其梯度都是減小的，正好起到了壓制衰減函數在外邊界處增速過快的效果?；谶@樣的特性，將修正函數與傳統的指數型衰減函數或者余弦衰減函數相加，適當地對衰減因子做出細微的改變，以達到更好的吸收效果。

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

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

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

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

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

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

3 數值模擬分析

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

(17)

依次改變γ和δ的值以及完全匹配層層數，分析其對邊界吸收效果的影響。同時將修正函數加入余弦衰減函數

(18)

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

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

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

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

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

圖9 應用不同衰減因子PML邊界條件的0.5s時刻波場vx分量模擬快照

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

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

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

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

4 結束語

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