基于CLG光學(xué)流和波場(chǎng)分解的逆時(shí)偏移角度道集提取方法研究

吳成梁， 王華忠， 馮波， 盛燊

同濟(jì)大學(xué)海洋與地球科學(xué)學(xué)院波現(xiàn)象與智能反演成像研究組， 上海 200092

0 引言

角度域共成像點(diǎn)道集(ADCIGs)包含背景速度和地下角度反射系數(shù)的信息，對(duì)于速度建模和儲(chǔ)層描述具有重要意義．與克?；舴蚱?、Beam偏移和單向波偏移相比，逆時(shí)偏移(Baysal et al.，1983；Whitmore，1983)能夠精確地處理橫向變速情況以及多波至現(xiàn)象，克服陡傾角限制，成為了復(fù)雜介質(zhì)的首選成像方法.逆時(shí)偏移角度道集可以采用擴(kuò)展成像條件方法(Sava and Fomel，2003，2006；Fomel，2004；Duveneck，2013)或局部平面波分解方法(Xie and Wu，2002；Xu et al.，2011；Yan and Xie，2011，Tang et al.，2013；Hu et al.，2015)實(shí)現(xiàn)．?dāng)U展成像條件方法首先生成地下局部偏移距道集，然后借用傾斜疊加等方法轉(zhuǎn)換為地下角度道集．局部平面波分解方法需要將波場(chǎng)分解為局部平面波分量，應(yīng)用角度域成像條件提取角度道集，計(jì)算代價(jià)較大．在逆時(shí)偏移角度道集中，算法的高效性和并行性是制約高效提取角度道集的重要因素．波矢量方向估計(jì)方法(Yoon et al.，2004；Yoon and Marfurt，2006；Dickens and Winbow，2011；Vyas et al.，2011a，b；Zhang and McMechan，2011a，b；王保利等，2013；Jin et al.，2014；Yoon，2017)在逆時(shí)偏移角度道集生成方法中是一類比較高效的方法．

在波矢量方向估計(jì)方法中，可以通過任意空間和任意時(shí)刻的振幅和相位梯度來計(jì)算傳播方向．一旦估計(jì)出震源波場(chǎng)和檢波點(diǎn)波場(chǎng)的方向矢量，就可以計(jì)算出反射張角和方位角，并利用互相關(guān)成像條件提取成像道集．在波矢量估計(jì)方法中，可以采用坡印廷矢量方法計(jì)算波場(chǎng)傳播方向．由于坡印廷矢量方法計(jì)算高效，具有較高的角度分辨率，在逆時(shí)偏移角度道集生成方法中被廣泛地應(yīng)用．但是由于在檢波點(diǎn)端波場(chǎng)比較復(fù)雜，此時(shí)計(jì)算的坡印廷矢量通常會(huì)出現(xiàn)不穩(wěn)定情況，Yoon等(2011)和Vyas等(2011a)提出利用坡印廷矢量計(jì)算震源端波場(chǎng)傳播方向和地下反射界面傾角來估計(jì)反射角；Zhao等(2012)提出通過偏移剖面預(yù)測(cè)反射層法向向量，利用穩(wěn)定的震源端方向波場(chǎng)結(jié)合反射層法向向量來計(jì)算角度道集．雖然坡印廷矢量方法比較高效，但是該方法無法處理波前重疊問題(Patrikeeva and Sava，2013)．

在生成逆時(shí)偏移角度道集中，外推波場(chǎng)的方向也可由光學(xué)流方法計(jì)算．光學(xué)流最先由Gibson于1950年提出，被廣泛應(yīng)用于解決連續(xù)圖像幀之間的視運(yùn)動(dòng)問題．目前有許多種計(jì)算光學(xué)流場(chǎng)的方法(Sobey and Srinivasan，1991；Barron et al.，1994；Brox et al.，2004；Bruhn et al.，2005；Szeliski，2010)．其中Lucas-Kanade方法(LK，或稱局部方法)(Lucas and Kanade，1981)和Horn-Schunck方法(HS方法，或稱全局方法)(Horn and Schunck，1981)是兩類主要的方法，其他方法大都是基于上述方法的補(bǔ)充．Vyas等(2011a)應(yīng)用LK光學(xué)流方法生成角度道集，Zhang(2014)采用HS光學(xué)流方法估計(jì)波場(chǎng)方向，提取角度道集．然而，LK光學(xué)流方法容易出現(xiàn)奇異現(xiàn)象，依賴于局部點(diǎn)的梯度值，而HS光學(xué)流方法依賴于整體的能量約束，距離較遠(yuǎn)的波場(chǎng)也會(huì)參與到波場(chǎng)方向計(jì)算中，導(dǎo)致計(jì)算的波場(chǎng)方向精度較低．

本文提出采用CLG(Combining Local and Global)光學(xué)流方法計(jì)算波場(chǎng)方向，提高波場(chǎng)方向的精度．CLG光學(xué)流方法結(jié)合局部LK方法和全局HS方法的優(yōu)點(diǎn)，在施加全局約束的基礎(chǔ)上，考慮局部點(diǎn)的加權(quán)作用，在保持波場(chǎng)方向估計(jì)穩(wěn)定的基礎(chǔ)上，提高估計(jì)的波場(chǎng)方向精度．另外，為了弱化光學(xué)流方法無法處理波前重疊問題，本文提出利用解析波場(chǎng)對(duì)波場(chǎng)進(jìn)行方向分解，將分解后的行波分別成像，提取角度道集．進(jìn)一步地，在提取角度道集過程中，需要將光學(xué)流場(chǎng)計(jì)算的波場(chǎng)方向轉(zhuǎn)換為反射張角和方位角，本文提出一種有效的歸一化方法和改進(jìn)的最小二乘除法，提高角度計(jì)算的準(zhǔn)確性，避免帶限子波引起的角度計(jì)算的不穩(wěn)定．最后，采用角度道集面元化和規(guī)整化來優(yōu)化角度道集．

1 方法理論

1.1 波場(chǎng)方向估計(jì)

光學(xué)流問題的基本假設(shè)可以描述為后續(xù)時(shí)刻外推的波場(chǎng)值不隨時(shí)間變化：

u(x+Δx,y+Δy,z+Δz,t+Δt)=u(x,y,z,t),

(1)

其中(x,y,z)是地下空間坐標(biāo)，t代表時(shí)間，u是波場(chǎng)．基于小位移假設(shè)，采用一階泰勒展開可得到光學(xué)流方程．

ux·vx+uy·vy+uz·vz+ut=0,

(2)

(ux,v)+ut=0,

(3)

其中x=(x,y,z)，v=(vx,vy,vz)，ux=(ux,uy,uz)．由于在方程(3)中未知量的個(gè)數(shù)大于方程的個(gè)數(shù)，該問題是一個(gè)欠定問題，需要引入額外的約束求解該問題．

在局部LK光學(xué)流方法中，通過假設(shè)未知的光學(xué)流矢量在局部鄰域Ω內(nèi)是恒定不變的，構(gòu)建如下加權(quán)最小二乘泛函求解上述光流問題：

(4)

其中WΩ為在局部鄰域Ω內(nèi)的加權(quán)系數(shù)，最小化ELK(v)得到如下線性方程組：

(5)

在每個(gè)外推時(shí)刻逐個(gè)空間點(diǎn)求解上述3×3的線性方程組即可得到光學(xué)流場(chǎng)．LK光學(xué)流方法計(jì)算量小，只受局部范圍的影響，距離當(dāng)前點(diǎn)較遠(yuǎn)區(qū)域的誤差不會(huì)影響當(dāng)前點(diǎn)的計(jì)算，誤差不具有傳播性．但是LK方法不穩(wěn)定，當(dāng)梯度消失時(shí)，線性方程組會(huì)出現(xiàn)奇異現(xiàn)象．因此，LK光學(xué)流方法通常在稀疏場(chǎng)中實(shí)現(xiàn)，然后采用不同的插值方法插密．此外，光學(xué)流的精度依賴于方程組中梯度張量的準(zhǔn)確性．

在全局HS光學(xué)流方法中，通過引入額外的全局能量的信息約束，求解上述欠定問題，構(gòu)建如下的誤差泛函：

(6)

(7)

其中Δ為拉普拉斯算子．全局HS光學(xué)流方法通過引入整體的能量約束，避免出現(xiàn)奇異現(xiàn)象，即使在某些點(diǎn)上沒有梯度值，但是由于全局能量約束，也能比較穩(wěn)定地求解光學(xué)流場(chǎng)，因此該方法能產(chǎn)生較密的光學(xué)流場(chǎng)．然而波場(chǎng)的傳播具有局部方向性，整體的能量約束帶來的問題是，距離計(jì)算點(diǎn)較遠(yuǎn)區(qū)域的波場(chǎng)也會(huì)參與到當(dāng)前點(diǎn)的計(jì)算中，導(dǎo)致計(jì)算的波場(chǎng)方向精度有所下降．

結(jié)合局部鄰域變化特征和全局的能量約束，本文提出構(gòu)建如下的CLG光學(xué)流誤差泛函：

(8)

(9)

CLG光學(xué)流方法不僅克服了梯度為零時(shí)無法估計(jì)光學(xué)流問題，而且實(shí)現(xiàn)了局部約束以提高抗噪性．采用CLG光學(xué)流方法估計(jì)的波場(chǎng)方向更加精確和穩(wěn)定．對(duì)比HS光學(xué)流方法(方程(7))，CLG光學(xué)流方法增加的計(jì)算量并不大．方程(9)和方程(7)的迭代格式是一致的，比較容易實(shí)施并行化處理，方便融合在逆時(shí)偏移的波場(chǎng)外推過程中．

1.2 波場(chǎng)方向分解

采用CLG光學(xué)流方法可以計(jì)算波場(chǎng)方向．但是光學(xué)流方法的假設(shè)前提是在每個(gè)成像點(diǎn)每個(gè)成像時(shí)刻只有一個(gè)傳播方向．因此，該方法不能處理波前重疊的情況．本文采用解析波場(chǎng)和空間傅里葉變換對(duì)外推的波場(chǎng)進(jìn)行行波方向分解，將分解后的行波分別計(jì)算光學(xué)流場(chǎng)，提取波場(chǎng)傳播方向構(gòu)建角度道集．

在解析波場(chǎng)中，僅包含正頻率信息，僅需對(duì)波場(chǎng)實(shí)施空間域的傅里葉變換即可實(shí)現(xiàn)波場(chǎng)分解，從而可以避免時(shí)間維的傅里葉變換．解析波場(chǎng)是復(fù)數(shù)波場(chǎng)，其實(shí)部可通過波動(dòng)方程外推獲得，其虛部可通過對(duì)地下波場(chǎng)進(jìn)行Hilbert變換得到．但是由于Hilbert變換是關(guān)于時(shí)間方向的卷積，理論上需要將整個(gè)波場(chǎng)存盤并進(jìn)行卷積，對(duì)內(nèi)存和存儲(chǔ)要求比較大．借助于波動(dòng)方程的源項(xiàng)和波場(chǎng)具有線性關(guān)系(Liu et al.，2006)，把對(duì)波場(chǎng)的Hilbert變換轉(zhuǎn)化為震源項(xiàng)的Hilbert變換，利用新的“震源”進(jìn)行傳播，從而構(gòu)建出解析波場(chǎng)的虛部(Hu et al.，2015)．以震源端為例，震源端的解析波場(chǎng)如下所示：

(10)

(11)

(12)

其中F算子代表傅里葉變換，ks為波數(shù)．因此，可以根據(jù)波數(shù)關(guān)系提取任意方向、任意角度寬度傳播的波場(chǎng)．

(13)

雖然理論上可以提取任意角度的波場(chǎng)．但是在逆時(shí)偏移過程中，波場(chǎng)方向分解需要付出較多的計(jì)算代價(jià)．為了平衡計(jì)算代價(jià)和方向估計(jì)精度，我們可以根據(jù)地下介質(zhì)的特征進(jìn)行適應(yīng)性的選擇．若把地下介質(zhì)抽象為在空間上廣泛分布的緩變的層狀沉積層，加上火山活動(dòng)，構(gòu)造運(yùn)動(dòng)等引起的大、小尺度的速度異常體(或波阻抗變化的異常體)(王華忠等，2015)．則只需要采用上下行波分解即可滿足大部分的成像需要，能夠在一定程度上弱化波前交叉問題．此時(shí)行波分解可以簡(jiǎn)單的寫成下式：

(14)

(15)

1.3 反射張角和方位角估計(jì)

光學(xué)流方法計(jì)算的波場(chǎng)方向矢量代表了波場(chǎng)的運(yùn)動(dòng)差異．在估計(jì)反射張角和方位角之前，需要實(shí)施歸一化處理．本文在歸一化過程中通過采用對(duì)光學(xué)流和波場(chǎng)同時(shí)約束，避免噪音干擾和弱值的影響．

(16)

反射張角θ和方位角φ一般可通過下式計(jì)算：

(17)

(18)

其中nx=(1,0,0)和nz=(0,0,1)分別為單位矢量．但是由于地震波是振蕩的．在計(jì)算反射張角和方位角時(shí)，涉及到除法，分母較小的位置會(huì)使得除法變得非常不穩(wěn)定．為了提高角度估計(jì)的穩(wěn)定性和精度，以及保護(hù)弱反射信號(hào)，我們采用最小二乘反演方法估計(jì)反射張角和方位角：

(19)

(20)

其中Ω={x,t}是關(guān)于時(shí)間和空間的局部平滑窗，在本文的測(cè)試中，僅采用了關(guān)于空間的平滑窗．

當(dāng)方向矢量值較小時(shí)會(huì)導(dǎo)致反射張角計(jì)算錯(cuò)誤，并且計(jì)算錯(cuò)誤的反射張角幾乎為90°．對(duì)于偏移成像，較大的入射角通常對(duì)應(yīng)需要濾除的低波數(shù)噪音．因此，本文采用改進(jìn)的最小二乘法來估計(jì)反射張角，把因方向矢量過小計(jì)算錯(cuò)誤的角度值轉(zhuǎn)移到更大的入射角度中，避免污染成像道集．最終計(jì)算反射張角的法方程為

(21)

1.4 角度域成像條件

(22)

1.5 角度道集面元化和規(guī)整化

由于采集的不規(guī)則地震數(shù)據(jù)，復(fù)雜的地下傳播路徑，在地下角度域中成像道集是不規(guī)則的．在成像點(diǎn)處，不同角度面元和地下不同位置的角度道集覆蓋情況是不同的．隨著地下空間位置的不同，角度道集會(huì)出現(xiàn)欠采樣和過采樣問題，實(shí)際的地下角度采樣是散亂不規(guī)則的．另外，由于數(shù)值計(jì)算的離散化問題，實(shí)際的角度道集不能放在正確的放置，不同的角度面元間隔會(huì)使得最終的不同角度的疊加次數(shù)不同，影響AVA分析．離散的、不規(guī)則的角度道集需要在角度域?qū)嵤┟嬖幚恚诒疚闹?，我們采用角度插值方法?duì)角度道集進(jìn)行面元化和平滑．進(jìn)一步地，還需要考慮角度照明和覆蓋問題．

另外，由于偏移速度的不準(zhǔn)確、地震子波的空間變化以及地下介質(zhì)的復(fù)雜性，角道集可能會(huì)出現(xiàn)剩余深度時(shí)差和子波拉伸等現(xiàn)象．成像道集后處理技術(shù)，如動(dòng)態(tài)時(shí)間規(guī)劃和局部相似疊加(Wu et al.，2019)，可用于優(yōu)化角度道集，提高成像道集質(zhì)量．

2 數(shù)值實(shí)驗(yàn)

基于CLG光學(xué)流和波場(chǎng)分解提取角道集流程如圖1所示．考慮實(shí)際的多炮地震數(shù)據(jù)特征和當(dāng)前的計(jì)算機(jī)集群的結(jié)構(gòu)特點(diǎn)，本文采用MPI主從模式對(duì)炮數(shù)據(jù)并行，在每個(gè)單炮偏移過程中，采用OpenMp多線程并行策略．在單炮偏移過程中，首先將震源子波和該炮觀測(cè)的地震數(shù)據(jù)進(jìn)行希爾伯特變換，然后分別進(jìn)行正演傳播，得到炮點(diǎn)端和檢波點(diǎn)端的解析波場(chǎng)．對(duì)解析波場(chǎng)進(jìn)行波場(chǎng)分解，得到分解后的上下行波(在本文測(cè)試中，波場(chǎng)分解僅實(shí)施上下行波分解)．然后采用CLG光學(xué)流方法分別計(jì)算上下行波場(chǎng)中的傳播方向，并將得到的光學(xué)流矢量場(chǎng)實(shí)施歸一化處理．然后采用公式(21)和(20)計(jì)算反射張角和方位角．最后采用角度域成像條件提取對(duì)應(yīng)的角度道集．然后將生成的角度道集進(jìn)行角度插值，避免角度面元化和數(shù)值化引起的誤差．每個(gè)單炮成像道集結(jié)果存放在每個(gè)處理器附帶的局部盤上．把所有炮成像道集數(shù)據(jù)體收集規(guī)約在一起，放到全局盤上，就可以產(chǎn)生多炮角度道集．最后，對(duì)角度道集進(jìn)行必要的后處理，形成最終的結(jié)果．

圖1 本文提出的基于CLG光學(xué)流和波場(chǎng)分解的逆時(shí)偏移角度道集方法流程圖Fig.1 The flow chart for extracting angle gathers in RTM based on the proposed CLG optical flow method with the wavefield decomposition

圖2 水平層狀速度模型及相應(yīng)的波場(chǎng)快照(a) 層狀速度模型； (b) 1.0 s時(shí)刻的震源端的波場(chǎng)快照； (c) 1.0 s時(shí)刻的檢波點(diǎn)端的波場(chǎng)快照．Fig.2 The layered velocity model and corresponding wave field snapshot(a) The layered velocity model; (b) The snapshot of the source wavefield at 1.0 s; (c) The snapshot of the receiver wavefield at 1.0 s.

2.1 二維合成數(shù)據(jù)測(cè)試

首先采用層狀模型說明本文提出的波場(chǎng)方向估計(jì)方法的有效性．其中水平層狀速度模型如圖2a所示，采用有限差分正演模擬方法生成的震源端波場(chǎng)和檢波點(diǎn)端波場(chǎng)快照如圖2b和2c所示．在該模型中，反射層深度為1500 m，分析可知只有反射界面附近的波前信息才對(duì)反射角道集有效．因此，該區(qū)域的波場(chǎng)方向必須被正確的估計(jì)．采用不同方法計(jì)算的反射界面附近的局部放大的檢波點(diǎn)端波場(chǎng)估計(jì)的方向如圖3所示．其中圖3a是坡印廷矢量方法的結(jié)果，采用HS和CLG光學(xué)流方法估計(jì)的波場(chǎng)方向分別如圖3b和3c所示．在光學(xué)流方法中，考慮到LK方法可能出現(xiàn)奇異現(xiàn)象，因此在本文中，僅對(duì)比HS光學(xué)流方法和與本文提出的CLG光學(xué)流方法．從圖3a中可以看出，在大多數(shù)區(qū)域，采用坡印廷矢量估計(jì)的波場(chǎng)方向都是不正確的，波場(chǎng)方向不能被有效的識(shí)別，特別是圖3a中的紅圈所示，該部分的波矢量方向雜亂無章，無法獲得正確的波場(chǎng)角度值．坡印廷矢量方法計(jì)算的波方向矢量中，有效的波矢量較少，存在較多幅值較小的波方向矢量(圖3中，箭頭之后線段的長(zhǎng)短代表了波矢量的幅值大小)，該部分波矢量無法用于波場(chǎng)方向計(jì)算，如圖3a中綠圈所示．而HS光學(xué)流方法的計(jì)算結(jié)果則優(yōu)于坡印廷矢量方法，其中HS光學(xué)流方法計(jì)算的波矢量幅值比坡印廷矢量方法更加平衡，波矢量方向更加準(zhǔn)確．這是因?yàn)槠掠⊥⑹噶康慕Y(jié)果只是對(duì)光學(xué)流場(chǎng)的初始估計(jì)(Zhang，2014)．但是在圖3b綠圈所示中，仍然存在較多幅值不均的波方向矢量，導(dǎo)致波場(chǎng)方向估計(jì)的精度降低．另外，在圖3b中紅圈中，雜亂無章的波方向矢量仍然存在，在該區(qū)域，計(jì)算的波場(chǎng)方向是錯(cuò)誤的．對(duì)比HS光學(xué)流方法和CLG光學(xué)流方法，CLG光學(xué)流方法能更準(zhǔn)確地估計(jì)波場(chǎng)方向，特別在圖3c中紅圈位置，差異更為明顯．CLG光學(xué)流方法計(jì)算的波矢量幅值比HS光學(xué)流方法更加平衡(如圖3c綠圈所示)，CLG光學(xué)流方法可以同時(shí)改善估計(jì)的方向和幅度．地震波場(chǎng)的方向在一定的范圍具有相似性(子波的帶限作用)，CLG光學(xué)流方法估計(jì)的波場(chǎng)方向更具有連續(xù)性，更符合地下介質(zhì)波場(chǎng)傳播的特征．其中0～60°角道集疊加的成像結(jié)果如圖4所示．由于坡印廷矢量方法波場(chǎng)方向不準(zhǔn)確，導(dǎo)致在對(duì)應(yīng)的成像剖面中出現(xiàn)了不該出現(xiàn)的大角度的低頻噪音．而在HS光流法方法中，低頻噪音被有效的去除．但是在大偏移距位置(如圖4中白圈所示)，成像結(jié)果不聚焦，彌散在成像剖面中．這是由于在該區(qū)域，波場(chǎng)的方向計(jì)算不準(zhǔn)確導(dǎo)致(如圖3中所示)．而在本文提出的方法中，成像效果較好，優(yōu)于HS方法，從而說明了方向估計(jì)的有效性．

接下來測(cè)試波場(chǎng)方向分解的有效性．圖5為采用上下行波波場(chǎng)分解方法得到的檢波點(diǎn)端波場(chǎng)快照．圖6a為在上行波場(chǎng)(圖5a)中應(yīng)用CLG光學(xué)流方法估計(jì)的局部波場(chǎng)方向(如圖5a中白框所示)．圖6b為在下行波場(chǎng)(圖5b)中采用CLG光學(xué)流方法估計(jì)的局部波場(chǎng)方向(如圖5b中白框所示)．圖6c為在未采用波場(chǎng)分解的檢波點(diǎn)端波場(chǎng)中應(yīng)用CLG光學(xué)流方法估計(jì)的波場(chǎng)方向(如圖2c中白框所示)．可以看到，在圖6c中，波形交叉部分(圖6c中紅圈所示)，仍然存在一些計(jì)算錯(cuò)誤的波矢量方向．而在波場(chǎng)分解之后，該區(qū)域的波場(chǎng)方向得到正確的估計(jì)(圖6b中紅圈所示)，有效提高了波矢量計(jì)算的精度，從而說明了波場(chǎng)分解的有效性．另外，圖6a中所示的波場(chǎng)方向與炮點(diǎn)端波場(chǎng)方向是相反的，互相關(guān)成像結(jié)果為大角度的低頻噪音，在偏移成像中是需要去除的．因此，采用上下行波分解，可以有效地減少這些噪音的影響．

接下來測(cè)試本文提出的角度估計(jì)的有效性．圖7a為常規(guī)除結(jié)果．本文提出的最小二乘反演除的結(jié)果如圖7b所示(圖7a和圖7b均采用了常規(guī)的HS光學(xué)流方法計(jì)算的波矢量結(jié)果來估計(jì)反射張角)．由于帶限子波的震蕩作用，常規(guī)方法計(jì)算的角度間斷、不連續(xù)．采用反演的方法估計(jì)的角度值更加穩(wěn)健精確．子波的帶限作用引起的角度誤差得到有效的消除．

然后采用Sigsbee 2A理論模型說明本文提出方法的有效性．圖8a為采用常規(guī)的HS光學(xué)流方法生成的角度道集，采用本文提出的CLG光學(xué)流方法+波場(chǎng)上下行波分解生成的角度道集如圖8b所示．可以看出，常規(guī)的HS光學(xué)流方法的角度道集中含有較多噪聲，模糊現(xiàn)象嚴(yán)重，角度道集彌散在比較廣的角度中，特別是在鹽丘的左側(cè)和上部，成像道集質(zhì)量比較差．然而在采用本文提出方法的成像道集中，角度道集比較連續(xù)、聚焦．對(duì)應(yīng)地0～60°角度道集疊加結(jié)果如圖9所示．可以明顯地看出，采用CLG光學(xué)流方法+波場(chǎng)上下行波分解之后的疊加成像中，假象更少，反射軸更加清楚．

在Sigsbee 2A模型上分析不同角度道集生成方法的單炮計(jì)算時(shí)間，如圖10所示．可以看出，坡印廷矢量方法計(jì)算效率最高，常規(guī)的HS光學(xué)流方法次之．而本文提出的CLG光學(xué)流方法+波場(chǎng)分解方法需要的計(jì)算時(shí)間則是最長(zhǎng)．這是由于在本文提出的方法中，需要采用解析波場(chǎng)和空間傅里葉變換對(duì)外推的波場(chǎng)進(jìn)行行波方向分解．波場(chǎng)方向分解依賴于解析波場(chǎng)，解析波場(chǎng)需要額外一倍的外推時(shí)間．本文的空間傅里葉變換是對(duì)整個(gè)波場(chǎng)進(jìn)行分解，不同于局部的波場(chǎng)分解方法(有關(guān)局部波場(chǎng)分解方法和常規(guī)的波矢量方法的效率對(duì)比可參考文獻(xiàn)(吳成梁等，2018))，需要的計(jì)算量并不大．另外，CLG光學(xué)流的計(jì)算，最小二乘角度估計(jì)，也都需要一定的計(jì)算量．對(duì)比常規(guī)的HS光學(xué)流方法和本文提出的方法，本文提出的CLG光學(xué)流+波場(chǎng)分解方法的計(jì)算效率是常規(guī)HS光學(xué)流方法的3倍左右．對(duì)比成像效果的提升和當(dāng)前計(jì)算機(jī)的快速發(fā)展，本文認(rèn)為這些計(jì)算代價(jià)是值得的．

圖3 采用不用方法計(jì)算的檢波點(diǎn)端的波場(chǎng)方向(a) 坡印廷矢量方法； (b) HS 光學(xué)流方法； (c) 本文提出的CLG光學(xué)流方法．Fig.3 The direction of the receiver wavefield estimated by different methods(a) Poynting vector method; (b) The HS optical flow method; (c) The proposed CLG optical flow method.

圖4 0～60°角度道集疊加結(jié)果(a) 坡印廷矢量方法； (b) HS 光學(xué)流方法； (c) 本文提出的CLG光學(xué)流方法．Fig.4 The imaging results stacking by the 0～60°angle gathers(a) Poynting vector method; (b) The HS optical flow method; (c) The proposed CLG optical flow method.

圖5 采用上下行波分解方法之后的檢波點(diǎn)端的波場(chǎng)快照(1.0 s時(shí)刻)(a) 上行波波場(chǎng)快照； (b) 下行波波場(chǎng)快照．Fig.5 The receiver wavefield with the wavefield decomposition method at 1.0 s(a) The upgoing wavefield; (b) The downgoing wavefield.

圖6 采用CLG光學(xué)流方法計(jì)算的檢波點(diǎn)端的波場(chǎng)方向(a) 采用上行波波場(chǎng)計(jì)算的波場(chǎng)方向； (b) 采用下行波波場(chǎng)計(jì)算的波場(chǎng)方向； (c) 未采用波場(chǎng)分解方法計(jì)算的波場(chǎng)方向．Fig.6 The direction of the receiver wavefield estimated by the proposed CLG optical flow method(a) With the upgoing wavefield; (b) With the downgoing wavefield; (c) Without the wavefield decomposition.

圖7 層狀模型中估計(jì)的反射張角結(jié)果(a) 常規(guī)的方法； (b) 本文提出的最小二乘除方法．Fig.7 The estimated reflection angle in the layered model by (a) conventional method and (b) the proposed least square division method

圖8 Sigsbee 2A模型每隔100個(gè)CDP點(diǎn)提取的角度道集(角度范圍0～60°，角度間隔1°)(a) 常規(guī)的HS光學(xué)流方法； (b) 本文提出的CLG光學(xué)流+波場(chǎng)分解方法．Fig.8 The ADCIGs of Sigsbee 2A model. The ADCIGs is selected with a 100-CDP interval and the angle range is 0～60 degrees with a 1 degree interval(a) The conventional HS optical flow method; (b) The proposed CLG optical flow method with the wavefield decomposition.

圖9 Sigsbee 2A模型0～60°角度道集疊加結(jié)果(a) 常規(guī)的HS光學(xué)流方法； (b) 本文提出的CLG光學(xué)流+波場(chǎng)分解方法．Fig.9 The stacked image of Sigsbee 2A model by ADCIGs. (The angle range is 0～60 degrees)(a) The conventional HS optical flow method; (b) The proposed CLG optical flow method with the wavefield decomposition.

圖10 采用不同方法計(jì)算的單炮時(shí)間對(duì)比Fig.10 Comparison of single shot time calculated by different methods

2.2 二維實(shí)際資料測(cè)試

接下來，采用某地區(qū)陸上實(shí)際資料說明本文提出方法的有效性．其中圖11a為常規(guī)的HS光學(xué)流方法生成的角度道集，圖11b為本文提出的方法生成的角度道集．對(duì)比可以看出，采用本文提出方法生成的角度道集假象比較少，成像道集比較聚焦．另外，從疊加成像結(jié)果(圖12)可以看出，常規(guī)的光學(xué)流方法結(jié)果中，存在較多的偏移假象，疊加剖面比較模糊，反射波同相軸不連續(xù)．在本文提出的方法中，成像結(jié)果連續(xù)性較好，剖面噪聲和偏移假象較少，成像質(zhì)量結(jié)果明顯好于常規(guī)的光學(xué)流方法．

圖11 陸上某地區(qū)2D實(shí)際資料，每隔100個(gè)CDP點(diǎn)提取的角度道集(角度范圍0～60°，角度間隔1°)(a) 常規(guī)的HS光學(xué)流方法； (b) 本文提出的CLG光學(xué)流+波場(chǎng)分解方法．Fig.11 The ADCIGs of real data. The ADCIGs is selected with a 100-CDP interval and the angle range is 0～60 degrees with a 1 degree interval (a) The conventional HS optical flow method; (b) The proposed CLG optical flow method with the wavefield decomposition.

3 討論與結(jié)論

利用波動(dòng)方程的時(shí)間演化進(jìn)行波矢量方向計(jì)算在逆時(shí)偏移角度道集中是一個(gè)核心內(nèi)容．局部波場(chǎng)分解方法能夠處理復(fù)雜波場(chǎng)，傳播方向計(jì)算準(zhǔn)確，但是效率低下，無法在大規(guī)模的3D逆時(shí)偏移角度道集中實(shí)用化．坡印廷矢量方法可以用于計(jì)算波場(chǎng)傳播方向，該方法計(jì)算簡(jiǎn)單，計(jì)算效率高，適合并行處理，但是坡印廷矢量通常會(huì)出現(xiàn)不穩(wěn)定情況，在波場(chǎng)復(fù)雜情況時(shí)，波場(chǎng)方向估計(jì)不準(zhǔn)確．

光學(xué)流方法可以提高波場(chǎng)方向估計(jì)的精度．光學(xué)流方程本身是欠定的，需要引入額外的約束計(jì)算光學(xué)流場(chǎng)．通過假設(shè)在局部鄰域Ω內(nèi)光學(xué)流矢量是不變的，LK光學(xué)流方法構(gòu)建局部的線性方程組計(jì)算方向矢量．然而，LK光學(xué)流方法依賴于局部點(diǎn)的梯度值，容易出現(xiàn)奇異現(xiàn)象．全局HS光學(xué)流方法通過引入全局的能量約束，能夠比較穩(wěn)定地求解光學(xué)流場(chǎng)，但是在該方法中，距離計(jì)算點(diǎn)較遠(yuǎn)區(qū)域的波場(chǎng)也會(huì)參與到當(dāng)前點(diǎn)的計(jì)算中，導(dǎo)致計(jì)算的波場(chǎng)方向精度有所下降．本文采用CLG光學(xué)流方法計(jì)算波場(chǎng)方向，在施加全局約束的基礎(chǔ)上，考慮局部點(diǎn)的加權(quán)作用，在保持波場(chǎng)方向估計(jì)穩(wěn)定的基礎(chǔ)上，提高波場(chǎng)方向計(jì)算精度．該方法簡(jiǎn)單高效，便于并行處理，其迭代格式與HS光學(xué)流方法相同，幾乎沒有太多的額外計(jì)算量，非常適合大規(guī)模地震數(shù)據(jù)處理．

為了克服光學(xué)流方法無法處理波前重疊問題，本文利用解析波場(chǎng)對(duì)波場(chǎng)進(jìn)行方向分解來減弱這種影響．通過構(gòu)建解析波場(chǎng)，僅需波場(chǎng)的空間傅里葉變換即可實(shí)現(xiàn)任意波場(chǎng)方向分解，有效地避免了逆時(shí)偏移外推過程存儲(chǔ)波場(chǎng)問題，減少了內(nèi)存和硬盤使用．理論上可以通過波數(shù)域?yàn)V波提取任意方向的波場(chǎng)，但是考慮到計(jì)算量問題以及地下介質(zhì)的特征主要以層狀沉積層加上大、小尺度的速度異常體，本文的數(shù)值實(shí)驗(yàn)僅考慮了上下行波分解情況．針對(duì)不同的地質(zhì)構(gòu)造特征，可以任選幾個(gè)方向角度進(jìn)行分解，波場(chǎng)分解的角度需要根據(jù)地質(zhì)特征確定．方向分解個(gè)數(shù)越多，波前交叉分的更仔細(xì)，但是需要的計(jì)算量越大．波場(chǎng)分解的個(gè)數(shù)需要在計(jì)算量和效果中達(dá)到平衡．

圖12 陸上某地區(qū)2D實(shí)際資料，0～60°角度道集疊加結(jié)果(a) 常規(guī)的HS光學(xué)流方法； (b) 本文提出的CLG光學(xué)流+波場(chǎng)分解方法．Fig.12 The stacked image of real data by ADCIGs. (The angle range is 0～60 degrees)(a) The conventional HS optical flow method; (b) The proposed CLG optical flow method with the wavefield decomposition.

進(jìn)一步地，本文提出采用有效的歸一化方法和最小二乘除方法，提高反射角的估計(jì)精度，有效避免帶限波場(chǎng)引起的震蕩作用．并對(duì)角度道集進(jìn)行了面元化和規(guī)整化處理，以獲得高質(zhì)量的角度道集．理論和實(shí)際資料證明了本文提出方法的有效性．

致謝感謝中石油勘探開發(fā)研究院及西北分院、中海油研究院和湛江分公司、中石化物探技術(shù)研究院和勝利油田分公司對(duì)波現(xiàn)象與智能反演成像研究組(WPI)研究工作的資助與支持．感謝審稿人提供的寶貴建議．

附錄A：CLG光學(xué)流控制方程推導(dǎo)

(A1)

(A2)

重新整理，誤差泛函ECLG(ω)可表示為

(A3)

采用變分法，上式誤差泛函對(duì)應(yīng)的Euler-Lagrange方程的解為

(A4)

最后，得到CLG光學(xué)流控制方程為

(A5)

其中Δ為拉普拉斯算子．