葉月明， 吳如山， 范國章， 丁仁偉

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

0 引言

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

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

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

圖1 起伏地表條件下的觀測系統(tǒng)Fig.1 Acquisition system in area with irregular surface

1 方法原理

1.1 起伏地表條件下的波場延拓

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

(1)

(2)

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

(3)

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

(4)

1.2 基于局部余弦基的小波束波場延拓

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

在頻率-空間域，標量聲波方程可以表示為

(5)

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

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

(6)

(7)

(8)

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

(9)

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

(10)

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

(11)

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

(12)

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

(13)

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

1.3 局部傾角域成像與照明補償

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

(14)

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

(15)

(16)

(17)

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

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

(18)

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

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

(19)

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

(20)

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

2 模型測試

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

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

圖3 (a) 起伏地表高程； (b) 起伏地表模型速度場； (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

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

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

圖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) 0°傾角照明補償因子; (c) 30°傾角照明補償因子; (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°傾角成像照明補償后結果; (e) 0°傾角成像照明補償后結果; (f) 30°傾角成像照明補償后結果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小波束疊前深度偏移后照明補償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小波束偏移照明補償前后層位拾取振幅值對比Fig.9 (a) Local velocity model; (b) Before and after illumination compensation of LCB beamlets migration

圖10 (a) 照明補償前不同位置局部傾角拾取的振幅; (b) 照明補償后不同位置局部傾角拾取的振幅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 結論

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