基于光流矢量行波分離的相關加權逆時偏移成像

解闖， 宋鵬，2，3*， 鄒志輝，2，3， 譚軍，2，3， 王紹文， 趙波，2，3

1 中國海洋大學海洋地球科學學院， 青島 266100 2 青島海洋科學與技術試點國家實驗室， 海洋礦產(chǎn)資源評價與探測技術功能實驗室， 青島 266100 3 中國海洋大學海底科學與探測技術教育部重點實驗室， 青島 266100

0 引言

逆時偏移起源于20世紀80年代(Whitmore，1983；Baysal et al.，1983；McMechan，1983)，其基于雙程波波動方程成像，能適應任意復雜速度模型，無成像角度限制，且理論上還可對回折波、棱柱波、繞射波、多次波等進行成像(楊仁虎，2021a)，是當前公認的高精度成像算法之一.然而，由于逆時偏移需要基于雙程波方程進行波場延拓，其在地震波傳播到反射界面時會產(chǎn)生背向反射，而常規(guī)互相關成像條件不加以區(qū)分地將所有的震源正傳波場與檢波點反傳波場直接應用于成像(杜啟振等，2013；Chen and He，2014；楊仁虎，2021b)，這會不可避免地產(chǎn)生大量的低波數(shù)偏移噪聲.

目前主要發(fā)展了三類方法用于壓制逆時偏移中的偏移噪聲.第一類是壓制背向反射方法，其通常采用無反射聲波方程進行成像.Baysal等(1984)首先提出了基于常波阻抗假設的無反射聲波方程，可顯著壓制垂向附近入射的地震波背向反射；宋鵬(2005)和Zhang等(2010)改進了無反射聲波方程，提升了背向反射的壓制效果.然而總體而言，這類方法的背向反射壓制效果尚不夠理想，難以達到有效消除偏移噪聲的目的.第二類是濾波方法，Mulder和Plessix(2004)直接采用高通濾波方法對成像剖面進行去噪處理，該方法易于實現(xiàn)但合適的濾波器閾值范圍較難選?。籞hang和Sun(2009)采用Laplace濾波方法對常規(guī)逆時偏移結果進行濾波，該方法對于低波數(shù)噪聲通常有很好的壓制效果，但其在去噪的同時也易損失部分低波數(shù)有效信息，且會引入高頻噪聲(杜啟振等，2013).第三類是選擇性成像方法，主要包括角道集疊加壓噪方法與行波分離方法.角道集疊加壓噪方法需對地下每個成像點進行角道集的求取，且需選擇適當?shù)慕嵌乳撝挡π∮谶@一角度的所有成像值進行疊加，因此該方法計算代價較高(杜啟振等，2013).相比角道集疊加壓噪方法，行波分離方法無需進行角道集的求取，而是直接通過修改成像條件的方式(將震源波場和檢波點波場分離成不同方向的行波，然后提取有效波場分量參與成像)來壓制低波數(shù)噪聲，其可高效地實現(xiàn)地下構造的精確成像，是一種較為理想的低波數(shù)噪聲壓制算法.Liu等(2011)、Fei等(2015)、胡江濤和王華忠(2015)、王一博等(2016)應用Hilbert變換實現(xiàn)了行波分離逆時偏移成像，有效地壓制了偏移噪聲.Chen和He(2014)采用Poynting矢量實現(xiàn)了震源波場與檢波點波場上、下、左、右四個方向的行波分離，其可在付出極小的計算代價的前提下有效提高低波數(shù)噪聲的壓制效果，因此該方法受到了專家學者們的廣泛關注(王鵬飛和何兵壽，2017； 王曉毅等，2021).

然而，當前基于Poynting矢量行波分離的逆時偏移算法依然存在如下問題：首先，Poynting矢量難以精確指示波場的傳播方向信息，且其計算過程容易出現(xiàn)不穩(wěn)定的現(xiàn)象，因此，基于Poynting矢量分離得到的各個方向的行波并不完全精確(Du et al.，2012；Zhang，2014；Duan and Sava，2015；Li and He，2020)；此外，對于行波分離后得到的多個成像剖面，諸多專家學者多是進行人為的選擇性組合(Chen and He，2014；王一博等，2016；Xue and Liu，2018)，且通常賦予各行波分離成像剖面一個經(jīng)驗權重，未能實現(xiàn)各個剖面的高精度疊加成像.

光流法最早應用于解決相鄰幀之間物體的運動信息問題(Horn and Schunck，1981；Lucas and Kanade，1981)，隨后被引入到逆時偏移角道集提取中(Zhang，2014；Gong et al.，2016；吳成梁等，2021).相比于Poynting矢量，光流矢量是一個經(jīng)過多次迭代運算后得到的更逼近于真實波場傳播方向信息的精確矢量，且在光流矢量的計算過程中還額外加入了正則化項，這也有效避免了由于計算得到的光流矢量值為0而出現(xiàn)的不穩(wěn)定現(xiàn)象.本文首先將光流矢量引入到基于行波分離的逆時偏移壓噪方法中，并基于光流矢量精確穩(wěn)定地實現(xiàn)了震源波場和檢波點波場的行波分離，由此分別得到了震源波場和檢波點波場的多方向行波，進而由每兩個行波互相關得到16個成像剖面，在此基礎上，將這16個成像剖面分別與參考剖面做相關計算，并將相關值作為權重分別賦予各成像剖面，最終實現(xiàn)了基于光流矢量行波分離的相關加權逆時偏移成像，顯著提升了逆時偏移的成像效果.

1 逆時偏移的波場延拓計算

在二維各向同性介質中，一階應力—速度聲波方程的表達式為：

(1)

式中，x、z分別為空間坐標，vx、vz分別為質點在x和z方向的振動速度，p為應力，t為時間，ρ為密度，v為聲波速度.

采用交錯網(wǎng)格對式(1)進行有限差分離散，從而實現(xiàn)正向延拓和逆時延拓.以正向延拓為例，式(1)差分離散后的高階差分格式為：

(2)

式中，k為時刻，i、j分別為x、z方向的空間離散點序號，N為差分階數(shù)的一半，Δt為時間采樣間隔，Δx、Δz為空間采樣間隔，Cm為差分系數(shù).

采用有限差分法進行波場延拓必須引入人工邊界來界定計算區(qū)域，為較好地壓制人工邊界反射，這里采用完全匹配層(PML)方法.PML邊界算法已有諸多學者進行研究(Berenger，1994；Collino and Tsogka，2001；Zhang and Shen，2010)，本文在此不再贅述.

2 基于光流矢量的行波分離

Poynting矢量又叫能流密度矢量，最早應用于電磁計算領域(Poynting，1884)，現(xiàn)在已成為地震波場計算中用于指示波場傳播方向的常用算法(Yoon and Marfurt，2006).

一階應力—速度聲波方程的Poynting矢量為：

(3)

(4)

式中，Sl(x,z,t)、Sr(x,z,t)、Su(x,z,t)、Sd(x,z,t)分別為震源波場的左、右、上、下行波.由式(3)可知，Poynting矢量的計算是由波場的時間導數(shù)與空間導數(shù)的乘積組成，當時間導數(shù)或空間導數(shù)為0時，Poynting矢量值即為0，會出現(xiàn)不穩(wěn)定的現(xiàn)象；此外，Zhang(2014)指出，Poynting矢量本身難以高精度地指示出波場的傳播方向信息.

光流矢量是一個經(jīng)過多次迭代運算后得到的，可穩(wěn)定且精確地指示波場傳播方向的矢量，本文將光流矢量引入逆時偏移的行波分離處理中.在二維逆時偏移中，光流問題的基本假設可以描述為鄰近點的相鄰時刻的波場是連續(xù)變化的，即：

u(x+dx,z+dz,t+dt)=u(x,z,t)，

(5)

式中，u為波場，x、z分別為空間坐標，將u(x+dx,z+dz,t+dt)進行泰勒展開并舍去二階以上的高階項，可得：

uxq+uzw+ut=0，

(6)

式中，ux、uz均為波場的空間導數(shù)，ut為波場的時間導數(shù),q和w分別為光流矢量的水平分量和垂直分量.

由于式(6)為欠定方程，需額外引入約束條件進行求解.本文引入全局平滑約束的正則化項(Horn and Schunck，1981)，并構建如下的誤差泛函：

E=?[(uxq+uzw+ut)2+α2C]dxdz，

(7)

式(7)可通過梯度法進行極小化求解，可得：

(8)

(9)

(10)

求解式(9)，可得：

(11)

然后根據(jù)第n次迭代得到的光流矢量(qn,wn)可計算出第n+1次迭代的光流矢量(qn+1,wn+1)，即為：

(12)

式中，n為計算光流矢量的迭代次數(shù)，其在本文模型實驗中的值為10.

由式(12)可知，光流矢量不再是簡單地以波場的時間導數(shù)與空間導數(shù)相乘的方式來直接求取波場傳播方向信息，而是通過給定一個初始光流矢量，經(jīng)過多次迭代運算后得到的一個更逼近于真實波場傳播方向信息的精確矢量，且由于在光流矢量的計算過程中還額外加入了正則化項，只有當時間導數(shù)與空間導數(shù)均為0時，光流矢量值才為0，其有效避免了光流矢量計算過程中的不穩(wěn)定情況.因此，從理論上來說基于光流矢量分離得到的各個方向的行波更為精確.

本文采用簡單速度模型來驗證基于光流矢量實現(xiàn)行波分離的精確性.圖1為雙層速度模型，第一層速度為2500 m·s-1，第二層速度為3000 m·s-1，該模型橫縱向均為1500 m，空間采樣間隔為5 m，震源位于(750 m，0 m)處，利用時間二階、空間十階交錯網(wǎng)格有限差分進行波場延拓，震源采用主頻30 Hz的雷克子波，時間采樣步長為0.5 ms，記錄時長為0.8 s.圖2為0.3 s時刻的震源波場快照.圖3a、b分別為基于Poynting矢量、光流矢量計算的反射界面附近(圖2紅框處)的波場方向.圖4a、b分別為該時刻Poynting矢量以及光流矢量的水平分量.圖5a、b分別為基于Poynting矢量、光流矢量實現(xiàn)波場分離后的左下行波.

圖1 雙層速度模型Fig.1 A two-layer velocity model

圖2 0.3 s的震源波場快照Fig.2 The source wavefield snapshot at 0.3 s

由圖3a、b(紅圈處)與圖4a、b(紅色箭頭處)可知，Poynting矢量難以精確指示波場的傳播方向信息，且易出現(xiàn)奇異值，而光流矢量更為光滑，且有效避免了不穩(wěn)定現(xiàn)象.對比圖5a、b可知，基于Poynting矢量雖可實現(xiàn)行波分離，但由于Poynting矢量計算不精確、不穩(wěn)定，其在實現(xiàn)過程中會引入如箭頭所指處的其他行波分量，而基于光流矢量實現(xiàn)行波分離不會引入其他行波分量，分離得到的行波更為精確.

圖3 波場方向(a) Poynting矢量； (b) 光流矢量.Fig.3 Wavefield direction(a) Poynting vector； (b) Optical flow vector.

圖4 水平分量(a) Poynting矢量； (b) 光流矢量.Fig.4 Horizontal component(a) Poynting vector； (b) Optical flow vector.

圖5 波場分離后的左下行波(a) Poynting矢量； (b) 光流矢量.Fig.5 Left-down going wavefield after decomposition(a) Poynting vector； (b) Optical flow vector.

3 基于參考剖面的相關加權逆時偏移成像

本文基于光流矢量，應用式(13)實現(xiàn)左上、左下、右上、右下四個方向的行波分離：

(13)

式中，Slu(x,z,t)、Sld(x,z,t)、Sru(x,z,t)、Srd(x,z,t)分別為震源波場的左上、左下、右上、右下行波.Rlu(x,z,t)、Rld(x,z,t)、Rru(x,z,t)、Rrd(x,z,t)分別為檢波點波場的左上、左下、右上、右下行波.根據(jù)行波分離結果，根據(jù)式(14)可進一步得到相應的16個歸一化互相關成像剖面：

(14)

式中，I1(x,z)、I2(x,z),…,I16(x,z)分別為相應的16個行波分離成像剖面，分別如圖6—圖9所示.

圖6 左上類型的震源波場與各檢波點波場互相關成像結果(a) SluRld； (b) SluRlu； (c) SluRrd； (d) SluRru.Fig.6 The RTM images constructed by the cross-correlation of the left-upgoing type source wavefield with all decomposed receiver wavefields

圖7 左下類型的震源波場與各檢波點波場互相關成像結果(a) SldRld； (b) SldRlu； (c) SldRrd； (d) SldRru.Fig.7 The RTM images constructed by the cross-correlation of the left-downgoing type source wavefield with all decomposed receiver wavefields

圖8 右上類型的震源波場與各檢波點波場互相關成像結果(a) SruRld； (b) SruRlu； (c) SruRrd； (d) SruRru.Fig.8 The RTM images constructed by the cross-correlation of theright-upgoing type source wavefield with all decomposed receiver wavefields

圖9 右下類型的震源波場與各檢波點波場互相關成像結果(a) SrdRld； (b) SrdRlu； (c) SrdRrd； (d) SrdRru.Fig.9 The RTM images constructed by the cross-correlation of the right-downgoing type source wavefield with all decomposed receiver wavefields

當前對于各個方向的正時、逆時波場相關成像結果，通常舍去同方向震源波場與檢波點波場的互相關成像剖面(即I2(x,z)、I5(x,z)、I12(x,z)、I15(x,z))，并對保留的行波分離成像剖面分別賦予一個經(jīng)驗權重.事實上，除I2(x,z)、I5(x,z)、I12(x,z)以及I15(x,z)這四個成像剖面外，其余的剖面上依然存在低波數(shù)噪聲，因此簡單舍去這四個剖面的做法不能最大限度地剔除偏移噪聲，且存在丟失有效信息的風險.為解決該問題，本文首先計算得到一個參考剖面Iref(x,z)(即將保留的成像剖面直接疊加)，其可表示為：

Iref(x,z)=I1(x,z)+I3(x,z)+I4(x,z)

+I6(x,z)+I7(x,z)+I8(x,z)

+I9(x,z)+I10(x,z)+I11(x,z)

+I13(x,z)+I14(x,z)+I16(x,z)，(15)

然后將式(14)得到的16個成像剖面分別與參考剖面做相關計算：

(16)

式中l(wèi)為1至16的序號，w為各波場分離的成像剖面與參考剖面的互相關值，求取這16個互相關值中的最大值wmax，可得：

wmax=max(wl)，

(17)

對互相關值w做歸一化處理：

(18)

式中，W為最終的相關值，并將所得的相關值W作為權重分別賦予各波場分離的成像剖面，得到最終的基于光流矢量行波分離的相關加權逆時偏移成像結果：

(19)

這里依然采用第2節(jié)的簡單速度模型來測試加權成像方法的效果.圖10a為常規(guī)逆時偏移的結果，圖10b為基于Poynting矢量行波分離的逆時偏移結果，圖10c為基于光流矢量行波分離的逆時偏移結果.本文以圖10c作為參考剖面，基于該參考剖面得到各個行波分離成像剖面的相關加權值分別為：0.0、0.0、0.0、0.0、0.0、0.99、0.01、0.002、0.0、0.0、0.0、0.0、0.013、0.002、0.0、1.0.基于光流矢量行波分離的相關加權逆時偏移結果見圖10d.圖10e為Laplace濾波方法逆時偏移結果.圖11為相應的波數(shù)譜.

由圖10、圖11可知，常規(guī)逆時偏移成像結果中存在明顯的低波數(shù)偏移噪聲.在基于Poynting矢量行波分離的逆時偏移結果中，大部分的低波數(shù)偏移噪聲得到壓制，但由于行波分離不精確，其仍有部分噪聲殘留；在基于光流矢量行波分離的逆時偏移結果中，低波數(shù)偏移噪聲得到了進一步的消除；通過對參考剖面相關加權的處理，低波數(shù)偏移噪聲基本完全消除，而有效的構造成像得以凸顯；從成像剖面上看，Laplace濾波方法雖然對低波數(shù)噪聲的壓制效果很好，但從波數(shù)譜中可以發(fā)現(xiàn)其處理結果中低波數(shù)信息損傷嚴重(如圖11e中的紅圈所示)，存在損傷低波數(shù)有效信息的風險；此外，該方法也會引入部分的高波數(shù)噪聲(如圖11e中的紅色箭頭所示).

圖10 逆時偏移的結果(a) 常規(guī)逆時偏移； (b) 基于Poynting矢量行波分離的逆時偏移； (c) 基于光流矢量行波分離的逆時偏移； (d) 基于光流矢量行波分離的相關加權逆時偏移； (e) Laplace濾波方法逆時偏移結果.Fig.10 The results of reverse time migration(a) Conventional reverse time migration； (b) Reverse time migration with wavefield decomposition based on Poynting vector； (c) Reverse time migration with wavefield decomposition based on optical flow vector； (d) Cross-correlation weighted reverse time migration with wavefield decomposition based on optical flow vector; (e) Reverse time migration using Laplace filtering.

圖11 波數(shù)譜(a) 圖10a對應的波數(shù)譜； (b) 圖10b對應的波數(shù)譜； (c) 圖10c對應的波數(shù)譜； (d) 圖10d對應的波數(shù)譜； (e) 圖10e對應的波數(shù)譜.Fig.11 Wavenumber spectra(a) Wavenumber spectrum corresponding to Fig.10a; (b) Wavenumber spectrum corresponding to Fig.10b; (c) Wavenumber spectrum corresponding to Fig.10c; (d) Wavenumber spectrum corresponding to Fig.10d; (e) Wavenumber spectrum corresponding to Fig.10e.

4 Marmousi模型實驗

本文采用Marmousi-II局部速度模型進一步進行成像效果測試.速度模型橫縱向分別為6500 m、3500 m(如圖12所示)，空間采樣間隔為5 m.成像所用數(shù)據(jù)共101炮，從模型左端開始放炮，炮點兩側各260道接收，炮間隔為65 m，道間隔為5 m，炮點和檢波點均在地表，深度為0 m.利用時間二階、空間十階交錯網(wǎng)格有限差分進行波場延拓，震源采用主頻30 Hz的雷克子波，時間采樣步長為0.4 ms，記錄時長為4 s.圖13a為常規(guī)逆時偏移的結果，圖13b為基于Poynting矢量行波分離的逆時偏移結果，圖13c為基于光流矢量行波分離的逆時偏移結果，圖13d為基于光流矢量行波分離的相關加權逆時偏移結果，圖13e為Laplace濾波方法逆時偏移結果.圖14為相應的波數(shù)譜.

圖12 Marmousi-II局部速度模型Fig.12 The local velocity model of Marmousi-II model

由圖13a可知，偏移噪聲嚴重影響了地層的成像質量，尤其在淺部地層和強波阻抗界面處，真實的成像構造被完全掩蓋，成像精度較低.由圖13b、c、d可知，三種行波分離方法均可起到顯著的噪聲壓制效果，且相比常規(guī)行波分離方法，本文方法的噪聲壓制效果更佳，地下構造更清晰.由圖13e可知，Laplace濾波方法雖可顯著壓制偏移噪聲，然而從波數(shù)譜(圖14)中可以發(fā)現(xiàn)，該方法在壓制低波數(shù)噪聲的同時極易損傷低波數(shù)有效信息(如圖14e中的紅圈所示)，且會引入部分高波數(shù)噪聲，因此從保護有效信息的角度來講，Laplace濾波方法需謹慎使用或可作為其他成像算法的一個后續(xù)去噪處理手段.

此外，本文對Poynting矢量法和光流矢量法進行單炮逆時偏移計算效率測試，其結果如表1所示(本實驗應用的GPU型號為Tesla V100).

由表1可知，本文方法的計算時間約為Poynting矢量方法的1.85倍，隨著以GPU集群為代表的高性能計算機集群在地球物理領域的應用與普及，逆時偏移等成像算法通常都會有幾十倍的計算效率提升，因此當前計算效率問題已不會成為限制本文算法在實際應用中的瓶頸.

5 實際資料試算

以東海某區(qū)域地震數(shù)據(jù)為例，測試本文算法的成像效果.該測線采用海洋拖攬單邊放炮的采集方式，炮點位于檢波點的右側.成像所用數(shù)據(jù)共1612炮，每炮648道接收.炮間隔為37.5 m，道間隔為12.5 m，最小偏移距為187.5 m.炮點和檢波點的深度為12.5 m.利用時間二階、空間八階交錯網(wǎng)格有限差分進行波場延拓，時間采樣步長為1 ms，記錄時長為8 s.圖15為該區(qū)域基于全波形反演方法得到的速度模型，橫向為58.75 km，縱向為3.5 km.圖16a為常規(guī)逆時偏移的結果，圖16b為基于Poynting矢量行波分離的逆時偏移結果，圖16c為基于光流矢量行波分離的逆時偏移結果，圖16d為基于光流矢量行波分離的相關加權逆時偏移結果.

圖15 實際資料的速度模型Fig.15 Velocity model of field data

圖16 逆時偏移的結果(a) 常規(guī)逆時偏移； (b) 基于Poynting矢量行波分離的逆時偏移； (c) 基于光流矢量行波分離的逆時偏移； (d) 基于光流矢量行波分離的相關加權逆時偏移.Fig.16 The results of reverse time migration(a) Conventional reverse time migration； (b) Reverse time migration with wavefield decomposition based on Poynting vector； (c) Reverse time migration with wavefield decomposition based on optical flow vector； (d) Cross-correlation weighted reverse time migration with wavefield decomposition based on optical flow vector.

由圖16a可知，在淺部地層處，低波數(shù)偏移噪聲掩蓋了真實的成像構造(如圖中的黑圈所示).在圖16b中，大部分的低波數(shù)偏移噪聲得到壓制.相比圖16b，圖16c中的偏移噪聲更少.而在圖16d中，低波數(shù)偏移噪聲最少.由此可知，本文方法同樣可以較好地實現(xiàn)實際資料的高精度成像.

6 結論與展望

本文首先基于光流矢量實現(xiàn)了震源波場和檢波點波場的行波分離，并分別得到分離后的各方向行波及相應的16個成像剖面，在此基礎上，將這16個成像剖面分別與參考剖面做相關計算，并將相關值作為權重分別賦予各成像剖面，最終實現(xiàn)了基于光流矢量行波分離的相關加權逆時偏移成像.理論模型數(shù)據(jù)和實際資料試算表明：

(1)光流矢量可精確、穩(wěn)定地實現(xiàn)行波分離，基于光流矢量行波分離的逆時偏移可有效壓制偏移噪聲.

(2)相關加權處理避免了過多的人為干預，且可進一步壓制偏移噪聲、凸顯有效的地下構造成像，從而顯著提升逆時偏移的成像精度.

需要注意的是，與Poynting矢量相同，本文引入的光流矢量同樣難以精確指示多波疊加情況下的行波傳播方向.Tang和McMechan(2016)、Tang 等(2017)提出了一種多方向矢量方法，較好地實現(xiàn)了多行波傳播方向的精確指示，但其需要付出較大的計算消耗，因此，進一步發(fā)展高效高精度的多傳播方向指示算法將是下一步的主要研究工作.