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

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

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

?

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

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

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

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

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

1 引言

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

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

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

2 方法原理

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

二維聲波波動方程為

(1)

其中，p為標量波場，v為速度，t為時間，x、z為空間坐標.

VSP逆時偏移的震源正向延拓和地震記錄逆時延拓的表達式分別為

(2)

(3)

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

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

(4)

(5)

(6)

(7)

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

(8)

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

(9)

(10)

當δmax越趨向于1，頻散越小.

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

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

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

(11)

(12)

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

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

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

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

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

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

2.3 混合吸收邊界條件

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

(13a)

(13b)

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

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

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

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

(14)

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

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

(15)

(16)

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

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

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

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

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

3.2 VSP逆時偏移

3.2.1 VSP全波逆時偏移分析

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

3.2.2 模型試算

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

(I) 層狀模型

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

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

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

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

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

圖10 (a) 層狀速度模型; (b) 自適應(yīng)泰勒級數(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.觀測系統(tǒng)為地面放炮，炮點x坐標范圍為0～2000 m，炮間距為50 m，共41個炮點.檢波點為x=1950 m,z=0～3000 m，檢波點間距為10 m，共301個檢波點.圖12a—c分別為x=1950 m、z=0 m的炮點波場分離后上行波記錄、下行波記錄和切除直達波的全波記錄.圖13為逆時偏移成像結(jié)果，可得：

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

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

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

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

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

圖13 VSP逆時偏移結(jié)果(a)上行波作為有效波; (b) 下行波作為有效波; (c) 全波波場作為有效波.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，記錄總時間為2 s，震源是主頻為30 Hz的Ricker子波，混合吸收邊界過渡區(qū)網(wǎng)格數(shù)為10，自適應(yīng)變空間算子長度方案參數(shù)為fmax=50 Hz，η=10-5.觀測系統(tǒng)參數(shù)為：地面放炮，炮點x坐標范圍為0～2000 m，炮間距為100 m，共21個炮；檢波點位于x=2000 m，檢波點深度范圍為0～2000 m，檢波點間距為10 m，共201個檢波點.圖15為VSP逆時偏移成像結(jié)果，可見：

(1) 上行波、下行波、全波波場分別作為有效波時，VSP逆時偏移成像位置都準確；

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

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

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

圖15 斷層模型VSP逆時偏移結(jié)果(a) 上行波作為有效波; (b) 下行波作為有效波; (c)全波波場作為有效波.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逆時偏移成像結(jié)果(圖15c)比僅用上行波或者下行波逆時偏移成像界面形態(tài)更清晰.

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

(III) Marmousi模型

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

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

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

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

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

波場分離后的上行波和下行波，都有部分波場信息的丟失，同時在波場分離中難免會損害部分有效信號，而全波波場信息較全，所以采用全波進行VSP逆時偏移效果更好.

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

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

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

圖16 Marmousi模型的VSP地震記錄，炮點位于(0 m, 0 m)(a)上行波場; (b)下行波場; (c)全波波場.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逆時偏移成像結(jié)果(a) 以上行波場為有效波; (b) 以下行波場為有效波; (c) 以全波波場為有效波.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模型波場快照(a) 固定算子長度的泰勒級數(shù)展開有限差分法; (b) 自適應(yīng)泰勒級數(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)化有限差分方法與泰勒級數(shù)展開有限差分方法的算子長度參數(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表明了該方法的效率性.綜上所述，相同精度前提下，在固定算子長度的泰勒級數(shù)展開有限差分方法中引入自適應(yīng)變空間差分算子長度方案可以提高計算效率；再用優(yōu)化有限差分方法替代泰勒級數(shù)展開有限差分方法可以進一步提高計算效率.

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

3.2.3 實際資料處理

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

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

4 討論

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

圖22 (a) VSP逆時偏移結(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) 地面逆時偏移; (b) VSP逆時偏移; (c) 井地聯(lián)合逆時偏移Fig.23 (a) Surface RTM result; (b) VSP RTM;(c) Surface combining with VSP RTM

5 結(jié)論

本文采用自適應(yīng)優(yōu)化有限差分方法求解聲波波動方程，該方法與泰勒級數(shù)展開有限差分方法的對比表明：在算子長度相同前提下，優(yōu)化有限差分方法模擬精度和VSP逆時偏移成像結(jié)果精度更高；在精度相同的條件下，優(yōu)化有限差分方法數(shù)值模擬及VSP逆時偏移效率更高；引入自適應(yīng)變差分算子長度方案可以進一步提高VSP逆時偏移計算效率.VSP逆時偏移中，全波逆時偏移比上行波逆時偏移成像結(jié)果波場信息更全、構(gòu)造形態(tài)更連續(xù)、成像范圍更廣，因為多次波在增加照明度以及拓展成像范圍起到了較大的作用.本文通過模型試算與實際資料VSP逆時偏移處理有效驗證了自適應(yīng)優(yōu)化有限差分方法聲波VSP逆時偏移的有效性.

致謝 感謝國家自然科學(xué)基金項目(41074100，41474110)、教育部新世紀優(yōu)秀人才支持計劃(NCET-10-0812)和殼牌地球物理優(yōu)秀博士生獎學(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.

附中文參考文獻

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

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

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

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

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

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

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

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

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

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

朱金明，董敏煜，李承楚．1989. VSP的雙程無反射波動方程逆時偏移. 石油地球物理勘探，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逆時偏移.地球物理學(xué)報，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收修定稿

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

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

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