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

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

0 引言

角度域共成像點道集(ADCIGs)包含背景速度和地下角度反射系數(shù)的信息，對于速度建模和儲層描述具有重要意義．與克?；舴蚱?、Beam偏移和單向波偏移相比，逆時偏移(Baysal et al.，1983；Whitmore，1983)能夠精確地處理橫向變速情況以及多波至現(xiàn)象，克服陡傾角限制，成為了復雜介質(zhì)的首選成像方法.逆時偏移角度道集可以采用擴展成像條件方法(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)實現(xiàn)．擴展成像條件方法首先生成地下局部偏移距道集，然后借用傾斜疊加等方法轉(zhuǎn)換為地下角度道集．局部平面波分解方法需要將波場分解為局部平面波分量，應用角度域成像條件提取角度道集，計算代價較大．在逆時偏移角度道集中，算法的高效性和并行性是制約高效提取角度道集的重要因素．波矢量方向估計方法(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)在逆時偏移角度道集生成方法中是一類比較高效的方法．

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

在生成逆時偏移角度道集中，外推波場的方向也可由光學流方法計算．光學流最先由Gibson于1950年提出，被廣泛應用于解決連續(xù)圖像幀之間的視運動問題．目前有許多種計算光學流場的方法(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)是兩類主要的方法，其他方法大都是基于上述方法的補充．Vyas等(2011a)應用LK光學流方法生成角度道集，Zhang(2014)采用HS光學流方法估計波場方向，提取角度道集．然而，LK光學流方法容易出現(xiàn)奇異現(xiàn)象，依賴于局部點的梯度值，而HS光學流方法依賴于整體的能量約束，距離較遠的波場也會參與到波場方向計算中，導致計算的波場方向精度較低．

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

1 方法理論

1.1 波場方向估計

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

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

(1)

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

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)中未知量的個數(shù)大于方程的個數(shù)，該問題是一個欠定問題，需要引入額外的約束求解該問題．

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

(4)

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

(5)

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

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

(6)

(7)

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

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

(8)

(9)

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

1.2 波場方向分解

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

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

(10)

(11)

(12)

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

(13)

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

(14)

(15)

1.3 反射張角和方位角估計

光學流方法計算的波場方向矢量代表了波場的運動差異．在估計反射張角和方位角之前，需要實施歸一化處理．本文在歸一化過程中通過采用對光學流和波場同時約束，避免噪音干擾和弱值的影響．

(16)

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

(17)

(18)

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

(19)

(20)

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

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

(21)

1.4 角度域成像條件

(22)

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

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

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

2 數(shù)值實驗

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

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

圖2 水平層狀速度模型及相應的波場快照(a) 層狀速度模型； (b) 1.0 s時刻的震源端的波場快照； (c) 1.0 s時刻的檢波點端的波場快照．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ù)測試

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

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

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

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

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

圖3 采用不用方法計算的檢波點端的波場方向(a) 坡印廷矢量方法； (b) HS 光學流方法； (c) 本文提出的CLG光學流方法．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 光學流方法； (c) 本文提出的CLG光學流方法．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 采用上下行波分解方法之后的檢波點端的波場快照(1.0 s時刻)(a) 上行波波場快照； (b) 下行波波場快照．Fig.5 The receiver wavefield with the wavefield decomposition method at 1.0 s(a) The upgoing wavefield; (b) The downgoing wavefield.

圖6 采用CLG光學流方法計算的檢波點端的波場方向(a) 采用上行波波場計算的波場方向； (b) 采用下行波波場計算的波場方向； (c) 未采用波場分解方法計算的波場方向．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 層狀模型中估計的反射張角結(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個CDP點提取的角度道集(角度范圍0～60°，角度間隔1°)(a) 常規(guī)的HS光學流方法； (b) 本文提出的CLG光學流+波場分解方法．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光學流方法； (b) 本文提出的CLG光學流+波場分解方法．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 采用不同方法計算的單炮時間對比Fig.10 Comparison of single shot time calculated by different methods

2.2 二維實際資料測試

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

圖11 陸上某地區(qū)2D實際資料，每隔100個CDP點提取的角度道集(角度范圍0～60°，角度間隔1°)(a) 常規(guī)的HS光學流方法； (b) 本文提出的CLG光學流+波場分解方法．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é)論

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

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

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

圖12 陸上某地區(qū)2D實際資料，0～60°角度道集疊加結(jié)果(a) 常規(guī)的HS光學流方法； (b) 本文提出的CLG光學流+波場分解方法．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.

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

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

附錄A：CLG光學流控制方程推導

(A1)

(A2)

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

(A3)

采用變分法，上式誤差泛函對應的Euler-Lagrange方程的解為

(A4)

最后，得到CLG光學流控制方程為

(A5)

其中Δ為拉普拉斯算子．