高效逆時偏移角度道集生成方法研究

吳成梁,周 陽,胡江濤,王華忠,孫 榮

(1.波現(xiàn)象與反演成像研究組(WPI),同濟大學海洋與地球科學學院,上海200092;2.油氣藏與開發(fā)工程國家重點實驗室,成都理工大學,四川成都610059;3.中國石油勘探開發(fā)研究院,北京100083)

角度域共成像點道集(ADCIGs)既包含背景速度的影響,也包含方位角度反射系數(shù)的信息,因此可用于偏移速度分析和更新速度模型[1-3]、照明分析[4]、提取AVA信息[5],還可用于介質(zhì)各向異性分析和提高構(gòu)造成像質(zhì)量等。

ADCIGs可以由射線類疊前偏移方法生成(比如克希霍夫偏移方法[6-8]或者Beam類偏移方法[9-10]),也可由波動方程類偏移方法(比如單程波方程偏移[11-12]或者逆時偏移(RTM)[13-15])生成。與克?；舴蚱啤eam偏移和單程波方程偏移相比,逆時偏移能夠精確地處理多波至現(xiàn)象和突破高陡傾角成像限制,成為了復雜區(qū)域的首選成像方法。所以,逆時偏移的角度道集生成方法研究非常重要[16-17]。

目前生成逆時偏移角度道集的方法主要有兩類:一類是成像后方法,包括擴展成像條件方法[18-21]和逆散射成像條件方法[22-24];另一類是成像前方法,包括波矢量方向估計方法[25-34]和波場局部平面波分解方法[17,35-39]。

擴展成像條件方法首先生成地下局部偏移距道集,然后將其轉(zhuǎn)換為角度域成像道集,該轉(zhuǎn)換過程需要進行高維傅里葉變換和局部傾斜疊加計算,存在內(nèi)存和I/O問題。在進行局部傾斜疊加時,時窗大小對生成道集的質(zhì)量影響很大,選擇的局部時窗小時,角度道集噪聲少,但分辨率低;選擇的時窗大時分辨率高,但是傾斜疊加產(chǎn)生的假象多,從而降低了角度域成像道集的質(zhì)量,影響AVA特征提取。

逆散射成像條件方法是利用方位角度濾波器對空間求導后的波場進行相關(guān)成像,得到單方位角度的成像結(jié)果,疊加后產(chǎn)生角度道集。盡管該方法可以自動地壓制低頻噪聲,并且增加的計算量不大,但是其成像結(jié)果的振幅(包括子波相位)保真性及波場方向計算的準確性等方面還存在問題。

波矢量方向估計方法是指在波場傳播過程中計算波矢量方向,利用估計的波矢量方向判斷波場的傳播方向,并利用角度域成像條件生成角度道集。由于坡印廷矢量方法計算效率高,且具有較高的角度分辨率,因此被廣泛應用于逆時偏移角度道集的生成。在角度道集估計過程中,DICKENS等[27]利用炮點端波場和檢波點端波場的傳播方向提取角度道集,并分析了角度道集的振幅特征。但是當檢波點端波場比較復雜時,計算的坡印廷矢量通常會出現(xiàn)不穩(wěn)定情況,因此YOON等[40]和VYAS等[28]利用坡印廷矢量計算的炮點端波場傳播方向和地下反射界面傾角來估計反射角;ZHAO等[41]利用偏移剖面預測反射層法向向量,并將其與穩(wěn)定的炮點端方向波場結(jié)合起來計算角度道集。王保利等[32]采用一階聲波波動方程求解坡印廷矢量,減少了坡印廷矢量計算時間,并利用高斯采樣函數(shù)解決了角度不等間隔問題。雖然坡印廷矢量方法的計算效率比較高,但是該方法無法解決波前重疊的問題[42]。VYAS等[29]和TANG[43]先將波場進行方向分解,然后利用分解后的波場來計算坡印廷矢量。在波矢量方向估計方法中,最重要的是波傳播方向的計算,ZHANG等[30-31]利用P波的位移極化矢量表示波的傳播方向;ZHANG[44]采用光學流計算波場傳播方向。

局部平面波分解方法是先將炮點端波場和檢波點端波場分解為各自的局部平面波分量,然后應用角度域成像條件來提取地下局部角度域成像道集,通常在時間-空間域或頻率-空間域中實現(xiàn)。XIE等[35]利用局部傾斜疊加方法提取偏移波場中的角度域信息。YAN等[37]提出在時間-空間域利用局部慢度方向提取炮點端和檢波點端的波場方向來構(gòu)建角度道集的方法,這種方法可以解決波前重疊的問題,但是存在角度分辨率低的問題,而且需要儲存足夠的時間片來計算波傳播方向。XU等[17]給出了一種在頻率-波數(shù)域直接構(gòu)建角度道集的方法,并采用抗泄露傅里葉變換[45-46]解決了角度分辨率問題,但是該方法需要在逆時偏移中將整個波場全部儲存下來進行關(guān)于時間維的傅里葉變換,計算代價很大。HU等[47]提出基于解析時間波場外推與局部波場分解的逆時偏移角度道集生成方法,解析時間波場只包含正頻率的波場,不需要進行時間維的傅里葉變換即可確定波場的傳播方向,解決了儲存波場的問題,但是該方法需要對震源和波場記錄實施Hilbert變換,分別進行正反傳計算,導致一次逆時偏移過程至少需要兩次正傳和兩次反傳,計算代價還是較大。

本文結(jié)合波矢量方向估計方法和局部平面波分解方法的優(yōu)點,提出一種混合構(gòu)建角度道集的方法。首先在常規(guī)逆時偏移基礎(chǔ)上利用局部平面波分解類方法提取局部平面波,然后利用波矢量方向估計方法計算的波場方向來確定最終的波前傳播方向,最后采用角度域成像條件生成角度道集。同時,利用結(jié)構(gòu)張量特征值快速判斷波前是否交叉,只在波前交叉的情況下使用本文提出的混合方法,而在單一的波前面采用高效的波矢量估計方法,進一步提高了計算效率。

1 波矢量方向估計方法

波矢量方向估計方法是指利用波場的振幅梯度和相位梯度等信息來計算地下任意時刻、任意點的波傳播方向,其最主要的難點是波前方向矢量的計算。計算波前方向矢量的方法通常有:坡印廷矢量方法、極化矢量方法、瞬時波數(shù)矢量方法以及光學流方法等。

坡印廷矢量原本是指電磁場中的能流密度矢量,用于計算單位時間內(nèi)穿過垂直于矢量方向的單位面積的電磁場能量。2006年,YOON等[26]首次將坡印廷矢量引入到逆時偏移中,用于確定波的傳播方向。坡印廷矢量可由下式計算:

(1)

(2)

其中,ps和pr分別代表炮點端波場和檢波點端波場在地下x點處t時刻的波場傳播方向;θ為地下反射角度,θ=arccos(ps-pr)。反射角的計算方法有多種,如DICKENS等[27]利用炮點端波場和檢波點端波場傳播方向計算反射角;YOON等[40]和VYAS等[28]利用炮點端波場傳播方向和反射截面傾角計算反射角。

坡印廷矢量方法只需要計算波場的空間導數(shù)和時間導數(shù),計算效率高,易實現(xiàn),具有較高的角度分辨率,但易受噪聲影響。特別是通過地表數(shù)據(jù)逆時延拓得到的檢波點端波場,由于只有地表的上邊界數(shù)據(jù)作為波動方程的邊值條件,在進行逆時外推時,波前無法干涉閉合,加上存在噪聲,因而坡印廷矢量方法計算的方向矢量連續(xù)性差,穩(wěn)定性和準確性難以保證。另外,坡印廷矢量方法假設(shè)地下任一點任一時刻只有一個傳播方向,當波場出現(xiàn)交叉或存在界面時,波前方向的計算不準確,不能得到正確的角度道集。

2 局部平面波分解方法

局部波場分解方法主要是利用局部傅里葉變換或局部傾斜疊加等方法將局部波場變換到其它域(空間-時間域、頻率-波數(shù)域或時間-波數(shù)域)中實現(xiàn)角度分解。在各種局部波場分解方法中,比較高效的是HU等[47]提出的基于解析時間波場外推與局部波場分解的角度道集生成方法。由于解析時間波場中只包含正頻率的波場信息,該方法只需要在時間域逆時外推中進行空間維的傅里葉變換即可實現(xiàn)波場的角度分解。

假設(shè)波場方向分解借助傅里葉變換來實現(xiàn),將波場實施時間維和空間維的傅里葉變換后,根據(jù)波場的頻率和波數(shù)即可定義波場的上下行波[48],從而實現(xiàn)波場的方向分解。在時間域中實現(xiàn)這種波場分離時,需要存儲地下整個波場,進行時間維的傅里葉變換,其I/O量和計算量都很大。如果能夠在外推過程中固定波場的頻率,那么波場的方向分解只需要利用波場的空間波數(shù)信息即可,而空間波數(shù)卻很容易獲得。在信號處理中有一類時間信號叫做解析時間信號,其實部為信號本身,而虛部為信號的Hilbert變換,它只包含正頻率的信號。HU等[47]將其引入到波場計算中,提出解析時間波場的概念。解析時間波場是復波場,其實部可以通過波動方程的傳播過程得到,其虛部可通過對地下波場進行Hilbert變換得到。但是由于Hilbert變換是關(guān)于時間方向的卷積,因此需要存儲整個波場并進行卷積,不利于計算。Hilbert變換為一個90°相移算子,借助于波動方程的震源項與波場具有線性關(guān)系[49],因此可以將波場的Hilbert變換轉(zhuǎn)化為震源項的Hilbert變換,利用新的‘震源’進行傳播,從而可以獲得解析時間波場的虛部。

以炮點端為例,解析時間波場的外推過程可表示為:

(3)

由于在局部的空間窗內(nèi)僅包含幾個局部的平面波,因此可以采用局部波場的稀疏分解方法來高效提取波場的局部角度分量[47]。在實現(xiàn)局部平面波稀疏分解時,定義角度分解的誤差泛函為:

(4)

其中,uwin(t,x)代表局部解析空間波場,uwin(t,x,θ)代表提取的局部平面波,θ為局部平面波的傳播角度。實現(xiàn)波場稀疏角度分解的過程為:首先對局部解析時間波場實施空間維的傅里葉變換,并投影到波數(shù)域,然后在波數(shù)域中將同一傳播角度的波數(shù)進行反傅里葉變換,實現(xiàn)局部平面波的分解。通過使目標泛函(4)最小化,可將局部窗內(nèi)的全部平面波提取出來。

對稀疏角度分解后提取的炮點端波場和檢波點端波場的角度分量應用角度域成像條件,提取成像結(jié)果,得到角度道集。采用解析時間波場外推與局部波場稀疏分解的方法不需要進行關(guān)于時間維的傅里葉變換,避免了整個外推波場的存儲,因此在逆時偏移中可以很方便地實現(xiàn)。同時該方法可以處理復雜的波前交叉情況,得到的角度道集物理意義更加明確,能夠比較好地反映地下界面的AVA特征。但是該方法需要對解析波場虛部進行波場外推,額外增加的波場外推計算量非常大,特別是在各向異性介質(zhì)中,實現(xiàn)波場外推需要付出巨大的計算代價。

3 混合方法

波矢量方向估計方法可以高效地實現(xiàn),但是在檢波點端存在方向估計不穩(wěn)定的問題,并且無法處理波場交叉的情況。利用局部平面波分解方法可以準確地提取波傳播方向,即使在波場交叉、檢波點波場比較復雜的情況下,仍然能夠很好地提取波傳播方向。但是該方法需要進行高維傅里葉變換,需要大量內(nèi)存和I/O,效率低下,即使采用解析時間波場,也需要對震源和波場記錄實施Hilbert變換,分別進行波場外推,一次逆時偏移至少需要做4次正演,計算量還是比較大。本文結(jié)合局部平面波場分解的準確性和波矢量方向估計方法的高效性,提出了一種混合生成RTM-ADCIGs的方法,其實現(xiàn)過程如下:

1) 在正演和逆時波場外推的每個時間片上,選擇局部波場進行空間維傅里葉變換,將局部波場從空間域變換到波數(shù)域。

2) 在波數(shù)域采用HU等[47]提出的稀疏波場分解方法,按照能量譜大小提取局部平面波(由于在波數(shù)域中能量譜是對稱的,因此只需提取正波數(shù)方向),此時可以得到該平面波的兩個波前方向φ或-φ(φ與-φ相差180°)。由于缺乏時間維的信息,因此無法確定該波前傳播方向是φ還是-φ,但是提取的局部波前矢量具有穩(wěn)健性,其方向是正確的。

3) 利用坡印廷矢量方法計算局部波場矢量方向,此時計算的波場方向在波場交叉和復雜情況下是不準確的。

4) 在選取的局部平面波空間窗內(nèi),沿著垂直于該平面波的波前矢量方向選擇帶狀的空間窗。在這個帶狀的空間窗內(nèi)計算的波矢量方向和平面波的波前方向矢量大體一致。因此,在選擇的帶狀空間窗內(nèi)統(tǒng)計波矢量方向落在計算的兩個波前方向φ或-φ的個數(shù),統(tǒng)計個數(shù)最多的方向就代表了波前的方向(可選擇一定的角度統(tǒng)計范圍),利用統(tǒng)計得到的波前方向來約束第2步中提取的波矢量方向的符號,從而得到該平面波準確的波前擴散方向。空間窗的寬度以能覆蓋有效波矢量方向為準則,一般為局部平面波寬度的1/3～1/2除以局部平面波個數(shù)。角度統(tǒng)計范圍根據(jù)具體的平面波個數(shù)限定,本文選用的范圍為30°除以局部平面波個數(shù)。

5) 利用提取的炮點端局部平面波和檢波點端局部平面波,應用角度域成像條件提取該時刻的ADCIGs。

6) 重復上述2)～5)步可提取下一個局部平面波,直到空間窗內(nèi)沒有平面波為止。

下面通過一個簡單的模型來論述角度道集的生成過程。如圖1a所示,該模型從上往下速度值依次為3000,3200,3400m/s。對應1s時刻的波場快照如圖1b所示。選取該時刻某一局部空間波場如圖2a 所示。采用HU等[47]提出的解析時間波場方法計算的波場虛部如圖2b所示。將圖2a所示的局部空間波場進行傅里葉變換,投影到波數(shù)域,其能量譜如圖3a所示。將圖2a和圖2b組合形成的解析時間波場變換到波數(shù)域,其能量譜如圖3b所示。由于解析時間波場中只包含正頻率的信息,兩個局部平面波只含有兩個波數(shù)域能量團,依次提取對應的最大能量團即可獲得該空間窗內(nèi)的局部平面波方向。而采用常規(guī)波場得到的波數(shù)域能量譜中兩個局部平面波包含4個能量團,無法直接根據(jù)波數(shù)域能量團確定局部平面波的傳播方向。分析發(fā)現(xiàn),圖3a中的波數(shù)域能量團是關(guān)于中心點對稱的,同一局部平面波的兩個波數(shù)能量團的方向相反,其夾角為180°,如圖3a中紅色箭頭和藍色箭頭所示。利用對稱性,只在正波數(shù)區(qū)域內(nèi)采用稀疏局部平面波分解方法提取最大波數(shù)能量團對應的局部平面波,如圖4a所示。由于缺乏時間維信息,一個局部平面波在波數(shù)域中對應兩個能量團,無法確定該局部平面波的傳播方向,只能得到兩個可能的傳播方向,如圖4a中紅色實箭頭和黑色虛箭頭所示。重復上述過程提取該區(qū)域內(nèi)的第二個局部平面波如圖4b所示,兩個可能的傳播方向分別為圖4b 中藍色實箭頭和黑色虛箭頭所示。圖5為采用坡印廷矢量方法提取的局部波場方向,可以看到,在波場交叉的區(qū)域,波矢量方向的估計是錯誤的,如圖5 中橙色圓圈所示。根據(jù)圖4a中提取的局部平面波確定帶狀的空間窗如圖5中紅色線框所示,根據(jù)圖4b 中提取的局部平面波確定帶狀的空間窗如圖5中藍色線框所示。在紅色線框內(nèi)統(tǒng)計波矢量方向,落在圖4a紅色實箭頭和黑色虛箭頭所指方向的個數(shù)多的方向即為該局部平面波的傳播方向(圖4a中紅色實箭頭所指的方向)。同樣,統(tǒng)計藍色方框內(nèi)波矢量方向落在圖4b中藍色實箭頭和黑色虛箭頭所指方向的個數(shù),根據(jù)統(tǒng)計結(jié)果判斷圖4b中藍色實箭頭方向為其局部平面波的傳播方向。

圖1 簡單模型(a)及其1s時刻波場快照(b)

圖2 局部的空間波場a 實波場; b 解析時間波場虛部

圖3 波數(shù)域能量譜a 常規(guī)波場; b 解析時間波場

圖4 提取的局部平面波a 第一個局部平面波分量(紅色實箭頭為提取的波前矢量方向); b 第二個局部平面波分量(藍色實箭頭為提取的波前矢量方向)

接下來采用圖6所示的層狀模型測試本文混合方法提取波場方向的有效性。正演采用的觀測系統(tǒng)參數(shù)如下:震源位置為(3000m,0),震源為30Hz的雷克子波,檢波點間距10m,偏移距范圍為(-3000m,3000m),時間采樣點數(shù)nt=5000,時間采樣間隔Δt=0.5ms。檢波點端波場在0.505s時刻采用不同方法獲得的波場方向如圖7所示,其中圖7a 為采用坡印廷矢量方法計算的波場方向(為了與局部平面波顯示一致,只選取了局部空間窗中心點的矢量方向代表該空間窗內(nèi)波場方向);圖7b為對常規(guī)波場進行空間維傅里葉變換,根據(jù)波數(shù)域能量譜大小獲得的局部平面波波場方向;圖7c為在圖7b基礎(chǔ)上施加坡印廷矢量約束后計算的波場方向(即本文提出的混合方法);圖7d為采用解析時間波場方法獲得的波場方向。對比圖7a至圖7d可見,只采用坡印廷矢量方法會在波場分解面處出現(xiàn)方向計算不準確的情況(圖7a);而對常規(guī)波場進行空間維傅里葉變換方法由于缺乏時間維的信息,也無法得到準確的波場方向(圖7b);施加坡印廷矢量方向約束后,可以獲得準確的波場方向(圖7c),其精度也與解析時間波場方法提取的波場方向相當(圖7d)。

圖5 坡印廷矢量方法計算的局部波場方向

圖6 層狀模型

圖7 檢波點端波場0.505s時刻采用不同方法計算的波場方向a 坡印廷矢量方法; b 對常規(guī)波場只采用空間維傅里葉變換的方法; c 在空間維傅里葉變換基礎(chǔ)上施加坡印廷矢量約束的方法; d 解析時間波場方法

混合方法的角度道集估計過程還可以擴展到不同類型的方法組合,比如局部平面波的提取可以選擇傅里葉變換、局部傾斜疊加和τ-p變換等方法,波矢量方向的計算方法也可擴展到極化矢量方法、瞬時波數(shù)矢量方法以及光學流方法等。同時,該方法很容易擴展到各向異性介質(zhì)中,與解析時間波場方法相比,增加的計算量更少。

4 改進的混合方法

分析地下介質(zhì)不同時刻的波場快照可知,不同時刻的波場復雜程度不同,在波場簡單情況下,大部分區(qū)域由單一的波前面組成,只在局部區(qū)域存在波前交叉情況;隨著波場復雜程度的增加,地下波場的交叉情況也變得越來越復雜,即使如此,在局部空間范圍內(nèi)也存在單一的波前面。因此,可以根據(jù)波前的復雜情況,選擇合適的角度道集提取方法,在波前簡單的情況下采用快速的波矢量估計方法,而在波前出現(xiàn)重疊、交叉的情況下采用局部波場分解方法。

考慮到波場在不同時空范圍內(nèi)的特征,本文借助于結(jié)構(gòu)張量特征值,快速判斷波場交叉的情況,從而只在波場能量交叉的情況下使用混合方法,估計一個或多個局部波場的方向,而在其它區(qū)域采用快速的波矢量估計方法提取波場方向,以進一步提高角度道集的計算效率。

結(jié)構(gòu)張量,又稱二階矩矩陣,在圖像處理中非常流行[50-51],也常用于地震數(shù)據(jù)分析。BAKKER[52]對結(jié)構(gòu)張量在地震數(shù)據(jù)濾波方面的應用做了詳細的描述,HALE[53]利用結(jié)構(gòu)張量實現(xiàn)了地震剖面沿構(gòu)造方向的平滑濾波處理,增強了地下結(jié)構(gòu)特征,保持了斷層結(jié)構(gòu)特征。

以二維為例,給定一波場u,其結(jié)構(gòu)張量矩陣可表示為:

(5)

(6)

其中,λe和λv分別是結(jié)構(gòu)張量矩陣G的特征矢量e和v的特征值。結(jié)構(gòu)張量矩陣特征值可以展開為:

(7a)

(7b)

結(jié)構(gòu)張量矩陣特征矢量可以展開為:

(8a)

(8b)

圖8 層狀模型的波場快照及對應的結(jié)構(gòu)張量的特征值a 層狀模型在1s時刻的波場快照; b 歸一化后垂直于局部波前的特征向量對應的特征值; c 歸一化后平行于局部波前的特征向量對應的特征值

另外,分析公式(1)和公式(7)可以發(fā)現(xiàn),利用公式(7)計算特征值時需要的ux和uz已在波矢量方向估計中進行了計算,額外增加的計算量非常少,因此可以在波矢量估計方法的框架中高效地實現(xiàn)波場交叉判斷。

圖9 局部Sigsbee模型

圖10 局部Sigsbee模型1.6s時刻的波場快照以及對應的結(jié)構(gòu)張量特征值a 1.6s時刻波場快照; b 歸一化后垂直于局部波前的特征向量對應的特征值; c 歸一化后平行局部波前的特征向量對應的特征值

圖11 局部Sigsbee模型3.2s時刻的波場快照以及對應的結(jié)構(gòu)張量特征值a 3.2s時刻波場快照; b 歸一化后垂直于局部波前的特征向量對應的特征值; c 歸一化后平行局部波前的特征向量對應的特征值

5 數(shù)值實驗

首先采用局部Sigsbee 2A模型(地表范圍0～6000m)測試了本文混合方法的有效性(圖12)。檢波點間隔10m,偏移距范圍為(-3000m,3000m)。圖13a 為正演計算的原始單炮炮集,其加入信噪比(S/N)為5的高斯隨機噪聲后如圖13b所示。采用圖13b所示含噪聲炮集計算角度道集,提取其中的3道如圖14所示,提取位置如圖12中藍色豎線所示。圖14a為采用坡印廷矢量估計方法計算的角度道集,圖14b為采用本文局部平面波分解+坡印廷矢量混合的方法得到的角度道集,圖14c 是采用結(jié)構(gòu)張量特征值約束改進后的混合方法得到的角度道集,圖14d是采用解析時間波場得到的角度道集。上述4種方法對應的0～60°角度道集疊加結(jié)果如圖15所示。由圖14a可以看到,由于隨機噪聲的加入,坡印廷矢量方法估計的方向不準確,角度道集中含有較多噪聲,模糊情況嚴重,對應的疊加成像結(jié)果(圖15a)也非常差,整個剖面淹沒在噪聲中,有效同相軸不易識別。而采用本文混合方法計算的角度道集(圖14b和圖14c)質(zhì)量與采用解析波場計算的角度道集(圖14d)質(zhì)量相當,改進后的混合方法產(chǎn)生的角度道集與混合方法產(chǎn)生的角度道集質(zhì)量也基本一致,只在局部區(qū)域質(zhì)量略低于混合方法。0～60°角度道集疊加結(jié)果也說明本文混合方法及改進方法具有很好的抗噪性。

圖12 Sigsbee 2A模型

圖13 正演計算得到的原始單炮炮集(a)及加入高斯隨機噪聲(S/N=5)的結(jié)果(b)

圖14 對含有高斯噪聲(S/N=5)的炮集采用不同方法生成的角度道集a 坡印廷矢量方法; b 局部平面波分解+坡印廷矢量估計的混合方法; c 結(jié)構(gòu)張量特征值約束的混合方法; d 解析時間波場方法(角度范圍:0～60°)

圖15 含有高斯噪聲(S/N=5)的炮集采用不同方法生成的角度道集疊加結(jié)果a 坡印廷矢量方法; b 局部平面波分解+坡印廷矢量估計的混合方法; c 結(jié)構(gòu)張量特征值約束的混合方法; d 解析時間波場方法

分析本文混合方法及改進方法計算的角度道集(圖14b和圖14c)可知,位于炮點下方的角度道集主要集中在小角度范圍,說明該炮下方的照明能量主要集中在小角度位置;隨著地表偏移距的增大,該炮對地下層位照明的入射角越來越大。說明本文混合方法生成的角度道集可以準確地反映地下層位的角度照明情況,為后續(xù)的角度照明分析等奠定基礎(chǔ)。

分析不同角度道集生成方法的單炮計算時間,如圖16所示?？梢钥闯?坡印廷矢量方法計算效率最高,解析時間波場方法效率最差,本文提出的混合方法與解析時間波場方法相比,計算效率有所提升,改進的混合方法在此基礎(chǔ)上有進一步提升。綜合效率和質(zhì)量看,本文改進的混合方法能在滿足一定精度的基礎(chǔ)上,高效地實現(xiàn)逆時偏移角道集的生成。

接下來采用不同方法計算了圖17所示Sigsbee 2A模型右側(cè)(地表范圍6000～10000m)的角度道集。該模型的觀測系統(tǒng)如下:炮點間距10m,檢波點間距10m,炮點范圍0～4000m;每一炮的檢波器接收排列方式為震源兩側(cè)各2000m。圖18a～圖18d分別是坡印廷矢量方法、本文混合方法、改進的混合方法和解析時間波場方法每隔10個CDP點提取的角度道集,其中角度范圍為0～60°。由于坡印廷矢量方向在復雜區(qū)域計算不準確,無法獲得正確的角度值,導致圖18a中角度道集出現(xiàn)較多噪聲,且分布間斷、不連續(xù)。而本文提出的混合方法及改進方法(圖18b 和圖18c)計算的波場方向準確,角度道集連續(xù)性較好,剖面噪聲也少,成像質(zhì)量基本與解析時間波場方法一致(圖18d)。

圖16 采用不同方法計算的每炮時間對比

圖17 Sigsbee 2A模型右側(cè)(地表范圍6000～10000m)

圖18 每隔10個CDP點提取的不同方法的角度道集(Sigsbee 2A模型右側(cè))a 坡印廷矢量方法; b本文混合方法; c 結(jié)構(gòu)張量特征值約束的混合方法; d 解析時間波場方法(角度范圍:0～60°)

對圖6所示水平層狀模型不同方法生成的角度道集進行了校正[27],得到的AVA關(guān)系(淺層)如圖19 所示。其中藍色實線為利用佐普利茲方程近似得到的角度道集隨角度變化的振幅值,紅色十字線為采用坡印廷矢量方法計算的角度道集隨角度變化的振幅值;藍色圓圈線為采用解析時間波場計算的角度道集隨角度變化的振幅值;綠色五星線為采用本文混合方法計算的角度道集隨角度變化的振幅值。可以看到,隨著角度的增大,振幅值也不斷增加,采用混合方法得到的AVA關(guān)系基本上符合實際情況。這說明簡單介質(zhì)和簡單波場情形下計算的角度道集振幅保真性是有保障的。在復雜介質(zhì)情況下,需結(jié)合保真RTM成像才能產(chǎn)生保真的角度成像道集。

圖19 不同方法計算的角度道集振幅值隨角度的變化關(guān)系

6 結(jié)論

本文結(jié)合波矢量方向估計方法和局部平面波分解方法的特點,提出一種混合構(gòu)建角度道集的方法。首先利用局部波場分解方法獲得波前方向,然后利用坡印廷矢量方法確定波的傳播方向,最后采用角度域成像條件提取成像道集,從而高效地獲取高質(zhì)量角度道集。該方法對波場方向的計算比較準確,具有一定的抗噪性,可解決波場交叉等問題,并且很容易擴展到各向異性介質(zhì)。在此基礎(chǔ)上,利用結(jié)構(gòu)張量的特征值快速判斷波場交叉情況,只在波場能量交叉情況下使用混合方法,而在其它區(qū)域采用快速的波矢量估計方法提取波場方向,進一步提高了角度道集的計算效率。然而,由于地下介質(zhì)的復雜性、疊前地震數(shù)據(jù)采集時觀測系統(tǒng)的不完備以及地下照明角度的非均勻性等因素影響,真正產(chǎn)生保真的角度道集還非常困難,需要進行深入的研究。

致謝:感謝中國石油天然氣股份有限公司勘探開發(fā)研究院及西北分院、中海油研究總院和湛江分公司、中國石油化工股份有限公司石油物探技術(shù)研究院和勝利油田分公司對波現(xiàn)象與反演成像研究組(WPI)研究工作的資助與支持。

參 考 文 獻

[1] SYMES W W.A differential semblance criterion for inversion of multioffset seismic reflection data[J].Journal of Geophysical Research:Solid Earth,1993,98(B2):2061-2073

[2] LIU Z,BLEISTEIN N.Migration velocity analysis:theory and an iterative algorithm[J].Geophysics,1995,60(1):142-153

[3] BIONDI B,SYMES W W.Angle-domain common-image gathers for migration velocity analysis by wavefield-continuation imaging[J].Geophysics,2004,69(5):1283-1298

[4] YANG H,XIE X B,LUO M,et al.Target oriented full-wave equation based illumination analysis[J].Expanded Abstracts of 78thAnnual Internat SEG Mtg,2008:2216-2220

[5] YAN R,XIE X B.AVA analysis based on RTM angle-domain common image gather[J].Expanded Abstracts of 82ndAnnual Internat SEG Mtg,2012:1190-1195

[6] FRENCH W S.Two-dimensional and three-dimensional migration of model-experiment reflection profiles[J].Geophysics,1974,39(3):265-277

[7] GARDNER G H F,FRENCH W S,MATZUK T.Elements of migration and velocity analysis[J].Geophysics,1974,39(6):811-825

[8] SCHNEIDER W A.Integral formulation for migration in two and three dimensions[J].Geophysics,1978,43(1):49-76

[9] GRAY S H.Gaussian beam migration of common-shot records[J].Geophysics,2005,70(4):S71-S77

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

[11] BIONDI B,PALACHARLA G.3-D prestack migration of common-azimuth data[J].Geophysics,1996,61(6):1822-1832

[12] ZHANG Y,ZHANG G,BLEISTEIN N.True amplitude wave equation migration arising from true amplitude one-way wave equations[J].Inverse Problems,2003,19(5):1113-1138

[13] BAYSAL E,KOSLOFF D D,SHERWOOD J W C.Reverse time migration[J].Geophysics,1983,48(11):1514-1524

[14] MCMECHAN G A.Migration by extrapolation of time-dependent boundary values[J].Geophysical Prospecting,1983,31(3):413-420

[15] WHITMORE D.Iterative depth migration by backward time propagation[J].Expanded Abstracts of 53rdAnnual Internat SEG Mtg,1983:382-385

[16] ZHANG Y,XU S,TANG B,et al.Angle gathers from reverse time migration[J].The Leading Edge,2010,29(11):1364-1371

[17] XU S,ZHANG Y,TANG B.3D angle gathers from reverse time migration[J].Geophysics,2011,76(2):S77-S92

[18] SAVA P C,FOMEL S.Angle-domain common-image gathers by wavefield continuation methods[J].Geophysics,2003,68(3):1065-1074

[19] SAVA P C,FOMEL S.Time-shift imaging condition in seismic migration[J].Geophysics,2006,71(6):S209-S217

[20] FOMEL S.Theory of 3-D angle gathers in wave-equation imaging[J].Expanded Abstracts of 74thAnnual Internat SEG Mtg,2004:1053-1056

[21] DUVENECK E.A pragmatic approach for computing full-volume RTM reflection angle/azimuth gathers[J].Extended Abstract of 75thEAGE Conference & Exhibition,2013:Tu-11-01

[22] WHITMORE N D,CRAWLEY S.Applications of RTM inverse scattering imaging conditions[J].Expanded Abstracts of 82ndAnnual Internat SEG Mtg,2012:1-6

[23] GUAN H,WILLIAMSON P,DENEL B,et al.Angle-domain common-image gathers extracted from pre-stack RTM images[J].Expanded Abstracts of 83rdAnnual Internat SEG Mtg,2013:3767-3772

[24] WHITMORE N D,CRAWLEY S,ZHU C,et al.Dynamic angle and azimuth decomposition of RTM image[J].Expanded Abstracts of 84thAnnual Internat SEG Mtg,2014:3801-3805

[25] YOON K,MARFURT K J,STARR E W.Challenges in reverse-time migration[J].Expanded Abstracts of 74thAnnual Internat SEG Mtg,2004:1057-1060

[26] YOON K,MARFURT K J.Reverse-time migration using the Poynting vector[J].Exploration Geophysics,2006,37(1):102-107

[27] DICKENS A,WINBOW G A.RTM angle gathers using Poynting vectors[J].Expanded Abstracts of 81stAnnual Internat SEG Mtg,2011:3109-3113

[28] VYAS M,NICHOLS D,MOBLEY E.Efficient RTM angle gathers using source directions[J].Expanded Abstracts of 81stAnnual Internat SEG Mtg,2011:3104-3108

[29] VYAS M,DU X,MOBLEY E,et al.Methods for computing angle gathers using RTM[J].Extended Abstract of 73rdEAGE Conference & Exhibition,2011:F020

[30] ZHANG Q,MCMECHAN G A.Common-image gathers in the incident phase-angle domain from reverse time migration in 2D elastic VTI media[J].Geophysics,2011,76(6):S197-S206

[31] ZHANG Q,MCMECHAN G A.Direct vector-field method to obtain angle-domain common-image gathers from isotropic acoustic and elastic reverse-time migration[J].Geophysics,2011,76(5):WB135-WB149

[32] 王保利,高靜懷,陳文超,等.逆時偏移中用 Poynting 矢量高效地提取角道集[J].地球物理學報,2013,56(1):262-268

WANG B L,GAO J H,CHEN W C,et al.Extracting efficiently angle gathers using Poynting vector during reverse time migration[J].Chinese Journal of Geophysics,2013,56(1):262-268

[33] JIN H,MCMECHAN G A,GUAN H.Comparison of methods for extracting ADCIGs from RTM[J].Geophysics,2014,79(3):S89-S103

[34] JIN H,MCMECHAN G A.Removing smearing-effect artifacts in angle-domain common-image gathers from reverse time migration[J].Geophysics,2015,80(3):U13-U24

[35] XIE X,WU R S.Extracting angle domain information from migrated wavefields[J].Expanded Abstracts of 72ndAnnual Internat SEG Mtg,2002:1360-1363

[36] YAN R,XIE X.Angle gather extraction for acoustic and isotropic elastic RTM[J].Expanded Abstracts of 81stAnnual Internat SEG Mtg,2011:3141-3146

[37] YAN R,XIE X B.An angle-domain imaging condition for elastic reverse time migration and its application to angle gather extraction[J].Geophysics,2012,77(5):S105-S115

[38] TANG B,XU S,ZHANG Y.3D angle gathers with plane-wave reverse-time migration[J].Geophysics,2013,78(2):S117-S123

[39] JIN H,MCMECHAN G A,NGUYEN B.Improving input/output performance in 2D and 3D angle-domain common-image gathers from reverse time migration[J].Geophysics,2015,80(2):S65-S77

[40] YOON K,GUO M,CAI J,et al.3D RTM angle gathers from source wave propagation direction and dip of reflector[J].Expanded Abstracts of 81stAnnual Internat SEG Mtg,2011:3136-3140

[41] ZHAO L,FENG B,WANG H.Efficient RTM angle gathers using source propagation direction and reflector normal direction[J].Expanded Abstracts of 82ndAnnual Internat SEG Mtg,2012:1-5

[42] PATRIKEEVA N,SAVA P.Comparison of angle decomposition methods for wave equation migration[J].Expanded Abstracts of 83rdAnnual Internat SEG Mtg,2013:3773-3778

[43] TANG C.Combining source direction vectors with wavefield decomposition to calculate angle gathers from reverse time migration[J].Expanded Abstracts of 84thAnnual Internat SEG Mtg,2014:3831-3836

[44] ZHANG Q.RTM angle gathers and specular filter(SF) RTM using optical flow[J].Expanded Abstracts of 84thAnnual Internat SEG Mtg,2014:3816-3820

[45] XU S,ZHANG Y,LAMBARé G.Antileakage Fourier transform for seismic data regularization in higher dimensions[J].Geophysics,2010,75(6):WB113-WB120

[46] XU S,ZHANG Y,PHAM D,et al.Antileakage Fourier transform for seismic data regularization[J].Geophysics,2005,70(4):V87-V95

[47] HU J,WANG H,WANG X.Angle gathers from reverse time migration using analytic wavefield propagation and decomposition in the time domain[J].Geophysics,2016,81(1):S1-S9

[48] LIU F,ZHANG G,MORTON S A,et al.An effective imaging condition for reverse-time migration using wavefield decomposition[J].Geophysics,2011,76(1):S29-S39

[49] LIU F,HANSON D W,WHITMORE N D,et al.Toward a unified analysis for source plane-wave migration[J].Geophysics,2006,71(4):S129-S139

[50] FEHMERS G C,H?CKER C F W.Fast structural interpretation with structure-oriented filtering[J].Geophysics,2003,68(4):1286-1293

[51] BROX T,WEICKERT J,BURGETH B,et al.Nonlinear structure tensors[J].Image and Vision Computing,2006,24(1):41-55

[52] BAKKER P.Image structure analysis for seismic interpretation[D].Delft:Delft University of Technology,2002

[53] HALE D.Structure-oriented smoothing and semblance[R].Golden:CWP report,2009