一種起伏地表?xiàng)l件下的照明補(bǔ)償方法

葉月明， 吳如山， 范國(guó)章， 丁仁偉

1 中國(guó)石油天然氣股份有限公司杭州地質(zhì)研究院, 杭州 310023 2 Modeling and Imaging Laboratory, University of California, Santa Cruz, U.S.A. 95064 3 山東科技大學(xué)地球科學(xué)與工程學(xué)院, 青島 266590

0 引言

隨著勘探的不斷深入，我國(guó)油氣戰(zhàn)略重點(diǎn)轉(zhuǎn)向了具有復(fù)雜地表和地下條件的雙復(fù)雜地區(qū).像中國(guó)西部山區(qū)這樣的起伏地表?xiàng)l件限制了采集孔徑范圍，造成了深層構(gòu)造照明不足、成像振幅不均衡，從而影響了成像質(zhì)量和后續(xù)的屬性分析.因此，開(kāi)展起伏地表?xiàng)l件下的照明補(bǔ)償方法研究對(duì)于提高起伏地表地區(qū)的成像品質(zhì)有著重要的現(xiàn)實(shí)意義.

傳統(tǒng)地震資料處理方法基于水平地表假設(shè)，不能直接用于地表起伏條件下的地震成像.起伏地表面的影響和深層弱照明問(wèn)題是地表起伏地區(qū)成像面臨的兩個(gè)主要問(wèn)題.目前工業(yè)生產(chǎn)用的各種靜校正方法(高程、折射波、層析等)都是基于地表一致性的假設(shè)條件，地震射線垂直入射和出射，但是實(shí)際受高程和近地表速度的變化射線并不是垂直入射和出射，從而引起靜校不“靜”，波場(chǎng)延拓的方法能夠解決這種靜校不“靜”的問(wèn)題.針對(duì)高程基準(zhǔn)面校正所帶來(lái)的問(wèn)題，Berryhill(1979)提出了波動(dòng)方程基準(zhǔn)面校正的方法，將已知任意形狀觀測(cè)面上的波場(chǎng)延拓到一個(gè)特定的基準(zhǔn)面上.零速層法(Beasley and Lynn, 1992)是在地表面與基準(zhǔn)面之間插入一個(gè)速度值很小的虛擬層，利用高程基準(zhǔn)面靜校正將野外數(shù)據(jù)校正到水平基準(zhǔn)面上，地震波在該層中橫向傳播可以忽略不計(jì)，即用波動(dòng)方程的方式抵消了高程校正的時(shí)移量.直接下延法可以避免偏移前的高程靜校正(Reshef，1991；Yang et al.，1999；葉月明等，2008).Al-Saleh等(2009)直接從起伏面開(kāi)始構(gòu)建橫向可變深度步長(zhǎng)的顯式波場(chǎng)延拓算子，相對(duì)于零速層技術(shù)有明顯的成像優(yōu)勢(shì).Shragge和Sava(2005)提出了基于保角變換的方法建立黎曼坐標(biāo)系統(tǒng)以適應(yīng)地表起伏面，再按照常規(guī)偏移方法成像.張新彥等(2017)利用數(shù)學(xué)變換手段將笛卡爾坐標(biāo)系的不規(guī)則模型映射到曲線坐標(biāo)系的規(guī)則模型，提出了起伏地形下的高精度反射波走時(shí)層析成像方法.隨著計(jì)算機(jī)技術(shù)的發(fā)展，相對(duì)耗時(shí)的起伏地表逆時(shí)偏移和最小二乘偏移方法也在被逐步應(yīng)用(曲英銘等，2015).侯爵等(2018)引入一種精度無(wú)損的處理起伏邊界的模型參數(shù)化方法，采用零延遲歸一化互相關(guān)成像條件實(shí)現(xiàn)了起伏地表?xiàng)l件下的彈性波場(chǎng)逆時(shí)偏移成像.周華敏等(2014)以地震波場(chǎng)在地下的照明強(qiáng)度作為迭代反演的預(yù)條件算子，加快了最小二乘偏移的收斂速度.近年來(lái)，基于高斯束的復(fù)雜地表偏移方法以其較高的偏移效率而被應(yīng)用，岳玉波等(2012)提出了復(fù)雜地表?xiàng)l件下保幅高斯束偏移，黃建平等(2016)研發(fā)了基于有效鄰域波場(chǎng)近似的起伏地表保幅高斯束偏移.有限的采集孔徑造成地下照明不足，為提高地下照明度，Wu等(2004)提出了水平地表情況下的局部角度域照明補(bǔ)償方法.基于傾斜疊加的局部角度域道集求取非常耗時(shí)(Xie and Wu, 2002; Xie et al., 2004)，而基于小波束分解的方法可以較快地得到局部角度域成像和照明補(bǔ)償因子(Cao and Wu，2009; 毛劍等，2010).陳生昌等(2007)、陳生昌和王漢闖(2010)提出基于平面波照明的補(bǔ)償方法，相對(duì)于局部角度域的照明補(bǔ)償方法具有較高的計(jì)算效率.

有限偏移孔徑造成的照明不足問(wèn)題在起伏地表?xiàng)l件下更為嚴(yán)重，為了能夠改善起伏地表地區(qū)的成像品質(zhì)，本文提出了一種基于小波束偏移算子的起伏地表區(qū)照明補(bǔ)償方法.在基準(zhǔn)面與起伏面之間插入虛擬層，通過(guò)波場(chǎng)逐步累加的延拓方法解決了起伏面對(duì)波場(chǎng)延拓的影響， 并引入空間濾波函數(shù)壓制了虛擬層內(nèi)偏移噪音.應(yīng)用了小波束波場(chǎng)延拓算子，其局部窗化的參考速度減小了背景速度擾動(dòng)量，提高了波場(chǎng)延拓算子的精度.局部指數(shù)標(biāo)架小波束波場(chǎng)分解求得了局部波數(shù)域成像矩陣和照明補(bǔ)償因子，最后在局部?jī)A角域?qū)崿F(xiàn)起伏地表?xiàng)l件下的照明補(bǔ)償.SEG起伏地表模型補(bǔ)償測(cè)試和振幅分析驗(yàn)證了本方法的有效性，深層照明度得到加強(qiáng)的同時(shí)也提高了成像振幅均衡度.

圖1 起伏地表?xiàng)l件下的觀測(cè)系統(tǒng)Fig.1 Acquisition system in area with irregular surface

1 方法原理

1.1 起伏地表?xiàng)l件下的波場(chǎng)延拓

常規(guī)偏移方法都是基于地表面水平的假設(shè)，無(wú)法在地表面起伏的條件下直接應(yīng)用.本文采用虛擬層內(nèi)波場(chǎng)零速傳播的方法來(lái)解決起伏面對(duì)波場(chǎng)延拓的影響.在起伏面與基準(zhǔn)面間插入虛擬層(圖1)，通過(guò)對(duì)虛擬層內(nèi)波場(chǎng)延拓的控制，使得虛擬層內(nèi)無(wú)波場(chǎng)傳播.為避免波場(chǎng)在起伏面處發(fā)生透射或反射，在虛擬層內(nèi)填充了近地表速度.以下行波場(chǎng)ud(x,z,w)為例，從基準(zhǔn)面開(kāi)始波場(chǎng)沿深度方向正向延拓，而在整個(gè)虛擬層內(nèi)是以零速傳播，當(dāng)遇到起伏地表面時(shí)，將對(duì)應(yīng)的震源波場(chǎng)累加到延拓波場(chǎng)中并繼續(xù)正向延拓.對(duì)于上行波場(chǎng)uu(x,z,w)，從基準(zhǔn)面開(kāi)始波場(chǎng)沿深度方向逆向延拓，在整個(gè)虛擬層內(nèi)同樣也是以零速傳播，當(dāng)遇到起伏面時(shí)，將對(duì)應(yīng)的地震記錄累加到延拓波場(chǎng)中，并繼續(xù)沿深度方向逆向延拓，如式(1)和式(2)所示：

(1)

(2)

其中，ω是頻率，(x,z)是空間位置，udc(x,z,ω)和uuc(x,z,ω)分別上一層延拓至(x,z)處的下行波場(chǎng)和上行波場(chǎng)，us(xs,zs,ω)和ur(xr,zr,ω)別是起伏地表面上對(duì)應(yīng)震源位置(xs,zs)與檢波點(diǎn)位置(xr,zr)的震源波場(chǎng)和地震記錄波場(chǎng)，在每個(gè)延拓深度層都會(huì)加權(quán)空間濾波因子f(x,z),以壓制虛擬層內(nèi)的噪音，如(3)式所示：

(3)

其中，udf(x,z,ω)和uuf(x,z,ω)分別是加權(quán)濾波因子后的下行和上行波場(chǎng)，空間濾波因子f(x,z)如式(4)所示：

(4)

1.2 基于局部余弦基的小波束波場(chǎng)延拓

傳統(tǒng)單程波偏移算子中的背景速度一般為偏移孔徑內(nèi)的平均速度或最小速度，在地表高程起伏較大的地區(qū)，虛擬層內(nèi)的速度與實(shí)際地層間的速度可能會(huì)存在較大差異，導(dǎo)致真實(shí)速度與背景速度間擾動(dòng)量偏大，影響算子精度.小波束波場(chǎng)延拓算子在偏移孔徑范圍內(nèi)劃分多個(gè)窗，背景速度為每個(gè)窗內(nèi)的平均速度，局部窗化的參考速度降低了速度擾動(dòng)量，也提高了波場(chǎng)延拓算子精度.本文應(yīng)用了基于局部余弦基(Local Cosine Basis，簡(jiǎn)稱LCB)的小波束波場(chǎng)延拓算子(Wu et al., 2008)，該算子得益于快速LCB變換而具有較高的計(jì)算效率.

在頻率-空間域，標(biāo)量聲波方程可以表示為

(5)

其中，v(x,z)是速度，u(x,z,ω)代表頻率-空間域的波場(chǎng)，為了方便起見(jiàn)，后文用u(x,z)來(lái)代替頻率-空間域的波場(chǎng).

將深度z上的波場(chǎng)沿空間x方向窗化后，按照式(6)分解波場(chǎng)如下：

(6)

(7)

(8)

am n(x)=e±iAnΔzbm n(x)，

(9)

其中，am n(x)是小波束bm n(x)經(jīng)過(guò)介質(zhì)傳播后的波束，An是平方根算子.經(jīng)過(guò)介質(zhì)的傳播后am n(x)不再是一個(gè)小波束，應(yīng)用同樣的基函數(shù)對(duì)am n(x)再進(jìn)行小波束分解表示為

(10)

其中，傳播算子矩陣在小波束域可以表示為

(11)

那么，延拓至深度z+Δz的波場(chǎng)在小波束域可以表示為

(12)

Pjl,m n是小波束傳播算子矩陣，決定著小波束的傳播和交互偶合，通過(guò)式(12)可以使波場(chǎng)在小波域從深度z延拓到z+Δz層.波場(chǎng)重建是波場(chǎng)分解的反過(guò)程，小波域的波場(chǎng)在頻率空間域重建后可以表示為式(13)：

(13)

小波束波場(chǎng)延拓算子的特點(diǎn)在于局部窗化的參考速度，將每一炮的偏移孔徑范圍分成若干個(gè)窗，在每個(gè)窗內(nèi)計(jì)算出該窗內(nèi)的背景速度和每個(gè)網(wǎng)格點(diǎn)上的擾動(dòng)速度，能夠確定波場(chǎng)在空間局部位置處的波數(shù)，每個(gè)小窗內(nèi)的波場(chǎng)就是一個(gè)小波束.實(shí)際資料處理中的每個(gè)小波束寬度取值在120 m左右，就能夠滿足小波束波場(chǎng)延拓算子的精度需要.

1.3 局部?jī)A角域成像與照明補(bǔ)償

本文的照明補(bǔ)償是在角度域完成，通過(guò)波場(chǎng)分解將頻率空間域波場(chǎng)變換到局部?jī)A角域.在波場(chǎng)延拓過(guò)程中采用的是基于局部余弦基的小波束波場(chǎng)延拓算子，其中局部余弦基函數(shù)(式(7))還可以表示為式(14)，兩個(gè)指數(shù)項(xiàng)表明該類(lèi)小波束有兩個(gè)對(duì)稱的波瓣，即沿深度軸對(duì)稱的兩個(gè)方向傳播，缺乏單一定義的方向性，不利于照明補(bǔ)償.

(14)

事實(shí)上，如式(14)所示的基函數(shù)都是由局部指數(shù)標(biāo)架線性組合而來(lái)的，局部指數(shù)標(biāo)架定義為

(15)

(16)

(17)

其中，θs和θg分別表示地下局部位置入射角和散射角，地層的局部?jī)A角θn更能反映出地下的局部真實(shí)情況，各角度之間的關(guān)系如圖2所示，傾角與地層反射角之間的關(guān)系如式(18)所示：

θn=(θs+θg)/2.

(18)

圖2 反射體傾角和反射角間的關(guān)系Fig.2 A sketch of dip angle and reflection angle

通過(guò)式(18)將局部角度域成像L(x,z,θs,θg)轉(zhuǎn)換到局部?jī)A角域L(x,z,θn)，局部角度域成像能夠反映出地下局部反射點(diǎn)在不同方向的入射波場(chǎng)和反射波場(chǎng)對(duì)成像的貢獻(xiàn)，稱之為局部成像矩陣(Local Image Matrix, 簡(jiǎn)稱LIM).局部散射矩陣(Local Scattering Matrix，簡(jiǎn)稱LSM)則是反映了地下固有特征，描述振幅與散射能量間的一定比例關(guān)系.由于采集孔徑限制以及傳播路徑等影響，LIM是被扭曲了的LSM，照明補(bǔ)償?shù)哪康木褪菫榱俗孡IM更加逼近LSM，從而最終反映出地下的真實(shí)情況(Wu et al., 2004).與水平地表情況類(lèi)似，局部波數(shù)域照明補(bǔ)償因子可以表示為(Cao and Wu, 2009)

(19)

其中，GF是格林函數(shù)正演過(guò)程，A(xg,xs)是空間采集孔徑.根據(jù)式(17)和式(18)，將局部波數(shù)域照明補(bǔ)償因子變換到局部?jī)A角域F(x,z,θn)，于是，局部?jī)A角域照明補(bǔ)償計(jì)算如式(20)：

(20)

其中，θmin和θmax分別是最小和最大地層傾角，ε是阻尼因子，以保證計(jì)算的穩(wěn)定性，目前本方法的阻尼因子是基于經(jīng)驗(yàn)值和測(cè)試后選定，選取的基本原則是比傾角域補(bǔ)償因子小兩個(gè)數(shù)量級(jí)，這樣既保證了計(jì)算的穩(wěn)定性，也不改變相對(duì)關(guān)系.

2 模型測(cè)試

為了驗(yàn)證本文方法的有效性，對(duì)加拿大逆掩斷層模型(Gray and Marfurt, 1995)進(jìn)行試算分析.該模型地表高程較大(圖3a)，最大高程近2 km，而且地表起伏劇烈，局部范圍2.5 km內(nèi)高程差能夠達(dá)到1.5 km，并且發(fā)育多套高陡斷層和逆掩構(gòu)造(圖3b)，起伏地表面和上覆復(fù)雜構(gòu)造影響了深層兩套傾斜地層的照明度，圖3c是原始單炮記錄，地表起伏造成的同相軸扭曲現(xiàn)象嚴(yán)重.模擬的激發(fā)方式是中間放炮，共277炮、480道接收、采樣率2000、4 ms采樣、道間距15 m、最大炮檢距3600 m、最小炮檢局為15 m、炮間距為90 m.

波場(chǎng)在介質(zhì)中是以小波束的形式傳播，圖4a和4b分別是15 Hz和40 Hz情況下的小波束傳播，二者都是下行波場(chǎng)，以垂直地表的方向傳播.對(duì)比可以看到，較高頻率的小波束在地下傳播的角度擴(kuò)散相對(duì)更大一些.圖5a、圖5b和圖5c分別是橫向位置相同、不同深度三個(gè)位置的局部角度域成像矩陣，橫、縱軸分別代表了入射角和散射角.從圖中可以看出，深度8 km位置處的局部成像角度貢獻(xiàn)更為均勻(圖5a)，局部成像矩陣中心在入射角和散射角都是0°附近.深度6 km位置處，局部成像矩陣中心在入射角-30°，散射角0°附近，而在深度4 km處，局部成像矩陣成橢圓形分布.局部?jī)A角域能夠更加客觀地反演局部地層傾向，因此在局部?jī)A角域進(jìn)行照明

圖3 (a) 起伏地表高程； (b) 起伏地表模型速度場(chǎng)； (c) 單炮記錄Fig.3 (a) Elevation of irregular surface; (b) Velocity model with irregular surface; (c) Single shot gather

圖4 (a) 15 Hz小波束在起伏地表模型中的傳播； (b) 40 Hz小波束在起伏地表模型中的傳播Fig.4 (a) Beamlet propagation in irregular surface model with 15 Hz frequency; (b) Beamlet propagation in irregular surface model with 40 Hz frequency

補(bǔ)償能夠根據(jù)不同角度的照明情況來(lái)針對(duì)性地補(bǔ)償.圖6a、6b、6c和圖6d分別是局部?jī)A角為-30°、0°、30°和傳統(tǒng)自動(dòng)增益控制(Auto Gain Control，簡(jiǎn)稱AGC)的補(bǔ)償因子，不同局部?jī)A角域補(bǔ)償因子也表現(xiàn)出不同補(bǔ)償?shù)姆较蛐?圖7a、7b和7c分別是局部?jī)A角為-30°、0°、30°情況下的局部?jī)A角域成像，圖7d、7e和7f分別是局部?jī)A角為-30°、0°、30°情況下的局部?jī)A角域成像照明補(bǔ)償后結(jié)果，從圖中可以明顯看出，照明補(bǔ)償后，不同角度的成像都有改善，尤其是解決了由于深層照明不足對(duì)構(gòu)造成像帶來(lái)的影響.

相對(duì)于傳統(tǒng)的分步傅里葉疊前深度偏移(圖8a)，LCB小波束疊前深度偏移算子局部窗化的參考速度降低了背景速度擾動(dòng)量，尤其是在起伏地表面附近的淺層成像效果有了顯著的改善(圖8b).將所有傾角域補(bǔ)償后的像疊加得到最終局部?jī)A角域照明補(bǔ)償結(jié)果，如圖8d所示，與補(bǔ)償前(圖8b)相比，深層構(gòu)造成像振幅的能量有了顯著的提升，而且補(bǔ)償效果上更優(yōu)于傳統(tǒng)AGC增益處理(圖8c).為了進(jìn)一步分析照明補(bǔ)償后的振幅，沿圖9a箭頭所指層位，拾取了照明補(bǔ)償前(圖8b)、AGC補(bǔ)償后(8c)和照明補(bǔ)償后(圖8d)的成像振幅如圖9b所示，AGC增益處理只是增強(qiáng)了深層振幅值，在層位上拾取的振幅值和照明補(bǔ)償前的LCB偏移結(jié)果相近，不同反射角對(duì)成像貢獻(xiàn)的相對(duì)關(guān)系并沒(méi)有變，并沒(méi)有達(dá)到照明補(bǔ)償?shù)哪康?，?yīng)用本文方法照明補(bǔ)償后，不同反射角的成像貢獻(xiàn)更加的均衡，如圖9b矩形框內(nèi)振幅曲線可以看出，照明補(bǔ)償后的振幅值局部范圍內(nèi)最高值和最低值的差最小，相對(duì)更加均衡.在圖9a中選了6個(gè)點(diǎn)進(jìn)行局部?jī)A角域照明補(bǔ)償前后的振幅隨角度(AVA)分析，圖10a和圖10b分別是這6個(gè)點(diǎn)補(bǔ)償前后的局部?jī)A角域成像振幅值，照明補(bǔ)償后成像振幅值在各個(gè)角度上的均衡性相對(duì)于補(bǔ)償前有著顯著的改善.

圖5 局部角度域成像矩陣. (a) x=5 km, z=8 km; (b) x=5 km, z=6 km; (c) x=5 km, z=4 kmFig.5 Local image matrix with location

圖6 (a) -30°傾角照明補(bǔ)償因子; (b) 0°傾角照明補(bǔ)償因子; (c) 30°傾角照明補(bǔ)償因子; (d) 常規(guī)AGC因子Fig.6 (a) Illumination compensation factor of -30 degree dip angle; (b) Illumination compensation factor of 0 degree dip angle; (c) Illumination compensation factor of 30 degree dip angle; (d) Conventional AGC factor

圖7 (a) -30°傾角成像; (b) 0°傾角成像; (c) 30°傾角成像; (d) -30°傾角成像照明補(bǔ)償后結(jié)果; (e) 0°傾角成像照明補(bǔ)償后結(jié)果; (f) 30°傾角成像照明補(bǔ)償后結(jié)果Fig.7 (a) Image of -30 degree dip angle; (b) Image of 0 degree dip angle; (c) Image of 30 degree dip angle; (d) Image of -30 degree dip angle after compensation; (e) Image of 0 degree dip angle after compensation; (f) Image of 30 degree dip angle after compensation

圖8 (a) 分步傅里葉疊前深度偏移; (b) LCB小波束疊前深度偏移; (c) LCB小波束疊前深度偏移后AGC; (d) LCB小波束疊前深度偏移后照明補(bǔ)償Fig.8 (a) Prestack depth migration based on split-step Fourier; (b) LCB beamlets pre-stack depth migration; (c) AGC after LCB beamlets pre-stack depth migration; (d) Illumination compensation after LCB beamlets pre-stack depth migration

圖9 (a) 局部速模型; (b) LCB小波束偏移照明補(bǔ)償前后層位拾取振幅值對(duì)比Fig.9 (a) Local velocity model; (b) Before and after illumination compensation of LCB beamlets migration

圖10 (a) 照明補(bǔ)償前不同位置局部?jī)A角拾取的振幅; (b) 照明補(bǔ)償后不同位置局部?jī)A角拾取的振幅Fig.10 (a) Local dip angle amplitude of different location before illumination compensation; (b)Local dip angle amplitude of different location after illumination compensation

3 結(jié)論

針對(duì)起伏地表?xiàng)l件下的有限采集孔徑和深層弱照明問(wèn)題，本文提出了一種基于小波束的起伏地表照明補(bǔ)償方法，解決了起伏面對(duì)波場(chǎng)延拓的影響，也補(bǔ)償由于采集孔徑不足造成的深層成像能量不均衡問(wèn)題，為后續(xù)的AVA分析提供更為可靠的偏移數(shù)據(jù).在偏移過(guò)程中，基于波場(chǎng)逐步累加的思路，使得波場(chǎng)延拓能夠適應(yīng)起伏的地表?xiàng)l件，引入LCB小波束偏移算子，在每一層上將背景速度窗化，降低了背景速度擾動(dòng)量，提高了成像精度，克服了傳統(tǒng)單程波算子全局參考速度導(dǎo)致的起伏面附近強(qiáng)成像噪音問(wèn)題.局部?jī)A角域的照明補(bǔ)償不但提升了深層照明度，而且成像振幅均衡度也得到了有效改善，有助于后期的屬性分析和反演.對(duì)復(fù)雜的SEG起伏地表模型測(cè)試說(shuō)明了該方法的正確性.在近地速度精度較高情況下，本方法能夠適用于地表高程差較大和地表橫向起伏變化劇烈地區(qū)的成像，并有望在起伏地表地區(qū)的實(shí)際資料處理中取得較好的效果，為山地、盆山邊界等復(fù)雜地表地區(qū)地震勘探發(fā)揮重要作用.