王紹文 宋 鵬*②③ 譚 軍②③ 解 闖 毛士博 王倩倩

(①中國(guó)海洋大學(xué)海洋地球科學(xué)學(xué)院，山東青島 266100；②青島國(guó)家海洋科學(xué)與技術(shù)實(shí)驗(yàn)室，山東青島 266100；③中國(guó)海洋大學(xué)海底科學(xué)與探測(cè)技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室，山東青島 266100)

0 引言

地震波正演模擬技術(shù)是研究地震波傳播規(guī)律的有效手段，并在地震勘探的采集、處理以及反演等各個(gè)環(huán)節(jié)中都有著廣泛的應(yīng)用。地震波正演模擬必須引入人工邊界界定計(jì)算區(qū)域，而人工邊界的處理不當(dāng)，往往會(huì)產(chǎn)生虛假反射污染中心波場(chǎng)，從而降低波場(chǎng)模擬的精度，因而人工邊界的處理一直是波動(dòng)方程數(shù)值計(jì)算的重要研究?jī)?nèi)容。

當(dāng)前的人工邊界處理方法主要有三類。第一類是衰減邊界方法[1-3]。該類方法是在靠近邊界處引入一個(gè)衰減區(qū)，在衰減區(qū)內(nèi)地震波傳播遵循阻尼波動(dòng)方程，能量按指數(shù)衰減。由于該類方法存在阻尼因子難以確定、計(jì)算量大和內(nèi)存消耗較大等問(wèn)題，因此應(yīng)用較少。第二類是完全匹配層(PML)方法。該類方法由Berenger[4]在電磁波數(shù)值模擬中首先提出，通過(guò)在中心波場(chǎng)計(jì)算區(qū)域外加上一系列的吸收層、層內(nèi)引入衰減因子以達(dá)到消除邊界反射的目的。王守東[5]、Hastings等[6]、Collino等[7]將PML法成功應(yīng)用于聲波和彈性波正演模擬。隨后， Komatitsch等[8-9]、熊章強(qiáng)等[10]實(shí)現(xiàn)了非分裂PML法；羅玉欽等[11]、楊茜娜等[12]發(fā)展了復(fù)頻移及多軸復(fù)頻移PML法；陳可洋[13]和Wang等[14]提出了基于余弦型和基于高斯型衰減函數(shù)的PML方法。PML法現(xiàn)已在地震波數(shù)值模擬中獲得廣泛應(yīng)用，然而為達(dá)到更好的吸收效果，PML法往往需要設(shè)計(jì)幾十甚至上百個(gè)匹配層，意味著巨大的計(jì)算量和內(nèi)存消耗，嚴(yán)重限制了該方法在大規(guī)模波場(chǎng)模擬中的應(yīng)用。第三類為基于單程波近似的吸收邊界條件。該類方法由Enguist等[15-16]首先提出；Higdon[17-19]提出了一種與其等價(jià)的漸進(jìn)吸收邊界條件，可以設(shè)計(jì)吸收任意角度的入射波，顯著提升了吸收精度和適應(yīng)性。然而由于單程波的高階近似會(huì)導(dǎo)致邊界條件方程中出現(xiàn)高階的混合時(shí)間、空間偏導(dǎo)數(shù)項(xiàng)，給數(shù)值求解帶來(lái)困難，因此常規(guī)方法多采用二階邊界條件，但對(duì)大角度入射波的吸收效果較差，其整體吸收效果也難以滿足高精度數(shù)值模擬的要求。為進(jìn)一步提升基于單程波近似的吸收邊界條件的效果，宋鵬等[20-21]發(fā)展了優(yōu)化系數(shù)的四階吸收邊界條件算法。Baffet等[22]、Hagstrom等[23]以及Potter等[24]也通過(guò)引入輔助變量實(shí)現(xiàn)了高階吸收邊界條件，但這些算法同樣存在計(jì)算過(guò)程復(fù)雜且計(jì)算效率低等問(wèn)題。

為提升常規(guī)基于單程波近似的吸收邊界條件的效果，在不提高階數(shù)的前提下，Liu等[25]提出了一種混合吸收邊界條件。該方法在中心波場(chǎng)與邊界之間增加一個(gè)過(guò)渡區(qū)域，通過(guò)線性加權(quán)使過(guò)渡區(qū)域內(nèi)雙程波波場(chǎng)與單程波波場(chǎng)平滑過(guò)渡，最終實(shí)現(xiàn)了人工邊界的高效吸收。隨后，Liu等[26-27]、任志明等[28]以及Liu等[29]實(shí)現(xiàn)了三維聲波及彈性波正演模擬中的混合吸收邊界，Ren等[30]將混合吸收邊界推廣到頻率域數(shù)值模擬中，取得了較好的吸收效果。劉學(xué)義等[31]將混合吸收邊界用于各向異性介質(zhì)，Wang等[32]、薛浩等[33]分別將混合吸收邊界應(yīng)用于逆時(shí)偏移和最小二乘逆時(shí)偏移。

當(dāng)前混合吸收邊界已在地震波場(chǎng)正反演延拓中得到廣泛應(yīng)用，然而，該方法的加權(quán)系數(shù)仍以線性和指數(shù)型[34]為主，未能兼顧內(nèi)邊界和外邊界的綜合吸收效果，在吸收效率上并沒(méi)有達(dá)到最優(yōu)；此外，當(dāng)前大規(guī)模的波場(chǎng)模擬往往借助于GPU并行提高計(jì)算效率，但混合吸收邊界差分格式的GPU加速效率遠(yuǎn)遠(yuǎn)低于中心波場(chǎng)，嚴(yán)重影響了整個(gè)波場(chǎng)模擬的計(jì)算效率。本文提出一種基于高斯型加權(quán)系數(shù)的混合吸收邊界條件方法，其更能兼顧內(nèi)、外吸收邊界的吸收效果，具有更高效的整體邊界吸收效率；同時(shí)本文還推導(dǎo)了一種更適用于GPU加速的一階Higdon混合吸收邊界條件差分格式，可顯著提升邊界計(jì)算的GPU并行效率。

1 基于高斯型加權(quán)系數(shù)的混合吸收邊界條件

1.1 彈性波數(shù)值模擬中混合吸收邊界原理

二維空間中的一階速度—應(yīng)力彈性波方程[35]可以表示為

(1)

式中：vx、vz分別為水平方向和垂直方向波場(chǎng)速度分量；σxx、σzz、σxz為三個(gè)應(yīng)力分量；ρ為密度，λ和μ為拉梅系數(shù)，與縱波速度VP和橫波速度VS關(guān)系為

(2)

(3)

式(1)的交錯(cuò)網(wǎng)格有限差分格式為[36-37]

(4)

(5)

(6)

(7)

(8)

式中：U和W分別表示速度分量vx、vz的離散形式；P、Q和R分別表示應(yīng)力分量σxx、σzz、σxz的離散形式；Δt、Δx、Δz分別為離散的時(shí)間步長(zhǎng)、空間水平和垂直網(wǎng)格間隔；i、j、k為空間和時(shí)間離散序號(hào)；2N為差分階次；am(1≤m≤N)為差分系數(shù)，可通過(guò)求解下式方程獲取[38]

(9)

應(yīng)用混合吸收邊界方法進(jìn)行數(shù)值模擬時(shí)，整個(gè)計(jì)算區(qū)域被劃分為中心波場(chǎng)區(qū)域Ⅰ、邊界過(guò)渡區(qū)域 Ⅱ 以及邊界區(qū)Ⅲ三部分(圖1)，其波場(chǎng)計(jì)算流程如下：

圖1 二維混合吸收邊界條件示意圖

(3)對(duì)于M層區(qū)域Ⅱ，對(duì)雙程波場(chǎng)值和單程波場(chǎng)值進(jìn)行加權(quán)求和得到混合波場(chǎng)

ul=(1-wl)ut+wlu°l=1，2，…，M

(10)

式中wl為加權(quán)系數(shù)。

而對(duì)于區(qū)域Ⅱ和區(qū)域Ⅲ的單程波場(chǎng)計(jì)算，由于二階Higdon吸收邊界條件穩(wěn)定性較差，通常采用改進(jìn)的一階Higdon吸收邊界條件[28]，其左邊界方程為

(11)

(12)

則可將一階Higdon左邊界算子表示為[19]

(13)

其差分格式為

(14)

式中與空間位置和網(wǎng)格剖分相關(guān)的常量為[19]

(15)

(16)

(17)

其中b為常數(shù)，通常取值為0.5。

混合邊界的角點(diǎn)區(qū)域需做特殊處理，以左上角速度水平分量為例，其角點(diǎn)差分格式為

(18)

1.2 基于高斯型加權(quán)系數(shù)的混合吸收邊界

混合邊界通過(guò)引入過(guò)渡區(qū)實(shí)現(xiàn)波場(chǎng)混疊壓制邊界反射，過(guò)渡區(qū)與中心波場(chǎng)相鄰的邊界可稱為內(nèi)邊界，過(guò)渡區(qū)的最外層邊界可稱為外邊界，混合邊界的殘余反射通常來(lái)自于內(nèi)邊界和外邊界反射。內(nèi)、外邊界的殘余反射能量可通過(guò)混合加權(quán)系數(shù)的選取進(jìn)行調(diào)整，常用的線性加權(quán)系數(shù)和指數(shù)型加權(quán)系數(shù)[34]分別為

(19)

(20)

線性加權(quán)系數(shù)(圖2黑線所示)在內(nèi)、外邊界處梯度變化相同，而由于內(nèi)邊界更靠近內(nèi)部波場(chǎng)，因此應(yīng)用線性加權(quán)系數(shù)的混合邊界其通常會(huì)產(chǎn)生較強(qiáng)的內(nèi)邊界反射；指數(shù)型加權(quán)系數(shù)(圖2綠線所示)在內(nèi)邊界處變化更加緩慢，使內(nèi)邊界與中心波場(chǎng)區(qū)域得到更好的耦合，因此其內(nèi)邊界殘余反射得到有效壓制，但另一方面，由于其在內(nèi)邊界處的變化過(guò)于緩慢導(dǎo)致外邊界處的變化過(guò)于劇烈，從而易產(chǎn)生較強(qiáng)的外邊界反射。因此一個(gè)更為理想的混合邊界條件必須全面考慮地震波在內(nèi)、外邊界處的綜合邊界反射壓制效果。

鑒于以上分析，本文提出了一種高斯型混合加權(quán)系數(shù)

(21)

式中α為高斯型權(quán)系數(shù)的階數(shù)，控制內(nèi)、外邊界處的變化率。大量的實(shí)驗(yàn)證明，當(dāng)α=2.5時(shí)，高斯型混合邊界既能保持穩(wěn)定性，又可達(dá)到良好的邊界吸收效果。

高斯型加權(quán)系數(shù)(圖2紅線)在內(nèi)邊界處梯度變化雖然比指數(shù)型加權(quán)系數(shù)稍顯劇烈，但其較線性加權(quán)系數(shù)已更為平緩；而外邊界處，高斯型加權(quán)系數(shù)又比指數(shù)型加權(quán)系數(shù)更平緩，較線性加權(quán)系數(shù)稍劇烈，因此理論上其應(yīng)有更優(yōu)的內(nèi)、外邊界殘余反射綜合壓制效果。

當(dāng)然，采用混合吸收邊界時(shí)，整個(gè)邊界的吸收效果還與過(guò)渡區(qū)層數(shù)有關(guān)，Liu等[27]證明了當(dāng)過(guò)渡區(qū)厚度較小時(shí)，混合邊界易不穩(wěn)定，并且考慮到為追求更高精度的波場(chǎng)模擬效果，當(dāng)前數(shù)值模擬時(shí)通常大多采用20層以上的混合吸收邊界。因此本文將以20層和30層為例討論高斯型混合加權(quán)系數(shù)的吸收效果。

圖2 線性、指數(shù)型、高斯型加權(quán)系數(shù)變化曲線

1.3 常速模型數(shù)值模擬實(shí)驗(yàn)

為驗(yàn)證本文提出的高斯型混合吸收邊界條件的吸收效果，本文設(shè)計(jì)3200m×3200m、VP=4000m/s、VS=2000m/s、ρ= 2100kg/m3的均勻介質(zhì)模型。縱、橫向空間網(wǎng)格間距均為8m，時(shí)間步長(zhǎng)為1ms，差分格式精度為空間10階。震源位于模型中心，采用主頻為25Hz的Ricker子波加載于速度垂直分量上，記錄長(zhǎng)度為2s。

應(yīng)用20層混合吸收邊界得到1s時(shí)刻的vx、vz分量波前快照如圖3所示，波前快照的水平和垂直切片(位置見(jiàn)圖3紅色虛線)如圖4所示。應(yīng)用30層混合吸收邊界得到1s時(shí)刻的vx、vz分量波前快照如圖5所示，波前快照的水平和垂直切片(位置見(jiàn)圖5紅色虛線)如圖6所示。

由圖3～圖6可見(jiàn)，無(wú)論是20層還是30層吸收邊界，對(duì)于vx、vz分量，線性加權(quán)系數(shù)其波前快照中仍存在較強(qiáng)的內(nèi)邊界殘余反射(圖3左、圖5左箭頭所示)，指數(shù)型加權(quán)系數(shù)其波前快照中則存在著明顯的外邊界殘余反射(圖3中、圖5中箭頭所示)，而高斯型加權(quán)系數(shù)的波前快照中無(wú)明顯邊界殘余反射波。

圖3 20層混合吸收邊界不同加權(quán)系數(shù)的1s時(shí)刻vx分量(a)、vz分量(b)波場(chǎng)模擬快照對(duì)比 左：線性；中：指數(shù)型；右：高斯型

圖4 20層混合吸收邊界不同加權(quán)系數(shù)的1s時(shí)刻vx分量(a)、vz分量(b)波場(chǎng)模擬切片對(duì)比

圖5 30層混合吸收邊界不同加權(quán)系數(shù)的1s時(shí)刻vx分量(a)、vz分量(b)波場(chǎng)模擬快照對(duì)比 左：線性；中：指數(shù)型；右：高斯型

圖6 30層混合吸收邊界不同加權(quán)系數(shù)的1s時(shí)刻vx分量(a)、vz分量(b)波場(chǎng)模擬切片對(duì)比

2 適用于GPU加速的一階Higdon混合吸收邊界條件差分格式

當(dāng)前大規(guī)模的波場(chǎng)模擬往往借助于GPU并行提高計(jì)算效率。但應(yīng)用于大規(guī)模波場(chǎng)模擬的混合邊界差分格式的GPU加速比遠(yuǎn)遠(yuǎn)低于中心波場(chǎng)計(jì)算的GPU加速比，嚴(yán)重影響了整個(gè)波場(chǎng)模擬的計(jì)算效率。本文推導(dǎo)了一種新的混合邊界差分計(jì)算格式，可顯著提高邊界計(jì)算的GPU加速比。

2.1 適用于GPU并行的混合吸收邊界差分格式

由式(14)可知，基于常規(guī)混合吸收邊界差分格式進(jìn)行GPU并行時(shí)，由于每層邊界上k+1時(shí)刻的波場(chǎng)值計(jì)算都要用到該層內(nèi)側(cè)點(diǎn)k+1時(shí)刻的波場(chǎng)值，所以應(yīng)用GPU計(jì)算邊界時(shí)無(wú)法同時(shí)啟動(dòng)大量的線程實(shí)現(xiàn)邊界各點(diǎn)值的同時(shí)計(jì)算，無(wú)法實(shí)現(xiàn)高效的GPU并行。為了解決常規(guī)邊界差分格式中需要由內(nèi)向外依次計(jì)算、無(wú)法實(shí)現(xiàn)GPU高效并行的問(wèn)題，本文推導(dǎo)了一種新的邊界條件差分格式，具體推導(dǎo)過(guò)程如下(以左邊界和左上角點(diǎn)為例)。

為避免邊界區(qū)域計(jì)算時(shí)，差分格式的右側(cè)不出現(xiàn)k+1時(shí)刻項(xiàng)，將時(shí)間偏導(dǎo)數(shù)和水平方向空間偏導(dǎo)數(shù)分別表示為

(22)

(23)

則可將一階Higdon邊界算子表示為

(24)

(25)

對(duì)于左上角點(diǎn)，可通過(guò)上邊界和左邊界的加權(quán)平均得到

(26)

對(duì)比式(14)與式(25)可知，當(dāng)計(jì)算k+1時(shí)刻波場(chǎng)值時(shí)，本文推導(dǎo)的差分格式只需用到已計(jì)算出的k時(shí)刻波場(chǎng)值，因此同一時(shí)間切片內(nèi)所有網(wǎng)格點(diǎn)可同時(shí)計(jì)算，能實(shí)現(xiàn)GPU的高效并行加速。

在常規(guī)混合吸收邊界算法中，由于第l層邊界計(jì)算與第l+1層邊界計(jì)算存在計(jì)算先后順序關(guān)系，過(guò)渡區(qū)Ⅱ中的混合波場(chǎng)計(jì)算需要由最內(nèi)側(cè)開(kāi)始逐層向外進(jìn)行，因此同一時(shí)間層內(nèi)，前者需多次啟動(dòng)邊界計(jì)算核函數(shù)(與邊界層數(shù)相當(dāng))，且每次核函數(shù)只能計(jì)算本吸收層內(nèi)的網(wǎng)格點(diǎn)數(shù)，嚴(yán)重降低了GPU的加速效率?；诒疚牟罘指袷降幕旌衔者吔缢惴?，只需一次啟動(dòng)和計(jì)算邊界計(jì)算核函數(shù)，可實(shí)現(xiàn)所有吸收層內(nèi)的所有網(wǎng)格點(diǎn)的同時(shí)計(jì)算，更適合于大規(guī)模波場(chǎng)模擬的GPU高效并行計(jì)算。

2.2 常速模型數(shù)值實(shí)驗(yàn)

應(yīng)用上述均勻模型和模擬參數(shù)，驗(yàn)證本文推導(dǎo)的混合吸收邊界差分格式的正確性。20、30層常規(guī)差分格式與本文差分格式vx、vz分量波前快照如圖7、圖8所示。

由圖7、圖8可以看出，在邊界層數(shù)相同的條件下，兩種差分格式的混合吸收邊界均能有效壓制邊界反射。

圖7 不同邊界差分格式的20層混合吸收邊界數(shù)值模擬的0.6(左)、1.0(中)和1.4s(右)時(shí)刻波前快照對(duì)比(a)常規(guī)差分格式vx分量； (b)常規(guī)差分格式vz分量； (c)本文差分格式vx分量； (d)本文差分格式vz分量

為評(píng)估兩種差分格式的計(jì)算效率，本文基于RTX-2080Ti-11GB的GPU應(yīng)用上述均勻模型進(jìn)行測(cè)試，結(jié)果如表1所示。

由表1可見(jiàn)，若采用20層邊界，本文邊界差分格式計(jì)算效率為常規(guī)差分格式3.25倍；若采用30層邊界，本文邊界差分格式計(jì)算效率為常規(guī)差分格式4.34倍，因此本文差分格式更適于大規(guī)模波場(chǎng)模擬的GPU高效并行計(jì)算。

表1 兩種邊界差分格式正演模擬計(jì)算效率對(duì)比

圖8 不同邊界差分格式30層混合吸收邊界數(shù)值模擬的0.6(左)、1.0(中)和1.4s(右)時(shí)刻波前快照對(duì)比(a)常規(guī)差分格式vx分量； (b)常規(guī)差分格式vz分量； (c)本文差分格式vx分量； (d)本文差分格式vz分量

3 復(fù)雜模型實(shí)驗(yàn)

本文截取了Marmousi模型的一部分(圖9)進(jìn)行正演模擬測(cè)試。該模型尺寸為3000m×3000m，縱橫波速比為1.73，常密度為2100kg/m3。數(shù)值模擬中，橫向和縱向網(wǎng)格間隔分別為5m和4m，時(shí)間步長(zhǎng)為0.4ms。震源采用加載于vz分量上、主頻為25Hz的Ricker子波，位于(1500m，4m)處。

未加吸收邊界、常規(guī)差分格式的高斯型混合吸收邊界(30層)和本文差分格式的高斯型混合吸收邊界(30層)模擬的1.4s時(shí)vx、vz分量的波前快照如圖10所示。由圖可見(jiàn)，未加吸收邊界的波場(chǎng)中有明顯的強(qiáng)邊界反射，而無(wú)論采用常規(guī)邊界差分格式還是本文差分格式的30層混合吸收邊界，均可有效吸收邊界反射，表明基于高斯型混合加權(quán)吸收邊界能夠適應(yīng)于復(fù)雜構(gòu)造模型的邊界處理。

基于單卡GPU(型號(hào)：RTX-2080Ti-11GB)，常規(guī)差分格式的高斯型混合吸收邊界的部分Mar-mousi模型計(jì)算時(shí)間為16.6335s，本文差分格式計(jì)算時(shí)間為5.3510s，可見(jiàn)本文差分格式的計(jì)算效率是常規(guī)格式的3倍以上，說(shuō)明本文的邊界差分格式更適合大規(guī)模數(shù)值模擬的GPU高效并行。

圖9 部分Marmousi縱波速度模型

圖10 Marmousi模型不同邊界模擬的1.5s時(shí)刻的vx分量(左)、vz分量(右)波場(chǎng)快照(a)不加邊界； (b)常規(guī)差分格式的高斯型混合吸收邊界； (c)本文差分格式的高斯型混合吸收邊界

4 結(jié)論

本文在系統(tǒng)分析線性混合邊界和指數(shù)型混合邊界殘余反射形成機(jī)制的基礎(chǔ)上，提出了一種能兼顧內(nèi)外邊界吸收效果的高斯型混合吸收邊界條件。在此基礎(chǔ)上，推導(dǎo)了更適用于GPU加速的一階Higdon混合吸收邊界條件差分格式及相應(yīng)的角點(diǎn)計(jì)算差分格式。模型實(shí)驗(yàn)結(jié)果表明，高斯型混合吸收邊界能更有效地壓制內(nèi)邊界和外邊界反射，比線性和指數(shù)型混合吸收邊界具有更高的吸收效率，適用于復(fù)雜構(gòu)造模型的邊界處理；本文推導(dǎo)的邊界及其角點(diǎn)的邊界差分格式，GPU并行的計(jì)算效率更高，更適用于大規(guī)模彈性波場(chǎng)數(shù)值模擬。

將本文提出的混合吸收邊界算法推廣至三維彈性波數(shù)值模擬及逆時(shí)偏移和全波形反演，是今后的研究方向。

十分感謝中國(guó)海洋大學(xué)地球探測(cè)軟件技術(shù)研發(fā)團(tuán)隊(duì)提供了資金和軟、硬件支持，以及團(tuán)隊(duì)成員提供了建設(shè)性意見(jiàn)。