張衡， 劉洪， 唐祥德， 王洋， 張寶金

1 國土資源部海底礦產(chǎn)資源重點(diǎn)實(shí)驗(yàn)室，廣州海洋地質(zhì)調(diào)查局， 廣州 510075 2 中國科學(xué)院地質(zhì)與地球物理研究所，中國科學(xué)院油氣資源研究重點(diǎn)實(shí)驗(yàn)室， 北京 100029 3 中國科學(xué)院大學(xué)， 北京 100049

?

基于平均導(dǎo)數(shù)優(yōu)化方法的VTI介質(zhì)頻率空間域正演

張衡1， 劉洪2*， 唐祥德2,3， 王洋2,3， 張寶金1

1 國土資源部海底礦產(chǎn)資源重點(diǎn)實(shí)驗(yàn)室，廣州海洋地質(zhì)調(diào)查局， 廣州 510075 2 中國科學(xué)院地質(zhì)與地球物理研究所，中國科學(xué)院油氣資源研究重點(diǎn)實(shí)驗(yàn)室， 北京 100029 3 中國科學(xué)院大學(xué)， 北京 100049

本文提出了一種新的基于平均導(dǎo)數(shù)優(yōu)化方法(average-derivative optimal method，簡稱ADM)的二維VTI介質(zhì)qP波波動(dòng)方程頻率空間域二階9點(diǎn)格式，這種新算法將二維VTI介質(zhì)qP波波動(dòng)方程中中心空間導(dǎo)數(shù)項(xiàng)的差分近似表示為正交方向上3個(gè)網(wǎng)格點(diǎn)的加權(quán)平均形式.通過最小二乘優(yōu)化方法求取空間導(dǎo)數(shù)項(xiàng)和加速度項(xiàng)的加權(quán)優(yōu)化系數(shù)從而使數(shù)值頻散達(dá)到極小化，每個(gè)波長所需要的網(wǎng)格點(diǎn)數(shù)在1%的誤差范圍內(nèi)僅為3.57個(gè)網(wǎng)格點(diǎn)數(shù)，而VTI介質(zhì)常規(guī)9點(diǎn)差分格式在相同的誤差范圍內(nèi)則需要約12個(gè)網(wǎng)格點(diǎn)數(shù)，新方法的計(jì)算精度明顯提高.復(fù)雜BP2007 2D VTI海洋標(biāo)準(zhǔn)模型數(shù)值模擬結(jié)果也驗(yàn)證了本文VTI介質(zhì)9點(diǎn)ADM算法的有效性和準(zhǔn)確性.

VTI介質(zhì)； 平均導(dǎo)數(shù)優(yōu)化方法； 頻率空間域

1 引言

地震各向異性在地球介質(zhì)中是普遍存在的(Tsvankin, 2001; Tsvankin et al．，2010)，如墨西哥灣、西非洲海岸和中國南海都有大量的各向異性地層存在，而橫向各向同性(TI)介質(zhì)是最常見的一種各向異性介質(zhì)(Thomsen，1986).

Thomsen(1986)提出了著名的弱各向異性理論，Tsvankin(1996)基于Thomsen(1986)的弱各向異性理論推導(dǎo)得到P-SV波VTI介質(zhì)相速度頻散方程，這是后來很多學(xué)者研究TI介質(zhì)正演模擬和偏移成像的出發(fā)點(diǎn).

TI介質(zhì)正演數(shù)值模擬是TI介質(zhì)逆時(shí)偏移(張巖和吳國忱，2013)和TI介質(zhì)全波形反演(Warner et al．，2013；Alkhalifah，2014)的基礎(chǔ).TI介質(zhì)正演數(shù)值模擬一般在時(shí)間域進(jìn)行，也可以在頻率空間域進(jìn)行，頻率空間域正演的優(yōu)勢在于可以多炮同時(shí)模擬且不存在累積誤差，此時(shí)頻率空間域波動(dòng)方程的求解問題變成了一個(gè)大型稀疏線性方程組的求解問題.近年來若干學(xué)者對(duì)TI介質(zhì)頻率空間域正演進(jìn)行了詳細(xì)的研究(Grini et al．，2007；Operto et al．，2007，2009；吳國忱等，2007；李桂花等，2011).吳國忱和梁鍇(2005，2007)通過引入25點(diǎn)加權(quán)優(yōu)化差分算子，實(shí)現(xiàn)了VTI介質(zhì)的頻率空間域正演模擬.吳國忱等(2007)將這種方法推廣到TTI介質(zhì)，推導(dǎo)出二維TTI介質(zhì)頻率空間域彈性波動(dòng)方程，通過引入25點(diǎn)加權(quán)優(yōu)化差分算子，采用高斯牛頓方法(Lines and Treitel，1984；Min et al．，2000)求取優(yōu)化系數(shù)，從而在頻率空間域模擬了彈性波在TTI介質(zhì)中的傳播過程.杜向東等(2009)將吳國忱和梁鍇(2005，2007)的方法推廣到三維，通過引入125點(diǎn)優(yōu)化差分算子，實(shí)現(xiàn)了三維VTI介質(zhì)的頻率空間域正演模擬.李桂花等(2011)通過引入25點(diǎn)優(yōu)化差分算子，實(shí)現(xiàn)了黏彈性VTI介質(zhì)的頻率空間域正演模擬.Grini等(2007)和Operto等(2007，2009)研究了黏聲TTI介質(zhì)一階方程形式的混合網(wǎng)格頻率空間域正演方法.Operto等(2014)將VTI介質(zhì)四階方程分解為一個(gè)橢圓型二階方程和非橢圓校正項(xiàng)的形式，實(shí)現(xiàn)了三維VTI介質(zhì)的頻率空間域正演模擬.

以相速度誤差1%為標(biāo)準(zhǔn)，VTI介質(zhì)qP波波動(dòng)方程常規(guī)9點(diǎn)頻率空間域差分格式需要每個(gè)波長12個(gè)網(wǎng)格點(diǎn)才能正確模擬，計(jì)算精度差.Chen(2012)提出了一種平均導(dǎo)數(shù)優(yōu)化方法的新思路，能有效提高計(jì)算精度，而且具有很好的靈活性，能適用于橫縱向網(wǎng)格不等間距的情形.張衡等(2014)基于Chen(2012)利用平均導(dǎo)數(shù)優(yōu)化方法構(gòu)建差分格式的新思路發(fā)展了一種新的基于平均導(dǎo)數(shù)優(yōu)化方法的聲波方程頻率空間域四階25點(diǎn)格式，能將每個(gè)波長所需要的網(wǎng)格點(diǎn)數(shù)降低到2.78個(gè)，相比四階常規(guī)9點(diǎn)方法明顯提高了計(jì)算精度，也克服了Shin和Sohn(1998)提出的四階25點(diǎn)混合格式只能適用于橫縱向網(wǎng)格等間距情形的局限性.本文將平均導(dǎo)數(shù)優(yōu)化方法進(jìn)一步推廣到VTI介質(zhì)qP波波動(dòng)方程，研究了基于平均導(dǎo)數(shù)優(yōu)化方法的VTI介質(zhì)頻率空間域二階9點(diǎn)算法.

本文首先推導(dǎo)了二維VTI介質(zhì)頻率空間域qP波波動(dòng)方程，然后提出了基于平均導(dǎo)數(shù)優(yōu)化方法(ADM)的二維VTI介質(zhì)qP波波動(dòng)方程頻率空間域二階9點(diǎn)格式.隨后通過最小二乘優(yōu)化方法求取優(yōu)化系數(shù)并做了頻散分析，表明VTI ADM 9點(diǎn)算法能夠?qū)⒚總€(gè)波長所需的網(wǎng)格點(diǎn)數(shù)降低到3.57個(gè)，相比VTI常規(guī)9點(diǎn)格式每個(gè)波長需要12個(gè)網(wǎng)格點(diǎn)數(shù)明顯提高了計(jì)算精度.最后進(jìn)行數(shù)值模擬，測試表明VTI ADM 9點(diǎn)算法能實(shí)現(xiàn)高精度模擬.

2 二維VTI介質(zhì)頻率空間域qP波波動(dòng)方程

利用VTI介質(zhì)彈性張量矩陣，結(jié)合彈性動(dòng)力學(xué)的本構(gòu)方程(廣義虎克定律)、運(yùn)動(dòng)微分方程(牛頓第二定律的微分形式)和幾何方程(描述位移與應(yīng)變關(guān)系)可得到VTI介質(zhì)中的彈性波動(dòng)方程.將平面波方程U=Aei(kxx+kyy+kzz-ω t)代入VTI介質(zhì)中的彈性波動(dòng)方程并忽略體力項(xiàng)，推導(dǎo)得到Christoffel方程(吳國忱，2006；程玖兵等，2013).

由Christoffel方程出發(fā)，并利用TI介質(zhì)中的Thomsen各向異性參數(shù)ε和δ推導(dǎo)得到P-SV波VTI介質(zhì)相速度頻散方程(Tsvankin，1996；Alkhalifah，1998，2000)：

(1)

式(1)中Vp(θ)表示相速度，θ表示相速度傳播角度，f中VSz表示橫波速度，VPz表示縱波速度，ε和δ表示Thomsen各向異性參數(shù)(Thomsen，1986).

采用Alkhalifah(1998)提出的TI聲波近似思想(即假設(shè)橫波速度VSz為零)基于式(1)即可推導(dǎo)得到VTI介質(zhì)頻率空間域qP波波動(dòng)方程.

(2)

(3)

式(3)從波數(shù)域變換到空間域即可得到二維VTI介質(zhì)頻率空間域qP波波動(dòng)方程：

(4)

其中，P為地震波場，ω為角頻率，Vp0為VTI介質(zhì)qP波速度.

當(dāng)Thomsen各向異性參數(shù)ε=δ時(shí)式(4)退化為二維橢圓各向異性介質(zhì)方程：

(5)

當(dāng)Thomsen各向異性參數(shù)ε=δ=0時(shí)式(4)退化為二維各向同性介質(zhì)方程：

(6)

因此二維橢圓各向異性介質(zhì)方程和二維各向同性介質(zhì)方程可看作二維VTI介質(zhì)頻率空間域qP波波動(dòng)方程的特例.

3 基于平均導(dǎo)數(shù)優(yōu)化方法的二維VTI介質(zhì)qP波波動(dòng)方程頻率空間域二階9點(diǎn)格式

式(4)表示的二維VTI介質(zhì)頻率空間域qP波波動(dòng)方程可表示成VTI常規(guī)9點(diǎn)格式的形式(圖1a)，其差分格式為：

×(4Pm,n-2(Pm,n+1+Pm,n-1+Pm+1,n+Pm-1,n)+Pm+1,n+1+Pm-1,n+1+Pm+1,n-1+Pm-1,n-1)=0．

(7)

但是VTI常規(guī)9點(diǎn)格式正演精度低，數(shù)值頻散嚴(yán)重(圖2)，需要發(fā)展新的優(yōu)化算法.本文基于Chen(2012)和張衡等(2014)采用平均導(dǎo)數(shù)優(yōu)化方法構(gòu)建聲波方程頻率空間域差分格式的思路，提出一種新的基于平均導(dǎo)數(shù)優(yōu)化方法的二維VTI介質(zhì)qP波波動(dòng)方程頻率空間域二階9點(diǎn)有限差分格式(圖1b).

×(4Pm,n-2(Pm,n+1+Pm,n-1+Pm+1,n+Pm-1,n)+Pm+1,n+1+Pm-1,n+1+Pm+1,n-1+Pm-1,n-1)

(8)

其中：

式(8)中α、β、c、d是加權(quán)系數(shù)，當(dāng)α=β=1且c=1,d=0時(shí)VTIADM9點(diǎn)差分格式(8)可寫成公式(7)的形式，即VTI常規(guī)9點(diǎn)差分是VTIADM9點(diǎn)格式的特例.本文通過引入VTIADM優(yōu)化方法，改善了數(shù)值頻散，提高了計(jì)算精度(圖2).

圖1 (a) VTI常規(guī)9點(diǎn)差分格式示意圖；(b) VTI ADM 9點(diǎn)差分格式示意圖Fig.1 (a) VTI conventional 9-point scheme；(b) VTI ADM 9-point scheme

4 系數(shù)優(yōu)化與頻散分析

采用經(jīng)典的頻散分析方法，引入一個(gè)平面波表示P(x,z,ω)=P0e-i(kxx+kzz)來進(jìn)行頻散分析的研究.

先考慮Δx≥Δz的情形，令r=Δx/Δz，將平面波P(x,z,ω)=P0e-i(kxx+kzz)代入VTIADM9點(diǎn)差分格式即式(8)可求取差分方程的相速度關(guān)系式，同時(shí)各向異性頻率空間域波動(dòng)方程頻散關(guān)系的求取還需要特別考慮求取波動(dòng)方程的相速度關(guān)系式(吳國忱，2006)，最后求得Δx≥Δz情形下的相速度頻散關(guān)系式如下：

(9)

其中：

式(9)中G表示每個(gè)波長所需的網(wǎng)格點(diǎn)數(shù). 令α=β=1且c=1,d=0時(shí)式(9)退化得VTI常規(guī)9點(diǎn)差分的相速度頻散公式：

(10)

其中：

系數(shù)α,β,c,d由求E(α,β,c,d)的極小得到：

(11)

優(yōu)化系數(shù)r=1r=1．5r=2r=2．5r=3r=3．5r=4α0．8055883170．8174548050．8205739330．8122823620．7971065140．7812683010．766020269β0．7728979590．8052414840．8143188900．8107244480．8062612950．8044679060．804287094c0．6283325460．6479670130．6413999530．6207767680．5969135810．5725920340．546171078d0．0942945230．0875583760．0922639810．1013067880．1115861500．1228480810．135721278

優(yōu)化系數(shù)r=1．5r=2r=2．5r=3r=3．5r=4α0．8052414840．8143188900．8107244480．8062612950．8044679060．804287094β0．8174548050．8205739330．8122823620．7971065140．7812683010．766020269c0．6479670130．6413999530．6207767680．5969135810．5725920340．546171078d0．0875583760．0922639810．1013067880．1115861500．1228480810．135721278

圖2 當(dāng)Δx≥Δz時(shí)對(duì)于不同常規(guī)9點(diǎn)格式和VTI ADM 9點(diǎn)格式相速度頻散曲線圖Fig.2 Phase velocity curves of the VTI conventional 9-point scheme and VTI ADM 9-point scheme

圖3 截取的BP2007 2D VTI海洋標(biāo)準(zhǔn)模型參數(shù)(a) VP； (b) ε； (c) δ.Fig.3 Parameters of truncated BP2007 2D VTI ocean standard model

圖4 (a)VTI ADM 9點(diǎn)格式計(jì)算的35 Hz單頻波場； (b)VTI常規(guī)9點(diǎn)格式計(jì)算得到的時(shí)間域地震記錄； (c)VTI ADM 9點(diǎn)格式計(jì)算得到的時(shí)間域地震記錄； (d)VTI介質(zhì)時(shí)間域12階高階有限差分方法計(jì)算得到的時(shí)間域地震記錄Fig.4 (a) 35 Hz monochromatic wavefield computed by VTI ADM 9-point scheme； (b) Time-domain seismograms computed with the VTI conventional 9-point scheme； (c) Time-domain seismograms computed with the VTI ADM 9-point scheme； (d) Time-domain seismograms computed with the time domain 12-order high order finite difference scheme for VTI media

5 數(shù)值模擬與分析

引入頻率空間域完全匹配層(Perfectly Matched Layer，PML)吸收邊界條件(Berenger，1994；吳國忱等，2007；張衡等，2014；Moreira et al．，2014)，構(gòu)建得到VTI ADM 9點(diǎn)格式PML波動(dòng)方程(附錄A)，并結(jié)合根據(jù)VTI模型不同的縱橫向采樣間隔求得的優(yōu)化系數(shù)，即可實(shí)現(xiàn)VTI ADM 9點(diǎn)格式頻率空間域數(shù)值模擬.下面通過一個(gè)復(fù)雜的BP2007 VTI海洋標(biāo)準(zhǔn)模型來分析VTI ADM 9點(diǎn)格式的正演精度和正演效率，VTI常規(guī)9點(diǎn)算法作為對(duì)比方法進(jìn)行比較.

本文數(shù)值模擬時(shí)采用截取的部分BP2007 2D VTI標(biāo)準(zhǔn)模型(圖3).截取的BP2007 2D VTI模型模擬的主要是近海岸的地質(zhì)構(gòu)造，此模型特點(diǎn)是左部有一個(gè)被各向異性地層包圍的高速各向同性鹽丘，模型參數(shù)包括縱波速度VP、Thomsen縱波各向異性參數(shù)ε、Thomsen變異系數(shù)δ，模型區(qū)域大小為橫向網(wǎng)格數(shù)Nx=401，縱向網(wǎng)格數(shù)Nz=361，模型橫向網(wǎng)格間距為10 m，縱向網(wǎng)格間距為4 m，此時(shí)橫縱向網(wǎng)格間距比達(dá)到了2.5，其中海水層為各向同性介質(zhì)(圖3).

數(shù)值模擬時(shí)，雷克震源子波主頻取15 Hz，震源位置設(shè)置在縱向第11層的橫向300個(gè)網(wǎng)格點(diǎn)處，坐標(biāo)為(3000 m，40 m)，檢波點(diǎn)排列置于海平面以下40 m(即電纜沉放深度)處，為測線所有點(diǎn)接收，每個(gè)網(wǎng)格點(diǎn)設(shè)置一個(gè)檢波器，全排列共401個(gè)檢波器.采用頻率空間域PML吸收邊界條件，PML層數(shù)設(shè)為40層，傳播時(shí)間為2 s，時(shí)間采樣間隔為0.004 s.

值得一提的是，實(shí)際測試表明VTI頻率空間域波場傳播過程中經(jīng)常會(huì)出現(xiàn)較強(qiáng)的偽橫波(SV波)噪聲，SV波干擾產(chǎn)生的原因是由于TI聲波近似導(dǎo)致的.式(12)表示TI聲波近似下SV波的相速度表達(dá)式，容易得出TI聲波近似(直接令VSz=0)雖然使得沿對(duì)稱軸方向(θ=0°)的橫波速度為零，但是其他傳播方向上的SV波相速度并不為零，從而產(chǎn)生了偽SV波噪聲 (Grechka et al．，2004).這種偽SV波噪聲相比qP波而言雖然體現(xiàn)出低速低能量的特征，但是在實(shí)際偏移和反演中仍會(huì)產(chǎn)生比較嚴(yán)重的干擾，而且可能會(huì)導(dǎo)致算法不穩(wěn)定，因此需要加以壓制和消除.對(duì)式(12)的分析表明若把震源激發(fā)的位置設(shè)置為橢圓各向異性介質(zhì)(ε=δ)或各向同性介質(zhì)(ε=δ=0)即可消除偽SV波噪聲.如通過采取在震源處添加一個(gè)各向同性薄層的策略即可有效實(shí)現(xiàn)頻率空間域qP波波場傳播過程中的偽SV波噪聲消除.

(12)

因?yàn)楸纠鹪粗糜诟飨蛲缘暮Ｋ畬又校琕TI數(shù)值模擬時(shí)不會(huì)出現(xiàn)影響qP波場傳播的偽橫波(SV波)噪聲，此種情況下不需要采取添加震源各向同性薄層消除偽SV波噪聲的策略.

分別采用VTI常規(guī)9點(diǎn)格式和VTI ADM 9點(diǎn)格式進(jìn)行頻率空間域正演模擬，橫縱向網(wǎng)格間距比為2.5情況下的VTI ADM 9點(diǎn)差分格式優(yōu)化系數(shù)如表1所示.從VTI ADM 9點(diǎn)格式頻率空間域模擬35 Hz的單頻波場切片(圖4a)可以看出頻率空間域qP波場特征明顯，其運(yùn)動(dòng)學(xué)特征得到了準(zhǔn)確的刻畫，在模型右側(cè)速度漸變帶，頻率空間域波場變化不大，而在左側(cè)陡峭鹽丘處由于速度和各向異性參數(shù)變化復(fù)雜，入射波和反射波的相互干涉導(dǎo)致了反射界面上頻率空間域波場的劇烈擾動(dòng).

對(duì)所有頻率的頻率空間域波場進(jìn)行反傅里葉變換到時(shí)間域就可得到所有時(shí)刻的時(shí)間域波場記錄，提取每一時(shí)刻海面測線所有檢波器采集得到的時(shí)間域波場記錄即可得到時(shí)間域地震記錄(圖4b，4c).從VTI常規(guī)9點(diǎn)格式的地震記錄(圖4b)和VTI ADM 9點(diǎn)格式的地震記錄(圖4c)對(duì)比可以清晰看出VTI ADM 9點(diǎn)的地震記錄幾乎沒有頻散，能精確模擬，而VTI常規(guī)9點(diǎn)的地震記錄則出現(xiàn)了嚴(yán)重的數(shù)值頻散，這是由于兩種方法數(shù)值模擬精度差異較大而導(dǎo)致VTI常規(guī)9點(diǎn)格式在大網(wǎng)格間距下數(shù)值頻散誤差較大，而VTI ADM 9點(diǎn)格式仍能保持高精度.具體而言由于VTI常規(guī)9點(diǎn)格式達(dá)到無頻散模擬所需的每個(gè)波長網(wǎng)格點(diǎn)數(shù)為12個(gè)，由于最大頻率約為40 Hz且最小速度(即海水速度)為1492 m·s-1，則無頻散模擬最大允許網(wǎng)格間距大概為1492/40/12=3.1 m，而VTI ADM 9點(diǎn)格式達(dá)到無頻散模擬所需的每個(gè)波長網(wǎng)格點(diǎn)數(shù)僅為3.57個(gè)，則無頻散模擬最大允許網(wǎng)格間距可達(dá)到1492/40/3.57=10.4 m，因此對(duì)于本例橫縱向網(wǎng)格間距分別為10 m和4 m的情況來說，橫向網(wǎng)格間距10 m已經(jīng)明顯超過了VTI常規(guī)9點(diǎn)格式精確模擬所允許的3.1 m間距，因此產(chǎn)生了比較嚴(yán)重的數(shù)值頻散，而VTI ADM 9格式仍能精確模擬.另外從地震記錄上來看頻率空間域PML吸收邊界條件的引入很好地消除了人工邊界反射(圖4b，4c).通過與圖4d中VTI介質(zhì)高精度時(shí)間域12階高階有限差分方法(Duveneck and Bakker，2011)的結(jié)果進(jìn)行對(duì)比，驗(yàn)證了本文方法與VTI介質(zhì)時(shí)間域高階方法精度相當(dāng).

對(duì)BP2007 2D VTI海洋標(biāo)準(zhǔn)模型相同計(jì)算區(qū)域同等計(jì)算精度下做了一個(gè)計(jì)算效率測試(本文同等計(jì)算精度即以恰恰不產(chǎn)生數(shù)值頻散為限)，計(jì)算平臺(tái)為聯(lián)想ThinkCentre M8500t臺(tái)式機(jī)電腦(酷睿i7八核，3.4 GHz).測試結(jié)果為VTI常規(guī)9點(diǎn)格式頻率空間域單炮模擬耗時(shí)10262 s，而VTI ADM 9點(diǎn)格式頻率空間域單炮模擬耗時(shí)僅為2730 s.因此在相同的精度要求下，VTI ADM 9點(diǎn)格式的計(jì)算效率要明顯高于VTI常規(guī)9點(diǎn)格式，這是由于VTI ADM 9點(diǎn)格式相比VTI常規(guī)9點(diǎn)格式計(jì)算精度提升明顯，因此可以通過采用更大的網(wǎng)格間距來提高計(jì)算效率.

6 結(jié)論

VTI介質(zhì)頻率空間域常規(guī)9點(diǎn)格式精度低，本文提出了基于平均導(dǎo)數(shù)優(yōu)化方法的二維VTI介質(zhì)qP波波動(dòng)方程頻率空間域二階9點(diǎn)格式，通過優(yōu)化方法使數(shù)值頻散達(dá)到極小化，每個(gè)波長所需要的網(wǎng)格點(diǎn)數(shù)在1%的誤差范圍內(nèi)僅為3.57個(gè)網(wǎng)格點(diǎn)數(shù)，而VTI常規(guī)9點(diǎn)差分格式在相同的誤差范圍內(nèi)則需要約12個(gè)網(wǎng)格點(diǎn)數(shù)，因此VTI常規(guī)9點(diǎn)格式波場模擬精度差，相較本文方法更容易出現(xiàn)數(shù)值頻散誤差，本文方法體現(xiàn)出在克服數(shù)值頻散誤差方面的優(yōu)越性.數(shù)值模擬結(jié)果表明本文的基于平均導(dǎo)數(shù)優(yōu)化方法的二維VTI介質(zhì)頻率空間域正演方法不僅具有很好的計(jì)算精度和計(jì)算效率，而且平均導(dǎo)數(shù)優(yōu)化方法的特性也使其同時(shí)具有很好的適用性和靈活性，能適用于不同的縱橫向采樣間隔情況.從新方法的應(yīng)用前景上來看，由于各向異性頻率空間域正演是各向異性頻率空間域全波形反演的重要基礎(chǔ)和關(guān)鍵環(huán)節(jié)，本文研究的基于平均導(dǎo)數(shù)優(yōu)化方法的高精度VTI介質(zhì)頻率空間域正演方法可進(jìn)一步應(yīng)用于VTI介質(zhì)頻率空間域全波形反演的研究中，具有良好的應(yīng)用前景.

致謝 衷心感謝兩位外審專家提出的很好的修改意見和本文編輯細(xì)致而負(fù)責(zé)的工作，從而使得本文更加完善.

附錄A 構(gòu)建VTI ADM 9點(diǎn)格式PML波動(dòng)方程

引入PML技術(shù)(Berenger，1994；Moreira et al．，2014)，將頻率域PML二維拉伸函數(shù)ex和ez與二維VTI介質(zhì)頻率空間域qP波波動(dòng)方程(式(4))相結(jié)合推導(dǎo)得到二維VTI介質(zhì)頻率空間域qP波PML波動(dòng)方程：

(A1)

其中：

其中fpeak是震源主頻，xi和zi分別為PML區(qū)域內(nèi)部點(diǎn)與邊界的x方向和z方向的距離，npml為PML層厚度，a0取經(jīng)驗(yàn)值1.79 (吳國忱等，2007).

式(A1)結(jié)合二階VTIADM9點(diǎn)有限差分格式(式(8))推導(dǎo)得到VTIADM9點(diǎn)格式的PML波動(dòng)方程如下：

αPm-1,n-1+βzPm,n-1+αPm+1,n-1+βxPm-1,n+γPm,n+βxPm+1,n+αPm-1,n+1+βzPm,n+1+αPm+1,n+1=0，

(A2)

其中：

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

Alkhalifah T. 1998. Acoustic approximations for processing in transversely isotropic media.Geophysics, 63(2): 623-631, doi: 10.1190/1.1444361.

Alkhalifah T. 2000. An acoustic wave equation for anisotropic media.Geophysics, 65(4): 1239-1250, doi: 10.1190/1.1444815. Alkhalifah T. 2014. Full Waveform Inversion in an Anisotropic World: Where are the Parameters Hiding? Houten, The Netherlands: EAGE Publications bv.

Berenger J P. 1994. A perfectly matched layer for the absorption of electromagnetic waves.JournalofComputationalPhysics, 114(2): 185-200, doi: 10.1006/jcph.1994.1159. Chen J B. 2012. An average-derivative optimal scheme for frequency-domain scalar wave equation.Geophysics, 77(6): T201-T210, doi: 10.1190/GEO2011-0389.1. Chen J B. 2014. Dispersion analysis of an average-derivative optimal scheme for Laplace-domain scalar wave equation.Geophysics, 79(2): T37-T42, doi: 10.1190/GEO2013-0230.1.

Cheng J B, Kang W, Wang T F. 2013. Description of qP-wave propagation in anisotropic media, Part 1: Pseudo-pure-mode wave equations.ChineseJ.Geophys. (in Chinese), 56(10): 3474-3486, doi: 10.6038/cjg20131022.

Du X D, Liu J R, Qi Y P, et al. 2009. Finite difference high-order weighted-averaging operator for frequency-space domain qP wave in 3-D VTI media.ProgressinGeophys. (in Chinese), 24(1): 211-222.

Duveneck E, Bakker P M. 2011. Stable P-wave modeling for reverse-time migration in tilted TI media.Geophysics, 76(2): S65-S75, doi: 10.1190/1. 3533964.

Grechka V, Zhang L B, Rector J W. 2004. Shear waves in acoustic anisotropic media.Geophysics, 69(2): 576-582, doi: 10.1190/1.1707077.

Grini M, Ribodetti A, Virieux J, et al. 2007. Finite-difference frequency-domain modeling of acoustic wave propagation in 2D TTI anisotropic media. ∥69th Annual EAGE Conference and Exhibition, EAGE Extended Abstracts, 323.

Li G H, Feng J G, Zhu G M. 2011. Quasi-P wave forward modeling in viscoelastic VTI media in frequency-space domain.ChineseJ.Geophys. (in Chinese), 54(1): 200-207, doi: 10.3969/j.issn.0001-5733.2011.01.021.

Lines L R, Treitel S. 1984. A review of least-squares inversion and its application to geophysical problems.GeophysicalProspecting, 32(2): 159-186, doi: 10.1111/j.1365-2478.1984.tb00726.x. Min D J, Shin C, Kwon B D, et al. 2000. Improved frequency-domain elastic wave modeling using weighted-averaging difference operators.Geophysics, 65(3): 884-895, doi: 10.1190/1.1444785.Moreira R M, Santos M A, Martins J L, et al. 2014. Frequency-domain acoustic-wave modeling with hybrid absorbing boundary conditions.Geophysics, 79(5): A39-A44, doi: 10.1190/GEO2014-0085.1. Operto S, Ribodetti A, Grini M, et al. 2007. Mixed-grid finite-difference frequency-domain viscoacoustic modeling in 2D TTI anisotropic media. ∥ 77th Annual International Meeting, SEG Expanded Abstracts, 2099-2103.

Operto S, Virieux J, Ribodetti A, et al. 2009. Finite-difference frequency-domain modeling of viscoacoustic wave propagation in 2D tilted transversely isotropic (TTI) media.Geophysics, 74(5): T75-T95, doi: 10.1190/1.3157243.

Operto S, Brossier R, Combe L, et al. 2014. Computationally efficient three-dimensional acoustic finite-difference frequency-domain seismic modeling in vertical transversely isotropic media with sparse direct solver.Geophysics, 79(5): T257-T275, doi: 10.1190/GEO2013-0478.1.

Shin C S, Sohn H J. 1998. A frequency-space 2-D scalar wave extrapolator using extended 25-point finite-difference operator.Geophysics, 63(1): 289-296, doi: 10.1190/1.1444323.

Thomsen L. 1986. Weak elastic anisotropy.Geophysics, 51(10): 1954-1966, doi: 10.1190/1.1442051.

Tsvankin I. 1996. P-wave signatures and notation for transversely isotropic media: An overview.Geophysics, 61(2): 467-483, doi: 10.1190/1.1443974.

Tsvankin I. 2001. Seismic Signatures and Analysis of Reflection Data in Anisotropic Media. United Kingdom: Elsevier Science.

Tsvankin I, Gaiser J, Grechka V, et al. 2010. Seismic anisotropy in exploration and reservoir characterization: An overview.Geophysics, 75(5): 75A15-75A29, doi: 10.1190/1.3481775.

Warner M, Ratclife A, Nangoo T, et al. 2013. Anisotropic 3D full-waveform inversion.Geophysics, 78(2): R59-R80, doi: 10.1190/GEO2012-0338.1. Wu G C, Liang K. 2005. Quasi P-wave forward modeling in frequency-space domain in VTI media.OilGeophysicalProspecting(in Chinese), 40(5): 535-545. Wu G C. 2006. Seismic Wave Propagation and Imaging in Anisotropic Media (in Chinese). Dongying: China University of Petroleum Press. Wu G C, Luo C M, Liang K. 2007. Frequency-space domain finite difference numerical simulation of elastic wave in TTI media.JournalofJilinUniversity(EarthScienceEdition) (in Chinese), 37(5): 1023-1033. Wu G C, Liang K. 2007. High precision finite difference operators for qP wave equation in VTI media.ProgressinGeophys. (in Chinese), 22(3): 896-904.

Zhang H, Liu H, Liu L, et al. 2014. Frequency domain acoustic equation high-order modeling based on an average-derivative method.ChineseJ.Geophys. (in Chinese), 57(5): 1599-1611, doi: 10.6038/cjg20140523. Zhang Y, Wu G C. 2013. Review of prestack reverse-time migration in TTI media.ProgressinGeophys. (in Chinese), 28(1): 409-420, doi: 10.6038/pg20130146.

附中文參考文獻(xiàn)

程玖兵, 康瑋, 王騰飛. 2013. 各向異性介質(zhì)qP波傳播描述1: 偽純模式波動(dòng)方程. 地球物理學(xué)報(bào), 56(10): 3474-3486, doi: 10.6038/cjg20131022.

杜向東, 劉軍榮, 戚艷平等. 2009. 三維VTI介質(zhì)qP波方程頻率空間域有限差分高階加權(quán)算子. 地球物理學(xué)進(jìn)展, 24(1): 211-222.

李桂花, 馮建國, 朱光明. 2011. 黏彈性VTI介質(zhì)頻率空間域準(zhǔn)P波正演模擬. 地球物理學(xué)報(bào), 54(1): 200-207, doi: 10.3969/j.issn.0001-5733.2011.01.021.

吳國忱, 梁鍇. 2005. VTI介質(zhì)頻率-空間域準(zhǔn)P波正演模擬. 石油地球物理勘探, 40(5): 535-545.

吳國忱. 2006. 各向異性介質(zhì)地震波傳播與成像. 東營: 中國石油大學(xué)出版社.

吳國忱, 羅彩明, 梁鍇. 2007. TTI介質(zhì)彈性波頻率-空間域有限差分?jǐn)?shù)值模擬. 吉林大學(xué)學(xué)報(bào)(地球科學(xué)版), 37(5): 1023-1033.

吳國忱, 梁鍇. 2007. VTI介質(zhì)qP波方程高精度有限差分算子. 地球物理學(xué)進(jìn)展, 22(3): 896-904.

張衡, 劉洪, 劉璐等. 2014. 基于平均導(dǎo)數(shù)優(yōu)化方法的聲波方程頻率空間域高階正演. 地球物理學(xué)報(bào), 57(5): 1599-1611, doi: 10.6038/cjg20140523.

張巖, 吳國忱. 2013. TTI介質(zhì)疊前逆時(shí)偏移成像研究綜述. 地球物理學(xué)進(jìn)展, 28(1): 409-420, doi: 10.6038/pg20130146.

(本文編輯 何燕)

Forward modeling of VTI media in the frequency-space domain based on an average-derivative optimal method

ZHANG Heng1，LIU Hong2*，TANG Xiang-De2,3，WANG Yang2,3，ZHANG Bao-Jin1

1MLRKeyLaboratoryofMarineMineralResources,GuangzhouMarineGeologicalSurvey,Guangzhou510075,China2KeyLaboratoryofPetroleumResourcesResearch,InstituteofGeologyandGeophysics,ChineseAcademyofSciences,Beijing100029,China3UniversityofChineseAcademyofSciences,Beijing100049,China

Seismic anisotropy exists extensively in the earth, whereas the TI medium is the most common anisotropic one. The conventional finite difference modeling (FDFDM) of 2D VTI media frequency-space domain with 9-point scheme can cause serious numerical dispersion because of its poor computational accuracy. The average-derivative optimal method (ADM) not only is suitable for equal and unequal directional sampling intervals, but also retains high precision. We introduce this idea into the qP wave equation of 2D VTI media to improve the computational accuracy.We firstly deal with the P-SV wave phase velocity dispersion equation for VTI media and derive the frequency-space domain qP wave equation for 2D VTI media following Alkhalifah′s TI acoustic approximation theory. Then we propose a new kind of qP wave equation for 2D VTI media with a second-order 9-point FDFDM scheme based on ADM. Specifically, we represent the finite-difference approximations of the second-order centered spatial-derivative terms as the weighted average of 3 grid points in orthogonal directions and the acceleration term as the weighted average of all 9 grid points. Afterwards we use the least-square optimal method to resolve the optimized coefficients of the VTI media with the ADM 9-point scheme and perform numerical dispersion analysis towards this new scheme. The results show that the ADM 9-point scheme for VTI media can decrease the required number of grid points per wavelength from 12 to 3.57 bounded by a phase velocity error range of 1%. Therefore this scheme significantly increases the computational accuracy of VTI media FDFDM compared to the VTI media conventional 9-point scheme. So the VTI media ADM 9-point scheme allows using a larger grid interval in the modeling and obviously improves the computational efficiency. Finally we derive the 2D VTI media ADM 9-point frequency-space domain with the perfectly matched layer (PML) wave equation and perform seismic wave modeling in the frequency-space domain.We use a complex BP2007 2D VTI ocean standard model to verify the validity and precision of the VTI media ADM 9-point scheme. As the VTI media conventional 9-point scheme is also suitable for unequal directional sampling intervals, we use this scheme for comparison. The horizontal sampling interval and the vertical sampling interval are 10 m and 4 m, respectively, so the ratio of the horizontal sampling interval to the vertical sampling interval comes to 2.5. The VTI media parameters include P-wave velocity and Thomsen anisotropic parameter. The frequency-space domain wavefield characteristic of qP wave is very clear and depicts well from the 35 Hz monochromatic wavefield computed by the VTI media ADM 9-point scheme. The simulation result with this scheme is accurate while the result with VTI media conventional 9-point scheme exhibits errors due to numerical dispersion. The simulation results also demonstrate the VTI media ADM 9-point scheme is in good agreement with the high-precision time domain 12-order finite difference scheme for VTI media. It is worth mentioning that the frequency-space domain PML absorbing boundary condition eliminates the artificial boundary reflection very well. The numerical example proves the precision and validity of the VTI media ADM 9-point scheme.We propose a new VTI ADM 9-point scheme which obviously increases the computational precision of VTI FDFDM compared to the VTI media conventional 9-point scheme. We obtain the optimized coefficients by the least-square method and decrease the required grid points per wavelength from 12 to 3.57. The numerical example demonstrates that the VTI ADM 9-point scheme can not only possess high computational precision and computational efficiency, but also possess applicability and flexibility. The VTI ADM 9-point scheme can be developed further and has many applications. Taking FWI for example, the VTI ADM 9-point scheme can be applied to the VTI FWI as a fast and accurate modeling engine.

VTI media； Average-derivative optimal method； Frequency-space domain

張衡，劉洪，唐祥德等. 2015. 基于平均導(dǎo)數(shù)優(yōu)化方法的VTI介質(zhì)頻率空間域正演.地球物理學(xué)報(bào)，58(9):3306-3316,

10.6038/cjg20150924.

Zhang H, Liu H, Tang X D, et al. 2015. Forward modeling of VTI media in the frequency-space domain based on an average-derivative optimal method.ChineseJ.Geophys. (in Chinese),58(9):3306-3316,doi:10.6038/cjg20150924.

10.6038/cjg20150924

P631

2014-11-25，2015-04-10收修定稿

國家油氣重大專項(xiàng)(2011ZX05003-003)，國家自然科學(xué)基金(41176056)和國家高技術(shù)研究發(fā)展計(jì)劃(863計(jì)劃)主題項(xiàng)目(2012AA061202)聯(lián)合資助．

張衡，男，1987年生，2014年博士畢業(yè)于中國科學(xué)院地質(zhì)與地球物理研究所，工程師，主要從事地震波模擬、偏移成像和全波形反演研究．E-mail: zhangheng5037@126.com

*通訊作者 劉洪，男，1959年生，中國科學(xué)院地質(zhì)與地球物理研究所研究員，主要從事地震波成像和油儲(chǔ)地球物理研究． E-mail: liuhong@mail.iggcas.ac.cn