劉少勇 韓冰凱* 顧漢明 宋桂橋

(①中國地質大學地球物理與空間信息學院，湖北武漢 430074; ②地球內部多尺度成像湖北省重點實驗室，湖北武漢 430074； ③中國石化油田勘探開發(fā)事業(yè)部，北京 100728)

帶限局部平面波傳播算子與偏移方法

劉少勇①②韓冰凱*①②顧漢明①②宋桂橋③

(①中國地質大學地球物理與空間信息學院，湖北武漢 430074; ②地球內部多尺度成像湖北省重點實驗室，湖北武漢 430074； ③中國石化油田勘探開發(fā)事業(yè)部，北京 100728)

劉少勇，韓冰凱，顧漢明，宋桂橋.帶限局部平面波傳播算子與偏移方法.石油地球物理勘探，2017,52(5):948-955.

地震信號的帶限特點限制了高頻近似下的射線類算子在地震波模擬和偏移中的應用。本文基于射線理論構建帶限信號局部平面波傳播算子，兼顧射線類方法具有時空局部性和實現(xiàn)靈活高效的特點，有效進行帶限地震波場的模擬和偏移。即分析傳統(tǒng)高頻射線在速度界面的傳播特征，利用帶限信號頻帶與菲涅耳帶的關系，構建帶限斯奈爾定律進行中心射線傳播；在局部平面波假設下，基于帶限中心射線做旁軸近似展開，構建帶限局部平面波傳播算子；在此基礎上發(fā)展基于該算子的帶限射線束偏移方法。數(shù)值實驗結果表明，帶限射線束較傳統(tǒng)基于高頻近似的射線束具有更強穿透能力，改善射線陰影區(qū)照明，從而提高復雜構造區(qū)成像質量。

射線束偏移 帶限信號 傳播算子 局部平面波

1 引言

基于高頻近似的射線類傳播算子被廣泛應用于描述地震波傳播和偏移成像，但實際采集的地震數(shù)據(jù)是具有一定帶寬的數(shù)字信號。在高頻漸進射線理論下，地震波旅行時通常被表示為慢度沿著射線路徑的積分[1]。Woodward[2]提出，在一定頻帶范圍內的地震波旅行時通常是由慢度在炮檢點之間射線路徑周圍一定范圍內的空間體積分決定的，該旅行時與地震波頻帶相關。

Biondi[3]通過對赫姆霍茲方程進行近似，推導出頻率相關的程函方程進行旅行時計算，該方法在計算低頻部分旅行時時受限于高頻端的頻帶范圍。Hogan等[4]討論了求解頻率相關程函方程計算低頻波場的條件，即介質的平滑程度。另一類頻率相關的旅行時計算方法可通過射線追蹤實現(xiàn)，F(xiàn)oreman[5]從赫姆霍茲方程出發(fā)推導了其對應的頻率相關的射線追蹤系統(tǒng)。這種方法相比于傳統(tǒng)的射線追蹤，計算復雜度顯著增加。Protasov等[6]分析該頻率相關的射線追蹤算法的特點，將其與射線追蹤、時間域有限差分正演的波場進行對比分析，指出此方法得到帶限射線的射線路徑與震源主頻有關，能較好地改善射線陰影區(qū)問題。Lomax[7]利用波長平滑進行帶限波場傳播，具體通過對垂直于射線平面的速度加權平均來實現(xiàn)，其加權函數(shù)的寬度正比于波長。Lomax等[8]分析了波長平滑算法的原理，指出通過特定的加權函數(shù)可控制帶限波場的運動學特征，但這種加權平均的方法在速度差異較大的不規(guī)則界面處難以準確地刻畫射線的傳播。Protasov等[9]研究了射線在復雜速度界面?zhèn)鞑ヌ卣?，其核心思想是在復雜界面利用Kirchhoff邊界積分對積分區(qū)域內射線進行先加權后傳播。Protasov等[10]指出利用該方法構造出的帶限射線方向對應的是等效波場透過界面的最大能量方向，并提出通過構建界面處菲涅耳帶透射波場的最大值問題來求解帶限射線路徑。Yarman等[11]對Kirchhoff邊界積分引入射線級數(shù)解，通過積分表達式中相位部分的近似導出積分區(qū)間為第一菲涅耳帶，并構建該積分區(qū)間和波場頻帶的關系，將此方法稱為帶限射線追蹤。

本文結合帶限射線和射線束的特點，從分析傳統(tǒng)高頻射線與帶限射線的傳播特征出發(fā)，研究適應于復雜介質的帶限射線追蹤。在此基礎上，通過旁軸近似構建帶限局部平面波傳播算子模擬帶限波場，改善傳統(tǒng)射線類方法在射線陰影區(qū)和中深層的照明及偏移效果。

2 方法原理

2.1 帶限局部平面波傳播算子

高頻射線在速度界面的傳播可由斯奈爾定律控制，帶限射線則對應于帶限波場傳播的最大能量方向[9]。帶限斯奈爾定律可控制帶限射線的傳播，但沒有傳統(tǒng)斯奈爾定律的簡潔數(shù)學表達式。帶限射線通過界面時的影響范圍不再是一個點，而是由菲涅耳帶控制的一個范圍，在數(shù)學上可表達為求解透射波場能量最大的優(yōu)化問題對應的信任區(qū)間。在速度界面處，中心射線透過界面上一點x0的透射波場可表示為

(1)

式中：x為地下成像點；A(x)為與格林函數(shù)有關的振幅項；T(x)為透射系數(shù)； 積分區(qū)間FZx0為中心射線路徑坐標x0附近的菲涅耳帶的范圍，由頻率f和與旅行時有關的相位項Δτ(x,x0)決定。如圖1所示，中心射線穿過速度界面時，在其附近的菲涅耳帶FZx0內構建局部平面Γx0。局部平面內，以高頻射線描述等效局部平面波的傳播，入射射線參數(shù)pin(x)和出射射線參數(shù)pout(x)遵循運動學射線方程組

(2)

式中：p表示射線的慢度矢量；v表示射線路徑上的速度；s和τ分別表示弧長和射線的旅行時；x∈Γx0。由高頻射線的出射射線參數(shù)pout(x)加權得到等效局部平面波的出射方向pe(x0)，對應透射波場u(x0)的最大能量方向，即帶限中心射線在速度界面處以帶限斯奈爾定律傳播

|u(x0)|

(3)

基于旁軸近似理論， 由帶限中心射線可構建帶限局部平面波傳播算子。局部平面波寬度范圍的旅行時由中心射線旁軸近似展開[24]

圖1 帶限射線在速度界面的傳播示意圖

(4)

式中： Δs為x在中心射線上投影點與x0間距離；r為x到中心射線的距離(圖2a)；mr表征波前曲率

(5)

(6)

式中：v0表示中心射線初射點的速度； Δθ0表示中心射線初射角度間隔(圖2b)。

圖2 笛卡爾坐標系下中心射線旅行時旁軸近似展開(a)與射線束寬度(b)

2.2 帶限射線束偏移

射線束偏移區(qū)別于Kirchhoff積分偏移單道輸入的特點，其對應的輸入數(shù)據(jù)也是由地表地震記錄分解的局部平面波。傳統(tǒng)τ-p變換可將地表地震記錄分解成不同方向的局部平面波。這種線性Randon變換方法由于其基函數(shù)不完備，往往存在能量泄露而降低局部平面波的分辨率。平面波分解通?？山y(tǒng)一在反演理論框架下進行，通過構建反問題求解以獲得更高的分辨率。本文通過在反問題中引入L0范數(shù)的約束，在壓縮感知框架下求解局部平面波問題?；诜纸夂蟮木植科矫娌ㄟM行帶限射線束的成像。

在頻率域中，構建平面波分解的正問題描述為：局部平面波源經(jīng)過傳播矩陣得到地表局部射線束范圍內地震數(shù)據(jù)

d(x,ω)=A(x,pe,ω)D(pe,ω)

(7)

式中：d(x,ω)表示地表局部地震數(shù)據(jù)；D(pe,ω)表示待求解的局部平面波；A(x,pe,ω)表示局部平面波源和地表地震數(shù)據(jù)之間的投影矩陣。通過求解對應反問題獲得射線束控制范圍內的局部平面波源，即給定相移矩陣和空間局部信號求解局部平面波。該反問題的實現(xiàn)通常借助于求解如下最小二乘泛函實現(xiàn)

J=‖A(x,pe,ω)D(pe,ω)-d(x,ω)‖2

(8)

矩陣A(x,pe,ω)通常為欠定矩陣，造成上述反演問題不穩(wěn)定，可通過加入模型約束求最小二乘解，但是其結果通常受到能量泄露的影響。疊前地震數(shù)據(jù)可以通過在變換域中實現(xiàn)數(shù)據(jù)的稀疏表達，求解欠定問題比較有效的方法是假設模型具有一定的稀疏性[28]，對應在局部平面波分解中，該稀疏性的物理解釋為：在局部范圍內僅存在少數(shù)幾個平面波源。這種假設在局部空間窗數(shù)據(jù)中通常成立，可借助于對模型的L1和L0范數(shù)約束來實現(xiàn)[29]。更進一步，為了得到更稀疏的局部平面波解，可將模型估計改為L0范數(shù)約束

min ‖D(pe,ω)‖0使得

‖A(x,pe,ω)D(pe,ω)-d(x,ω)‖2<ε

(9)

該模型的L0范數(shù)約束代表解的個數(shù)極小[28]，在射線束合成中代表局部平面波的個數(shù)較少。模型L0范數(shù)約束是模型稀疏約束比較極端的情況，通常該泛函借助一些貪婪算法實現(xiàn)，比如匹配追蹤等方法。

基于以上局部平面波數(shù)據(jù)，利用帶限局部平面波傳播算子進行炮檢端的波傳播加上成像條件就可以進行成像處理。射線束成像條件類似于傳統(tǒng)Kirchhoff積分偏移的等時成像條件，只是在成像過程前已經(jīng)對數(shù)據(jù)進行了方向篩選。傳統(tǒng)的Kirchhoff積分偏移成像公式可寫成[30]

eiωτdωdΩ

(10)

式中：Ω表示包含炮點xs和中心檢波點xr0的成像孔徑；ps和pr分別表示炮點、檢波點射線參數(shù)；W為加權函數(shù)；D表示頻率域局部平面波；τ為旅行時。在聲學介質中，旅行時表達是常數(shù)，D在頻率域的積分可以通過與δ函數(shù)卷積得到時間域表達的局部平面波d。由此，通過引入激勵時間成像條件去除式(10)的頻率域積分，得到時空域成像公式

pr(xr0),τ]dΩ|τ=τs+τr

(11)

3 帶限射線束偏移的實現(xiàn)

地震數(shù)據(jù)偏移的實現(xiàn)通常有兩個核心步驟：炮檢波場的傳播和成像條件的時間。射線束偏移一般采用局部平面波傳播算子，本節(jié)從方法的具體實現(xiàn)出發(fā)，對帶限局部平面波傳播算子的構建和帶限射線束的偏移流程進行說明。

如圖3所示，時空域共炮集記錄的帶限射線束偏移流程中帶限局部平面波傳播算子的構建分為三個步驟： ①構建局部平面波，當中心射線經(jīng)過模型空間的非均勻區(qū)域時，在射線束寬度范圍內的局部平面內以一定密度高頻射線近似局部平面波的傳播； ②以帶限射線追蹤計算射線路徑和旅行時，在局部平面內加權計算等效射線參數(shù)作為中心射線的參數(shù)，以帶限斯奈爾定律描述中心射線的傳播； ③帶限中心射線旅行時和振幅信息由旁軸近似展開分配到射線束寬度范圍內的成像網(wǎng)格點上。

圖3 帶限射線束偏移實現(xiàn)流程圖

結合當前計算集群多節(jié)點、多核心共享存儲的特點，為了提高帶限射線束偏移的效率，其實現(xiàn)過程需遵循如下原則： ①炮點端各方向格林函數(shù)依次計算并保存到內存中，避免后續(xù)炮檢點射線束相交成像過程中對炮點正問題的重復計算； ②對炮循環(huán)利用消息傳遞接口(MPI)進行并行，內部空間方向成像利用OpenMP線程并行，利用計算機單節(jié)點多核的特點節(jié)省內存； ③一次輸入單炮數(shù)據(jù)，成像結果記錄于當前計算節(jié)點的本地盤中，最后通過并行規(guī)約(Reduce)統(tǒng)一抽取成像道集，提高節(jié)點本地存儲的I/O利用率。

4 數(shù)值算例

首先對水平海底模型進行帶限射線傳播的測試，其中海底以下為高速(4000m/s)介質，用水平海底模擬速度躍變界面。在該模型坐標(1200m，400m)處設置震源，從震源出發(fā)進行帶限射線束傳播(射線出射角度為25°)，并以有限差分正演模擬的波場快照作為參考(主頻為20Hz，時刻為1.0s)。如圖4所示，傳統(tǒng)的高頻射線無法穿過速度躍變海底，帶限射線卻能穿過速度躍變的海底，帶限射線寬度由綠色實線表示，帶限射線的1.0s旅行時對應坐標位置和局部波前分別由紅色小圓圈和紅色實線表示。帶限射線穿透能力較高頻射線強，且隨著帶限射線頻率的降低，其穿透能力增強。對比波場快照，帶限射線旅行時與有限差分計算的波前能準確地吻合。

圖4 水平海底模型高頻射線與帶限射線路徑比較(射線初射角度為25°)

(a)高頻射線；(b)、(c)、(d)依次代表高、中、低頻帶限射線。紅色實線表示中心射線，綠色實線標注射線束寬度， 紅色小圓圈及其對應紅色實線標注旅行時1．0s的局部波前位置，對應20Hz有限差分正演模擬1．0s波場快照

進一步，為驗證帶限射線束進行波傳播的有效性和準確性，同時驗證帶限射線束成像的效果，設計一個崎嶇海底模型進行測試。在圖5a中，以正弦函數(shù)型界面模擬速度躍變的崎嶇海底，對于高頻射線，崎嶇海底下方是典型的射線陰影區(qū)。如圖5b和圖5c所示，在模型坐標(4000m，0)處設置震源，高頻射線和帶限射線的射線路徑疊加在有限差分正演模擬波場快照(主頻為20Hz，時刻為1.2s)的背景上。對比崎嶇海底模型中的射線分布可以看出，帶限射線較傳統(tǒng)高頻射線具有更強的穿透能力，改善了崎嶇海底下方的照明，帶限射線也能更好地與波場快照吻合。圖5d和圖5e分別為傳統(tǒng)射線束和帶限射線束的成像結果，相比于基于高頻射線的射線束偏移結果，帶限射線束偏移結果提高了目標層的連續(xù)性，其中紅色箭頭標注的水平目標層的成像能量強且均勻。在偏移效率方面，帶限射線束和傳統(tǒng)射線束由于在相同的實現(xiàn)框架下，因此理論上有類似的計算效率。其執(zhí)行效率的差別源于正算子中計算格林函數(shù)效率不同，單炮測試中，帶限射線束傳播算子相比傳統(tǒng)射線束傳播算子增加計算量達373.7%。但由于該模塊并不是偏移過程中的熱點，因此其導致的偏移效率計算量只增加了不到18%，詳細數(shù)據(jù)如表1所示。

最后，選擇更一般的Salt Bag模型[31]進行鹽下成像的測試，圖6a所示模型包含不規(guī)則鹽丘體，鹽丘上界面為一崎嶇界面。如圖6b和圖6c所示，在模型坐標(5000m，0)處設置震源，高頻射線和帶限射線的射線路徑疊加在有限差分正演模擬波場快照(主頻為20Hz，時刻為1.2s)的背景上。對比高頻射線與帶限射線的射線路徑可知，帶限射線能更好地穿透鹽丘，對鹽下區(qū)域有更好的照明，這是利用改善射線類方法對鹽下成像的基礎。圖6d和圖6e分別為對應的傳統(tǒng)射線束和帶限射線束的偏移成像結果，對比紅色矩形框內的成像區(qū)域，帶限射線束偏移對鹽丘下界面和鹽丘下方的目標層位成像效果優(yōu)于傳統(tǒng)高頻射線束偏移。

圖5 崎嶇海底模型高頻射線與帶限射線傳播對比以及對應射線束成像結果

(a)崎嶇海底模型，海底下方為高速(4500m/s)鹽丘體； (b)高頻射線路徑疊加于20Hz有限差分正演模擬1.2s波場快照，紅色實線表示射線路徑，綠色實心圓點射線路徑對應1.2s旅行時的波前，圖c同； (c)帶限射線路徑疊加于20Hz有限差分正演模擬1.2s波場快照； (d)基于高頻射線的射線束成像剖面； (e)帶限射線束成像剖面

圖6 Salt Bag模型高頻射線與帶限射線傳播對比以及對應射線束成像結果

(a)Salt Bag速度模型； (b)高頻射線路徑疊加于20Hz有限差分正演模擬1.2s波場快照，紅色實線表示射線路徑，綠色實心圓點射線路徑對應1.2s旅行時的波前，圖c同； (c)帶限射線路徑疊加于20Hz有限差分正演模擬1.2s波場快照； (d)基于高頻射線的射線束成像剖面； (e)帶限射線束偏移成像剖面

表1 崎嶇海底模型傳統(tǒng)射線束與帶限射線束偏移效率對比

5 結論與討論

本文在射線理論框架下，基于帶限中心射線，利用旁軸近似構建帶限局部平面波傳播算子，進行帶限地震波場的傳播和偏移。該算子在理論上保留了傳統(tǒng)射線束傳播算子靈活高效和具有局部時空性的優(yōu)點，同時兼顧波動類算子能夠刻畫帶限信號的特點。在實際應用中帶限局部平面波傳播能提高射線束的穿透能力，改善傳統(tǒng)高頻射線理論對應的射線陰影區(qū)等問題，有利于改善復雜構造的照明?；谠搸迋鞑ニ阕拥某上穹椒ㄌ岣吡似閸绾５缀望}丘等復雜構造的偏移成像質量。數(shù)值實驗結果也表明本文提出的帶限局部平面波傳播算子能夠準確有效地描述帶限波場的傳播，改善射線陰影區(qū)和鹽下構造等復雜探區(qū)的成像質量。本文研究為基于帶限射線束的層析和反演等后續(xù)研究提供了基礎。

[1] Babich V M and Buldyrev V S.Asymptotic Methods in Problems of Diffraction of Short Waves.Nauka，Moscow,1972.

[2] Woodward M J.Wave-equation tomography.Geo-physics，1992，57(1)：15-26.

[3] Biondi B.Solving the frequency dependent eikonal equation.SEG Technical Program Expanded Abstracts，1992，11：1315-1319.

[4] Hogan C M，Margrave G F.Ray-tracing and eikonal solutions for low-frequency wavefields.CREWES Research Report，2007.

[5] Foreman T.A Frequency Dependent Ray Theory[D].The University of Texas at Austin，1987.

[6] Protasov M，Gadylshin K.Exact frequency dependent rays on the basis of Helmholtz solver.SEG Technical Program Expanded Abstracts，2015，34：3739-3743.

[7] Lomax A.The wavelength-smoothing method for approximating broad-band wave propagation through complicated velocity structures.Geophysical Journal International，1994，117(2)：313-334.

[8] Lomax A，Snieder R.Estimation of finite-frequency waveforms through wavelength-dependent averaging of velocity.Geophysical Journal International，1996，126(2)：369-381.

[9] Protasov M I，Yarman C E，Nichols D et al.Frequency-dependent ray-tracing through rugose interfaces.SEG Technical Program Expanded Abstracts，2011，30：2992-2996.

[10] Protasov M I，Osypov K S.Frequency dependent ray tracing for irregular boundaries.Seismic Technology,2014,11(3):1-11.

[11] Yarman C E，Cheng X，Osypov K et al.Band-limited ray tracing.Geophysical Prospecting，2013，61(6)：1194-1205.

[13] Hill N R.Gaussian beam migration.Geophysics，1990，55(11)：1416-1428.

[14] Hill N R.Prestack Gaussian-beam depth migration.Geophysics，2001，66(4)：1240-1250.

[15] Hale D.Migration by the Kirchhoff，slant stack，and Gaussian beam methods.CWP Annual Project Review Meeting，Colorado，1992.

[16] Gray S H.Gaussian beam migration of common-shot records.Geophysics，2005，70(1)：953-959.

[17] 鄧飛，劉超穎，趙波等.高斯射線束正演與偏移.石油地球物理勘探，2009，44(3)：265-269. Deng Fei，Liu Chaoying，Zhao Bo et al.Gaussian beams forward simulation and migration.OGP，2009，44(3)：265-269.

[18] 李振春，岳玉波，郭朝斌等.高斯波束共角度保幅深度偏移.石油地球物理勘探，2010，43(3)：360-365. Li Zhenchun，Yue Yubo，Guo Chaobin et al.Gaussian beam common angle preserved-amplitude migration.OGP，2010，45(3)：360-365.

[19] Zhu T，Gray S H，Wang D.Prestack Gaussian-beam depth migration in anisotropic media.Geophysics，2007，72(3)：S133-S138.

[20] 韓建光，王赟，張曉波等.VTI介質高斯束疊前深度偏移.石油地球物理勘探，2015，50(2)：267-273. Han Jianguang，Wang Yun，Zhang Xiaobo et al.Gaussian beam prestack depth migration in VTI media.OGP，2015，50(2)：267-273.

[21] 劉強，張敏，李振春等.各向異性介質共炮域高斯束偏移.石油地球物理勘探，2016，51(5)：930-937. Liu Qiang，Zhang Min，Li Zhenchun et al.Common-shot domain Gaussian beam migration in anisotropic media.OGP，51(5)：930-937.

[22] 代福材，黃建平，李振春等.角度域黏聲介質高斯束疊前深度偏移方法.石油地球物理勘探，2017，52(2)：283-293. Dai Fucai，Huang Jianping，Li Zhenchun et al.Angle domain prestack Gaussian beam migration for visco-acoustic media.OGP，2017，52(2)：283-293.

[23] Hu C S，Stoffa P L.Slowness-driven Gaussian-beam prestack depth migration for low-fold seismic data.Geophysics，2009，74(6)：WCA35-WCA45.

[24] Liu J，Palacharla G.Multiarrival Kirchhoff beam migration.Geophysics，2011，76(5)：WB109-WB118.

[25] Bube K P，Washbourne J K.Wave tracing：Ray tracing for the propagation of band-limited signals：Part 1 —Theory.Geophysics，2008，73(5)：VE377-VE384.

[26] Washbourne J K，Bube K P，Carillo P et al.Wave tra-cing：Ray tracing for the propagation of band-limited signals：Part 2 — Applications.Geophysics，2008，73(5)：VE385-VE393.

[27] Gray S H，Xie Y，Notfors C et al.Taking apart beam migration.The Leading Edge，2009，28(9)：1098-1108.

[28] Elad M.Sparse and Redundant Representations：from Theory to Applications in Signal and Image Processing.Springer Science & Business Media，2010.

[29] Liu S Y，Gu H and Feng B.Characteristic wave imaging in TTI media via compressed sensing.SEG Technical Program Expanded Abstracts，2016，35：4424-4429.

[30] Schneider W A.Integral formulation for migration in two and three dimensions.Geophysics，1978，43(1)：49-76.

[31] Popov M M，Semtchenok N M，Popov P M et al.Depth migration by the Gaussian beam summation method.Geophysics，2010，75(2)：S81-S93.

(本文編輯：朱漢東)

1000-7210(2017)05-0948-08

P631

A

10.13810/j.cnki.issn.1000-7210.2017.05.007

劉少勇 博士，1985年生； 2008年畢業(yè)于中國石油大學(華東)地球物理學專業(yè)，獲學士學位； 2015年畢業(yè)于同濟大學固體地球物理學專業(yè)，獲博士學位； 2015年11月至今，在中國地質大學(武漢)從事教學科研工作，主要研究方向為地震波傳播及成像、數(shù)字信號處理。

*湖北省武漢市洪山區(qū)魯磨路388號中國地質大學(武漢)地球物理與空間信息學院，430074。 Email：rhanbk@163.com

本文于2016年11月28日收到，最終修改稿于2017年8月14日收到。

本項研究受國家重大專項(2016ZX05026001-001)、國家自然科學基金(41604092)、中國博士后科學基金(2016M590725)及中國石化國家重點實驗室開放基金(33550006-16-FW0399-0021)等聯(lián)合資助。