李青陽,吳國忱,楊凌云,王玉梅

(1.中國石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院,山東青島266580;2.海洋國家實驗室海洋礦產(chǎn)資源評價與探測技術(shù)功能實驗室,山東青島266071;3中國石油化工股份有限公司勝利油田分公司物探研究院,山東東營257022)

基于雙程波動方程的逆時偏移方法可以處理不同種類的波動,如反射波、折射波、繞射波等,且不受成像傾角的限制,具備對高陡構(gòu)造進行高精度成像的能力[1-5]。逆時偏移在實際應(yīng)用中仍然受到一些因素的限制,由于波場延拓采用雙程波動方程,故對于存儲需求及I/O帶來很大的負擔。關(guān)于計算效率和內(nèi)存優(yōu)化的方法主要有兩大類,一類是高性能計算,另一類是改進算法。FOLTINEK等[6]給出了高階有限差分逆時偏移在GPU上的實現(xiàn)方法;SUH等[7]研究了基于集群的逆時偏移并行計算實施策略;MCGARRY等[8]通過計算域邊界上檢波點記錄的波場來近似重建震源波場。

正演模擬是逆時偏移的基礎(chǔ),正演算法的優(yōu)化能夠有效提高計算效率,采用波動方程方法進行地震波正演時,一般都需要增加人工邊界[9],將無限區(qū)域截斷成有限區(qū)域。目前有多種吸收邊界方法,例如基于傍軸近似的單程波動方程吸收邊界條件[9-10]、海綿吸收邊界條件[11-12]以及完全匹配層(perfectly matched layer,PML)吸收邊界條件[13]等。在以上吸收邊界條件中,完全匹配層吸收邊界條件具有更好的吸收效果。經(jīng)典PML是由復(fù)坐標伸展變換(complex coordinate stretching,CCS)推導(dǎo)給出的[14-15]。王守東[16]給出了聲波模擬中PML吸收邊界條件的基本原理。殷文等[17]將PML吸收邊界應(yīng)用于頻率域彈性波正演。王永剛等[18]研究推導(dǎo)了二維聲波方程PML吸收邊界條件,并與旁軸近似吸收邊界條件進行了比較,證明PML吸收邊界條件更好。陳可洋[19]對傳統(tǒng)的聲波完全匹配層進行改進,提出了改進的聲波分裂PML吸收邊界條件。COLLINO等[20]對一階方程實現(xiàn)PML的應(yīng)用,KOMATITSCH等[21]將其擴展到二階方程,并推導(dǎo)出分裂的PML公式(SPML)。然而CCS存在一定的缺陷,即無法吸收高角度的入射波和瞬逝波。KUZUOGLU等[22]提出復(fù)頻移(complex frequency shift,CFS)PML,能夠?qū)Ω呓嵌热肷洳ㄟM行衰減,BéRENGER[23-24]同樣驗證了CFS-PML吸收邊界條件對瞬逝波具有顯著的吸收效果。但是CFS-PML邊界條件需要對微分系統(tǒng)引入過多輔助變量,求解該系統(tǒng)需要計算大量卷積或者對波場進行分裂,計算代價巨大。秦臻等[25]在頻率域彈性波模擬中給出輔助變量推導(dǎo)了微分形式下的CFS-PML吸收邊界條件。熊章強等[26]提出了針對一階速度-應(yīng)力聲波方程數(shù)值模擬CFS-PML吸收邊界條件。RODEN等[27]提出卷積PML(CPML)吸收邊界條件并利用遞歸方法求解CFS-PML中的卷積項[28]。KOMATITSCH等[29]用CPML解決一階速度-應(yīng)力系統(tǒng)彈性波方程。楊凌云等[30]提出了一種新的二階系統(tǒng)下聲波方程的CPML(NCPML)吸收邊界條件,大大減少了內(nèi)存使用,提高了計算效率。張魯新等[31]將完全匹配層應(yīng)用于孔隙介質(zhì)中。PML吸收邊界條件被引入到多孔彈性介質(zhì)[32-33],各向異性介質(zhì)[34-37]以及黏彈性介質(zhì)[38-39]。PASALIC等[40]給出基于CPML吸收邊界條件的二階各向同性和各向異性聲波方程數(shù)值模擬方法,DROSSAERT等[41]提出遞歸積分法(RI-PML),以積分計算替代卷積計算。除了常用的有限差分法數(shù)值模擬,李青陽等[42]研究了高階偽譜法彈性波PML吸收邊界條件,提高正演計算精度。覃發(fā)兵等[43]研究了完全匹配層在時域有限元彈性波數(shù)值模擬中的應(yīng)用。此外,近似PML[44](nearly PML,NPML)與上述PML對坐標進行拉伸不同,NPML直接對波場進行拉伸變換無需卷積運算。劉守偉等[45]在三維逆時偏移GPU/CPU機群實現(xiàn)中引入NPML吸收邊界條件。羅玉欽等[46]提出多軸復(fù)頻移NPML吸收邊界條件,提升了NPML對大角度入射波的吸收能力與NPML的穩(wěn)定性。

CFS-PML吸收邊界條件在均勻各向同性介質(zhì)下的三維二階聲波方程應(yīng)用中,SPML吸收邊界條件易于實現(xiàn),但是波場分裂引入了9個波場變量[47](含3個中間變量),常規(guī)CPML吸收邊界條件通過引入6個輔助變量[40]大大減少了內(nèi)存消耗。趙茂強[48]通過忽略衰減因子空變特性消去了SPML的中間變量,將新增波場數(shù)量降為6個,但并未分析衰減因子忽略帶來的影響。楊凌云等[30]進一步將趙茂強的處理思路[48]引入CPML研究中,提出了NCPML吸收邊界條件,進一步減少了輔助變量的數(shù)量,同時文中分析指出犧牲衰減因子的空變特性后吸收效果完全可以接受。由于楊凌云等[30]提出的是二維聲波方程的NCPML吸收邊界條件,本文將二階系統(tǒng)下的NCPML吸收邊界推廣至三維并成功應(yīng)用于聲波逆時偏移。首先回顧了二階系統(tǒng)下的SPML吸收條件在聲波方程中的應(yīng)用,然后推導(dǎo)出含NCPML吸收邊界條件的三維聲波方程,最后進行三維均勻模型正演和SEG/EAGE推覆體模型和實際資料逆時偏移測試,通過對比NCPML吸收邊界條件與常規(guī)二階SPML吸收邊界條件在數(shù)值模擬中的計算效率與內(nèi)存占用情況,驗證了NCPML在計算效率和穩(wěn)定性上相比常規(guī)SPML所具有的優(yōu)勢。

1 方法原理

1.1 二階系統(tǒng)下的SPML吸收邊界條件

時間域二階常密度聲波方程表示為:

(1)

式中:x,y,z為笛卡爾坐標系中三個不同正交方向;P為標量波場;vP為速度參數(shù);f為震源項。

忽略震源項,將方程(1)變換到頻率域:

(2)

(3)

將方程(3)代入頻率域聲波方程(2)得到尺度變換后的聲波方程:

(4)

式中:sx,sy,sz是x,y,z3個方向的衰減因子。

通過傅里葉反變換得到的時間域方程右端會產(chǎn)生卷積項,以x方向為例:

(5)

式中:F-1表示傅里葉反變換算子;*表示卷積符號。

傳統(tǒng)SPML邊界條件是將波場分裂為笛卡爾坐標系下的3個正交波場,其目的是避免對方程(4)進行傅里葉反變換時出現(xiàn)卷積項。忽略衰減因子1/sn空變特性(即令方程(5)右端第二項為0),應(yīng)用于方程(4)和方程(5),經(jīng)傅里葉反變換后得到常規(guī)SPML吸收邊界條件:

(6)

如果不忽略衰減因子,將會出現(xiàn)波場對時間的三階導(dǎo)數(shù)項,并因此引入3個中間變量。本文中SPML吸收邊界條件采用趙茂強[48]的研究成果,不討論含中間變量的情況,并參考了劉友山等[49]對彈性波PML吸收邊界條件應(yīng)用的研究。

1.2 二階系統(tǒng)下的NCPML吸收邊界條件

雖然SPML方法算法簡單,計算效率高,但是該方法對于計算機內(nèi)存要求較高,分裂的波場提高了正演的內(nèi)存使用量。KOMATITSCH等[29]提出基于一階系統(tǒng)方程的CFS-PML(簡寫為:CPML)吸收邊界條件,通過直接計算卷積來節(jié)省內(nèi)存。前文提到SPML提出的初衷即是避免計算卷積,顯而易見,SPML和CPML本質(zhì)上是在計算復(fù)雜度和計算內(nèi)存取舍方面處于對立位置。PASALIC等[40]將CPML推廣到二階系統(tǒng),其中新增6個輔助變量。楊凌云等[30]在二階系統(tǒng)下忽略了衰減因子空變特性推導(dǎo)出NCPML吸收邊界條件,輔助變量的數(shù)量進一步減少。以上研究均在二維條件下展開,本文將NCPML推廣至三維。

楊凌云等[30]使用遞歸卷積方法有效地解決了方程(5)中的卷積如下:

(7)

(8)

式中:k表示當前時刻;其它衰減參數(shù)為:

(9)

將方程(7)至方程(9)代入方程(1)可以得到三維聲波方程NCPML吸收邊界條件:

(10)

(10)式中輔助變量僅為3個,進一步降低了二階聲波波動方程數(shù)值模擬中內(nèi)存的使用。表1展示了SPML和NCPML吸收邊界條件下一次波場延拓所需要的的內(nèi)存。

1.3 逆時偏移

在各向同性介質(zhì)中,逆時偏移的理論基礎(chǔ)是Claerbout提出的時間一致性原理[50],即反射面存在于地層內(nèi)下行波波至?xí)r間與某一上行波波至?xí)r間相一致之處。成像條件可以表示為震源波場(下行波場)和檢波點波場(上行波場)之間的零延遲互相關(guān)。

(11)

式中:R,U和D分別表示成像剖面、檢波點波場和震源波場;smax和tmax分別表示最大炮數(shù)和最大記錄時間;ti和sj分別表示初至?xí)r間和震源。

表1 內(nèi)存分配情況統(tǒng)計

逆時偏移包含3個步驟[51]:①震源波場正向延拓;②檢波點波場逆向延拓;③應(yīng)用成像條件?；ハ嚓P(guān)成像的量綱為波場振幅的平方,并且該方法對所有波進行成像,偏移剖面低頻噪聲較嚴重,KAELIN等[52]提出了震源歸一化互相關(guān)成像條件,利用震源波場的照明度歸一化互相關(guān)成像剖面,不僅可以消除部分低頻噪聲,而且可以消除震源強度對剖面的影響,獲得的偏移剖面數(shù)值接近于真實的反射系數(shù)。震源歸一化互相關(guān)成像條件表示為:

(12)

三維逆時偏移涉及大量的波場延拓,其對計算機的計算效率、硬件性能有很高的要求,而單一的PC機不能滿足三維RTM計算需求,需要借助集群并行計算和圖形處理器(graphic processing unit,GPU)來實現(xiàn)三維RTM。MPI是最常用的編程工具,為多節(jié)點并行計算提供了一組良好定義的庫函數(shù),CPU集群解決了粗粒度并行問題,可以用于不同炮的正、反傳運算。GPU則解決了細粒度并行問題,波場進行有限差分運算中,空間域中每個網(wǎng)格點都是獨立的,所有網(wǎng)格點可以細粒度并行計算。本文中逆時偏移利用MPI控制的CPU和GPU協(xié)同并行策略實現(xiàn),如圖1所示。第一級為MPI控制命令,控制第二級多個CPU節(jié)點并行運算,CPU節(jié)點控制第三級多個GPU運算。

圖1 MPI控制下CPU/GPU協(xié)同并行策略

綜上所述,本文采用的三維逆時偏移CPU/GPU集群實現(xiàn)方案如圖2所示。NCPML邊界條件應(yīng)用于波場的正向和反向延拓,由于二者延拓方向相反,互相關(guān)成像條件在同一時刻互相關(guān),因此本文采用FENG等[53]提出的逐階遞增模擬精度的源波場重建方法。其中,三維RTM框架最外層為炮循環(huán),采用MPI控制CPU集群進行多炮偏移并行運算,每一次循環(huán)所計算的炮數(shù)為GPU卡數(shù)總和。假設(shè)每個節(jié)點包含b塊GPU卡,共有a臺CPU機器,則每次循環(huán)共算出a×b個炮集數(shù)據(jù)的成像結(jié)果。因為GPU顯存較小,因此每一塊GPU卡只負責(zé)某一炮的波場正反延拓中的差分運算。如果存在GPU型號不同的情況,也需要保證GPU最小顯存能容納一次波場重構(gòu)和反向延拓所需要的的數(shù)組內(nèi)存之和。所有炮集數(shù)據(jù)偏移過程中的I/O均在一塊共享硬盤上完成,并在主節(jié)點進行最后多炮成像結(jié)果的疊加。

圖2 三維逆時偏移CPU/GPU機群實現(xiàn)方案

2 數(shù)值測試

2.1 3D均勻介質(zhì)模型

為了測試NCPML完全匹配層的吸收衰減效果,建立了一個均勻介質(zhì)模型。模型大小為4000m×4000m×4000m,地震波速度為3000m/s?？臻g步長為10m,震源采用主頻為15Hz的雷克子波,震源位于坐標(2000m,2000m,2000m)處,時間采樣間隔為0.8ms;PML吸收層數(shù)為60層;采用時間二階差分精度和空間十階差分精度,在擴展坐標變量中,方向的參數(shù)如下[54]:

(13)

其中,

(14)

式中:x0和d分別表示NCPML邊界的起始位置坐標和厚度;R0是理論上的反射系數(shù)[55],這里取1×10-6;n1和n2是控制NCPML層衰減變化的指數(shù)因子,通常2≤n1≤5且0≤n2≤1,n3是NCPML層的頻移變換尺度的指數(shù)因子,通常取n1=2,n2=0,n3=1。在y和z方向上參數(shù)選擇與x方向類似。

圖3a至圖3c分別表示未加吸收邊界、常規(guī)二階SPML吸收邊界和二階NCPML吸收邊界進行波場模擬計算得到的1100ms時刻的波場快照,從圖3a可以看出,不加吸收邊界時,外行波在人工邊界處完全反射,而加載SPML和NCPML吸收邊界后的波場(如圖3b和3c所示)邊界處外行波明顯被衰減和吸收。由此可見,常規(guī)SPML吸收邊界與NCPML吸收邊界對外行波都具有良好的吸收效果,在邊界處不會產(chǎn)生明顯的虛假反射。

圖3 1100ms時刻三維波場快照a 未加吸收邊界; b 常規(guī)SPML吸收邊界; c 二階NCPML吸收邊界

圖4a顯示了分別利用常規(guī)二階SPML吸收邊界和本文NCPML吸收邊界算法數(shù)值模擬檢波點(1000m,1000m,0)處地震記錄的數(shù)值解及其與理論值的對比,該理論值是通過大幅度擴大計算范圍以消除邊界影響得到數(shù)值解近似而來的,并不是通過給定一個具有一定頻帶范圍的子波,然后在三維均勻空間中求取的理論解。兩種吸收邊界模擬的結(jié)果與理論值很接近,誤差主要出現(xiàn)在980～1200ms的時間段,我們截取該時間段地震記錄單獨顯示(圖4b),NCPML和SPML吸收邊界模擬結(jié)果均存在誤差,二者波形振幅在理論值附近波動,而SPML的誤差比NCPML的誤差更大,為更加直觀地定量分析二者的誤差,將常規(guī)二階SPML和二階NCPML吸收邊界模擬的地震記錄與理論值求差,結(jié)果如圖4c所示。由圖4c可見,SPML吸收邊界模擬記錄的誤差最大值達到5.2×10-3,而NCPML吸收邊界模擬記錄的誤差最大值僅為2.7×10-3,大約是常規(guī)二階SPML的1/2。因此,NCPML吸收邊界對于邊界處的吸收效果要優(yōu)于常規(guī)二階SPML。表2是兩種不同邊界得到的地震記錄的累計絕對誤差和累計誤差的統(tǒng)計表,即對誤差值取平方求和,從表中可以看到,本文提出的NCPML吸收邊界比傳統(tǒng)的SPML吸收邊界的計算精度更高。

圖5a和圖5b顯示了理論值以及利用兩種吸收邊界條件模擬(對應(yīng)于圖3中紅色豎線)得到的1100ms和1200ms時刻的波場值。如圖所示,邊界附近仍然存在一定程度的反射波,但NCPML吸收邊界下的反射波形比SPML吸收邊界下的反射波形在振幅上與理論值更接近。因此本文推導(dǎo)的NCPML吸收邊界的計算精度比常規(guī)SPML更高。

圖4 SPML和NCPML吸收邊界條件在接收點(2000m,2000m,0)位置處的地震記錄(a)、地震記錄局部顯示(b)以及SPML和NCPML吸收邊界條件的地震記錄與理論值之差(c)

表2 誤差統(tǒng)計表

2.2 三維SEG/EAGE推覆體模型

利用三維SEG/EAGE推覆體模型進行逆時偏移測試,速度模型如圖6所示,圖6a為真實速度模型,用以合成地震記錄,圖6b為偏移速度模型。模型大小為16000m×16000m×4800m,x,y和z方向網(wǎng)格間距均為40m。模型測試整體共設(shè)置676炮(26×26),炮點均勻分布于地表,起始炮點坐標(0,0),相鄰兩炮640m間隔向地表兩個正交方向延拓。檢波點布置于地表,間距為40m,即道距和線距均為40m。觀測系統(tǒng)最大偏移距為±4000m,當炮點位置靠近模型邊緣時,例如炮點坐標為(0,0),該炮在偏移過程中僅100條接收線,每條線有100道,而炮點位于模型中央時則接收線數(shù)達到最大的200條,同時每條接收線包含了200道。模型數(shù)據(jù)震源子波為主頻8Hz的雷克子波。地震記錄的總時長為1.6s,采樣點數(shù)為1600,時間采樣間隔1ms。

圖5 SPML和NCPML邊界條件下1100ms(a)和1200ms(b)時刻波場快照抽道對比

圖6 三維SEG/EAGE推覆體模型a 真實速度模型; b 偏移速度模型

三維推覆體模型逆時偏移使用CPU集群共18個節(jié)點,其中共66塊GPU顯卡,GPU型號是Tesla K10C,GPU并行環(huán)境為CUDA,MPICH2版本。表3是三維推覆體模型進行676炮的逆時偏移內(nèi)存和計算耗時結(jié)果,包括所有的計算與I/O讀寫時間,內(nèi)存對比是指RTM過程中單個CPU節(jié)點和單個GPU卡的對比?？梢钥闯鍪褂肗CPML吸收邊界進行三維逆時偏移時在內(nèi)存和計算效率上均得到很大優(yōu)化。

圖7a和圖7b分別為SEG/EAGE推覆體模型應(yīng)用SPML和NCPML邊界條件的逆時偏移三維剖面。三維推覆體模型的目標層位于中深層的河道處,從水平切片中可以看出兩種邊界條件下的三維逆時偏移方法均能夠?qū)拥捞幍臉?gòu)造細節(jié)準確歸位,具有較高的精度。從速度縱向剖面與對應(yīng)的應(yīng)用SPML和NCPML邊界條件的逆時偏移縱向切片對比,可以看出中深層成像能量均衡,淺層由于震源因素影響成像結(jié)果,有稍許噪聲,中深層成像精度較高。

表3 三維推覆體模型逆時偏移所需內(nèi)存和耗時

綜上所述,在本節(jié)的三維SEG/EAGE推覆體模型逆時偏移模型測試中,應(yīng)用NCPML邊界條件與常規(guī)SPML邊界條件計算得到的偏移結(jié)果精度幾乎一致。但是NCPML在內(nèi)存占用上比SPML邊界條件減少了約35%,計算效率也有小幅提高。需指出該數(shù)據(jù)僅在本文所用測試模型中有效,對于大型三維逆時偏移,隨著正演模型增大,NCPML的優(yōu)勢更加明顯。

圖7 不同吸收邊界條件下三維推覆體模型逆時偏移剖面a SPML; b NCPML

2.3 實際應(yīng)用

選取中國東部某三維工區(qū)的疊前地震資料進行試算。數(shù)據(jù)覆蓋區(qū)域約120km2,滿覆蓋區(qū)域約為40km2。該區(qū)地表起伏平緩,地震資料信噪比較高。所選區(qū)塊共1944炮,最大偏移距5464m,面元網(wǎng)格25m×25m。炮集采樣間隔2ms,通過對炮集濾波處理,地震數(shù)據(jù)主頻集中在20Hz。采用偏移速度分析所得的速度模型如圖8所示,最大深度為4000m,網(wǎng)格間距為10m×10 m×10m。測試中所用的偏移距為0～2000m,即單一炮的網(wǎng)格大小為400m×400m×400m。時間采樣間隔0.8ms,針對炮集網(wǎng)格稀疏與偏移中模擬網(wǎng)格規(guī)則化的不符,我們采用多項式插值法進行觀測系統(tǒng)自適應(yīng)。偏移中使用了所有炮以增加覆蓋次數(shù),壓制噪聲。

數(shù)據(jù)預(yù)處理階段,由于三維逆時偏移只用到資料中反射波的信息,因此在實施之前去除資料中的面波和異常噪聲,在盡量保留有效信號的前提下壓制噪聲,提供高保真的道集數(shù)據(jù)。同時,我們對所選工區(qū)實際資料做能量補償以消除目的層縱、橫向能量差異,有利于目的層弱信號的成像。

圖8三維初始速度模型

三維疊前資料逆時偏移利用GPU集群8個節(jié)點,所用GPU型號是M2090,共12個GPU,GPU并行環(huán)境為CUDA,節(jié)點間的并行環(huán)境為MPICH2。表4為利用上述運行環(huán)境進行三維逆時偏移所需的內(nèi)存消耗及時間消耗,這里包括所有的I/O讀寫時間。由于計算設(shè)備數(shù)量和運算性能的限制,三維逆時偏移測試的運算效率較慢,因此只進行了NCPML邊界條件的應(yīng)用。

表4 三維實際資料逆時偏移所需內(nèi)存與耗時

圖9是應(yīng)用NCPML邊界條件的三維逆時偏移測試結(jié)果,分別為抽取In-line測線和Cross-line測線的逆時偏移結(jié)果,可以看出逆時偏移結(jié)果深淺層能量均衡,成像精度高,這是使用歸一化成像條件的原因,使得中、深層能量得到照明補償。邊界處覆蓋次數(shù)較低,因而成像結(jié)果中噪聲較多。淺層由于震源等因素的影響,有一定的噪聲。綜上所述,NCPML在實際資料應(yīng)用中計算效率和內(nèi)存占用都比較令人滿意。

圖9 不同測線的三維逆時偏移效果a In-Line測線;b Cross-Line測線

3 結(jié)論

本文將二階系統(tǒng)下的NCPML吸收邊界條件推廣至三維并成功應(yīng)用于聲波逆時偏移。均勻模型正演以及三維SEG/EAGE推覆體模型逆時偏移測試結(jié)果表明:NCPML應(yīng)用效果優(yōu)于常規(guī)SPML,同時計算效率和內(nèi)存使用方面NCPML比SPML更具優(yōu)勢。實際資料的逆時偏移測試證明了NCPML吸收邊界條件在實際應(yīng)用中的穩(wěn)定性。此外,針對二階聲波方程提出的NCPML可以推廣到更復(fù)雜的介質(zhì)波動方程模擬,例如聲學(xué)各向異性建模和彈性介質(zhì)。