蔡曉慧， 劉洋*， 王建民， 王維紅， 任志明

1 中國(guó)石油大學(xué)(北京)油氣資源與探測(cè)國(guó)家重點(diǎn)實(shí)驗(yàn)室， 北京 102249 2 中國(guó)石油大學(xué)(北京)CNPC物探重點(diǎn)實(shí)驗(yàn)室， 北京 102249 3 大慶油田有限責(zé)任公司勘探開發(fā)研究院， 大慶 163318

?

基于自適應(yīng)優(yōu)化有限差分方法的全波VSP逆時(shí)偏移

蔡曉慧1,2， 劉洋1,2*， 王建民3， 王維紅3， 任志明1,2

1 中國(guó)石油大學(xué)(北京)油氣資源與探測(cè)國(guó)家重點(diǎn)實(shí)驗(yàn)室， 北京 102249 2 中國(guó)石油大學(xué)(北京)CNPC物探重點(diǎn)實(shí)驗(yàn)室， 北京 102249 3 大慶油田有限責(zé)任公司勘探開發(fā)研究院， 大慶 163318

與地面地震資料相比，VSP資料具有分辨率高、環(huán)境噪聲小及能更好地反映井旁信息等優(yōu)點(diǎn).常規(guī)VSP偏移主要對(duì)上行反射波進(jìn)行成像，存在照明度低、成像范圍受限等問題.為了增加照明度、拓寬成像范圍、提高成像精度，本文采用直達(dá)波除外的所有聲波波場(chǎng)數(shù)據(jù)(全波)，包括一次反射波、多次反射波等進(jìn)行疊前逆時(shí)偏移成像.針對(duì)逆時(shí)偏移中的四個(gè)關(guān)鍵問題，即波場(chǎng)延拓、吸收邊界條件、成像條件及低頻噪聲的壓制，本文分別采用自適應(yīng)變空間差分算子長(zhǎng)度的優(yōu)化有限差分方法(自適應(yīng)優(yōu)化有限差分方法)求解二維聲波波動(dòng)方程以實(shí)現(xiàn)高精度、高效率的波場(chǎng)延拓，采用混合吸收邊界條件壓制因計(jì)算區(qū)域有限所引起的人工邊界反射，采用震源歸一化零延遲互相關(guān)成像條件進(jìn)行成像，采用拉普拉斯濾波方法壓制逆時(shí)偏移中產(chǎn)生的低頻噪聲.本文對(duì)VSP模型數(shù)據(jù)的逆時(shí)偏移成像進(jìn)行了分析，結(jié)果表明：自適應(yīng)優(yōu)化有限差分方法比傳統(tǒng)有限差分方法具有更高的模擬精度與計(jì)算效率，適用于VSP逆時(shí)偏移成像；全波場(chǎng)VSP逆時(shí)偏移成像比上行波VSP逆時(shí)偏移的成像范圍大、成像效果好；相對(duì)于反褶積成像條件，震源歸一化零延遲互相關(guān)成像條件具有穩(wěn)定性好、計(jì)算效率高等優(yōu)點(diǎn).將本文方法應(yīng)用于某實(shí)際VSP資料的逆時(shí)偏移成像，進(jìn)一步驗(yàn)證了本文方法的正確性和有效性.

VSP； 自適應(yīng)優(yōu)化有限差分方法； 全波； 逆時(shí)偏移

1 引言

近年來，隨著VSP技術(shù)的發(fā)展以及地震勘探對(duì)象的日趨復(fù)雜，可適應(yīng)起伏地表、高陡構(gòu)造、復(fù)雜速度區(qū)域的VSP成像在實(shí)際工業(yè)中越來越受重視.因?yàn)榈孛娴卣鹳Y料處理通常不能提供全面、可靠的目標(biāo)層巖丘下方的構(gòu)造信息，而VSP技術(shù)在查清井旁構(gòu)造，探明巖丘下方構(gòu)造，提取縱、橫波信息等方面有著重要的作用(Campbell等，2002; Chon等，2003; Hornby等，2006)，因此VSP偏移成像的研究顯得十分必要.VSP成像方法主要包括：Harwijanto等(1987)提出的單炮記錄反演成像、Nicoletis等(1995)的基于射線理論的偏移成像方法、Xiao和Schuster(2009)的基于格林函數(shù)的偏移成像、Chang和McMechan(1986)的基于波動(dòng)方程的逆時(shí)偏移成像方法等.VSP逆時(shí)偏移是一種基于全波波動(dòng)方程的深度域偏移方法，其基本原理是以炮點(diǎn)處震源進(jìn)行正向延拓，再利用記錄的地震波場(chǎng)作為該記錄點(diǎn)的邊界條件進(jìn)行逆時(shí)延拓，并利用成像條件實(shí)現(xiàn)對(duì)地下各點(diǎn)的成像.逆時(shí)偏移能適應(yīng)任意復(fù)雜速度模型，無傾角限制，理論上可以實(shí)現(xiàn)對(duì)所有類型波(反射波、繞射波、回轉(zhuǎn)波以及多次反射波等)的成像，而不需要進(jìn)行上、下行波的分離(Hou and He，1990).20世紀(jì)80年代初期，Whitemore(1983)首次提出逆時(shí)偏移思想后，疊后逆時(shí)偏移得到較快發(fā)展.Baysal等(1983)采用偽譜法，Chang等(1986)采用有限差分方法，Loewenthal和Mufti(1983)在空間-頻率域分別實(shí)現(xiàn)了疊后逆時(shí)偏移.隨后，疊前逆時(shí)偏移開始迅猛發(fā)展.Chang和McMechan(1987，1990，1993，1994)應(yīng)用有限差分方法分別實(shí)現(xiàn)了2D和3D的聲波、彈性波及各向異性逆時(shí)偏移；Causse和Ursin(2000)證明了基于黏彈性介質(zhì)逆時(shí)偏移的有效性；Wu等(1996)實(shí)現(xiàn)了3D高階有限差分逆時(shí)偏移.21世紀(jì)以來，隨著地震勘探對(duì)象的復(fù)雜化與并行計(jì)算機(jī)的發(fā)展，再一次掀起了逆時(shí)偏移研究與應(yīng)用的熱潮.Xu等(2011)實(shí)現(xiàn)了角道集逆時(shí)偏移；Yan和Liu(2013a)發(fā)展了時(shí)空域有限差分的聲波逆時(shí)偏移；Yan和Liu(2013b)、Zhao等(2014)實(shí)現(xiàn)了黏滯聲波逆時(shí)偏移成像；Li等(2014)將最小二乘逆時(shí)偏移應(yīng)用于黏彈介質(zhì).除此之外，不少學(xué)者對(duì)逆時(shí)偏移中的難點(diǎn)進(jìn)行了專門的研究，包括逆時(shí)偏移的噪聲壓制(Yoon et al., 2004; Zhang et al., 2009; Du et al., 2013)、成像方法(Chattopadhyay and Mcmechan，2008)、保幅性(Zhang et al., 2007)、計(jì)算量和存儲(chǔ)問題(Li et al., 2010; Liu et al., 2010a; Liu et al., 2010b; Wang et al., 2012; Hu et al., 2013)等.針對(duì)VSP資料的逆時(shí)偏移，Chang等(1986)率先將逆時(shí)偏移方法應(yīng)用到VSP資料；Zhu等(1989)提出了消除內(nèi)部界面反射波和層間多次波干擾的雙程無反射波動(dòng)方程逆時(shí)偏移；Ketil等(1998)和Hokstad等(1988)實(shí)現(xiàn)了基于Walkaway VSP資料的彈性波逆時(shí)偏移；Xiao和Leaney(2010)實(shí)現(xiàn)了基于VSP資料的透射波逆時(shí)偏移；Sun和McMechan(2010)提出了一種基于VSP資料的非線性彈性逆時(shí)反演；Sun和Sun(2010)應(yīng)用了偽譜法實(shí)現(xiàn)VSP逆時(shí)偏移；Sun等(2011)實(shí)現(xiàn)了相對(duì)保幅的角度域VSP逆時(shí)偏移.

VSP逆時(shí)偏移中四個(gè)關(guān)鍵問題為波場(chǎng)延拓、吸收邊界條件、成像條件以及低頻噪聲的壓制.波場(chǎng)延拓的實(shí)質(zhì)是求解波動(dòng)方程，由于有限差分方法具有計(jì)算效率高、精度較高、穩(wěn)定性較好的優(yōu)點(diǎn)，因而被廣泛應(yīng)用于求解波動(dòng)方程.目前主要有兩種基于頻散關(guān)系計(jì)算空間導(dǎo)數(shù)有限差分系數(shù)的方法，一種是基于泰勒級(jí)數(shù)展開的方法，另一種是基于優(yōu)化的方法.這兩種方法都有一定的優(yōu)越性和局限性，泰勒級(jí)數(shù)展開法求解有限差分系數(shù)的計(jì)算效率高，但是只能在較小的波數(shù)或者頻率范圍內(nèi)達(dá)到較高的精度.對(duì)于優(yōu)化的方法，局部?jī)?yōu)化方法得到的差分系數(shù)難以達(dá)到全局優(yōu)化，而全局優(yōu)化方法計(jì)算效率一般較低，且不一定能得到全局最優(yōu)解.Liu(2013)提出了一種基于最小二乘的全局優(yōu)化有限差分方法，這是一種線性求解有限差分系數(shù)的方法，計(jì)算效率高且能在給定的頻帶范圍內(nèi)達(dá)到較高的精度.本文分析對(duì)比了該優(yōu)化有限差分方法與泰勒級(jí)數(shù)展開有限差分方法的頻散、模擬精度以及VSP逆時(shí)偏移結(jié)果.模型試算結(jié)果表明：在差分算子長(zhǎng)度相同的情況下，優(yōu)化有限差分法比泰勒級(jí)數(shù)展開有限差分法的模擬精度高，優(yōu)化有限差分逆時(shí)偏移成像結(jié)果更準(zhǔn)確；在相同模擬精度的前提下，優(yōu)化有限差分法比泰勒級(jí)數(shù)展開有限差分法的空間差分算子長(zhǎng)度短，即波場(chǎng)延拓的效率高.因此，我們采用優(yōu)化有限差分法來提高波場(chǎng)延拓的計(jì)算效率.為進(jìn)一步節(jié)省計(jì)算時(shí)間，引入自適應(yīng)變空間算子長(zhǎng)度方案(Liu and Sen，2011)，通過在高速區(qū)域采用較短的空間差分算子來提高計(jì)算效率.模型試算結(jié)果驗(yàn)證了自適應(yīng)變空間算子長(zhǎng)度方案與優(yōu)化有限差分方法結(jié)合(自適應(yīng)優(yōu)化有限差分方法)的優(yōu)越性.本文采用混合吸收邊界條件(Liu and Sen，2010)壓制因計(jì)算區(qū)域有限所引起的人工邊界反射，對(duì)比分析震源歸一化零延遲互相關(guān)成像條件與反褶積成像條件(Alejandro and Biondo，2002)，采用拉普拉斯濾波(Zhang and Sun，2009)壓制逆時(shí)偏移成像產(chǎn)生的低頻噪聲.

目前，VSP逆時(shí)偏移主要以上行反射波作為有效波進(jìn)行逆時(shí)延拓，但僅用上行反射波做逆時(shí)延拓會(huì)導(dǎo)致偏移成像結(jié)果照明度低、成像范圍窄(Fleury，2013; Soni and Verschuur，2014)；并且當(dāng)上行反射波作為有效波時(shí)，需要先進(jìn)行波場(chǎng)分離，在波場(chǎng)分離中一般會(huì)對(duì)原始波場(chǎng)造成一定程度的損害，損失部分有效信息.本文以直達(dá)波除外的所有聲波波場(chǎng)(全波)作為有效波，即一次多次波、多次反射波等都被視為逆時(shí)偏移的有效波進(jìn)行逆時(shí)延拓.為驗(yàn)證全波VSP逆時(shí)偏移的有效性，本文采用上行波、下行波、全波分別作為逆時(shí)延拓的有效波，對(duì)層狀模型、斷層模型和Marmousi模型的VSP數(shù)據(jù)進(jìn)行逆時(shí)偏移.模型試算的結(jié)果表明：以全波作為有效波的成像范圍大、成像效果好.為驗(yàn)證自適應(yīng)優(yōu)化有限差分方法的效率與精度，本文將該方法與自適應(yīng)泰勒級(jí)數(shù)展開有限差分方法和固定算子長(zhǎng)度的泰勒級(jí)數(shù)展開有限差分法進(jìn)行了對(duì)比分析.分析結(jié)果表明：與固定算子長(zhǎng)度的泰勒級(jí)數(shù)展開有限差分法相比，自適應(yīng)優(yōu)化有限差分方法可以在不降低精度的前提下，較大地提高計(jì)算效率.最后，采用自適應(yīng)優(yōu)化有限差分方法、混合吸收邊界條件、震源歸一化零延遲互相關(guān)成像條件以及拉普拉斯濾波去噪方法實(shí)現(xiàn)了某地區(qū)實(shí)際VSP資料的逆時(shí)偏移，進(jìn)一步驗(yàn)證了本文方法的有效性.

2 方法原理

2.1 優(yōu)化有限差分方法

二維聲波波動(dòng)方程為

(1)

其中，p為標(biāo)量波場(chǎng)，v為速度，t為時(shí)間，x、z為空間坐標(biāo).

VSP逆時(shí)偏移的震源正向延拓和地震記錄逆時(shí)延拓的表達(dá)式分別為

(2)

(3)

其中，pF為震源延拓波場(chǎng)，pB為逆時(shí)延拓波場(chǎng)，f(t)為震源，δ為脈沖函數(shù)，Q(xr,zr;t)為地震記錄，(xr,zr)為檢波點(diǎn)坐標(biāo)，(xs,zs)為炮點(diǎn)坐標(biāo).

分別采用二階中心差分計(jì)算二階時(shí)間導(dǎo)數(shù)和2M階精度差分計(jì)算二階空間導(dǎo)數(shù)，公式(2)和(3)變形為

(4)

(5)

(6)

(7)

(n=1,2,…,M).

(8)

二維聲波波動(dòng)方程數(shù)值模擬的相速度頻散δ表示為

(9)

(10)

當(dāng)δmax越趨向于1，頻散越小.

圖1是優(yōu)化有限差分方法(optimal-basedFDmethod)、空間域泰勒級(jí)數(shù)展開有限差分方法(spaceTE-basedFDmethod)和時(shí)空域泰勒級(jí)數(shù)展開有限差分方法(time-spaceTE-basedFDmethod)的頻散曲線對(duì)比圖.由圖可見：優(yōu)化有限差分方法在區(qū)間[0，π]內(nèi)，總體比空間域泰勒級(jí)數(shù)展開有限差分方法和時(shí)空域泰勒級(jí)數(shù)展開有限差分方法的數(shù)值頻散都小，并且空間域泰勒級(jí)數(shù)展開有限差分方法和時(shí)空域泰勒級(jí)數(shù)展開有限差分方法的頻散相近.

2.2 自適應(yīng)變空間差分算子長(zhǎng)度方案

在求解波動(dòng)方程時(shí)，通常采用固定的算子長(zhǎng)度計(jì)算空間導(dǎo)數(shù)，而自適應(yīng)變空間差分算子長(zhǎng)度方案采用變化的算子長(zhǎng)度計(jì)算空間導(dǎo)數(shù).自適應(yīng)變差分算子長(zhǎng)度方案是基于算法的精度控制不同速度網(wǎng)格點(diǎn)的M值，其表達(dá)式為(Liu and Sen，2011)

(11)

(12)

其中，ε為單位網(wǎng)格內(nèi)傳播時(shí)間的相對(duì)誤差，fmax為子波最大頻率，η為最大誤差范圍.利用公式(12)可以自動(dòng)確定不同速度網(wǎng)格點(diǎn)的M值.一般在速度小的網(wǎng)格點(diǎn)用長(zhǎng)空間差分算子，在速度大的網(wǎng)格點(diǎn)用短空間差分算子.這種方法可以在不降低模擬精度的情況下，提高計(jì)算效率.

圖1 優(yōu)化有限差分、空間域泰勒級(jí)數(shù)展開有限差分和時(shí)空域泰勒級(jí)數(shù)展開有限差分方法的最大頻散值δmax頻散曲線(b=2.74，r=0.15)(a) M=4; (b) M=14.Fig.1 Maximum dispersion value δmax curves for the optimal-based FD method, the space TE-based FD method and time-space TE-based FD method (b=2.74， r=0.15)

為驗(yàn)證自適應(yīng)變空間差分算子長(zhǎng)度方案的有效性，設(shè)計(jì)一個(gè)速度模型，其速度的變化范圍為1500～5000 m·s-1，h=10 m，τ=1 ms，fmax=40 Hz，同時(shí)對(duì)比自適應(yīng)優(yōu)化有限差分方法(adaptive optimal-based FD method)、自適應(yīng)時(shí)空域泰勒級(jí)數(shù)展開有限差分方法(adaptive time-space TE-based FD method) (后文中的時(shí)空域泰勒級(jí)數(shù)展開有限差分方法(time-space TE-based FD method)簡(jiǎn)寫為泰勒級(jí)數(shù)展開有限差分法(TE-based FD method))在相應(yīng)的速度網(wǎng)格點(diǎn)的M值，最大誤差η分別為10-4、10-5.圖2為M隨速度的變化圖，由圖可知：

(1) 對(duì)于相同的η，速度越小，M值越大；

(2) 對(duì)于優(yōu)化有限差分方法，η越小，則精度越高，M越大；

(3) 當(dāng)η=10-5時(shí)，優(yōu)化有限差分方法的M比泰勒級(jí)數(shù)展開有限差分法的更小，即計(jì)算效率更高.

圖2 有限差分算子長(zhǎng)度參數(shù)M隨速度變化圖Fig.2 Variation of FD operator length parameter M with velocity

2.3 混合吸收邊界條件

在VSP逆時(shí)偏移中，一個(gè)不可避免的問題是如何有效壓制由計(jì)算區(qū)域邊界引起的邊界反射能量，本文采用混合吸收邊界條件進(jìn)行壓制.混合吸收邊界條件的基本原理是把計(jì)算區(qū)域分為內(nèi)部區(qū)、過渡區(qū)、邊界區(qū).第一步，采用雙程波波動(dòng)方程計(jì)算內(nèi)部區(qū)和過渡區(qū)波場(chǎng)值；第二步，采用單程波波動(dòng)方程計(jì)算邊界區(qū)與過渡區(qū)波場(chǎng)值；第三步，采用雙程波波場(chǎng)值與單程波波場(chǎng)值的線性加權(quán)得到過渡區(qū)最終波場(chǎng)值.過渡區(qū)起到了對(duì)波場(chǎng)平滑過渡的作用，避免由內(nèi)部區(qū)的雙程波波動(dòng)方程到邊界區(qū)的單程波波動(dòng)方程的急劇變化而導(dǎo)致的邊界干擾反射無法得到有效壓制的問題.以模型的上邊界和左上角為例，其單程波波動(dòng)方程表達(dá)式分別為(Clayton and Engquist，1997；Liu and Sen，2010)

(13a)

(13b)

當(dāng)過渡區(qū)的網(wǎng)格厚度為1時(shí)，混合吸收邊界條件等效于Clayton-Engquist吸收邊界條件；過渡區(qū)網(wǎng)格厚度為10時(shí)能達(dá)到較好的吸收效果.以重新采樣后的Marmousi速度模型(圖3)為例，速度模型大小為5000m×3500m，震源是位于(2500m，0m)的主頻為30Hz的Ricker子波，h=10 m，τ=1 ms，記錄時(shí)間為4 s.圖4為Marmousi速度模型的無邊界條件與邊界網(wǎng)格點(diǎn)數(shù)為10的混合吸收邊界條件模擬地震記錄，由圖可見混合吸收邊界條件能有效地壓制邊界反射.

圖3 Marmousi速度模型Fig.3 Marmousi velocity model

2.4 震源歸一化零延遲互相關(guān)成像條件與反褶積成像條件的對(duì)比分析

震源歸一化零延遲互相關(guān)成像條件的表達(dá)式為(Chattopadhyay and Mcmechan，2008)

(14)

其中，S(x,z,t)表示震源波場(chǎng)，R(x,z,t)表示檢波點(diǎn)波場(chǎng)，Tmax為地震記錄長(zhǎng)度，I(x,z)表示反射系數(shù)，μ是穩(wěn)定因子.

反褶積成像條件的表達(dá)式為(Alejandroetal., 2002)

(15)

(16)

為對(duì)比震源歸一化零延遲互相關(guān)成像條件與反褶積成像條件的成像效果，設(shè)計(jì)一個(gè)雙層層狀模型，模型大小為1000m×1000m，h=10 m，τ=1 ms，震源為主頻為30 Hz的Ricker子波，混合吸收邊界過渡區(qū)網(wǎng)格點(diǎn)數(shù)為10，自適應(yīng)變空間算子長(zhǎng)度方案參數(shù)為fmax=75 Hz、η=10-4，記錄長(zhǎng)度為1 s.觀測(cè)系統(tǒng)參數(shù)為：地面放炮，炮點(diǎn)x坐標(biāo)范圍為0～990 m，炮間距為100 m，共10炮；5個(gè)檢波點(diǎn)位于x=500 m，z=610～650 m處，檢波點(diǎn)間距為10 m.圖5為該層狀模型的VSP逆時(shí)偏移成像結(jié)果.圖5 a、d的對(duì)比表明：反褶積成像條件比互相關(guān)成像條件分辨率高.但是從圖5b—d可以看出：穩(wěn)定因子對(duì)反褶積成像的影響較大，如果μ取值不當(dāng)會(huì)造成嚴(yán)重成像噪聲.因此，本文采用震源歸一化零延遲互相關(guān)成像條件實(shí)現(xiàn)VSP逆時(shí)偏移.

圖4 Marmousi速度模型地震記錄(a) 無邊界條件; (b) 邊界網(wǎng)格數(shù)為10的混合吸收邊界條件.Fig.4 Seismic records for Marmousi velocity model(a) Non-ABC; (b) Hybrid ABC with the width of 10.

圖5 VSP逆時(shí)偏移結(jié)果(a) 震源歸一化零延遲互相關(guān)成像條件; (b) 反褶積成像條件，μ=0.01; (c) 反褶積成像條件，μ=0.6; (d) 反褶積成像條件，μ=1.Fig.5 Reverse VSP RTM results(a) Source-normalized cross-correlation imaging conditions; (b) Deconvolution imaging condition with μ=0.01; (c) Deconvolution imaging condition with μ=0.6; (d) Deconvolution imaging condition with μ=1.

圖6 (a)傾斜界面模型; (b) 泰勒級(jí)數(shù)展開有限差分法，M=8; (c)泰勒級(jí)數(shù)展開有限差分法，M=32; (d)優(yōu)化有限方法，M=8Fig.6 (a) Dipping reflector model; (b) TE-based FD method, M=8; (c) TE-based FD method, M=32; (d) Optimal-based FD method, M=8

3 模型算例

3.1 數(shù)值模擬

為了驗(yàn)證自適應(yīng)優(yōu)化有限差分方法的可行性，設(shè)計(jì)一個(gè)傾斜界面模型如圖6 a所示，模型大小為4000 m×4000 m，h=20 m，τ=1 ms，數(shù)值模擬的震源采用主頻為20 Hz的Ricker子波，且位于(2000 m，0 m)，混合吸收邊界過渡區(qū)網(wǎng)格點(diǎn)數(shù)為10.分別采用泰勒級(jí)數(shù)展開有限差分法與優(yōu)化有限差分方法實(shí)現(xiàn)波場(chǎng)模擬.圖6b—d為1.2 s時(shí)刻數(shù)值模擬的波場(chǎng)快照，可見圖6b的頻散比較大，圖6c與d頻散相近，其結(jié)果表明：

(1) 對(duì)于相同算子長(zhǎng)度(M=8)，優(yōu)化有限差分方法(圖6d)模擬精度高于泰勒級(jí)數(shù)展開有限差分方法 (圖6b)；

(2)M=8時(shí)的優(yōu)化有限差分方法(圖6d)與M=32時(shí)的泰勒級(jí)數(shù)展開有限差分方法(圖6c)模擬精度相近.

傾斜界面模型試算結(jié)果表明：在相同算子長(zhǎng)度情況下，優(yōu)化有限差分方法模擬精度高于泰勒級(jí)數(shù)展開有限差分方法模擬精度，表明了優(yōu)化有限差分方法適用于求解波動(dòng)方程空間差分系數(shù)，同時(shí)試算結(jié)果也表明混合吸收邊界條件能有效地壓制邊界反射.

3.2 VSP逆時(shí)偏移

3.2.1 VSP全波逆時(shí)偏移分析

設(shè)計(jì)一個(gè)層狀模型，其速度模型如圖7所示，模型大小為1000 m×1200 m，h=10 m，τ=1 ms，記錄時(shí)間為2 s，地表為自由邊界，震源是主頻為30 Hz 的Ricker子波.觀測(cè)系統(tǒng)為：炮點(diǎn)位于(500 m，0 m)處，檢波點(diǎn)位于(0 m，400 m)處.圖8為檢波點(diǎn)處地震記錄，地震記錄中可見直達(dá)波、一次反射波、一階多次波、二階多次波以及高階多次波.分別對(duì)這5種波進(jìn)行逆時(shí)偏移成像，偏移結(jié)果如圖9所示，其中圖9中深度z=650 m處的黑線為界面位置.圖9表明：直達(dá)波只會(huì)造成炮點(diǎn)到檢波點(diǎn)路徑上的噪聲(圖9 a)；一次反射波在R1點(diǎn)處成像(圖9 b)；一階多次波在R2點(diǎn)處成像，與R1點(diǎn)相比，R2點(diǎn)距離井的位置更遠(yuǎn)(圖9 c)，所以多次波可以拓寬成像區(qū)域；二階多次波可以貢獻(xiàn)兩個(gè)成像點(diǎn)(圖9 d)，一個(gè)比R1點(diǎn)距井源更近，一個(gè)比R2點(diǎn)距井源更加遠(yuǎn)，說明了多次波可以增加照明度及拓寬成像范圍；高階多次波的成像范圍比R1點(diǎn)到R2點(diǎn)相距范圍更廣(圖9e中箭頭所示).以上5種波中，直達(dá)波和一階多次波為下行波；一次反射波和二階多次波為上行波；高階多次波中既有上行也有下行波，如圖8所示，其能量比一次反射波能量弱.若只用上行反射波作為逆時(shí)偏移的有效波場(chǎng)則會(huì)損失部分有效信息.理論分析表明：除直達(dá)波以外的所有波場(chǎng)對(duì)界面的成像均有貢獻(xiàn)，所以多次波可以作為有效波進(jìn)行逆時(shí)延拓，以增加照明度、拓寬偏移成像范圍.

3.2.2 模型試算

首先，我們針對(duì)層狀介質(zhì)模型分別實(shí)現(xiàn)自適應(yīng)優(yōu)化有限差分方法和自適應(yīng)泰勒級(jí)數(shù)展開有限差分方法的VSP逆時(shí)偏移，以驗(yàn)證自適應(yīng)優(yōu)化有限差分方法的有效性和效率.其次，我們采用自適應(yīng)優(yōu)化有限差分方法完成三個(gè)模型試算，分別用分離后的上行波、分離后的下行波以及分離前的全波實(shí)現(xiàn)VSP逆時(shí)偏移，下文中逆時(shí)偏移采用的波場(chǎng)均不包括直達(dá)波.雖然上行波包括一次反射波以及多次波，但是直達(dá)波除外的下行波都是多次波.我們采用下行波偏移即多次波偏移，偏移效果可以證明多次波偏移的有效性，而上行波偏移結(jié)果與全波偏移結(jié)果的對(duì)比也可以驗(yàn)證全波偏移的優(yōu)越性.

(I) 層狀模型

模型如圖10 a所示，模型大小為4000 m×4000 m，h=20 m，τ=1 ms，震源為主頻為20 Hz的Ricker子波，記錄長(zhǎng)度為2.5 s.混合吸收邊界過渡區(qū)網(wǎng)格點(diǎn)數(shù)為10.自適應(yīng)變空間算子長(zhǎng)度方案參數(shù)為fmax=30 Hz，η=5×10-6.觀測(cè)系統(tǒng)參數(shù)為：地面放炮，炮點(diǎn)x坐標(biāo)范圍為0～3980 m，炮間距為20 m，共200個(gè)炮點(diǎn)；檢波點(diǎn)位于x=2000 m，z=0～1980 m，檢波點(diǎn)間距為20 m，共100個(gè)檢波點(diǎn).圖10b、c分別為

圖7 層狀介質(zhì)速度模型Fig.7 Layered velocity model

圖8 層狀介質(zhì)模型的地震記錄Fig.8 Seismic record for the layered model

圖9 VSP逆時(shí)偏移分別采用 (a) 直達(dá)波; (b) 一次反射波; (c) 一階多次波; (d) 二階多次波; (e) 高階多次波Fig.9 VSP RTM results with (a) Direct wave; (b) Primary wave; (c) First-order multiple; (d) Second-order multiple;(e) Higher-order multiples

自適應(yīng)泰勒級(jí)數(shù)展開有限差分方法和自適應(yīng)優(yōu)化有限差分方法的VSP逆時(shí)偏移成像結(jié)果，從圖10中可以看到自適應(yīng)泰勒級(jí)數(shù)展開有限差分方法的VSP逆時(shí)偏移成像結(jié)果中有連續(xù)的虛假界面出現(xiàn)，其假軸由數(shù)值模擬中的頻散產(chǎn)生；自適應(yīng)優(yōu)化有限差分VSP逆時(shí)偏移成像中構(gòu)造清晰，無假軸出現(xiàn).所以，在算子長(zhǎng)度相同的前提下，自適應(yīng)優(yōu)化有限差分VSP逆時(shí)偏移成像精度更高，以下VSP逆時(shí)偏移均采用自適應(yīng)優(yōu)化有限差分方法.

我們?cè)O(shè)計(jì)另一個(gè)層狀模型(圖11)分析多次波在VSP逆時(shí)偏移中的影響.模型大小為2000 m×3000 m，h=10 m，τ=1 ms，記錄總時(shí)長(zhǎng)為2.5 s，震源是主頻為30 Hz的Ricker子波，混合吸收邊界過渡區(qū)網(wǎng)格數(shù)為10.自適應(yīng)變空間算子方案參數(shù)為

圖10 (a) 層狀速度模型; (b) 自適應(yīng)泰勒級(jí)數(shù)展開有限差分法，M=2～3; (c) 自適應(yīng)優(yōu)化有限差分方法，M=2～3.Fig.10 (a) Layered velocity model; (b) Adaptive TE-based FD method FD method, M=2～3; (c) Adaptive optimal-based FD method, M=2～3.

圖11 層狀模型Fig.11 Layered velocity model

fmax=40 Hz，η=10-5.觀測(cè)系統(tǒng)為地面放炮，炮點(diǎn)x坐標(biāo)范圍為0～2000 m，炮間距為50 m，共41個(gè)炮點(diǎn).檢波點(diǎn)為x=1950 m,z=0～3000 m，檢波點(diǎn)間距為10 m，共301個(gè)檢波點(diǎn).圖12a—c分別為x=1950 m、z=0 m的炮點(diǎn)波場(chǎng)分離后上行波記錄、下行波記錄和切除直達(dá)波的全波記錄.圖13為逆時(shí)偏移成像結(jié)果，可得：

(1) 上行波作為有效波時(shí)(圖13 a)，成像位置準(zhǔn)確，井旁界面成像清晰；

(2) 下行波作為有效波時(shí)(圖13 b)，遠(yuǎn)井源距的界面成像比上行波偏移成像結(jié)果更清晰，而井旁界面成像不如上行波偏移好，成像結(jié)果有一定噪聲，但其噪聲能量較弱；

(3) 全波作為有效波(圖13c)時(shí)，井旁界面和遠(yuǎn)井源距界面均準(zhǔn)確成像.由于下行波成像的噪聲能量與上行波的成像能量差異大，在全波成像結(jié)果中幾乎看不到成像噪聲.

此層狀模型試算結(jié)果表明，全波作為有效波時(shí)成像效果比僅用上行波或者僅用下行波作為有效波的成像效果更好.

圖12 層狀模型的VSP地震記錄(a) 上行波場(chǎng); (b) 下行波場(chǎng); (c) 全波波場(chǎng).Fig.12 Seismic records for layered model(a) Up-going wavefield; (b) Down-going wavefield; (c) Full-wavefield.

圖13 VSP逆時(shí)偏移結(jié)果(a)上行波作為有效波; (b) 下行波作為有效波; (c) 全波波場(chǎng)作為有效波.Fig.13 VSP RTM results(a) Up-going wavefield as input; (b) Down-going wavefield as input; (c) Full-wavefield as input.

(II) 斷層模型

斷層速度模型(圖14)大小為2000 m×2000 m，h=10 m，τ=1 ms，記錄總時(shí)間為2 s，震源是主頻為30 Hz的Ricker子波，混合吸收邊界過渡區(qū)網(wǎng)格數(shù)為10，自適應(yīng)變空間算子長(zhǎng)度方案參數(shù)為fmax=50 Hz，η=10-5.觀測(cè)系統(tǒng)參數(shù)為：地面放炮，炮點(diǎn)x坐標(biāo)范圍為0～2000 m，炮間距為100 m，共21個(gè)炮；檢波點(diǎn)位于x=2000 m，檢波點(diǎn)深度范圍為0～2000 m，檢波點(diǎn)間距為10 m，共201個(gè)檢波點(diǎn).圖15為VSP逆時(shí)偏移成像結(jié)果，可見：

(1) 上行波、下行波、全波波場(chǎng)分別作為有效波時(shí)，VSP逆時(shí)偏移成像位置都準(zhǔn)確；

(2) 上行波作為有效波時(shí)(圖15a)，井旁界面成像較清晰；

(3) 下行波作為有效波時(shí)(圖15b)，噪聲比上行波偏移成像噪聲嚴(yán)重；

圖14 斷層速度模型Fig.14 Fault velocity model

圖15 斷層模型VSP逆時(shí)偏移結(jié)果(a) 上行波作為有效波; (b) 下行波作為有效波; (c)全波波場(chǎng)作為有效波.Fig.15 VSP RTM result for fault model(a) Up-going wavefield as input; (b) Down-going wavefield as input; (c) Full-wavefield as input.

(4) 全波VSP逆時(shí)偏移成像結(jié)果(圖15c)比僅用上行波或者下行波逆時(shí)偏移成像界面形態(tài)更清晰.

對(duì)于此斷層模型，全波VSP逆時(shí)偏移成像的結(jié)果比僅用上行波或下行波成像結(jié)果更理想.

(III) Marmousi模型

層狀模型與斷層模型試算表明了全波在簡(jiǎn)單構(gòu)造中VSP逆時(shí)偏移中的有效性，下面針對(duì)復(fù)雜構(gòu)造模型進(jìn)行試算.以Marmousi速度模型(圖3)為例，τ=1 ms，記錄長(zhǎng)度為4 s，震源為主頻為30 Hz的Ricker子波，混合吸收邊界條件的過渡帶點(diǎn)數(shù)為10，自適應(yīng)變空間算子長(zhǎng)度方案參數(shù)為fmax= 75 Hz，η=10-4，差分算子長(zhǎng)度參數(shù)M為2～6.觀測(cè)系統(tǒng)參數(shù)為：地面放炮，炮點(diǎn)x坐標(biāo)范圍為0～5000 m，炮間距為250 m，共21個(gè)炮點(diǎn)；檢波點(diǎn)位于x=0 m，檢波點(diǎn)深度范圍為0～3500 m，檢波點(diǎn)間距為10 m，共351個(gè)檢波點(diǎn).圖16a—c分別為炮點(diǎn)位于(0 m，0 m)的上行波、下行波、全波的VSP地震記錄.圖17為Marmousi速度模型VSP逆時(shí)偏移成像結(jié)果，由圖可知：

(1) 上行波、下行波及全波VSP逆時(shí)偏移成像位置準(zhǔn)確，構(gòu)造清晰；

(2) 上行波作為有效波時(shí)(圖17a)，井旁構(gòu)造成像較清晰，但遠(yuǎn)井源距處成像結(jié)果不清晰，部分構(gòu)造無法顯現(xiàn)；

(3) 下行波作為有效波時(shí)(圖17b)，遠(yuǎn)井源距構(gòu)造成像效果較好，但近井源距成像結(jié)果不如上行波成像結(jié)果連續(xù)；

(4) 全波作為有效波時(shí)(圖17c)，井旁構(gòu)造更為清晰、連續(xù)，并且遠(yuǎn)井源距構(gòu)造也能較好成像.

波場(chǎng)分離后的上行波和下行波，都有部分波場(chǎng)信息的丟失，同時(shí)在波場(chǎng)分離中難免會(huì)損害部分有效信號(hào)，而全波波場(chǎng)信息較全，所以采用全波進(jìn)行VSP逆時(shí)偏移效果更好.

為驗(yàn)證自適應(yīng)優(yōu)化有限差分算法的效率性，我們做了一個(gè)計(jì)算效率測(cè)試，CPU的型號(hào)是Intel(R) Xeon(R) CPU E5-26400 @ 2.5GHz，分別采用固定算子長(zhǎng)度的泰勒級(jí)數(shù)展開有限差分方法(fixed-length TE-based FD method)、自適應(yīng)泰勒級(jí)數(shù)展開有限差分方法(adaptive TE-based FD method)、自適應(yīng)優(yōu)化有限差分方法(adaptive optimal-based FD method)三種有限差分方法實(shí)現(xiàn)Marmousi模型全波VSP逆時(shí)偏移.以上三種有限差分方法的M值如圖18所示.圖19a—c分別為以上三種有限差分方法在時(shí)刻t=1200 ms時(shí)的波場(chǎng)快照.由圖19可見，這三種有限差分方法的波場(chǎng)快照模擬精度相近，圖19d為固定算子長(zhǎng)度的泰勒級(jí)數(shù)展開有限差分方法M=7(圖19a)與自適應(yīng)優(yōu)化有限差分方法M=2～4(圖19c)波場(chǎng)快照的差值，由圖19d可知，兩者的誤差很小，即兩者的精度相近.本文對(duì)比了以上三種方法的效率(表1)，可得如下結(jié)論：

(1) 自適應(yīng)泰勒級(jí)數(shù)展開有限差分方法與固定算子長(zhǎng)度的泰勒級(jí)數(shù)展開有限差分方法對(duì)比，效率提高了15%左右；

(2) 自適應(yīng)優(yōu)化有限差分方法與固定算子長(zhǎng)度的泰勒級(jí)數(shù)展開有限差分方法對(duì)比，計(jì)算量約節(jié)省了28%.

圖16 Marmousi模型的VSP地震記錄，炮點(diǎn)位于(0 m, 0 m)(a)上行波場(chǎng); (b)下行波場(chǎng); (c)全波波場(chǎng).Fig.16 VSP seismic records for Marmousi model with source located at (0 m, 0 m)(a) Up-going wavefield; (b) Down-going wavefield; (c) Full-wavefield.

圖17 VSP逆時(shí)偏移成像結(jié)果(a) 以上行波場(chǎng)為有效波; (b) 以下行波場(chǎng)為有效波; (c) 以全波波場(chǎng)為有效波.Fig.17 VSP RTM imaging result(a) Up-going wavefield as input; (b) Down-going wavefield as input; (c) Full-wavefield as input.

圖19 Marmousi模型波場(chǎng)快照(a) 固定算子長(zhǎng)度的泰勒級(jí)數(shù)展開有限差分法; (b) 自適應(yīng)泰勒級(jí)數(shù)展開有限差分法; (c) 自適應(yīng)優(yōu)化有限差分方法; (d) 為(a)與(c)的誤差圖.Fig.19 Snapshots for Marmousi model(a) Fixed-length TE-based FD method FD method； (b) Adaptive TE-based FD method FD method； (c) Adaptive optimal-based FD method； (d) Error of (a) and (c).

圖18 優(yōu)化有限差分方法與泰勒級(jí)數(shù)展開有限差分方法的算子長(zhǎng)度參數(shù)M隨速度變化圖Fig.18 Variation of FD operator length parameter M with velocity by optimal-based FD method and TE-based FD method

圖19表明了自適應(yīng)優(yōu)化有限差分方法的精度，圖18和表1表明了該方法的效率性.綜上所述，相同精度前提下，在固定算子長(zhǎng)度的泰勒級(jí)數(shù)展開有限差分方法中引入自適應(yīng)變空間差分算子長(zhǎng)度方案可以提高計(jì)算效率；再用優(yōu)化有限差分方法替代泰勒級(jí)數(shù)展開有限差分方法可以進(jìn)一步提高計(jì)算效率.

表1 自適應(yīng)優(yōu)化有限差分方法逆時(shí)偏移效率Table 1 The efficiency of the adaptive optimal-based FD RTM

3.2.3 實(shí)際資料處理

為進(jìn)一步驗(yàn)證自適應(yīng)優(yōu)化有限差分方法的VSP逆時(shí)偏移的有效性與正確性，本文將該有限差分方法應(yīng)用于某地區(qū)的實(shí)際VSP資料逆時(shí)偏移中.原始地震記錄為395炮，炮間距不固定，單個(gè)共檢波點(diǎn)道集如圖20a所示，其炮點(diǎn)沒有位于整網(wǎng)格點(diǎn)上.本文采用抗泄露Fourier變換規(guī)則化重建方法(Gao et al., 2011)對(duì)圖20a處理，規(guī)則化后共411炮，炮間距為10 m，如圖20b所示.檢波器深度為610～1000 m，共40道，道間距為10 m，每道記錄長(zhǎng)度為3.5 s，h=10 m，τ=1 ms.速度模型如圖21所示.我們采用自適應(yīng)優(yōu)化有限差分方法、混合吸收邊界條件、震源歸一化零延遲互相關(guān)成像條件及拉普拉斯濾波去噪方法實(shí)現(xiàn)該VSP實(shí)際資料的逆時(shí)偏移，其中自適應(yīng)變空間算子方案參數(shù)為fmax=75 Hz，η=10-4，其差分算子長(zhǎng)度參數(shù)M=2～6，混合吸收邊界過渡區(qū)網(wǎng)格點(diǎn)數(shù)為10.我們采用共檢波點(diǎn)道集實(shí)現(xiàn)VSP逆時(shí)偏移.有效成像區(qū)域的偏移結(jié)果如圖22 a所示，圖22 b為實(shí)際資料的地面常規(guī)偏移的結(jié)果，圖22 c為VSP逆時(shí)偏移剖面嵌入地面地震偏移剖面圖.由圖22可見，VSP逆時(shí)偏移成像同相軸清晰可辨，與地面地震成像結(jié)果構(gòu)造的走向相同，VSP剖面與地面地震剖面部分層位的略有錯(cuò)位，這種錯(cuò)位主要是速度模型不準(zhǔn)導(dǎo)致的，準(zhǔn)確的速度模型的建立還有待進(jìn)一步的研究.

圖21 實(shí)際資料速度模型Fig.21 Field data velocity model

4 討論

相對(duì)于地面地震資料而言，VSP資料能更好地反映井旁信息，所以理論上對(duì)井旁的成像處理有著天然的優(yōu)勢(shì).如果能結(jié)合VSP資料與地面資料進(jìn)行井地聯(lián)合逆時(shí)偏移，理論上可以更加清晰地描述地下構(gòu)造.以圖3所示的Marmousi速度模型為例，分別實(shí)現(xiàn)地面逆時(shí)偏移和VSP逆時(shí)偏移，地面勘探的觀測(cè)系統(tǒng)為地面放炮，炮點(diǎn)x坐標(biāo)范圍為0～4900 m，炮間距為100 m，共50炮，地面接收，接收點(diǎn)x坐標(biāo)范圍為0～4990 m，檢波點(diǎn)間距為10 m，共500個(gè)檢波點(diǎn).VSP資料的炮點(diǎn)觀測(cè)系統(tǒng)與地面勘探相同，其檢波點(diǎn)坐標(biāo)位于x=0 m，井中接收點(diǎn)深度為0～3490 m，檢波點(diǎn)間距為10 m，共有350個(gè)檢波點(diǎn).地面資料與VSP資料逆時(shí)偏移的參數(shù)均為：τ=1 ms，總記錄時(shí)間為4 s，混合吸收邊界過渡區(qū)網(wǎng)格點(diǎn)數(shù)為10，震源為主頻為30 Hz的Ricker子波，自適應(yīng)變空間算子長(zhǎng)度方案參數(shù)為fmax=75 Hz，η=10-4，M=2～6.地面逆時(shí)偏移(圖23 a)在深度為0～1000 m的水平構(gòu)造和小傾角構(gòu)造清晰，深度為2000～3000 m的角度不整合、地層隆起都得到了較好成像.VSP逆時(shí)偏移(圖23 b)是以全波作為逆時(shí)延拓的有效波，VSP逆時(shí)偏移成像的井旁構(gòu)造和中淺部大傾角構(gòu)造(如箭頭處)的成像效果比地面逆時(shí)偏移的效果好.圖23c為兩口井分別位于x=0 m和x=5000 m的VSP記錄以及地面地震記錄進(jìn)行聯(lián)合逆時(shí)偏移成像(井地聯(lián)合逆時(shí)偏移)結(jié)果.圖23a—c結(jié)果表明：井地聯(lián)合成像結(jié)果(圖23c)結(jié)合了VSP逆時(shí)偏移與地面逆時(shí)偏移的優(yōu)點(diǎn)，其成像在井旁構(gòu)造、大和小傾角構(gòu)造、水平層狀構(gòu)造、角度不整合和地層隆起構(gòu)造區(qū)域成像結(jié)果都較好.針對(duì)于同時(shí)進(jìn)行地面地震勘探與VSP地震勘探的區(qū)域，可以采用井地聯(lián)合偏移方法實(shí)現(xiàn)地下構(gòu)造的成像以進(jìn)一步提高成像精度.

圖22 (a) VSP逆時(shí)偏移結(jié)果; (b) 地面常規(guī)偏移結(jié)果; (c) VSP成像結(jié)果與地面成像結(jié)果的嵌入圖Fig.22 (a) VSP RTM result; (b) Surface conventional migration result; (c) VSP RTM result merged into surface migration result

圖23 (a) 地面逆時(shí)偏移; (b) VSP逆時(shí)偏移; (c) 井地聯(lián)合逆時(shí)偏移Fig.23 (a) Surface RTM result; (b) VSP RTM;(c) Surface combining with VSP RTM

5 結(jié)論

本文采用自適應(yīng)優(yōu)化有限差分方法求解聲波波動(dòng)方程，該方法與泰勒級(jí)數(shù)展開有限差分方法的對(duì)比表明：在算子長(zhǎng)度相同前提下，優(yōu)化有限差分方法模擬精度和VSP逆時(shí)偏移成像結(jié)果精度更高；在精度相同的條件下，優(yōu)化有限差分方法數(shù)值模擬及VSP逆時(shí)偏移效率更高；引入自適應(yīng)變差分算子長(zhǎng)度方案可以進(jìn)一步提高VSP逆時(shí)偏移計(jì)算效率.VSP逆時(shí)偏移中，全波逆時(shí)偏移比上行波逆時(shí)偏移成像結(jié)果波場(chǎng)信息更全、構(gòu)造形態(tài)更連續(xù)、成像范圍更廣，因?yàn)槎啻尾ㄔ谠黾诱彰鞫纫约巴卣钩上穹秶鸬搅溯^大的作用.本文通過模型試算與實(shí)際資料VSP逆時(shí)偏移處理有效驗(yàn)證了自適應(yīng)優(yōu)化有限差分方法聲波VSP逆時(shí)偏移的有效性.

致謝 感謝國(guó)家自然科學(xué)基金項(xiàng)目(41074100，41474110)、教育部新世紀(jì)優(yōu)秀人才支持計(jì)劃(NCET-10-0812)和殼牌地球物理優(yōu)秀博士生獎(jiǎng)學(xué)金的資助.

Alejandro A V, Biondo B. 2002. Deconvolution imaging condition for reverse-time migration, Stanford Exploration Project, 112: 83-96.

Baysal E, Kosloff D D, Sherwood J W C. 1983. Reverse time migration.Geophysics, 48(11): 1514-1524.

Campbell A, Gross J, Walker B, et al. 2002. Evaluation of a near-salt reservoir with an offset VSP (a case study in the Gulf of Mexico). 72nd SEG meeting, Salt Lake City, Utah, USA, Expanded Abstracts, 2353-2356.

Causse E, Ursin B. 2000. Viscoacoustic reverse time migration.JournalofSeismicExploration, 9(1): 165-184.

Chang W, McMechan G A. 1986. Reverse-time migration of offset vertical seismic profiling data using the excitation-time imaging condition.Geophysics, 51(1): 67-84.

Chang W, McMechan G A. 1987. Elastic reverse-time migration.Geophysics, 52(10): 1365-1375.

Chang W, McMechan G A. 1990. 3D acoustic prestack reverse-time migration.GeophysicalProspecting, 38(7): 737-756.Chang W, McMechan G A. 1993. 3D prestack migration in anisotropic media.Geophysics, 58(1): 79-90.Chang W, McMechan G A. 1994. 3D elastic prestack reverse-time depth migration.Geophysics, 59(4): 597-609.

Chattopadhyay S, Mcmechan G A. 2008. Imaging conditions for prestack reverse time migration.Geophysics, 73(3): S81-S89.

Chon Y, Souza J, Planchart. 2003. Offset VSP imaging with elastic reverse-time migration. 73rd SEG meeting, Dallas, Texas, USA, Expanded Abstracts, 1079-1053.

Clayton R W, Engquist B. 1977. Absorbing boundary conditions for acoustic and elastic wave equations.BulletinoftheSeismologicalSocietyofAmerica, 6: 1529-1540.Dablain M A. 1986. The application of high-order differencing to the scalar wave equation.Geophysics, 51(1): 54-66

Du Q, Zhu Y, Zhang M Q, et al. 2013. A study on the strategy of low wavenumber noise suppression for prestack reverse-time depth migration.ChineseJournalofGeophysics, 56(7): 2391-2401, doi:10.6038/cjg20130725.

Fleury C. 2013. Increasing illumination and sensitivity of reverse-time migration with internal multiples.GeophysicalProspecting, 2013, 61: 891-906.Gao J, Chen X, Wang F F, et al. 2011. Study on regularized reconstruction of uneven seismic traces.ProgressinGeophysics, 26(3): 983-991, doi:10.3969/j.issn.1004-2903.2011.03.025.Harwijanto J A, Wapenaar C P A, Berkhout A J. 1987. VSP migration by shot record inversion.FirstBreak, 5(7): 247-255.

Hokstad K, Mittet R, Landral M. 1988. Elastic reverse time migration of marine walkaway vertical seismic profiling data.Geophysics, 63(5): 1685-1695.

Hornby B, Yu J, Sharp J A, et al. 2006. VSP: Beyond time-to-depth.TheLeadingEdge, 25: 446-452.

Hou A, He Q. 1990. VSP reverse time migration.JournalofChangchunUniversityofEarthScience(in Chinese), 20(2): 227-233.

Hu H, Liu Y, Chang X, et al. 2013. Analysis and application of boundary treatment for the computation of reverse-time migration.ChineseJ.Geophys. (in Chinese), 56(6): 2033-2042, doi:10.6038/cjg20130624.Ketil H, Rune M, Martin L. 1998. Elastic reverse time migration of marine walkaway vertical seismic profiling data.Geophysics, 63(5): 1685-1695.

Li B, Liu H, Liu G F. 2010. Computational strategy of seismic pre-stack reverse time migration on CPU/GPU.ChineseJ.Geophys. (in Chinese), 53(12): 2938-2943, doi:10.3969/j.issn.0001-5733.2010.12.017..

Li Z, Guo Z, Tian K. 2014. Least-square reverse time migration in visco-acoustic medium.ChineseJ.Geophys. (in Chinese), 57(1): 214-228, doi:10.6038/cjg20140118.

Liu H W, Li B, Liu H, et al. 2010a. The algorithm of high order finite difference pre-stack reverse time migration and GPU implementation.ChineseJ.Geophys. (in Chinese), 53(7):1725-1733.

Liu H W, Liu H, Zou Z. 2010b. The problems of denoise and storage in seismic reverse time migration.ChineseJ.Geophys. (in Chinese), 53(9): 2171-2180, doi:10.3969/j.issn.0001-5733.2010.07.024.

Liu Y, Sen M K. 2009. A new time-space domain high-order finite-different method for the acoustic wave equation.JournalofComputationalPhysics, 228: 8779-8806.

Liu Y, Sen M K. 2010. A hybrid scheme for absorbing edge reflections in numerical modeling of wave propagation.Geophysics, 75(2):A1-A6.Liu Y, Sen M K. 2011. Finite-difference modeling with adaptive variable-length spatial operators.Geophysics, 76(4): T79-T89.

Liu Y. 2013. Globally optimal finite-difference schemes based on least squares.Geophysics, 78(4): T113-T132.

Loewenthal D, Mufti I R. 1983. Reverse time migration in the spatial frequency domain.Geophysics, 48(5): 627-635.

Nicoletis L, Mendes M, Compte P, et al. 1995. Quantitative ray-born elastic migration of VSP’s data. 57th Annual International Meeting, European Association of Geoscientists and Engineers, Expanded Abstracts.

Soni A K, Verschuur D J. 2014. Full-wavefield migration of vertical seismic profiling data: using all multiples to extend the illumination area.GeophysicalProspecting, 62(4): 740-759.

Sun R, McMechan G A. 2010. Nonlinear reverse-time inversion of elastic offset vertical seismic profile data.Geophysics, 53(10): 1295-1302.

Sun W B, Sun Z D. 2010. VSP reverse time migration based on the pseudo-spectral method and its applications.ChineseJ.Geophys. (in Chinese), 53(9): 2196-2203, doi:10.3969/j.issn.0001-5733.2010.09.020.

Sun W, Sun Z, Zhu X. 2011. Amplitude preserved VSP reverse time migration for angle-domain CIGs extraction.AppliedGeophysics, 8(2): 141-149.

Wang B L, Gao J H, Chen W S, et al. 2012. Efficient boundary storage for seismic reverse time migration.ChineseJ.Geophys. (in Chinese), 55(7): 2412-2421, doi:10.6038/j.issn.0001-5733.2012.07.025.Whitmore N D. 1983. Iterative depth migration by backward time propagation. 53rd Annual International Meeting, SEG, 827-830.

Wu W, Lines L R, Lu H. 1996. Analysis of higher-order, finite-difference schemes in 3-D reverse time migration.Geophysics, 61(3): 845-856.

Xiao X, Leaney W S. 2010. Local vertical seismic profiling (VSP) elastic reverse-time migration and migration resolution: Salt-flank imaging with transmitted P-to-S waves.Geophysics, 75(2): S35-S47.

Xiao X, Schuster G T. 2009. Local migration with extrapolated VSP Green′s functions.Geophysics, 74(1): SI15-SI26.

Xu S, Zhang Y, Tang B. 2011. 3D angle gathers from reverse time migration.Geophysics, 76(2): S77-S92.

Yan H, Liu Y. 2013a. Acoustic prestack reverse time migration using the adaptive high-order finite-difference method in time-space domain.ChineseJ.Geophys. (in Chinese), 56(3): 971-984.Yan H, Liu Y. 2013b. Visco-acoustic prestack reverse-time migration based on the time-space domain adaptive high-order finite-difference method.GeophysicalProspecting, 2013, 61: 941-954.

Yoon K, Marfurt K J, Houston U, et al. 2004. Challenges in reverse time migration. 74th Annual International Meeting, SEG, Expanded Abstracts, 1057-1060.Zhang Y, Sun J, Gray S H. 2007. Reverse-time migration: Amplitude and implementation issues. 77th Annual International Meeting, SEG, Expanded Abstracts, 2145-2149.

Zhang Y, Sun J. 2009. Practical issues in reverse time migration:

true amplitude gathers, noise removal and harmonic-source encoding.FirstBreak, 26(1): 29-35.

Zhao Y, Liu Y, Ren Z. 2014. Viscoacoustic prestack reverse time migration based on the optimal time-space domain high-order finite-difference method.AppliedGeophysics, 11(1): 50-62.

Zhu J M, Dong M Y, Li C C. 1989. VSP reverse-time migration using two-way nonreflection wave equation.OilGeophysicalProspecting, 24(3): 256-270.

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

杜啟振，朱釔同，張明強(qiáng)等. 2013. 疊前逆時(shí)深度偏移低頻噪聲壓制策略研究. 地球物理學(xué)報(bào)，56(7)：2391-2401, doi:10.6038/cjg20130725.

高建軍，陳小宏，王芳芳等. 2011. 不規(guī)則地震道數(shù)據(jù)規(guī)則化重建方法研究．地球物理學(xué)進(jìn)展，26(3)：983-991, doi:10.3969/j.issn.1004-2903.2011.03.025.

侯安寧，何樵登. 1990. VSP數(shù)據(jù)的逆時(shí)偏移. 長(zhǎng)春地質(zhì)學(xué)院學(xué)報(bào)，20(2)：227-233.

胡昊，劉伊克，常旭等. 2013. 逆時(shí)偏移計(jì)算中的邊界處理分析及應(yīng)用. 地球物理學(xué)報(bào)，56(6)：2033-2042, doi:10.6038/cjg20130624.

李博，劉紅偉，劉國(guó)峰等. 2010. 地震疊前逆時(shí)偏移算法的CPU/GPU實(shí)施對(duì)策. 地球物理學(xué)報(bào)，53(12)：2938-2943, doi:10.3969/j.issn.0001-5733.2010.12.017.

李振春，郭振波，田坤. 2014. 黏聲介質(zhì)最小平方逆時(shí)偏移. 地球物理學(xué)報(bào)，57(1)：214-2288, doi:10.6038/cjg20140118.

劉紅偉，李博，劉洪等. 2010. 地震疊前逆時(shí)偏移高階有限差分算法及GPU實(shí)現(xiàn). 地球物理學(xué)報(bào)，53(7)：1725-1733, doi:10.3969/j.issn.0001-5733.2010.07.024.

劉紅偉，劉洪，鄒振等. 2010. 地震疊前逆時(shí)偏移中的去噪與存儲(chǔ). 地球物理學(xué)報(bào)，53(9)：2171-2180, doi:10.3969/j.issn.0001-5733.2010.09.017.

孫文博，孫贊東. 2010. 基于偽譜法的VSP逆時(shí)偏移及其應(yīng)用研究. 地球物理學(xué)報(bào)，53(9)：2196-2203, doi:10.3969/j.issn.0001-5733.2010.09.020.

王保利，高靜懷，陳文超等. 2012. 地震疊前逆時(shí)偏移的有效邊界存儲(chǔ)策略. 地球物理學(xué)報(bào)，55(7)：2412-2421, doi:10.6038/j.issn.0001-5733.2012.07.025.

朱金明，董敏煜，李承楚．1989. VSP的雙程無反射波動(dòng)方程逆時(shí)偏移. 石油地球物理勘探，24(3)：256-270.

(本文編輯 劉少華)

Full-wavefield VSP reverse-time migration based on the adaptive optimal finite-difference scheme

CAI Xiao-Hui1,2, LIU Yang1,2*, WANG Jian-Min3, WANG Wei-Hong3, REN Zhi-Ming1,2

1StateKeyLaboratoryofPetroleumResourcesandProspecting,ChinaUniversityofPetroleum,Beijing102249,China2CNPCKeyLaboratoryofGeophysicalProspecting,ChinaUniversityofPetroleum,Beijing102249,China3ExplorationandDevelopmentResearchInstituteofDaqingOilfieldCompanyLimited,Daqing163318,China

Vertical seismic profiling (VSP) data provides higher resolution information, lower environmental noise and better information beside wells than surface seismic data. Meanwhile, reverse-time migration (RTM) based on two-way wave equation has no dip limitations on the image. It is significant to implement the RTM on VSP data.The four key factors of RTM are wavefield extrapolation, absorbing boundary conditions (ABC), imaging condition and noise suppressing. For the first factor, the paper compares the finite-difference (FD) method based on optimal scheme (optimal-based FD method) with the FD method based on space Taylor series expansion (space TE-based FD method) and the FD method based on time-space Taylor series expansion (time-space TE-based FD method) to indicate the accuracy of the optimal-based FD method. Moreover, to improve the efficiency, the adaptive variable-length spatial operators method is introduced into the optimal-based FD method (adaptive optimal-based FD method). The adaptive optimal-based FD method is adopted to solve the acoustic wave equation and thus implements the wavefield extrapolation of RTM efficiently and accurately. Second, the hybrid ABC is used to suppress artificial reflections from the model boundaries caused by the limited computational boundaries. Then the paper compares normalized cross-correlation imaging condition with deconvolution imaging condition for imaging, which indicates the prior is more stable and efficient. Finally, the Laplace operator filtering is applied to remove the low frequency noises.The conventional migration from VSP data using only up-going primary wavefield often suffers from poor illumination and limited migration scope. To enhance the illumination, enlarge the migration scope and improve the accuracy, the analysis of VSP RTM imaging with direct wave, primary wave, first-order multiple, second-order multiple and higher-order multiples implies that the full-wavefield (complete acoustic wavefield data without direct wave) can be treated as effective wave to implement prestack RTM.RTM test on VSP model data shows that the adaptive optimal-based FD method is more efficient and accurate than the fixed-length TE-based FD method, and suitable to VSP RTM. VSP RTM result using full-wavefield is clearer and images a larger area than VSP RTM only using up-going field. Imaging test on VSP field data demonstrates the efficiency and accuracy of our method.

VSP; Adaptive optimal-based FD method; Full acoustic wavefield; RTM

蔡曉慧， 劉洋， 王建民等. 2015. 基于自適應(yīng)優(yōu)化有限差分方法的全波VSP逆時(shí)偏移.地球物理學(xué)報(bào)，58(9):3317-3334,

10.6038/cjg20150925.

Cai X H, Liu Y, Wang J M, et al. 2015. Full-wavefield VSP reverse-time migration based on the adaptive optimal finite-difference scheme.ChineseJ.Geophys. (in Chinese),58(9):3317-3334,doi:10.6038/cjg20150925.

10.6038/cjg20150925

P631

2014-08-14，2015-04-16收修定稿

國(guó)家自然科學(xué)基金項(xiàng)目(41074100，41474110)和教育部新世紀(jì)優(yōu)秀人才支持計(jì)劃(NCET-10-0812)聯(lián)合資助.

蔡曉慧，女，1990年生，中國(guó)石油大學(xué)(北京)博士研究生，主要從事VSP地震資料處理和逆時(shí)偏移成像方面的研究工作. E-mail:caixh_1990@sina.com

*通訊作者 劉洋，男，1972年生，博士，教授，主要從事地震波正演、成像、反演及多分量地震等方面的研究工作.E-mail：wliuyang@vip.sina.com