張 衡,劉 洪,李 博,丁仁偉,李麗青,李福元

(1.國(guó)土資源部海底礦產(chǎn)資源重點(diǎn)實(shí)驗(yàn)室,中國(guó)地質(zhì)調(diào)查局廣州海洋地質(zhì)調(diào)查局,廣東廣州510075;2.中國(guó)科學(xué)院地質(zhì)與地球物理研究所,中國(guó)科學(xué)院油氣資源研究重點(diǎn)實(shí)驗(yàn)室,北京100029;3.中國(guó)石油化工股份有限公司石油物探技術(shù)研究院,江蘇南京211103;4.山東科技大學(xué)地球科學(xué)與工程學(xué)院,山東青島266590)

?

VTI介質(zhì)聲波方程非分裂式PML吸收邊界條件研究

張 衡1,劉 洪2,李 博3,丁仁偉4,李麗青1,李福元1

(1.國(guó)土資源部海底礦產(chǎn)資源重點(diǎn)實(shí)驗(yàn)室,中國(guó)地質(zhì)調(diào)查局廣州海洋地質(zhì)調(diào)查局,廣東廣州510075;2.中國(guó)科學(xué)院地質(zhì)與地球物理研究所,中國(guó)科學(xué)院油氣資源研究重點(diǎn)實(shí)驗(yàn)室,北京100029;3.中國(guó)石油化工股份有限公司石油物探技術(shù)研究院,江蘇南京211103;4.山東科技大學(xué)地球科學(xué)與工程學(xué)院,山東青島266590)

針對(duì)垂直橫向各向同性(VTI)介質(zhì)聲波波動(dòng)方程的特點(diǎn),研究了VTI介質(zhì)非分裂式完全匹配層(Unsplit perfectly matched layer,UPML)吸收邊界條件。首先從VTI介質(zhì)波場(chǎng)傳播穩(wěn)定性和減少偽橫波噪聲的角度對(duì)常見(jiàn)的幾種VTI介質(zhì)聲波波動(dòng)方程進(jìn)行了歸納和對(duì)比,發(fā)現(xiàn)VTI介質(zhì)聲波近似方程更適用于VTI介質(zhì)情形。然后以VTI介質(zhì)聲波近似二階波動(dòng)方程為基礎(chǔ),推導(dǎo)出UPML波動(dòng)方程,并給出求解該波動(dòng)方程的具體數(shù)值實(shí)現(xiàn)方法。數(shù)值模擬結(jié)果表明,UPML吸收邊界條件能達(dá)到很好的人工邊界反射吸收效果,相比于優(yōu)化海綿吸收邊界條件,其人工邊界反射吸收效果明顯提高。

VTI介質(zhì);非分裂式完全匹配層;偽橫波噪聲;聲波近似方程;二階波動(dòng)方程

VTI介質(zhì)正演數(shù)值模擬是逆時(shí)偏移和全波形反演的基礎(chǔ)[1]。VTI介質(zhì)參數(shù)又稱(chēng)Thomsen各向異性參數(shù)(包括縱波各向異性參數(shù)ε和變異系數(shù)δ),由THOMSEN[2]于1986年引入地震偏移領(lǐng)域,因其能夠描述地球介質(zhì)中大多數(shù)各向異性情形而得到廣泛應(yīng)用。THOMSEN[2]在對(duì)全球范圍內(nèi)已知的地球介質(zhì)各向異性數(shù)據(jù)統(tǒng)計(jì)的基礎(chǔ)上,歸納出弱各向異性、強(qiáng)各向異性和極端各向異性3種情形,并根據(jù)統(tǒng)計(jì)結(jié)果得出大多數(shù)各向異性介質(zhì)為弱各向異性介質(zhì)即橫向各向同性介質(zhì)(TI介質(zhì))的結(jié)論。垂直橫向各向同性介質(zhì)(VTI介質(zhì))作為一種最基本的弱各向異性介質(zhì),其研究具有重要的理論和實(shí)際意義。VTI介質(zhì)可以看作是具有傾斜對(duì)稱(chēng)軸的橫向各向同性(TTI)介質(zhì)在對(duì)稱(chēng)軸角度為零時(shí)的一種特殊情形。

1996年,TSVANKIN[3]基于Thomsen弱各向異性理論推導(dǎo)得到P-SV波VTI介質(zhì)精確相速度頻散方程。2006年,ZHOU等[4]從TSVANKIN相速度頻散方程出發(fā),采用橫波速度為零的VTI介質(zhì)聲波近似,推導(dǎo)得到了P-SV波耦合的四階VTI介質(zhì)qP波波動(dòng)方程。因?yàn)樗碾A各向異性波動(dòng)方程難以求解,ZHOU等[4]引入輔助波場(chǎng)項(xiàng)將四階各向異性波動(dòng)方程轉(zhuǎn)換為兩個(gè)二階的各向異性波動(dòng)方程,但是其仍然是一個(gè)P-SV波耦合的波動(dòng)方程,因?yàn)楦飨虍愋圆▓?chǎng)傳播過(guò)程中偽橫波(SV波)仍然存在[5],在角度劇變情況下波場(chǎng)傳播存在嚴(yán)重的不穩(wěn)定性問(wèn)題[6]。為了解決TTI介質(zhì)聲波近似各向異性波場(chǎng)傳播存在的不穩(wěn)定性問(wèn)題,2009年,FLETCHER等[7]引入適當(dāng)?shù)臋M波分量來(lái)減小角度劇變時(shí)產(chǎn)生的不穩(wěn)定性,但是同時(shí)也帶來(lái)了較強(qiáng)的偽SV波假象。DUVENECK等[8]從Hooke定律和運(yùn)動(dòng)學(xué)方程出發(fā)推導(dǎo)得到新的VTI介質(zhì)方程,ZHANG等[9-10]以VTI介質(zhì)彈性波方程為基礎(chǔ),利用自共軛旋轉(zhuǎn)方式推導(dǎo)了自共軛3D VTI波動(dòng)方程,但這些方程仍然存在偽SV波假象且在角度劇變時(shí)不穩(wěn)定。HESTHOLM[11]將二階VTI方程變換成6個(gè)一階耦合偏微分方程的形式,研究了VTI介質(zhì)的高階有限差分模擬。LIU等[12]研究了一種新的基于時(shí)空域頻散關(guān)系優(yōu)化的VTI介質(zhì)數(shù)值模擬算法,相比常規(guī)算法,該算法明顯提高了計(jì)算精度。FOWLER等[13]從特征值矩陣分析的角度對(duì)各種VTI介質(zhì)耦合方程的精度和穩(wěn)定性進(jìn)行了綜合分析。

VTI介質(zhì)正演時(shí),吸收邊界條件至關(guān)重要,因?yàn)槲者吔鐥l件的效果會(huì)極大影響正演模擬效果。各向異性地震波場(chǎng)正演數(shù)值模擬時(shí)計(jì)算區(qū)域有限,必須研究有效的吸收邊界條件對(duì)強(qiáng)人工邊界反射進(jìn)行有效吸收。目前常用的吸收邊界條件包括單程波旁軸近似吸收邊界條件[14]和衰減吸收邊界條件[15-16]兩大類(lèi)。CLAYTON等[14]提出的單程波旁軸近似吸收邊界條件應(yīng)用廣泛,但是其對(duì)大角度入射的波吸收效果較差。CERJAN等[15]提出了添加阻尼層的海綿吸收邊界條件,通過(guò)添加阻尼層達(dá)到衰減效果,但是邊界吸收效果仍較差。BORDING[17]對(duì)CERJAN方法進(jìn)行改進(jìn),提出了優(yōu)化的海綿吸收邊界條件,但是計(jì)算效率相對(duì)CERJAN方法有所降低。完全匹配層(PML)吸收邊界條件最早于1994年由BERENGER[16]在電磁學(xué)領(lǐng)域提出,并于1996年由CHEW等[18]將其應(yīng)用于地震波數(shù)值模擬,取得了很好的應(yīng)用效果。BéCACHE等[19]指出了各向異性介質(zhì)應(yīng)用PML吸收邊界條件時(shí)的不穩(wěn)定性。PML吸收邊界條件是目前應(yīng)用效果最好的吸收邊界條件,但是其公式推導(dǎo)和編程實(shí)現(xiàn)復(fù)雜,需要對(duì)不同的波動(dòng)方程推導(dǎo)不同的PML方程形式。PML吸收邊界條件技術(shù)的研究代表吸收邊界技術(shù)研究的前沿發(fā)展方向。

作為PML吸收邊界條件技術(shù)研究的一個(gè)重要方面,UPML吸收邊界條件已經(jīng)得到了廣泛應(yīng)用。WANG等[20]在彈性波介質(zhì)模擬中提出了一種新的UPML算法。DROSSAERT等[21]提出了一種基于遞歸積分復(fù)頻移UPML算法。QIN等[22]推導(dǎo)得到一種簡(jiǎn)潔有效的非褶積實(shí)現(xiàn)的UPML算法,該算法既不需要對(duì)波場(chǎng)分量進(jìn)行分裂,也不需要做復(fù)雜的褶積運(yùn)算,且具有良好的吸收衰減性能。

本文分析了常見(jiàn)的幾種VTI介質(zhì)波動(dòng)方程,并從波場(chǎng)傳播穩(wěn)定性和減少偽橫波噪聲的角度對(duì)VTI介質(zhì)波動(dòng)方程進(jìn)行了對(duì)比,得出了VTI介質(zhì)聲波近似方程更適用于VTI介質(zhì)情形的結(jié)論。然后從VTI介質(zhì)聲波近似二階波動(dòng)方程出發(fā)推導(dǎo)得到2D VTI介質(zhì)時(shí)間域UPML波動(dòng)方程,并給出了2D VTI介質(zhì)時(shí)間域UPML波動(dòng)方程高階規(guī)則網(wǎng)格有限差分算法求解的具體實(shí)現(xiàn)方法。最后利用簡(jiǎn)單均勻2D VTI介質(zhì)模型和復(fù)雜BP2007 2D VTI介質(zhì)海洋標(biāo)準(zhǔn)模型對(duì)UPML波動(dòng)方程進(jìn)行了測(cè)試,結(jié)果表明本文研究的VTI介質(zhì)時(shí)間域UPML吸收邊界條件能很好地吸收邊界反射,能夠?qū)崿F(xiàn)VTI介質(zhì)高精度數(shù)值模擬。

1 方法原理

1.1 VTI介質(zhì)耦合qP波波動(dòng)方程

P-SV波VTI介質(zhì)相速度頻散方程[3]為:

(1)

式中:vP(θ)表示精確相速度;θ表示相速度傳播角度;vSZ表示橫波速度;vPZ表示縱波速度;ε和δ表示Thomsen各向異性參數(shù)[1]。

直接引入VTI介質(zhì)聲波近似(VTI介質(zhì)聲波近似情況下直接令橫波速度vSZ=0)[23],推導(dǎo)得到的2D VTI介質(zhì)波動(dòng)方程為[4]:

(2)

本文將采用直接令橫波為零的VTI介質(zhì)聲波近似思想推導(dǎo)得到的VTI介質(zhì)方程稱(chēng)為基于VTI介質(zhì)聲波近似的VTI介質(zhì)聲波近似方程。但是這種方程在實(shí)際計(jì)算中往往會(huì)產(chǎn)生數(shù)值不穩(wěn)定問(wèn)題,尤其是在角度劇變的情況下[6]。FLETCHER等[7]完全放棄了采用VTI介質(zhì)聲波近似的思路,提出在方程推導(dǎo)中保留橫波分量,通過(guò)引入橫波分量的方式能有效消除VTI介質(zhì)聲波近似產(chǎn)生的不穩(wěn)定問(wèn)題,FLETCHER等推導(dǎo)的2D VTI介質(zhì)縱橫波耦合qP波方程為:

(3)

一般情況下α=1,本文σ取經(jīng)驗(yàn)常數(shù)0.75[7]。

本文將FLETCHER等[7]和李博等[24]通過(guò)引入橫波分量的方式推導(dǎo)得到的VTI介質(zhì)方程稱(chēng)為VTI介質(zhì)縱橫波耦合方程。

VTI介質(zhì)波場(chǎng)模擬首先必須滿足物理穩(wěn)定性條件,即P-SV波VTI介質(zhì)精確相速度頻散方程((1)式)中相速度的平方不能小于0[25]:

(4)

由(4)式得到TI介質(zhì)波場(chǎng)模擬中的物理穩(wěn)定性條件:

(5)

特別地,在VTI介質(zhì)聲波近似情況下(vSZ=0),即f=1,由(5)式可推導(dǎo)出TI介質(zhì)聲波近似情況下的物理穩(wěn)定性條件為:

(6)

本文設(shè)計(jì)了一個(gè)楔形2D VTI介質(zhì)模型進(jìn)行算法測(cè)試,其主要特點(diǎn)是模型左部有一個(gè)楔形構(gòu)造,上層為各向同性介質(zhì)(圖1)。模型參數(shù)中橫向網(wǎng)格數(shù)和縱向網(wǎng)格數(shù)都為300,模型橫向網(wǎng)格間距和縱向網(wǎng)格間距都為10m。因?yàn)楸纠鹪粗糜诟飨蛲缘暮Ｋ畬又?VTI數(shù)值模擬時(shí)不會(huì)出現(xiàn)影響qP波場(chǎng)傳播的偽SV波噪聲,此種情況下不需要采取在震源處添加各向同性薄層消除偽SV波噪聲的策略[5,26]。

數(shù)值模擬時(shí),雷克震源子波主頻取25Hz,傳播時(shí)間為3s,時(shí)間采樣間隔為0.5ms,吸收邊界條件采用優(yōu)化海綿吸收型[17],分別采用VTI介質(zhì)聲波近似方程(ZHOU等的方程)和引入橫波分量的VTI介質(zhì)縱橫波耦合方程(FLETCHER方程)。

對(duì)波場(chǎng)進(jìn)行分析,從t=1.0s和t=1.5s時(shí)刻的波場(chǎng)快照(圖2,圖3)看出,ZHOU等的方程都能保持穩(wěn)定傳播而且沒(méi)有出現(xiàn)影響qP波場(chǎng)傳播的偽SV噪聲,而FLETCHER等的方程由于引入了橫波分量解決了地震波在各向異性介質(zhì)中傳播的不穩(wěn)定性問(wèn)題,但是也帶來(lái)了較多的偽SV噪聲(圖2b和圖3b中紅色箭頭處)。因此VTI介質(zhì)情況下更適合采用基于VTI介質(zhì)聲波近似方程,因此在二維和三維VTI介質(zhì)正演中,我們將采用穩(wěn)定的ZHOU等的方程。原因是VTI介質(zhì)的對(duì)稱(chēng)軸角度為零,在VTI介質(zhì)聲波近似假設(shè)下基本不存在穩(wěn)定性問(wèn)題,另外VTI介質(zhì)聲波近似方程相比VTI介質(zhì)縱橫波耦合方程而言偽SV噪聲也更少。

綜上所述,基于VTI介質(zhì)精確相速度頻散方程推導(dǎo)出的VTI介質(zhì)波動(dòng)方程可以分為兩類(lèi),一類(lèi)為基于VTI介質(zhì)聲波近似的VTI介質(zhì)聲波近似方程,另一類(lèi)為考慮橫波分量的VTI介質(zhì)縱橫波耦合方程,其中VTI介質(zhì)聲波近似方程可以視為VTI介質(zhì)縱橫波耦合方程在橫波速度為零時(shí)的特例。通過(guò)研究發(fā)現(xiàn),VTI介質(zhì)聲波近似方程更適用于VTI介質(zhì)。

圖1 楔形2D VTI介質(zhì)模型各向異性參數(shù)a 縱波速度vP(km/s); b 縱波各向異性參數(shù)ε; c 變異系數(shù)δ

圖2 楔形2D VTI介質(zhì)模型t=1.0s時(shí)的波場(chǎng)快照a ZHOU等的方程; b FLETCHER方程

圖3 楔形2D VTI介質(zhì)模型t=1.5s時(shí)的波場(chǎng)快照a ZHOU等的方程; b FLETCHER方程

1.2 2D VTI介質(zhì)時(shí)間域UPML波動(dòng)方程推導(dǎo)

VTI介質(zhì)波動(dòng)方程的特點(diǎn)是不含比較復(fù)雜的交叉導(dǎo)數(shù)項(xiàng),因此可以采用計(jì)算效率更高的UPML吸收邊界條件來(lái)實(shí)現(xiàn)。針對(duì)VTI介質(zhì)波動(dòng)方程上述特點(diǎn),我們以ZHOU等的二階2D VTI介質(zhì)聲波近似時(shí)間域波動(dòng)方程(公式(2))為基礎(chǔ),推導(dǎo)2D VTI介質(zhì)時(shí)間域UPML波動(dòng)方程。

首先對(duì)二階2D VTI介質(zhì)時(shí)間域波動(dòng)方程(公式(2))做傅里葉變換,得到其頻率域波動(dòng)方程的形式:

(7)

再引入X和Z方向的頻率域PML拉伸函數(shù)ξx,ξz:

(8)

其中,dx,dz為PML衰減因子,其表達(dá)式為[27]:

(9)

式中:LPML是PML邊界層的厚度;θ是計(jì)算點(diǎn)距PML邊界的距離;Rcoeff是理論反射系數(shù)值,本文Rcoeff取0.0001。

(10)

對(duì)公式(10)求導(dǎo)并結(jié)合公式(8)推導(dǎo)得到:

(11)

在正常計(jì)算區(qū)域內(nèi),加PML衰減項(xiàng)的聲波波動(dòng)方程與不加PML衰減項(xiàng)的聲波波動(dòng)方程是一致的,因?yàn)樵谡Ｓ?jì)算區(qū)域內(nèi)di=0。如果PML計(jì)算區(qū)域內(nèi)di>0,則在PML區(qū)域內(nèi)聲波波動(dòng)方程的平面波解為指數(shù)衰減的形式:

(12)

當(dāng)波從內(nèi)部區(qū)域過(guò)渡到PML邊界層后波場(chǎng)逐漸衰減。

將公式(11)代入2D VTI介質(zhì)頻率域波動(dòng)方程(公式(7)),并引入6個(gè)輔助波場(chǎng)項(xiàng)Epx,Epz,Eqx,Hpx,Hpz,Hqx進(jìn)行推導(dǎo),最后得到2D VTI介質(zhì)頻率域UPML波動(dòng)方程,將此方程反傅里葉變換到時(shí)間域即可得到2D VTI介質(zhì)時(shí)間域UPML波動(dòng)方程:

(13)

其中6個(gè)輔助波場(chǎng)項(xiàng)Epx,Epz,Eqx,Hpx,Hpz,Hqx對(duì)應(yīng)的控制方程為:

公式(13)的右端項(xiàng)由兩部分組成,一部分是正常項(xiàng)(不含輔助波場(chǎng)項(xiàng)),另一部分是衰減項(xiàng)(含輔助波場(chǎng)項(xiàng))。如果衰減項(xiàng)為零,則UPML波動(dòng)方程就退化為常規(guī)的2D VTI介質(zhì)波動(dòng)方程(公式(2))。

通過(guò)一系列公式推導(dǎo)可以發(fā)現(xiàn)不同2D VTI介質(zhì)方程的UPML方程的形式不一樣,差分離散形式和輔助波場(chǎng)項(xiàng)的表達(dá)形式也不一致,其具體形式根據(jù)不同的VTI介質(zhì)方程形式而定。本文只討論二維VTI介質(zhì)情形,但是這種方法也可以擴(kuò)展到三維VTI介質(zhì)情形(詳細(xì)推導(dǎo)見(jiàn)附錄A),在此不再展開(kāi)討論。

1.3 2D VTI介質(zhì)UPML波動(dòng)方程數(shù)值實(shí)現(xiàn)

從二維UPML波動(dòng)方程數(shù)值實(shí)現(xiàn)示意圖(圖4)可以看出,UPML具體實(shí)現(xiàn)的時(shí)候不再有如分裂式完全匹配層(SPML)波動(dòng)方程一樣的邊和角的概念,只有方向的概念,我們只需要考慮X和Z2個(gè)方向即上下左右4個(gè)PML區(qū)域波場(chǎng)的衰減,而且UPML不需要對(duì)波場(chǎng)項(xiàng)進(jìn)行方向分裂,相比二維SPML需要考慮8種情形并對(duì)波場(chǎng)項(xiàng)進(jìn)行方向分裂而言,UPML簡(jiǎn)化了計(jì)算量,也減少了存儲(chǔ)空間,編程工作量大大減少。

圖4 二維UPML數(shù)值實(shí)現(xiàn)示意

2D VTI介質(zhì)UPML波動(dòng)方程具體數(shù)值實(shí)現(xiàn)過(guò)程如下:

1) 在所有區(qū)域進(jìn)行波場(chǎng)計(jì)算。采用12階高階規(guī)則網(wǎng)格有限差分算法對(duì)公式(13)的正常項(xiàng)進(jìn)行波場(chǎng)計(jì)算。

2) 在PML區(qū)域內(nèi)進(jìn)行波場(chǎng)更新。在PML計(jì)算區(qū)域內(nèi)對(duì)公式(13)的衰減項(xiàng)進(jìn)行有限差分計(jì)算并更新PML計(jì)算區(qū)域內(nèi)的波場(chǎng),共考慮4個(gè)PML區(qū)域情形,其中4個(gè)角域需要計(jì)算2次。

3) 在PML區(qū)域內(nèi)更新輔助波場(chǎng)項(xiàng)Epx,Epz,Eqx,Hpx,Hpz,Hqx。在PML計(jì)算區(qū)域內(nèi)對(duì)公式(13)的6個(gè)輔助波場(chǎng)項(xiàng)的控制方程進(jìn)行有限差分計(jì)算并更新PML計(jì)算區(qū)域內(nèi)的輔助波場(chǎng)項(xiàng),共考慮4個(gè)PML區(qū)域情形,其中4個(gè)角域需要計(jì)算2次。

2 數(shù)值計(jì)算實(shí)例

2.1 簡(jiǎn)單均勻2D VTI介質(zhì)模型

為了驗(yàn)證本文提出的VTI介質(zhì)UPML吸收邊界條件的有效性,首先采用簡(jiǎn)單均勻2D VTI介質(zhì)模型進(jìn)行UPML點(diǎn)源響應(yīng)測(cè)試。模型參數(shù)中橫向網(wǎng)格數(shù)和縱向網(wǎng)格數(shù)都為256,模型橫向網(wǎng)格間距和縱向網(wǎng)格間距都為10m,速度為常速3000m/s。VTI介質(zhì)的各向異性參數(shù)ε=0.3,δ=0.1。數(shù)值模擬時(shí),雷克震源子波主頻取25Hz,震源位置設(shè)置在模型的中心,傳播時(shí)間為1.5s,采樣間隔為1.5ms,采用12階高階規(guī)則網(wǎng)格有限差分格式進(jìn)行模擬。采用UPML吸收邊界條件,PML邊界吸收層數(shù)分別設(shè)為10層、20層和30層。圖5為3種不同PML吸收層數(shù)情形下,t=0.54s時(shí)刻的正演波場(chǎng)快照,可以看出,PML吸收層數(shù)為10層時(shí)邊界反射吸收效果不理想(圖5a);PML吸收層數(shù)為20層時(shí)邊界反射吸收效果明顯改善,但是仍然存在著微弱的邊界反射(圖5b);當(dāng)PML吸收層數(shù)為30層時(shí)邊界反射已經(jīng)基本被消除(圖5c)。因此本例選定PML層數(shù)為30層。圖6a到圖6f分別是UPML波動(dòng)方程在各個(gè)時(shí)刻的波場(chǎng)快照,由圖可見(jiàn)UPML對(duì)邊界反射吸收效果明顯。

2.2 復(fù)雜BP2007 2D VTI介質(zhì)海洋標(biāo)準(zhǔn)模型

以復(fù)雜BP2007 2D VTI介質(zhì)海洋標(biāo)準(zhǔn)模型為例,并且與優(yōu)化海綿吸收邊界條件[17]進(jìn)行對(duì)比，進(jìn)一步說(shuō)明VTI介質(zhì)UPML吸收邊界條件的應(yīng)用效果。數(shù)值模擬時(shí)采用截取的部分BP2007 2D VTI介質(zhì)海洋標(biāo)準(zhǔn)模型(圖7),該部分主要模擬近海岸的地質(zhì)構(gòu)造,其特點(diǎn)是左邊有一個(gè)被VTI各向異性地層包圍的高速各向同性鹽丘。模型參數(shù)包括縱波速度vP,Thomsen縱波各向異性參數(shù)ε和Thomsen變異系數(shù)δ,模型區(qū)域大小為:橫向網(wǎng)格數(shù)nx=601,縱向網(wǎng)格數(shù)nz=361,模型橫向網(wǎng)格間距和縱向網(wǎng)格間距都為6.25m,其中海水層為各向同性介質(zhì)。數(shù)值模擬時(shí),雷克震源子波主頻取25Hz,采用UPML吸收邊界條件,PML邊界設(shè)為30層,傳播時(shí)間為2s,時(shí)間采樣間隔為0.4ms,采用12階高階規(guī)則網(wǎng)格有限差分格式進(jìn)行模擬。因?yàn)楸纠鹪粗糜诟飨蛲缘暮Ｋ畬又?VTI數(shù)值模擬時(shí)不會(huì)出現(xiàn)影響qP波場(chǎng)傳播的偽SV波噪聲,此種情況也不需要采取在震源處添加各向同性薄層來(lái)消除偽SV波噪聲的策略[5,26]。從BP2007 2D VTI介質(zhì)海洋標(biāo)準(zhǔn)模型的正演單炮記錄來(lái)看,UPML消除邊界反射效果良好(圖8a),而優(yōu)化海綿吸收邊界條件仍然有著比較強(qiáng)的邊界反射(圖8b 紅色箭頭標(biāo)示),UPML吸收邊界條件明顯優(yōu)于優(yōu)化海綿吸收邊界條件的邊界反射吸收效果。

圖5 二維UPML波動(dòng)方程t=0.54s時(shí)刻不同PML吸收層數(shù)情況下正演波場(chǎng)快照a PML吸收層數(shù)為10層; b PML吸收層數(shù)為20層; c PML吸收層數(shù)為30層

圖6 均勻2D VTI介質(zhì)模型采用二維UPML波動(dòng)方程不同時(shí)刻的正演波場(chǎng)快照a t=0.39s; b t=0.45s; c t=0.51s; d t=0.57s; e t=0.90s; f t=1.35s

圖7 截取的部分BP2007 2D VTI介質(zhì)海洋標(biāo)準(zhǔn)模型a 縱波速度vP; b 縱波各向異性參數(shù)ε; c 變異系數(shù)δ

圖8 BP2007 2D VTI介質(zhì)海洋標(biāo)準(zhǔn)模型正演單炮記錄a UPML吸收邊界條件方法; b 優(yōu)化海綿吸收邊界條件方法

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

本文對(duì)常見(jiàn)的幾種VTI介質(zhì)波動(dòng)方程進(jìn)行了研究,發(fā)現(xiàn)基于VTI介質(zhì)精確相速度頻散方程推導(dǎo)出的VTI介質(zhì)波動(dòng)方程可以分為基于VTI介質(zhì)聲波近似的VTI介質(zhì)聲波近似方程和考慮橫波分量的VTI介質(zhì)縱橫波耦合方程,就波場(chǎng)傳播穩(wěn)定性和減少偽橫波噪聲而言,VTI介質(zhì)聲波近似方程更適用于VTI介質(zhì)情形。從基于VTI介質(zhì)聲波近似的二階波動(dòng)方程出發(fā),推導(dǎo)得到了VTI介質(zhì)UPML波動(dòng)方程。研究發(fā)現(xiàn)不同2D VTI介質(zhì)UPML方程的形式不一樣,差分離散形式和輔助波場(chǎng)項(xiàng)的表達(dá)形式也不一致,其具體形式根據(jù)不同的VTI方程形式而定。利用簡(jiǎn)單均勻2D VTI介質(zhì)模型和復(fù)雜BP2007 2D VTI介質(zhì)海洋標(biāo)準(zhǔn)模型對(duì)2D VTI介質(zhì)UPML方程進(jìn)行了數(shù)值模擬,結(jié)果表明,本文研究的UPML吸收邊界條件能取得很好的人工邊界反射吸收效果,相比優(yōu)化海綿吸收邊界條件,其人工邊界反射吸收效果明顯提高。在附錄中,我們給出了3D VTI介質(zhì)時(shí)間域UPML波動(dòng)方程的推導(dǎo)結(jié)果,從而可以將本文方法進(jìn)一步推廣到三維VTI介質(zhì)的高精度模擬中。

[1] WARNER M,RATCLIFE A,NANGOO T,et al.Anisotropic 3D full-waveform inversion[J].Geophysics,2013,78(2):R59-R80

[2] THOMSEN L.Weak elastic anisotropy[J].Geophysics,1986,51(10):1954-1966

[3] TSVANKIN I.P-wave signatures and notation for transversely isotropic media:an overview[J].Geophysics,1996,61(2):467-483

[4] ZHOU H B,ZHANG G Q,BLOOR R.An anisotropic acoustic wave equation for VTI media[J].Expanded Abstracts of 68thEAGE Annual Conference,2006:H033

[5] GRECHKA V,ZHANG L,RECTOR J W.Shear waves in acoustic anisotropic media[J].Geophysics,2004,69(2):576-582

[6] BAKKER P M,DUVENECK E.Stability analysis for acoustic wave propagation in tilted TI media by finite differences[J].Geophysical Journal International,2011,185(2):911-921

[7] FLETCHER R P,DU X,FOWLER P J.Reverse time migration in tilted transversely isotropic (TTI) media[J].Geophysics,2009,74(6):179-187

[8] DUVENECK E,MILCIK P,BAKKER P M.Acoustic VTI wave equations and their application for anisotropic reverse-time migration[J].Expanded Abstracts of 78thAnnual Internat SEG Mtg,2008:2186-2190

[9] ZHANG Y,ZHANG H Z.A stable TTI reverse time migration and its implementation[J].Expanded Abstracts of 79thAnnual Internat SEG Mtg,2009:2794-2798

[10] ZHANG Y,ZHANG H Z,ZHANG G Q.A stable TTI reverse time migration and its implementation[J].Geophysics,2011,76(3):WA3-WA11

[11] HESTHOLM S.Acoustic VTI modeling using high-order finite differences[J].Geophysics,2009,74(5):T67-T73

[12] LIU Y,SEN M K.Acoustic VTI modeling with a time-space domain dispersion-relation-based finite-difference scheme[J].Geophysics,2010,75(3):A11-A17

[13] FOWLER P J,DU X,FLETCHER R P.Coupled equations for reverse time migration in transversely isotropic media[J].Geophysics,2010,75(1):S11-S22

[14] CLAYTON R,ENQUIST B.Absorbing boundary conditions for acoustic and elastic wave equations[J].Bulletin of the Seismological Society of America,1977,67(6):1529-1540

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

[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] BORDING R P.Finite difference modeling-nearly optimal sponge boundary conditions[J].Expanded Abstracts of 74thAnnual Internat SEG Mtg,2004:1921-1924

[18] 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

[20] WANG T,TANG X M.Finite-difference modeling of elastic wave propagation:A nonsplitting perfectly matched layer approach[J].Geophysics,2003,68(5):1749-1755

[21] DROSSAERT F H,GIANNOPOULOS A.A nonsplit complex frequency-shifted PML based on recursive integration for FDTD modeling of elastic waves[J].Geophysics,2007,72(2):T9-T17

[22] QIN Z,LU M H,ZHENG X D,et al.The implementation of an improved NPML absorbing boundary condition in elastic wave modeling[J].Applied Geophysics,2009,6(2):113-121

[23] ALKHALIFAH T.Acoustic approximations for processing in transversely isotropic media[J].Geophysics,1998,63(2):623-631

[24] 李博,李敏,劉紅偉,等.TTI介質(zhì)有限差分逆時(shí)偏移的穩(wěn)定性探討[J].地球物理學(xué)報(bào),2012,55(4):1366-1375 LI B,LI M,LIU H W,et al.Stability of reverse time migration in TTI media[J].Chinese Journal of Geophysics,2012,55(4):1366-1375

[25] 李博.地震偏移成像和GPU超算技術(shù)[D].北京:中國(guó)科學(xué)院大學(xué),2011 LI B.Seismic migration and GPU supercomputing technology[D].Beijing:University of Chinese Academy of Sciences,2011

[26] 張衡,劉洪,唐祥德,等.基于平均導(dǎo)數(shù)優(yōu)化方法的VTI介質(zhì)頻率空間域正演[J].地球物理學(xué)報(bào),2015,58(9):3306-3316 ZHANG H,LIU H,TANG X D,et al.Forward modeling of VTI media in the frequency-space domain based on an average-derivative optimal method[J].Chinese Journal of Geophysics,2015,58(9):3306-3316

[27] COLLINO F,TSOGKA C.Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media[J].Geophysics,2001,66(1):294-307

附錄A 3D VTI介質(zhì)時(shí)間域UPML波動(dòng)方程

3D VTI介質(zhì)時(shí)間域UPML波動(dòng)方程的推導(dǎo)基于ZHOU等的3D VTI介質(zhì)聲波近似波動(dòng)方程:

(A1)

從公式(A1)出發(fā),對(duì)方程做傅里葉變換推導(dǎo)得到其頻率域波動(dòng)方程的形式:

(A2)

引入X,Y和Z方向的頻率域PML拉伸函數(shù)ξx,ξy,ξz:

(A3)

其中:dx,dy,dz為PML衰減因子,其表達(dá)式為:

(A4)

式中:LPML是PML邊界層的厚度;θ是計(jì)算點(diǎn)距PML邊界的距離;Rcoeff是理論反射系數(shù)值,取0.0001。

(A5)

對(duì)公式(A5)求導(dǎo)并結(jié)合公式(A3)推導(dǎo)得到:

(A6)

將公式(A6)代入3D VTI介質(zhì)頻率域波動(dòng)方程(公式(A2)),并引入10個(gè)輔助波場(chǎng)項(xiàng)Epx,Epy,Epz,Eqx,Eqy,Hpx,Hpy,Hpz,Hqx,Hqy,最后得到3D VTI介質(zhì)頻率域UPML波動(dòng)方程,將此方程反傅里葉變換到時(shí)間域即可得到3D VTI介質(zhì)時(shí)間域UPML波動(dòng)方程:

(A7)

其中10個(gè)輔助波場(chǎng)項(xiàng)Epx,Epy,Epz,Eqx,Eqy,Hpx,Hpy,Hpz,Hqx,Hqy對(duì)應(yīng)的控制方程為:

分析公式(A7)的右端項(xiàng),其由兩部分組成,一部分是正常項(xiàng)(不含輔助波場(chǎng)項(xiàng)),另一部分是衰減項(xiàng)(含輔助波場(chǎng)項(xiàng))。如果衰減項(xiàng)為零,則UPML波動(dòng)方程就退化為常規(guī)的3D VTI介質(zhì)波動(dòng)方程(公式(A1))。

通過(guò)一系列公式推導(dǎo)研究可以發(fā)現(xiàn)不同3D VTI介質(zhì)方程的UPML方程的形式不一樣,差分離散形式和輔助波場(chǎng)項(xiàng)的表達(dá)形式也不一致,其具體形式根據(jù)不同的VTI介質(zhì)方程形式而定。

對(duì)于三維UPML數(shù)值實(shí)現(xiàn)而言,只需要考慮X,Y和Z3個(gè)方向即上下左右前后共6個(gè)PML區(qū)域波場(chǎng)的衰減,而且不需要對(duì)波場(chǎng)項(xiàng)進(jìn)行方向分裂,相比三維SPML需要考慮26種情形并對(duì)波場(chǎng)項(xiàng)進(jìn)行方向分裂而言,UPML簡(jiǎn)化了計(jì)算量,也減少了存儲(chǔ)空間,編程工作量大大減少。

(編輯：朱文杰)

The research on unsplit PML absorbing boundary conditions of acoustic equation for VTI media

ZHANG Heng1,LIU Hong2,LI Bo3,DING Renwei4,LI Liqing1,LI Fuyuan1

(1.MLRKeyLaboratoryofMarineMineralResources,GuangzhouMarineGeologicalSurvey,ChinaGeologicalSurvey,Guangzhou510075,China;2.KeyLaboratoryofPetroleumResourcesResearch,InstituteofGeologyandGeophysics,ChineseAcademyofSciences,Beijing100029,China;3.SinopecGeophysicalResearchInstitute,Nanjing211103,China;4.CollegeofEarthScienceandEngineering,ShandongUniversityofScienceandTechnology,Qingdao266590,China)

We study unsplit perfectly matched layer (UPML) absorbing boundary conditions for vertical transversely isotropic (VTI) media based on the characteristics of VTI media acoustic wave equation.We firstly summarize several common VTI media acoustic wave equations and compare these wave equations in detail from the perspective of wavefield propagation stability and pseudo SV-wave noise reduction for VTI media.We conclude that the VTI acoustic approximation equation is applicable for VTI media.Then we derive the UPML wave equation from the VTI acoustic approximation second-order wave equation.Afterwards we give the concrete numerical implementation method.The numerical modeling examples show that the UPML boundary condition can absorb the artificial boundary reflection very well.The boundary reflection absorbing effect of the proposed UPML algorithm is much better compared with the optimal sponge boundary condition.

VTI media,unsplit perfectly matched layer (UPML),pseudo SV-wave noise,acoustic approximation equation,second-order wave equation

2016-02-15;改回日期：2016-05-24。

張衡(1987—),男,博士,工程師,主要從事復(fù)雜各向異性介質(zhì)地震波場(chǎng)正演數(shù)值模擬和全波形反演研究工作。

國(guó)土資源部海底礦產(chǎn)資源重點(diǎn)實(shí)驗(yàn)室開(kāi)放基金(No.KLMMR-2015-A-14)、國(guó)家自然科學(xué)基金項(xiàng)目(41604110)、國(guó)家高技術(shù)研究發(fā)展計(jì)劃(863計(jì)劃)項(xiàng)目(2013AA092501)和國(guó)家科技重大專(zhuān)項(xiàng)(2011ZX05003-003)聯(lián)合資助。

This research is financially supported by the Open Fund of MLR Key Laboratory of Marine Mineral Resources (Grant No.KLMMR-2015-A-14),the National Natural Science Foundation of China (Grant No.41604110),the National High-tech R & D Program of China (863 Program) (Grant No.2013AA092501) and the National Science and Technology Major Project of China (Grant No.2011ZX05003-003).

P631

A

1000-1441(2016)06-0781-12

10.3969/j.issn.1000-1441.2016.06.002