柯 璇,石 穎,2,宋利偉,李松齡

(1.東北石油大學(xué),黑龍江大慶163318;2.黑龍江省普通高?？萍紕?chuàng)新團隊“斷層變形、封閉性及與流體運移”,黑龍江大慶163318)

基于褶積完全匹配吸收邊界的聲波方程數(shù)值模擬

柯 璇1,石 穎1,2,宋利偉1,李松齡1

(1.東北石油大學(xué),黑龍江大慶163318;2.黑龍江省普通高?？萍紕?chuàng)新團隊“斷層變形、封閉性及與流體運移”,黑龍江大慶163318)

完全匹配層(PML)吸收邊界條件已廣泛應(yīng)用于地震波場模擬中,但PML吸收邊界條件存在一些問題,如對低頻和地震波為掠射波時會產(chǎn)生虛假反射等。褶積完全匹配層(CPML)吸收邊界條件能夠有效地解決PML吸收邊界條件存在的問題。推導(dǎo)了帶有CPML條件的一階速度-應(yīng)力波動方程,并且在方程中引入記憶變量代替復(fù)雜褶積項的運算,將CPML與交錯網(wǎng)格有限差分算法相結(jié)合,發(fā)揮了該邊界條件節(jié)省計算存儲空間、易于編程實現(xiàn)等優(yōu)點。數(shù)值模擬結(jié)果表明,在不增加計算量的情況下,CPML吸收邊界條件有效地提高了波場模擬精度。

褶積完全匹配層;交錯網(wǎng)格;數(shù)值模擬;聲波方程;邊界條件

隨著計算機技術(shù)的飛速發(fā)展,越來越多的數(shù)值方法被應(yīng)用于基于波動方程的地震波場數(shù)值模擬中,主要包括有限差分法[1]、有限元法[2]、偽譜法[3]等,這些方法的應(yīng)用提升了逆時偏移成像[4-7]、全波形反演[8-9]等技術(shù)解決實際問題的能力。高精度地震波場的模擬成為這些新技術(shù)的基礎(chǔ),也是眾多地球物理學(xué)家研究的重點[10-11]。

實際情況下的地震波場是在半無限區(qū)域傳播的,而數(shù)值模擬只能對有限區(qū)域模擬計算,這樣在邊界內(nèi)側(cè)和外側(cè)形成了具有不同波阻抗的兩種介質(zhì),計算邊界相當于一個反射面,當波場能量傳播到計算邊界時會給計算區(qū)域帶來虛假反射,從而導(dǎo)致錯誤的結(jié)果,因此需要人為設(shè)置邊界介質(zhì)將能量吸收,避免能量反射。在過去幾十年,發(fā)展了海綿吸收[12]、最優(yōu)化條件[13]、特征值分解法[14]、球面精確吸收條件[15]等吸收邊界技術(shù),但是,這些技術(shù)在應(yīng)用中受到一定的限制,如在地震波為掠射波情況下會產(chǎn)生虛假反射、所有角度入射會產(chǎn)生低頻反射,而且實施過程中存在大量數(shù)值計算。20世紀90年代,BERENGER[16]為了對電磁波數(shù)據(jù)精確建模提出了一種新的邊界條件,稱之為完全匹配層(perfectly matched layer,PML)吸收邊界條件,其基本思想是:在吸收邊界區(qū)域匹配與計算區(qū)域相同波阻抗,使入射波無反射的進入吸收層,吸收層中是衰減介質(zhì),使得入射波迅速衰減直至消除。理論上PML邊界能夠吸收任何入射角度和頻率的入射波,實踐也證明PML吸收邊界條件十分有效并被成功應(yīng)用于聲波[17-20]和彈性波[21]的研究。

如果在波場延拓的過程中采用基于波動方程時域有限差分進行求解,那么其中的離散化處理會降低PML吸收邊界的效果,尤其在掠入射的情況下,虛假反射更加明顯[22],使得PML在對薄層介質(zhì)、震源位置靠近網(wǎng)格邊緣或者檢波器處于遠偏移距的模型的數(shù)值模擬效果并不理想。為了解決上述問題,KUZUOGLU等[2]提出了一種復(fù)頻移方法,這種方法被RODEN等[23]稱為褶積完全匹配層(convolutional perfectly matched layer,CPML)吸收邊界。

在前人的研究基礎(chǔ)上,本文對CPML吸收邊界下的一階速度-應(yīng)力方程進行了詳盡推導(dǎo),并將其與交錯網(wǎng)格有限差分相結(jié)合,進行波場數(shù)值模擬。對比分析波場模擬結(jié)果可知,CPML具有更好的吸收效果并且能夠解決以往吸收邊界在掠入射情況下吸收效果不理想的問題。最后用3個模型對本文方法進行了測試。

1 一階速度-應(yīng)力聲波方程CPML的推導(dǎo)

當?shù)卣鸩ㄒ月尤肷溥M入到經(jīng)典PML吸收層時,并沒有傳播至吸收層較深部位置,而是沿著平行于吸收層傳播,CPML的思想是在經(jīng)典PML吸收邊界的基礎(chǔ)上,引入更加廣譜的濾波器將這部分能量衰減并且阻止其離開吸收層。由于交錯網(wǎng)格在不增加計算量情況下能夠壓制數(shù)值頻散,為了便于將交錯網(wǎng)格與CPML吸收邊界相結(jié)合,本文利用一階速度應(yīng)力聲波方程進行推導(dǎo)。在直角坐標系中,二維一階速度-應(yīng)力方程的時域形式[11]如下:

式中:v和w分別為質(zhì)點在x和z方向上介質(zhì)震動速度;p代表聲壓;t為時間;ρ為介質(zhì)密度;k為介質(zhì)的體模量。

為了在波動方程中引入CPML吸收邊界條件,將聲波方程進行傅里葉變換得:

(7)

式中:dx(s)是在x方向的衰減系數(shù);衰減程度依賴于傳播距離;κx主要影響以掠射形式入射到衰減層地震波的衰減;αx主要影響入射到衰減層的低頻成分。如果κx=1,αx=0,CPML邊界條件就退化為經(jīng)典PML邊界條件。

(8)

式中:

(9)

(10)

公式(9)有如下的形式:

(11)

因此,

(12)

式中:δ是狄拉克函數(shù);H是階躍函數(shù)。令:

(13)

在時間域中公式(8)變?yōu)?

(14)

(15)

由于對時間的積分采用了交錯網(wǎng)格,所以對空間的導(dǎo)數(shù)可以定義在mΔt和(m+1)Δt中間的時刻上[24],因此(15)式變?yōu)?

(16)

式中:

對(16)式進一步推導(dǎo):

(17)

本文中αx,dx,κx的形式分別為:

2 基于CPML的交錯網(wǎng)格差分格式

求解帶有CPML吸收邊界的波動方程,需要將其離散化,交錯網(wǎng)格較常規(guī)網(wǎng)格更能夠有效壓制頻散[25],在不增加計算量的情況下可以提高數(shù)值模擬精度,因此本文采用交錯網(wǎng)格有限差分的方法對波場進行正演計算。波場及相關(guān)參數(shù)分布在相差半網(wǎng)格點的主體網(wǎng)格和交錯網(wǎng)格點上,以上定義能夠滿足遞推要求。網(wǎng)格布局如圖1所示,網(wǎng)格中介質(zhì)波場及參數(shù)見表1。

圖1 網(wǎng)格布局

網(wǎng)格波場及參數(shù)1p,k2v,ρ3w,ρ

方程(18)至方程(20)的2N階交錯網(wǎng)格差分格式為:

式中:an為差分系數(shù)。

3 數(shù)據(jù)測試

3.1 測試一

為了驗證CPML邊界條件的有效性,進行了各向同性介質(zhì)二維聲波波場數(shù)值模擬。數(shù)值模型計算采用時間二階、空間十階交錯網(wǎng)格有限差分法,模型的計算區(qū)域為水平與垂直方向均為200網(wǎng)格點,水平和垂直方向的空間網(wǎng)格大小均為10m,介質(zhì)速度為2500m/s,密度1000kg/m3,震源位于(100,100)網(wǎng)格點處,采用主頻為20Hz的雷克子波,時間采樣間隔為1ms,利用公式(21)至公式(23)對衰減系數(shù)進行設(shè)置,κx取1,正演波場快照如圖2所示。圖2a表示計算區(qū)域外有20層CPML吸收層的波場快照,圖上基本上觀察不到來自邊界的反射,計算區(qū)域的波場更加純凈,圖2b表示計算區(qū)域外有20層經(jīng)典PML吸收層的波場快照,雖然經(jīng)典PML層吸收了大部分的能量,但是圖上還是能看到來自邊界較弱能量的反射,這些反射能量對計算區(qū)域波場形成了噪聲干擾,不利于波場的精確模擬,圖2c為計算區(qū)域外有50層經(jīng)典PML吸收層的波場快照,可看出,當?shù)卣鸩▊鞑サ竭吔鐣r,沒有出現(xiàn)明顯的能量反射現(xiàn)象。對比分析可知,在相同衰減層數(shù)的情況下,CPML吸收邊界的吸收效果明顯優(yōu)于經(jīng)典PML吸收邊界;當增加吸收層數(shù)時,經(jīng)典PML吸收邊界也能夠達到理想效果。因此,為實現(xiàn)相同的吸收效果,CPML吸收邊界條件所需吸收層數(shù)更少,從而降低了計算所需的內(nèi)存空間。

圖2 不同邊界條件及不同吸收層吸收效果對比a CPML 20層; b 經(jīng)典PML 20層; c 經(jīng)典PML 50層

3.2 測試二

逆時偏移成像過程中,吸收邊界條件的選取直接關(guān)系到波場模擬的精度,最終影響到成像剖面的刻畫。為了說明在掠入射情況下,CPML吸收邊界較經(jīng)典PML具有更好的吸收效果,特設(shè)計均勻介質(zhì)模型如圖3所示。計算區(qū)域為水平方向440網(wǎng)格點,垂直方向240網(wǎng)格點,震源s位于(220,25)網(wǎng)格點處,在(60,30)網(wǎng)格點布置檢波器1,在(220,25)網(wǎng)格點布置檢波器2,在計算區(qū)域外均采用30層的吸收層,其它參數(shù)同測試一,正演模擬波場快照如圖4所示。

圖3 模型示意

從圖4a可以看到,地震波與吸收層入射角較小時,經(jīng)典PML吸收邊界對入射地震波的吸收效果較好,但是對掠入射地震波的吸收能力尚有不足,有來自邊界的反射現(xiàn)象,這對有效信號形成噪聲干擾,不利于波場的精確模擬;而圖4b采用CPML吸收邊界,不但具有更高的吸收能力,更重要的是解決了地震波與吸收層形成掠入射情況下的吸收問題,計算區(qū)域波場較為純凈。這是由于經(jīng)典PML吸收邊界條件的衰減系數(shù)依賴于入射波的方向,當入射地震波與吸收層入射角較小時吸收系數(shù)有較大值,當入射角較大直至發(fā)生掠入射的情況,吸收系數(shù)有極小值,因此能看到圖4a箭頭位置的邊界反射現(xiàn)象,而在CPML吸收邊界條件中引入的是復(fù)坐標變換,一方面可以將入射波校正到垂直入射的情況;另一方面可以削弱低頻反射現(xiàn)象,從而提高了正演模擬精度。

圖4 不同邊界條件吸收效果對比a PML; b CPML

為了更加清晰地觀察吸收邊界條件對波場正演模擬精度的影響,采集圖3所示模型中檢波器1和檢波器2的單道信號,并與理論解進行比較。本文所述的理論解為在波場正演模擬中,對模擬區(qū)域進行擴展,從而保證地震波場不含任何計算邊界影響時所得到的用于參考的信號記錄。圖5為分別采用CPML和經(jīng)典PML吸收邊界條件采集的檢波器1單道信號,由于檢波器1距激發(fā)點較遠,且位置接近上邊界,而地震波入射到上邊界屬于掠入射。圖中紅色曲線為理論解,黑色虛線為檢波器記錄的單道信號,由圖5a可見,采用CPML吸收邊界條件的單道記錄與理論單道記錄吻合較好,雖然屬于掠入射情況,但是CPML邊界能對掠射波進行很好的吸收;由圖5b可見,經(jīng)典PML吸收邊界條件對掠射波吸收的效果并不理想,正演波場出現(xiàn)誤差。圖6為分別采用CMPL和經(jīng)典PML吸收邊界采集的檢波器2單道記錄,兩種單道記錄與理論記錄相吻合,說明在地震波與吸收層成小角度入射情況下,兩種邊界條件吸收效果均較理想。

圖5 不同吸收邊界條件下采集的檢波器1單道記錄a CPML; b PML

圖6 不同吸收邊界條件下采集的檢波器2單道記錄a CPML; b PML

3.3 測試三

為考察不同吸收邊界對地震記錄的影響,設(shè)計水平層狀介質(zhì)模型,第1層厚度620m,速度為2300m/s;第2層厚度460m,速度為2600m/s;第3層厚度420m,速度為2800m/s。對該3層水平層狀介質(zhì)進行正演模擬,計算區(qū)域為水平方向300網(wǎng)格點,垂直方向150網(wǎng)格點,震源位于(150,2)網(wǎng)格點,地震記錄時長1.3s,兩種吸收邊界的層數(shù)均為20層,其它參數(shù)均相同,地震記錄如圖7所示。從圖7a箭頭處可以明顯觀察到來自邊界的反射,而圖7b中只有經(jīng)過2個界面反射的雙曲同相軸。對比可知,利用CPML吸收邊界條件可以得到更加純凈的地震記錄。

圖7 不同吸收邊界條件下水平層狀介質(zhì)地震記錄對比a PML; b CPML

4 結(jié)論與討論

本文詳細推導(dǎo)了CPML吸收邊界條件下的一階速度-應(yīng)力聲波方程,并將CPML技術(shù)與交錯網(wǎng)格有限差分相結(jié)合,對均勻介質(zhì)進行了數(shù)值模擬。對比經(jīng)典PML和CPML吸收邊界條件吸收效果可知,在相同吸收層數(shù)的情況下,CPML吸收邊界條件吸收效果明顯優(yōu)于經(jīng)典PML吸收邊界條件,并且CPML吸收邊界能夠解決低頻反射和掠入射的問題,從而提高地震波場正演模擬精度。若將CPML吸收邊界應(yīng)用于逆時偏移成像等技術(shù)中,可以提升波場傳播精度,而且較少的吸收層就能夠達到理想效果,節(jié)約了計算存儲。理論上,CPML吸收邊界也可以應(yīng)用到基于qP波動方程的各向異性介質(zhì)正演模擬中,相對于PML吸收邊界,它可以對qP波和qS波在遠偏移距的情況下形成更好的吸收效果。

[1] ALTERMAN Z.Propagation of elastic waves in layered media by finite difference methods[J].Bulletin of the Seismological Society of America,1968,58(1):367-398

[2] KUZUOGLU M,Mittra R.Frequency dependence of the constitutive parameters of causal perfectly matched anisotropic absorbers[J].IEEE Microwave and Guided Wave Letters,1996,12(6):447-449

[3] CARCIONE J M.The wave equation in generalized coordinates[J].Geophysics,1994,59(12):1911-1919

[4] 劉守偉,王華忠,陳生昌,等.VSP上下行反射波聯(lián)合成像方法研究[J].地球物理學(xué)報,2012,55(9):3126-3133 LIU S W,WANG H Z,CHEN S C,et al.Joint imaging method of VSP upgoing and downgoing reflection wave[J].Chinese Journal of Geophysics,2012,55(9):3126-3133

[5] 陳生昌,馬在田,陳林.三維VSP數(shù)據(jù)的波動方程偏移成像[J].天然氣工業(yè),2008,28(3):51-53 CHEN S C,MA Z T,CHEN L.Wave equation migration for 3-D VSP data[J].Natural Gas Industry,2008,28(3):51-53

[6] 王維紅,柯璇,郭雪豹,等.基于Qt的Walkaway VSP逆時偏移軟件開發(fā)研究[J].石油物探,2015,54(4):452-458 WANG W H,KE X,GUO X B,et al.Qt based software development of reverse time migration for walkaway VSP seismic data[J].Geophysical Prospecting for Petroleum,2015,54(4):452-458

[7] 郭書娟,馬方正,段心標,等.最小二乘逆時偏移成像方法的實現(xiàn)與應(yīng)用研究[J].石油物探,2015,54(3):301-308 GUO S J,MA F Z,DUAN X B,et al.Research of least-squares reverse-time migration imaging method and its application [J].Geophysical Prospecting for Petroleum,2015,54(3):301-308

[8] FICHTNER A,KENNETT B,IGEL H,et al.Full seismic waveform tomography for upper-mantle structure in the Australasian region using adjoint methods[J].Geophysical Journal International,2009,179(3):1703-1725

[9] 董良國,遲本鑫,陶紀霞,等.聲波全波形反演目標函數(shù)性態(tài)[J].地球物理學(xué)報,2013,56(10):3445-3460 DONG L G,CHI B X,TAO J X,et al.Objective function behavior in acoustic full-waveform inversion[J].Chinese Journal of Geophysics,2013,56(10):3345-3460

[10] 王維紅,柯璇,裴江云.完全匹配層吸收邊界條件應(yīng)用研究[J].地球物理學(xué)進展,2013,28(5):2508-2514 WANG W H,KE X,PEI J Y.Application investigation of perfectly matched layer absorbing boundary condition[J].Progress in Geophysics,2013,28(5):2508-2514

[11] 王保利,高靜懷,陳文超,等.地震疊前逆時偏移的有效邊界存儲策略[J].地球物理學(xué)報,2012,55(7):2412-2421 WANG B L,GAO J H,CHEN W C,et al.Efficient boundary storage strategies for seismic reverse time migration[J].Chinese Journal of Geophysics,2012,55(7):2412-2421

[12] CERJAN C,KOSLOFF D,KOSLOFF R,et al.A nonreflecting boundary condition for discrete acoustic and elasticwave equation[J].Geophysics,1985,50(4):705-708

[13] PENG C B,TOKSOZ M N.An optimal absorbing boundary condition for elastic wave modeling[J].Geophysics,1995,60(1):296-301

[14] DONG L D,SHE D P,GUAN L P,et al.An eigenvalue decomposition method to construct absorbing boundary conditions for acoustic and elastic wave equations[J].Journal of Geophysics and Engineering,2005,2(3):192-198

[15] GROTE M J.Nonreflecting boundary conditions for elastodynamic scattering[J].Journal of Computational Physics,2000,161(1):331-353

[16] BERENGER J P.A perfectly matched layer for the absorption of electromagnetic waves[J].Journal of Computational Physics,1994,114(2):185-200

[17] LIU Q H,TAO J P.The perfectly matched layer for acoustic waves in absorptive media[J].Journal of the Acoustical Society of America,1997,102(4):2072-2082

[18] 王守東.聲波方程完全匹配層吸收邊界[J].石油地球物理勘探,2003,38(1):31-34 WANG S D.Absorbing boundary condition for acoustic wave equation by perfectly matched layer[J].Oil Geophysical Prospecting,2003,38(1):31-34

[19] 張衡,劉洪,李博,等.VTI介質(zhì)聲波方程非分裂式PML吸收邊界條件研究[J].石油物探,2016,55(6):781-792 ZHANG H,LIU H,LI B,et al.The research on unsplit PML absorbing boundary conditions of acoustic for VTI media[J].Geophysical Prospecting for Petroleum,2016,55(6):781-792

[20] 吳國忱,梁鍇.VTI介質(zhì)準P波頻率空間域組合邊界條件研究[J].石油物探,2005,44(4):301-307 WU G C,LIANG K.Combined boundary conditions ofquasi-P wave within ferquency-space domain in VTI media [J].Geophysical Prospecting for Petroleum,2005,44(4):301-307

[21] CHEW W C,LIU Q H.Perfectly matched layers for elastodynamics:A new absorbing boundary condition[J].Journal of Computational Acoustics,1996,4(4):341-359

[22] WINTON S C,RAPPAPORT C M.Specifying PML conductivities by considering numerical reflection dependencies[J].IEEE Transactions on Antennas and Propagation,2000,48(7):1055-1063

[23] RODEN J A,GEDNEY S D.Convolution PML(CPML):an efficient FDTD implementation of the CFS-PML forarbitrary media[J].IEEE Microwave and Optical Technology Letters,2000,27(5):334-339

[24] KOMATITSCH D,MARTIN R.An unsplit convolutional perfectly matched layer improved at grazing incidence for the seismic wave equation[J].Geophysics,2007,72(5):SM155-SM167

[25] MADARIAGA R.Dynamics of an expanding circular fault[J].Bulletin of the Seismological Society of America,1976,66(3):639-666

(編輯：朱文杰)

Numericalmodelingofacousticwaveequationsbasedonconvolutionalperfectlymatchedlayerabsorbingboundarycondition

KE Xuan1,SHI Ying1,2,SONG Liwei1,LI Songling1

(1.NortheastPetroleumUniversity,Daqing163318,China；2.HeilongjiangProvincecollegescienceandtechnologyinnovationteam“faultsdeformation,sealingabilityandfluidmigration”,Daqing163318,China)

The perfectly matched layer (PML) absorbing boundary condition has been widely undertaken in numerical modeling of seismic wave fields; however,PML reflects a large amount of spurious energy at grazing incidence or a low frequency at all angles of incidence.The convolutional perfectly matched layer (CPML) absorbing boundary condition can effectively solve the above problem.We deduce a first-order velocity-stress wave equation with CPML,substitute the memory variables for complex convolution operation in the equation,and combine the CPML and a staggered-grid finite difference method in the forward modeling.CPML can save the memory of computers and programs easily from the implementation process.The result of the numerical simulation testing indicated that CPML can improve the precision effectively in the wave field simulation without increasing the computational cost.

convolutional perfectly matched layer (CPML),staggered-grid,numerical modeling,acoustic wave equation,boundary condition

P631

：A

1000-1441(2017)05-0637-07DOI:10.3969/j.issn.1000-1441.2017.05.003

柯璇,石穎,宋利偉,等.基于褶積完全匹配吸收邊界的聲波方程數(shù)值模擬[J].石油物探,2017,56(5):643

KE Xuan,SHI Ying,SONG Liwei,et al.Numerical modeling of acoustic wave equations based on convolutional perfectly matched layer absorbing boundary condition

[J].Geophysical Prospecting for Petroleum,2017,56(5):643

2016-11-04;改回日期：2017-03-30。

柯璇(1989—),男,博士在讀,主要從事地震波場正演模擬與偏移成像等方面的研究。

石穎(1976—),女,教授,博士生導(dǎo)師,主要從事地震資料處理方面的研究。

國家自然科學(xué)基金(41474118)、黑龍江省杰出青年科學(xué)基金項目(JC2016006)、東北石油大學(xué)研究生創(chuàng)新科研項目(YJSCX2016-002NEPU)以及大連理工大學(xué)海岸和近海工程國家重點實驗室開放基金項目(LP1509)聯(lián)合資助。

This research is financially supported by the Natural Science Foundation of China (Grant No.41474118),the Outstanding Youth Science Fund of Heilongjiang Province (Grant No.JC2016006),the Northeast Petroleum University Innovation Foundation For Postgraduate (Grant No.YJSCX2016-002NEPU)and the Open Foundation Project of Dalian University of Technology,State Coastal and Offshore Engineering Key Laboratory (Grant No.LP1509).