蔡杰雄

(中國石油化工股份有限公司石油物探技術(shù)研究院,江蘇南京211103)

地震勘探所利用的數(shù)據(jù)主要是反射波。反射波場由背景速度擾動(dòng)和高波數(shù)速度擾動(dòng)共同構(gòu)成[1-2],因此基于反射波的成像分為疊前深度偏移成像(確定高波數(shù)擾動(dòng))和反射波層析成像(確定背景速度擾動(dòng))兩個(gè)步驟。偏移成像用于定位地下反射位置,其實(shí)現(xiàn)了從數(shù)據(jù)空間到成像空間的映射;層析成像用于反演傳播路徑上的速度場,其實(shí)現(xiàn)了從成像空間到模型空間的映射。當(dāng)然,這兩個(gè)步驟也可以耦合在一起進(jìn)行數(shù)據(jù)域的迭代反演成像(如全波形反演)。相比于常規(guī)射線走時(shí)層析,結(jié)合偏移成像在成像域進(jìn)行波動(dòng)方程線性化走時(shí)層析速度反演是當(dāng)前比較實(shí)用有效且精度較高的技術(shù)組合。

基于傍軸近似的常規(guī)射線類成像方法是目前應(yīng)用于地震數(shù)據(jù)層析及偏移的重要方法。傳統(tǒng)的射線追蹤方法一般局限于對射線路徑及走時(shí)的描述,其實(shí)現(xiàn)靈活高效、沒有傾角限制且容易拓展到起伏地表。但是,高頻近似下的常規(guī)射線追蹤認(rèn)為中心射線代表著地震波的主能量,在實(shí)現(xiàn)上僅僅利用中心射線來描述地震波傳播,這樣的近似處理只能反映地震波的高頻運(yùn)動(dòng)學(xué)特征。且在數(shù)值計(jì)算上對于復(fù)雜介質(zhì)可能存在焦散及多到達(dá)時(shí)問題,因此常規(guī)射線類成像方法(包括偏移和層析)在復(fù)雜構(gòu)造情況下的應(yīng)用效果并不十分理想。

高斯束傳播算子是射線傳播算子的發(fā)展,HILL[3-4],HALE[5]以及POPOV等[6-7]將高斯束方法應(yīng)用到偏移處理中并取得了較好的成像效果。Hill于1990年首先提出高斯束疊后深度偏移方法,并于2001年將該方法推廣到疊前深度偏移,利用共偏移距道集中的對稱性解決了高斯束疊前深度偏移中的執(zhí)行效率問題,非常適合海上拖纜地震數(shù)據(jù)的成像處理。對于陸上地震數(shù)據(jù)的成像處理,GRAY[8]于2005年提出了快速精確的、適合于陸地起伏地表數(shù)據(jù)的炮域高斯束疊前深度偏移方法,解決了高斯束偏移方法中的局部平面波分解因地表高程變化及地表速度橫向變化帶來的問題。自HILL奠定了高斯束偏移方法的理論基礎(chǔ)以來,衍生出了一系列實(shí)用化的束偏移技術(shù)(如控制束偏移和快速束偏移及自適應(yīng)束偏移等)。高斯束疊前深度偏移技術(shù)的實(shí)現(xiàn)主要包括單個(gè)獨(dú)立的高斯束傳播的求解及所有高斯束疊加成像兩個(gè)步驟。單個(gè)獨(dú)立的高斯束傳播分兩步求得,即通過運(yùn)動(dòng)學(xué)射線追蹤求取中心射線的路徑及走時(shí),通過動(dòng)力學(xué)射線追蹤獲取中心射線附近的走時(shí)及高頻能量分布。利用相互獨(dú)立的高斯束描述波傳播,既保持了射線方法的高效性和靈活性,又考慮了波場的動(dòng)力學(xué)特征。高斯束偏移利用相互獨(dú)立的高斯束疊加并成像,解決了射線類方法中的多路徑問題,兼具了初至波到達(dá)時(shí)Kirchhoff積分偏移的靈活性和波動(dòng)方程偏移的精確性,是一種精致、精確且實(shí)現(xiàn)上靈活高效的深度偏移方法。

與成像方法類似,在層析速度反演中也有介于常規(guī)射線理論和波動(dòng)理論之間的方法,如胖射線層析[9]、菲涅爾體層析[10-12]、高斯束層析[13]、波路徑層析[14]等。胖射線層析從射線層析向波動(dòng)層析走近了一步,但其理論基礎(chǔ)不夠嚴(yán)格,僅僅是對波動(dòng)理論感性認(rèn)識(shí)的結(jié)果?；诜颇鶢栿w的波形層析方法是近年來發(fā)展較快的重要方法。它利用射線周圍的波動(dòng)性質(zhì)來同時(shí)擬合走時(shí)和波場的振幅,其本質(zhì)是有限頻層析的一種。SEMTCHENOK等[13]最早提出了高斯束層析方法,利用高斯束展布范圍內(nèi)的波場分布建立層析方程,對于一個(gè)炮檢對通過建立多個(gè)層析方程來減小層析反演的病態(tài)性,在一定程度上提高了反演穩(wěn)定性。但是SEMTCHENOK在建立層析敏感度核函數(shù)(以下簡稱為核函數(shù))時(shí)未給出詳細(xì)的推導(dǎo)過程,僅僅是根據(jù)直觀認(rèn)識(shí)利用高斯束的局部波場振幅作為加權(quán)系數(shù)進(jìn)行求解計(jì)算,其理論基礎(chǔ)不嚴(yán)格。上面所述的層析速度反演方法大多在數(shù)據(jù)域?qū)崿F(xiàn),在面對復(fù)雜構(gòu)造低信噪比數(shù)據(jù)時(shí)，受疊前數(shù)據(jù)復(fù)雜波場及低信噪比的影響,實(shí)際應(yīng)用較為困難。邵榮峰等[15]將高斯束層析方法應(yīng)用到成像域,并通過自動(dòng)拾取等技術(shù)實(shí)現(xiàn)了偏移速度分析,萬弘等[16]將構(gòu)造約束等實(shí)用化方法加入到高斯束層析中以提高反演的穩(wěn)健性,但兩者的理論基礎(chǔ)仍然與文獻(xiàn)[13]相似,同樣是簡單采用了高斯束振幅作為加權(quán)系數(shù)構(gòu)建層析核函數(shù)。李輝等[17]從高斯束作為正演工具的角度重新推導(dǎo)了成像域高斯束層析核函數(shù),推導(dǎo)得到的層析核函數(shù)仍然是以高斯束振幅作為加權(quán)系數(shù)進(jìn)行計(jì)算。為了更好地將偏移成像與層析反演結(jié)合起來,本文將直接從偏移成像條件出發(fā),推導(dǎo)建立成像走時(shí)擾動(dòng)與速度擾動(dòng)的線性關(guān)系,更直觀地體現(xiàn)層析與偏移的關(guān)系,給出新的成像域走時(shí)層析核函數(shù)。鑒于高斯束傳播算子在疊前深度偏移領(lǐng)域中具有較強(qiáng)的實(shí)用性[18-22],且能夠高效提取三維共方位-反射角度成像道集[23],本文從高斯束偏移角度道集出發(fā),在波動(dòng)方程的一階Born近似和Rytov近似下,推導(dǎo)成像域反射波走時(shí)層析方程及其核函數(shù),從而發(fā)展一種新的成像域的波動(dòng)方程線性化走時(shí)層析反演方法,并利用高斯束傳播算子計(jì)算該走時(shí)層析核函數(shù)(簡稱基于高斯束傳播算子的成像域波動(dòng)方程線性化走時(shí)層析方法為高斯束層析),從而形成結(jié)合高斯束偏移及高斯束層析的迭代反演建模及成像方法和流程。

結(jié)合疊前深度偏移成像提取的角度域共成像點(diǎn)道集,成像域的反射波層析成像可以將反投影路徑分解成兩個(gè)透射波分支,類似于高斯束或單程波偏移成像的處理方式。其基本思想是認(rèn)為成像點(diǎn)的走時(shí)擾動(dòng)由炮點(diǎn)到成像點(diǎn)的走時(shí)擾動(dòng)以及成像點(diǎn)到檢波點(diǎn)的走時(shí)擾動(dòng)之和決定。從而,成像域反射波走時(shí)層析成像可以建立起包含兩個(gè)核函數(shù)分支的走時(shí)擾動(dòng)對速度擾動(dòng)的線性關(guān)系。這兩個(gè)核函數(shù)分別代表從炮點(diǎn)到成像點(diǎn)以及從成像點(diǎn)到檢波點(diǎn)的透射波路徑?？紤]到陸上地震數(shù)據(jù)采樣時(shí)的地表非水平、空間采樣不規(guī)則和數(shù)據(jù)信噪比低,為了更方便高效地提取方位-反射角成像道集,將偏移成像與層析成像更好地配套起來,本文選擇高斯束算子作為波傳播算子用于計(jì)算格林函數(shù),發(fā)展起高斯束偏移和高斯束層析配套方法。相對于常規(guī)射線層析,由于高斯束層析核函數(shù)具有更嚴(yán)格的理論推導(dǎo)、更逼近波傳播實(shí)際而反演精度更高,同時(shí)由于構(gòu)建的層析矩陣稀疏性較低而使得反演穩(wěn)定性更強(qiáng)。基于高斯束傳播算子的成像域走時(shí)層析方法與高斯束偏移相結(jié)合,可形成具有典型特征波成像特點(diǎn)、適應(yīng)復(fù)雜構(gòu)造低信噪比數(shù)據(jù)的成像與建模工具。

1　高斯束疊前深度偏移及方位-反射角度道集提取方法

高斯束偏移從實(shí)現(xiàn)上講,主要包括運(yùn)動(dòng)學(xué)和動(dòng)力學(xué)射線追蹤及高斯束波場疊加兩個(gè)部分。前者給出了單個(gè)獨(dú)立高斯束的計(jì)算,后者則表達(dá)了空間中任意一點(diǎn)的格林函數(shù)可由其鄰域內(nèi)所有高斯束的疊加得到,如(1)式所示。該格林函數(shù)可用于實(shí)現(xiàn)高斯束疊前深度偏移,亦可用于求解高斯束層析核函數(shù)。

(1)

(1)式表示從震源點(diǎn)x=(x,y,z)出發(fā),在目標(biāo)點(diǎn)x′=(x′,y′,z′)處的格林函數(shù)由每個(gè)相互獨(dú)立的高斯束在該點(diǎn)的積分得到。其中,ω為圓頻率;p為慢度矢量,代表傳播方向;uGB表示從震源點(diǎn)x=(x,y,z)出發(fā),在目標(biāo)點(diǎn)x′=(x′,y′,z′)接收的單個(gè)高斯束波場,其表達(dá)為(2)式,表示高斯束的構(gòu)建通過運(yùn)動(dòng)學(xué)射線追蹤確定其中心射線的軌跡和走時(shí)τ(s),通過動(dòng)力學(xué)射線追蹤獲取中心射線附近的動(dòng)力學(xué)參數(shù)Q和M,從而得到單個(gè)高斯束波場。該表達(dá)式在文獻(xiàn)[3]中已經(jīng)給出并討論了其具體計(jì)算方法,本文僅給出表達(dá)式,不再贅述。

(2)

給出了高斯束計(jì)算格林函數(shù)的表達(dá)式之后,通過頻率域波場相關(guān)的成像條件即可實(shí)現(xiàn)高斯束疊前深度偏移。進(jìn)一步,根據(jù)高斯束疊前深度偏移成像條件,可以推導(dǎo)成像域高斯束走時(shí)層析核函數(shù),并通過公式(1) 利用高斯束積分表達(dá)格林函數(shù)實(shí)現(xiàn)該核函數(shù)的計(jì)算,從而將高斯束偏移與高斯束層析有機(jī)結(jié)合起來,實(shí)現(xiàn)復(fù)雜介質(zhì)的偏移與層析迭代。

對于炮域?qū)崿F(xiàn)的高斯束疊前深度偏移,需要從炮點(diǎn)和接收點(diǎn)分別進(jìn)行高斯束正演傳播,以求得從炮點(diǎn)出發(fā)的下行高斯束波場和從反射點(diǎn)出發(fā)到達(dá)檢波器的上行高斯束波場。然后類似于單程波波動(dòng)方程偏移,使用上行波場和下行波場的互相關(guān)成像條件:

(3)

(3)式給出的是單炮成像結(jié)果。其中,G(x,xs,ω)與G(x,xr,ω)分別表示從炮點(diǎn)xs和從檢波點(diǎn)xr到成像點(diǎn)x的格林函數(shù),“*”代表共軛;Ds(xs,xr,pr,ω)表示基于炮道集的加高斯窗局部傾斜疊加。將(1)式代入(3)式,進(jìn)一步表示成:

(4)

(4)式說明了任一單炮的成像結(jié)果是由從炮點(diǎn)出發(fā)的所有下行方向高斯束波場與從檢波點(diǎn)出發(fā)的所有上行方向高斯束波場的相關(guān)得到。當(dāng)速度存在誤差導(dǎo)致成像存在誤差時(shí),可以從(4)式出發(fā)推導(dǎo)建立速度誤差與成像誤差的線性關(guān)系,給出基于高斯束偏移成像道集的成像域高斯束走時(shí)層析核函數(shù),將在下一節(jié)詳細(xì)介紹。

在疊前深度偏移過程中提取方位-反射角度道集是進(jìn)行后續(xù)成像域走時(shí)層析速度反演的必備過程。計(jì)算地下成像點(diǎn)的方位-反射角的一種有效且高效的方法是分別估算從炮點(diǎn)和檢波點(diǎn)出發(fā)到達(dá)成像點(diǎn)處的波場傳播方向。一旦獲得兩者的傳播方向,即可通過簡單的向量代數(shù)運(yùn)算得到方位角和反射角(張角)。對于射線類偏移方法而言,這個(gè)過程相對容易,只需通過旅行時(shí)場的空間導(dǎo)數(shù)分別計(jì)算震源波場的射線慢度ps和檢波點(diǎn)波場射線慢度pr即可。對高斯束函數(shù)式((2)式)解析求解,即可方便高效地計(jì)算旅行時(shí)場的空間導(dǎo)數(shù)。

假設(shè)高斯束鄰域內(nèi)任意一點(diǎn)Q,其對應(yīng)的中心射線上的點(diǎn)R,根據(jù)(2)式所示的高斯束函數(shù),則點(diǎn)Q的旅行時(shí)可以由點(diǎn)R的旅行時(shí)及動(dòng)力學(xué)射線追蹤參量表示為:

(5)

得到高斯束的旅行時(shí)可進(jìn)一步通過旅行時(shí)場的空間導(dǎo)數(shù)得到點(diǎn)Q的高斯束波場慢度矢量P(Q)=?T(Q)/?x,因此可以得到慢度矢量的解析表達(dá)式:

(6)

k=1,2,3

分別得到炮點(diǎn)波場和檢波點(diǎn)波場的慢度矢量后,根據(jù)反射角和方位角的定義,即可給出具體求解方法,詳見參考文獻(xiàn)[21],此處不再贅述。

2　基于高斯束成像道集的成像域高斯束層析方法

層析反演方法可以統(tǒng)一表示成正問題的一階近似:

(7)

仍然基于層析方程式(7),本文從高斯束偏移成像條件出發(fā),推導(dǎo)出基于高斯束積分計(jì)算格林函數(shù)的成像域波動(dòng)方程線性化走時(shí)層析(簡稱為高斯束層析)核函數(shù)。從而實(shí)現(xiàn)基于高斯束偏移成像道集的成像域高斯束層析方法。

從公式(4)出發(fā),去掉所有的積分項(xiàng),僅考慮任意一個(gè)炮檢對、任一傳播方向的成像結(jié)果(所有炮檢對、所有方向高斯束的成像結(jié)果疊加即為單炮成像結(jié)果),將角度域高斯束偏移(GBM)成像條件重寫為(8)式。為了后續(xù)表達(dá)易于區(qū)分,將從震源出發(fā)的下行高斯束波場表示為S(x,ps,ω;xs,xr),將從檢波點(diǎn)出發(fā)的上行高斯束波場表示為R(x,pr,ω;xs,xr),有:

(8)

(8)式是對應(yīng)頻率域任一炮檢對某一傳播方向p的高斯束成像結(jié)果(共方位-反射角成像道集)。其中,x=(x,y,z)表示成像點(diǎn)坐標(biāo);θ,φ分別表示成像點(diǎn)的反射張角及方位角。

在波動(dòng)方程的一階Born近似下,波場U可以分解為背景波場U0和擾動(dòng)波場ΔU:

(9)

因此從成像條件(8)式可以近似得到擾動(dòng)像:

(10)

式中:ΔS,ΔR分別為一階Born近似散射場;S0,R0分別為從炮點(diǎn)和檢波點(diǎn)出發(fā)到成像點(diǎn)的背景波場。該成像條件表明,成像點(diǎn)x處像的擾動(dòng)來自炮點(diǎn)端和檢波點(diǎn)端兩個(gè)分支的影響。

進(jìn)一步地,根據(jù)文獻(xiàn)[24]和文獻(xiàn)[25],一階Born近似散射場ΔS與ΔR可以表示為:

式中:VS和VR分別表示從炮點(diǎn)和從檢波點(diǎn)到成像點(diǎn)的Born波路徑;k0=ω/v0表示背景模型波數(shù);Δv為待反演的速度擾動(dòng);G(x,p,ω;y)表示由高斯束積分計(jì)算的從點(diǎn)y到點(diǎn)x的格林函數(shù);S0(y,ps,ω;xs,xr)和R0(y,pr,ω;xs,xr)分別表示在背景速度模型中從炮點(diǎn)和檢波點(diǎn)傳播到空間任意一點(diǎn)y的波場。

如果將(11)式代入(10)式,形式上可以給出像的擾動(dòng)(左端項(xiàng))與速度擾動(dòng)(右端項(xiàng))的關(guān)系式。但實(shí)際操作時(shí),顯然不能直接利用該式進(jìn)行層析反演。對比數(shù)據(jù)域?qū)游龇囱菘梢苑治?數(shù)據(jù)域反演利用正演數(shù)據(jù)與實(shí)測數(shù)據(jù)的差在某種范數(shù)下(一般是二范數(shù))最小作為誤差泛函,其實(shí)測數(shù)據(jù)是客觀的,可直接用來做逼近標(biāo)準(zhǔn);如果直接利用(10)式進(jìn)行成像域反演,則由于客觀上無法得到真實(shí)像IGBM(x,θ,φ,ω)從而無法得到擾動(dòng)像ΔIGBM(x,θ,φ,ω),因此其本質(zhì)是利用一個(gè)未知的中間量來估計(jì)另一未知量Δv。直接利用像的擾動(dòng)這個(gè)概念來進(jìn)行反演缺乏嚴(yán)格的判斷標(biāo)準(zhǔn)。從另一個(gè)角度分析,像的擾動(dòng)是一個(gè)綜合概念,其實(shí)質(zhì)包括了走時(shí)(位置)擾動(dòng)和振幅擾動(dòng)。實(shí)際計(jì)算時(shí)需要將像的擾動(dòng)退化為走時(shí)擾動(dòng)(深度擾動(dòng))或振幅擾動(dòng),從而分別建立與速度擾動(dòng)的關(guān)系式。由于像域的振幅擾動(dòng)影響因素太復(fù)雜,實(shí)際層析反演一般退化為僅利用走時(shí)擾動(dòng)。

考慮退化到僅利用走時(shí)擾動(dòng)進(jìn)行層析反演,進(jìn)一步地,將(10)式重寫為:

(12)

整理得:

(13)

在波動(dòng)方程的Rytov近似下,波場可以表示為u=(A0+ΔA)ei(φ0+Δφ)。成像值是兩個(gè)波場相關(guān)得到,因此同樣可以表示I=(A0+ΔA)ei(φ0+Δφ),I0(x,θ,φ,w)=A0eiφ0,由于ΔA?A0,則(13)式可以進(jìn)一步表示成:

(14)

由于Δφ趨于零,則(14)式左端項(xiàng)近似等于iΔφ,兩邊取虛部得:

(15)

進(jìn)一步地,根據(jù)文獻(xiàn)[26],單頻相位擾動(dòng)與單頻走時(shí)擾動(dòng)有近似關(guān)系:

(16)

同時(shí)將(11)式及(16)式代入(15)式整理,最終得到成像域單頻走時(shí)擾動(dòng)與速度擾動(dòng)的關(guān)系式:

(17)

其中,KF為頻率域走時(shí)層析核函數(shù),其兩個(gè)分支分別表示為:

(18a)

(18b)

從(18)式可以看出,成像域走時(shí)層析核函數(shù)的表現(xiàn)形式與LIU等[12]給出的數(shù)據(jù)域菲涅爾體走時(shí)層析核函數(shù)的表現(xiàn)形式類似,其本質(zhì)都是有限頻核函數(shù),關(guān)鍵是背景速度下格林函數(shù)的計(jì)算。

需要說明的是,在背景模型中的波場可以表示成子波與格林函數(shù)的乘積,因此在推導(dǎo)得到(18)式的過程中約去了子波項(xiàng)(假設(shè)子波不變),僅剩下格林函數(shù)項(xiàng)。

由于實(shí)際操作時(shí),走時(shí)擾動(dòng)Δt在時(shí)空域(角度域成像道集)測量得到,與頻率無關(guān)。因此,定義帶限地震信號的走時(shí)擾動(dòng)可以用單頻走時(shí)擾動(dòng)加權(quán)疊加得到:

(19)

(20)

最終得到成像域帶限走時(shí)擾動(dòng)與速度擾動(dòng)關(guān)系式:

(21)

其中,KT為帶限走時(shí)層析核函數(shù),其兩個(gè)分支分別表示為:

(22a)

(22b)

至此,基于高斯束角度道集、在波動(dòng)方程的一階Born近似和Rytov近似下導(dǎo)出了成像域帶限走時(shí)層析方程(21)式及其對應(yīng)的核函數(shù)表達(dá)式(22)式。從(22)式可以看出,該核函數(shù)的本質(zhì)是有限頻核函數(shù),其求解主要是背景波場中格林函數(shù)的計(jì)算。利用(1)式的高斯束積分計(jì)算格林函數(shù)是一種精度較高且計(jì)算量較小的實(shí)用化計(jì)算方式。

需要說明的是,盡管最終得到的有限頻核函數(shù)在推導(dǎo)過程引入了多次線性化處理,但有限頻核函數(shù)相對于高頻近似的射線層析核函數(shù)仍然更接近于波動(dòng)傳播的實(shí)際情況,是對波動(dòng)層析的一種折中近似。用此方法替換常規(guī)射線層析能提高精度,且能借鑒常規(guī)射線層析流程實(shí)現(xiàn)技術(shù)的實(shí)用化。為進(jìn)一步提高對復(fù)雜地質(zhì)體(復(fù)雜構(gòu)造、強(qiáng)橫向變速等)的成像效果,結(jié)合高斯束成像的優(yōu)勢,有必要在后續(xù)研究中加強(qiáng)非線性項(xiàng)的合理近似和算法優(yōu)化的方法研究。

下面通過數(shù)值試驗(yàn)分析單頻和帶限走時(shí)層析核函數(shù)的特征。設(shè)計(jì)背景模型參數(shù):網(wǎng)格個(gè)數(shù)Nx=Ny=Nz=201,網(wǎng)格間距dx=dy=dz=10m,速度為3000m/s,震源和檢波點(diǎn)坐標(biāo)分別為(500,1000,1000)和(1500,1000,1000),帶限走時(shí)層析核函數(shù)的頻率范圍取0～40Hz,單頻核函數(shù)計(jì)算取中心頻率20Hz。圖1給出了對應(yīng)的單頻及帶限走時(shí)層析核函數(shù)。從圖1g和圖1h上可以看出,核函數(shù)在炮檢連線的中心線上的敏感度為零,這與常規(guī)射線層析的核函數(shù)僅分布在射線上的假設(shè)完全相反。劉玉柱[10]詳細(xì)分析了這種差異,并指出了常規(guī)射線層析之所以仍然能取得較好效果的原因。

圖1　單頻和帶限走時(shí)層析核函數(shù)a 對應(yīng)單頻10Hz的走時(shí)層析核函數(shù)(XZ切片); b 對應(yīng)單頻10Hz的走時(shí)層析核函數(shù)(YZ切片); c 對應(yīng)單頻20Hz的走時(shí)層析核函數(shù)(XZ切片); d 對應(yīng)單頻20Hz的走時(shí)層析核函數(shù)(YZ切片); e 對應(yīng)單頻30Hz的走時(shí)層析核函數(shù)(XZ切片); f 對應(yīng)單頻30Hz的走時(shí)層析核函數(shù)(YZ切片); g 帶限走時(shí)層析核函數(shù)(XZ切片); h 帶限走時(shí)層析核函數(shù)(YZ切片)

圖1反映的是單頻和帶限走時(shí)層析核函數(shù)。根據(jù)(22)式,該核函數(shù)僅僅是成像域走時(shí)層析核函數(shù)的一個(gè)分支,兩個(gè)分支(炮點(diǎn)到成像點(diǎn)及成像點(diǎn)到檢波點(diǎn))則構(gòu)成了成像域走時(shí)層析反演核函數(shù)的完整形態(tài),如圖2所示。用此帶有一定寬度的高斯束層析核函數(shù)替換常規(guī)射線層析核函數(shù)能更準(zhǔn)確逼近波實(shí)際傳播方式,提高層析反演精度和穩(wěn)定性。

圖2　成像域高斯束走時(shí)層析反演核函數(shù)

成像域高斯束層析反演的實(shí)際操作流程及具體計(jì)算方法可以完全借鑒常規(guī)射線層析反演的框架和算法,包括成像剖面層位及成像道集RMO的自動(dòng)拾取,矩陣求解,層位約束正則化方法等等。兩者的區(qū)別僅僅是在核函數(shù)的表達(dá)與計(jì)算上,因此該方法的實(shí)用性較強(qiáng)。

需要進(jìn)一步說明的是,成像域走時(shí)層析的具體實(shí)現(xiàn)與其走時(shí)誤差的具體計(jì)算形式有關(guān)。XIE等[27-29]給出了基于炮偏移后按照炮索引排列的成像道集的剩余時(shí)差(RMO)所對應(yīng)的核函數(shù)的表現(xiàn)形式,該核函數(shù)的兩個(gè)分支并不對稱,這與成像道集上每一個(gè)RMO代表一炮的成像走時(shí)誤差有關(guān)。本文所推導(dǎo)的核函數(shù)與XIE等給出的核函數(shù)表達(dá)形式一致,但

由于是針對高斯束偏移所提取的方位角度域成像道集,所提取的RMO僅與方位角及反射角有關(guān)[30],核函數(shù)則表現(xiàn)為關(guān)于剖面上反射軸法方向?qū)ΨQ的兩個(gè)分支,這一點(diǎn)并不會(huì)影響理論上的精度,但會(huì)使得層析的實(shí)現(xiàn)更加自然方便。

3　理論模型及實(shí)際數(shù)據(jù)應(yīng)用

3.1　理論模型數(shù)據(jù)

根據(jù)某地區(qū)實(shí)際地質(zhì)構(gòu)造設(shè)計(jì)的速度模型如圖3a 所示,斷層發(fā)育,地層高陡,速度橫向變化較大。橫縱向采樣點(diǎn)為731×550,橫向采樣間隔為10m,縱向采樣間隔為5m。利用聲波正演得到疊前炮記錄(圖3b),炮間距和道間距都是10m,中心放炮,每炮361道。利用等梯度速度模型作為層析初始模型(圖4a)。對初始模型進(jìn)行高斯束疊前深度偏移,輸出偏移剖面(圖4b)及角度道集(圖5),可以看出初始模型對應(yīng)的偏移剖面上的各層位成像深度并不準(zhǔn)確,繞射波也沒有完全收斂。圖5對比了CDP400位置處提取的初始模型角度道集和層析后角度道集,由于速度偏低,初始角度道集同相軸上翹;經(jīng)過層析反演更新速度模型,偏移后角度道集得到了拉平。圖6是高斯束層析經(jīng)過5次迭代后得到的速度模型及其相應(yīng)的高斯束偏移剖面。圖7是常規(guī)射線層析經(jīng)過8次迭代得到的速度模型及其高斯束偏移剖面。與常規(guī)射線層析對比可知,本文發(fā)展的高斯束層析方法能夠提供更加豐富的速度信息,迭代次數(shù)也小于常規(guī)射線層析。從抽取的單道速度曲線對比(圖8)可以看出,高斯束層析反演的速度值分辨率更高,更逼近真實(shí)速度。對比圖6b和圖7b的偏移深度可以看出,高斯束層析偏移結(jié)果更加逼近真實(shí)深度,反射界面歸位到正確的深度位置,繞射波收斂,斷層位置聚焦更好,更新后偏移提取的角度道集更加接近真實(shí)速度模型下偏移提取的角度道集(圖5),說明高斯束層析對復(fù)雜構(gòu)造模型的速度反演精度優(yōu)于常規(guī)射線層析方法。

圖3　某地區(qū)復(fù)雜速度模型(a)及其炮集記錄(b)

圖4　初始速度模型(a)及其高斯束疊前深度偏移剖面(b)

圖5　CDP400角度道集a 初始角度道集; b 正確速度模型角度道集; c 高斯束層析角度道集; d 常規(guī)射線層析角度道集

圖6　高斯束層析更新速度場(a)及其層析后高斯束偏移剖面(b)

3.2　實(shí)際地震資料

實(shí)際資料來自中國東北某復(fù)雜構(gòu)造工區(qū)。該區(qū)基底之下是上古生界變質(zhì)巖,變質(zhì)巖之上發(fā)育火山巖。目的層波場復(fù)雜,存在多組斷裂系統(tǒng),速度變化大,深度域建模困難,精確成像難度較大。處理要求縱向上能刻畫火山機(jī)構(gòu)特征和噴發(fā)期次,橫向上能分辨接觸關(guān)系。合理突出與火山巖接觸的地層反射。前期進(jìn)行了時(shí)間域速度建模,以此作為深度域初始模型進(jìn)行后續(xù)高斯束層析反演建模。

本文的研究工作是在前期處理人員進(jìn)行了預(yù)處理、疊前時(shí)間偏移及深度域初始建模及偏移成像的基礎(chǔ)上開展的,主要通過高斯束層析反演后的速度模型進(jìn)行高斯束疊前深度偏移來體現(xiàn)本文方法的實(shí)際應(yīng)用效果。從圖9所示的三維體模型上可以看出,高斯束層析反演的速度模型分辨率明顯高于常規(guī)射線層析反演結(jié)果。進(jìn)一步對比過井速度剖面及偏移剖面可以更清楚地反映該效果。

從圖10所示的層析反演的速度模型上看,高斯束層析結(jié)果體現(xiàn)了速度模型的更多細(xì)節(jié)信息,進(jìn)一步,通過與井曲線的對比可以看出,高斯束層析反演速度模型的精度更高(圖11)。

圖7　常規(guī)射線層析更新速度場(a)及其層析后高斯束偏移剖面(b)

圖8　真實(shí)速度模型、初始速度模型、射線層析反演及高斯束層析反演速度模型中抽取CDP400處的速度對比

圖9　某實(shí)際資料常規(guī)射線層析(a)與高斯束層析(b)速度模型體

圖10　某過井線射線層析速度模型(a)及高斯束層析速度模型(b)(綠線為井位置)

圖12為利用射線層析模型及高斯束層析模型對實(shí)際數(shù)據(jù)進(jìn)行高斯束偏移的結(jié)果,圖13為圖12中藍(lán)框的放大顯示。從圖12及圖13可以看出,同樣利用高斯束偏移方法,高斯束層析反演的速度模型對應(yīng)的偏移結(jié)果斷裂成像質(zhì)量更高,斷點(diǎn)更干脆,同相軸更連續(xù),整體成像質(zhì)量更高。從該實(shí)際資料處理可以認(rèn)識(shí)到,結(jié)合高斯束層析與高斯束偏移技術(shù)的處理方式,可以更好地適應(yīng)復(fù)雜構(gòu)造資料的深度域速度建模,提供更高質(zhì)量的成像結(jié)果。

圖11　射線層析模型、高斯束層析模型及實(shí)鉆井速度曲線對比

圖12　利用射線層析模型(a)及高斯束層析模型(b)對實(shí)際數(shù)據(jù)進(jìn)行了高斯束偏移的結(jié)果

圖13　利用射線層析模型(a)及高斯束層析模型(b)對實(shí)際數(shù)據(jù)進(jìn)行高斯束偏移的結(jié)果(局部放大顯示)

4　結(jié)論

本文從高斯束積分表達(dá)格林函數(shù)出發(fā),引出了高斯束偏移及提取方位反射角成像道集方法?；诮嵌扔蚋咚故瞥上駰l件,在波動(dòng)方程的一階Born近似和Rytov近似下,推導(dǎo)給出了成像域波動(dòng)方程線性化走時(shí)層析成像方程及其核函數(shù)表達(dá)式,并利用高斯束傳播算子計(jì)算該核函數(shù)。利用高斯束層析核函數(shù)替代常規(guī)射線層析核函數(shù)(該核函數(shù)為常數(shù)1),可以改進(jìn)層析反演精度,加快反演收斂,從而形成了基于高斯束算子的偏移成像與層析成像的聯(lián)合迭代速度建模與偏移處理方法,相對于常規(guī)方法,該方法具有更高的精度和更強(qiáng)的實(shí)用性。

本文發(fā)展的高斯束層析反演方法利用高斯束傳播算子計(jì)算成像域走時(shí)層析核函數(shù),提供了一種新的成像域波動(dòng)方程線性化近似走時(shí)層析反演方法。該方法主要利用反射數(shù)據(jù)的走時(shí)信息,對初始模型的依賴性較低。但該方法僅能反演速度模型的低波數(shù)成分,因而主要用于背景速度建模,為偏移成像服務(wù)及為更高精度的反演方法(如FWI)提供初始模型。高斯束層析方法與高斯束偏移技術(shù)相結(jié)合,可形成具有典型特征波成像特點(diǎn)的、適應(yīng)低信噪比數(shù)據(jù)的成像與建模工具,真正體現(xiàn)偏移成像與速度建模一體化的處理思想。

[1]王華忠,馮波,王雄文,等.特征波反演成像理論框架[J].石油物探,2017,56(1):38-49

WANG H Z,FENG B,WANG X W,et al.The theoretical framework of characteristic wave inversion imaging[J].Geophysical Prospecting for Petroleum,2017,56(1):38-49

[2]王華忠,馮波,王雄文,等.地震波反演成像方法與技術(shù)核心問題分析[J].石油物探,2015,54(2):115-125

WANG H Z,FENG B,WANG X W,et al.Analysis of seismic inversion imaging and its technical core issues[J].Geophysical Prospecting for Petroleum,2015,54(2):115-125

[3]HILL N R.Gaussian beam migration[J].Geophysics,1990,55(11):1416-1428

[4]HILL N R.Prestack Gaussian-beam depth migration[J].Geophysics,2001,66(4):1240-1250

[5]HALE D.Migration by the Kirchhoff,slant stack,and Gaussian beam methods[J].Center for Wave Phenomena,1992:CWP-126

[6]POPOV M M,SEMTCHENOK N M,VERDEL A R,et al.Reverse time migration with Gaussian beams and velocity analysis applications[J].Expanded Abstracts of 70thEAGE Annual Conference,2008:F048

[7]POPOV M M,SEMTCHENOK N M,VERDEL A R,et al.Depth migration by the Gaussian beam summation method[J].Geophysics,2010,75(2):S81-S93

[8]GRAY S H.Gaussian beam migration of common-shot records[J].Geophysics,2005,70(1):953-959

[9]VASCO D W,PETERSON J E,MAJER E L.Beond ray tomography:wavepaths and Fresnel volumes[J].Geophysics,1995,60(6):1790-1804

[10]劉玉柱.菲涅爾體地震層析成像理論與應(yīng)用研究[D].上海:同濟(jì)大學(xué),2011

LIU Y Z.Theory and applications of Fresnel volume seismic tomography[D].Shanghai:Tongji University,2011

[11]劉玉柱,謝春,楊積忠.基于Born波路徑的高斯束初至波波形反演[J].地球物理學(xué)報(bào),2014,57(9):2900-2909

LIU Y Z,XIE C,YANG J Z.Gaussian beam first-arrival waveform inversion based on Born wavepath[J].Chinese Journal of Geophysics,2014,57(9):2900-2909

[12]LIU Y Z,DONG L G,WANG Y M,et al.Sensitivity kernels for seismic Fresnel volume tomography[J].Geophysics,2009,74(5):U35-U46

[13]SEMTCHENOK N M,POPOV M M,VERDEL A R.Gaussian beam tomography[J].Expanded Abstracts of 71stEAGE Annual Conference,2009:U32

[14]BAKKER P,GERRITSEN S,CAO Q.3D RTM-based wave path tomography tested at a realistic scale[J].Expanded Abstracts of 77thEAGE Annual Conference,2015:WS05-C02

[15]邵榮峰,方伍寶,蔡杰雄,等.高斯束層析偏移速度建模方法及應(yīng)用[J].石油物探,2016,55(1):91-99

SHAO R F,FANG W B,CAI J X,et al.A method of migration velocity analysis based on Gaussian beam tomography and its application[J].Geophysical Prospecting for Petroleum,2016,55(1):91-99

[16]萬弘,楊勤勇,蔡杰雄,等.地質(zhì)構(gòu)造約束高斯束層析反演方法與應(yīng)用[J].石油物探,2017,56(5):707-717

WAN H,YANG Q Y,CAI J X,et al.A method of geological structure constrained tomographic inversion based on Gaussian beam and its application[J].Geophysical Prospecting for Petroleum,2017,56(5):707-717

[17]李輝,王華忠,劉守偉.基于高斯束的速度層析方法研究[J].石油物探,2017,56(1):116-125

LI H,WANG H Z,LIU S W.A velocity tomography algorithm based on Gaussian beam[J].Geophysical Prospecting for Petroleum,2017,56(1):116-125

[18]蔡杰雄,方伍寶,王華忠.高斯束深度偏移的實(shí)現(xiàn)與應(yīng)用研究[J].石油物探,2012,51(5):469-475

CAI J X,FANG W B,WANG H Z.Realization and application of Gaussian beam depth migration[J].Geophysical Prospecting for Petroleum,2012,51(5):469-475

[19]黃建平,楊繼東,李振春,等.基于有效鄰域波場近似的起伏地表保幅高斯束偏移[J].地球物理學(xué)報(bào),2016,59(6):2245-2256

HUANG J P,YANG J D,LI Z C,et al.An amplitude-preserved Gaussian beam migration based on wave field approximation in effective vicinity under irregular topographical conditions[J].Chinese Journal of Geophysics,2016,59(6):2245-2256

[20]岳玉波,李振春,錢忠平,等.復(fù)雜地表?xiàng)l件下保幅高斯束偏移[J].地球物理學(xué)報(bào),2012,55(4):1376-1383

YUE Y B,LI Z C,QIAN Z P,et al.Amplitude-preserved Gaussian beam migration under complex topographic conditions[J].Chinese Journal of Geophysics,2012,55(4):1376-1383

[21]劉強(qiáng),張敏,李振春,等.各向異性介質(zhì)共炮域高斯束偏移[J].石油地球物理勘探,2016,51(5):930-937

LIU Q,ZHANG M,LI Z C,et al.Common-shot domain Gaussian beam migration in anisotropic media[J].Oil Geophysical Prospecting,2016,51(5):930-937

[22]代福材,黃建平,李振春,等.角度域黏聲介質(zhì)高斯束疊前深度偏移方法[J].石油地球物理勘探,2017,52(2):283-293

DAI F C,HUANG J P,LI Z C,et al.Angle domain prestack Gaussian beam migration for visco-acoustic media[J].Oil Geophysical Prospecting,2017,52(2):283-293

[23]蔡杰雄,王華忠,王立歆.基于三維高斯束算子解析的方位反射角道集提取技術(shù)研究[J].石油物探,2016,55(1):76-83

CAI J X,WANG H Z,WANG L X.Azimuth-opening angle domain common-image gathers from 3D Gaussian beam migration[J].Geophysical Prospecting for Petroleum,2016,55(1):76-83

[24]WOODWARD M J.Wave-equation tomography[J].Geophysics,1992,57(1):15-26

[25]JEROEN J,SPETZLER J,SMEULDERS D,et al.Validation of first-order diffraction theory for the traveltimes and amplitudes of propagating waves[J].Geophysics,2006,71(6):167-177

[26]SPETZLER G,SNIEDER R.The fresnel volume and transmitted waves[J].Geophysics,2004,69(3):653-663

[27]XIE X B,YANG H.A migration velocity updating method based on the shot index common image gather and finite-frequency sensitivity kernel[J].Expanded Abstracts of 77thAnnual Internat SEG Mtg,2007:2767-2771

[28]XIE X B,YANG H.The finite-frequency sensitivity kernel for migration residual moveout and its applications in migration velocity analysis[J].Geophysics,2008,73(6):S241-S249

[29]XIE X B,YANG H.A wave-equation migration velocity analysis approach based on the finite-frequency sensitivity kernel[J].Expanded Abstracts of 78thAnnual Internat SEG Mtg,2008:3093-3097

[30]蔡杰雄,王華忠,陳進(jìn),等.基于高斯束傳播算子的成像域走時(shí)層析成像方法[J].地球物理學(xué)報(bào),2017,60(9):3539-3554

CAI J X,WANG H Z,CHEN J,et al.Traveltime tomography in the image domain based on the Gaussian-beam-propagator[J].Chinese Journal of Geophysics,2017,60(9):3539-3554