劉夢麗 黃建平* 李 闖 崔 超 任英俊(①中國石油大學(華東)地球科學與技術學院，山東青島 266580；②海洋國家實驗室海洋礦產資源評價與探測技術功能實驗室，山東青島 266071)

1 引言

近年來，地震勘探對構造成像精度的要求越來越高，傳統(tǒng)的偏移方法如Kirchhoff偏移、相移法等難以對橫向速度變化劇烈的復雜構造準確成像。逆時偏移(RTM)對各種地震波型(一次波、多次波)和高陡構造成像具有明顯優(yōu)勢[1]。成像條件是影響逆時偏移成像效果的主要因素，目前基于波場零延遲互相關成像條件得到廣泛使用[2]，但是該條件忽略了地震波場的多路徑和方向性特征。因此RTM直接使用該條件時會不可避免地引入低波數(shù)強振幅噪聲，在淺層尤為明顯[3-5]。業(yè)界為解決這一問題提出了許多方法，其中高通濾波是一種常用的去噪方法,但忽略了波場的多路徑特點，造成有效反射信號被破壞[6]。Guitton等[7]在系統(tǒng)對比導數(shù)濾波與拉普拉斯濾波方法的基礎上提出了最小平方濾波方法，在去噪的同時也能保存有效信息。此外，低頻成像噪聲主要由強反射界面反射波引起，在計算區(qū)域的噪聲位置采用方向衰減的邊界條件，利用無反射方程也能抑制噪聲[8]。

在成像過程中消除噪聲的思路成為業(yè)界關注的熱點問題。坡印廷矢量去噪方法首先由坡印廷矢量的定義得到波場傳播方向，然后對成像條件參數(shù)乘以一個方向衰減因子進行去噪，不但需要額外計算和存儲波場傳播矢量，而且在波場復雜區(qū)域提取的波場傳播方向誤差較大[9]； 此外，通過把波場分解為(波數(shù)域)上、下、左、右行波，運用成像條件時只保留方向相反的波場分量成像結果，這種方法可以有效去噪，但波場分解增加了額外計算量，且分解三維波場可能會損失有效信息[10]。Rocha[11]發(fā)展了能量范數(shù)成像條件，把成像條件看作震源波場與檢波點波場之間的內積，可壓制常規(guī)互相關造成的低頻噪聲。OP’t Root等[12]、Whitmore等[13]、Kiyashchenko等[14]根據(jù)Stolk提出的逆散射理論，結合入射波和反射波的時間導數(shù)成像條件和空間梯度成像條件壓制低頻噪聲。

20世紀80年代，Tarantola[15]首先提出構建最小二乘的誤差泛函方法求解反演問題。Nemeth等[16]首先使用Kirchhoff最小二乘偏移(LSM)方法對不規(guī)則反射地震數(shù)據(jù)進行偏移。LSM將偏移成像看作是最小二乘反演框架下的反演問題，能夠克服常規(guī)RTM對儲層成像分辨率不高且振幅不均衡等問題，且能在一定程度上消除偏移假象[17-21]。因此，人們將RTM與LSM的優(yōu)勢相結合發(fā)展了最小二乘逆時偏移(LSRTM)進行復雜構造的保幅成像處理[22,23]。針對LSRTM計算成本高的問題，基于編碼策略的疊前平面波LSM成像能極大地提高計算效率[24,25]。

本文在前人研究的基礎上，首先分析了RTM中低頻噪聲產生的原因，提出了一種滿足能量守恒的逆散射成像條件，并與LSRTM結合，實現(xiàn)了基于角度濾波成像的最小二乘逆時偏移算法(ALSRTM)。最后，通過對SEG/EAGE二維鹽丘模型的成像測試驗證了ALSRTM 的有效性，并對比、分析了常規(guī)LSRTM及ALSRTM對低信噪比炮數(shù)據(jù)、含誤差速度場的適應性。

2 ALSRTM

2.1 基于互相關成像條件的常規(guī)RTM

RTM通常分為震源波場正向延拓、記錄波場逆時傳播和互相關成像三部分。其中時間域的正傳和反傳過程通過對雙程波波動方程

(1)

進行差分求解得到。式中：P(x,t,xs)為震源波場；Q(x,t,xr)為檢波點波場；fs(xs,t)為震源函數(shù)；dobs(xs,xr,t)為檢波點xr處在時刻t采集到的震源點xs處的地震記錄；v(x)為地震波在介質中的傳播速度。

應用傳統(tǒng)的互相關成像條件可得偏移剖面

(2)

圖1為低頻噪聲成像機制示意圖。采用式(2)進行互相關成像條件的前提假設是：炮點上行波和檢波點下行波場(圖1中虛線)為零，僅僅利用炮點下行波和檢波點上行波(圖1中實線)信息進行互相關并成像，可得到理想的反射層成像結果。根據(jù)地震反射機理，若存在強反射界面，由于RTM以雙程波波動方程為理論基礎，震源波場在波阻抗界面處必然產生上行反射波，類似地，檢波點波場在逆時延拓時也必然產生下行反射波，即同時存在炮點上行波場和檢波點逆時延拓過程中的下行波場，在滿足ts+tr等于地震道走時的位置均可成像，其成像結果為構造假象，嚴重影響成像精度。

圖1 低頻噪聲成像機制示意圖

式(2)代表的成像過程把波場看作標量，僅僅利用波場的數(shù)值成像，沒有區(qū)分波場方向，因此成像結果中包含炮點上行波場和檢波點逆時延拓過程中的下行波場造成的低頻噪聲。另外，從低頻噪聲成像機制示意圖(圖1)可以看出，在反射點(圖1中的位置1)處的反射角較小，而在產生低頻噪聲的位置(圖1中的位置2)震源波場方向與檢波點波場方向之間的角度接近180°。

2.2 角度濾波成像條件

對于單個震源，局部震源波場P(x,t)和檢波點波場Q(x,t)的射線參數(shù)分別為[13,26]

(4)

(5)

化簡可得

(6)

(7)

由式(7)得到的偏移剖面相當于從成像結果中去除了反射角為θa對應的成像值。由低頻噪聲成像機制示意圖(圖1)可見，在低頻噪聲產生位置(圖1中的位置2)的反射角θ大多在接近90°的大角度范圍內，而在反射層成像位置(圖1中的位置1)的反射角較小。對于正向傳播的震源波場和反向傳播的檢波點波場來說，兩者夾角主要集中在較小的區(qū)間內，在此區(qū)間內成像能突出有效波的能量；在低頻噪聲產生位置的反射角θ大多在接近90°的大角度范圍內，則近似認為僅產生噪聲，因此在此區(qū)間內的互相關不參與成像過程。因此，與坡印廷矢量成像去除低頻噪聲的思想類似，僅在θ較大的區(qū)間進行互相關(文中為了節(jié)省計算量，近似認為θ為90°)。

對比式(2)與式(7)可見：傳統(tǒng)互相關成像條件(式(2))只包含了波場的數(shù)值大小，刻畫了主要構造信息；能量成像條件(式(7))是波場的空間梯度和時間導數(shù)的線性疊加，速度是連接時間導數(shù)與空間梯度的“橋梁”，作為加權項作用于空間梯度，波場的局部時間導數(shù)則蘊含著波場的動力學特征(波形、振幅、頻率、相位)以及地下地層結構、巖石結構及流體性質之間的關系，可使成像過程更準確。

為了驗證新成像條件(式(7))的可靠性，從能量守恒的角度出發(fā)，定義各向同性介質聲波方程的能量為[27]

(8)

(9)

(10)

2.3 照明補償預條件算子

由于常規(guī)偏移算子僅為正演算子的轉置而不是逆，在成像時會受到地表觀測系統(tǒng)的影響，造成對地下構造的不均勻照明，尤其是深部及高速鹽體之下的構造，由于照明太弱而無法成像。因此，為了對深部構造更好地照明，本文使用炮照明算子進行補償，該算子在時間域可表示為[6]

(12)

2.4 基于逆散射能量保幅成像條件的LSRTM

為了得到與記錄數(shù)據(jù)最佳匹配的偏移結果，把LSRTM引入反演，其誤差函數(shù)定義為

(13)

式中：L為Born近似下的反偏移算子；m為偏移剖面；d為觀測數(shù)據(jù)。求解式(13)，第k次迭代得到的梯度為

gk=LT[Lmk-1-d]

(14)

本文利用共軛梯度法迭代求解式(13)，得到模擬數(shù)據(jù)與觀測數(shù)據(jù)差別最小的最佳反射系數(shù)模型，迭代過程使成像結果具有更高的分辨率和更可靠的振幅保真度。實現(xiàn)流程如下：

(1)輸入初始模型m0(對m0賦值0，即第一次迭代結果為逆時偏移剖面)、觀測數(shù)據(jù)dobs、梯度誤差閾值δ；

(2)計算第k次迭代的梯度gk、共軛方向修正因子β和共軛梯度zk

(3)計算步長α并對成像結果進行更新，得到

(4)判斷是否滿足終止迭代條件，若滿足則退出并保存成像結果，否則退回到步驟(2)、步驟(3)繼續(xù)迭代，直到滿足迭代終止條件，最后一次迭代即為最終的成像結果。

3 模型試算

為了分析ALSRTM與LSRTM的地震成像效果，對SEG/EAGE二維鹽丘模型(圖2)進行成像測試。模型主要包括鹽上沉積層、鹽丘及鹽下構造。主要考察ALSRTM成像算法對鹽丘側翼高陡小構造、鹽上小斷裂帶及鹽下小斷裂帶的成像效果，其中鹽丘側翼的高陡小構造會引起嚴重的構造假象， 而且由于鹽丘的屏蔽作用，導致鹽下構造的振幅較弱，難以準確成像[28]。稀疏采樣引起成像結果中出現(xiàn)嚴重的偏移假象，因此能更好地測試、對比常規(guī)LSRTM和ALSRTM去除淺層低頻噪聲的效果。圖3為RTM成像結果。由圖可見：在拉普拉斯濾波前，RTM剖面中含有較強的低波數(shù)噪聲和震源效應，難以識別地下構造(圖3a)； 在拉普拉斯濾波后，在模型淺部仍然存在低頻噪聲和頂部的震源效應，成像能量不均衡、成像精度低，仍然難以識別鹽下構造，但成像效果好于未做拉普拉斯濾波的RTM結果(圖3b)。圖4為LSRTM和ALSRTM成像結果對比。由圖可見： ①當?shù)?次時，在常規(guī)LSRTM成像結果中，鹽上沉積層存在明顯的偏移假象(圖4a)； ALSRTM能較徹底、有效地壓制低頻噪聲，且震源采集效應更弱，另外，對于鹽丘上部小斷裂的成像分辨率更高(圖4b)。②當?shù)?0次時，常規(guī)LSRTM和ALSRTM均有效壓制了低頻噪聲，采集腳印減弱，成像效果明顯改善(圖4c、圖4d)，但常規(guī)LSRTM成像結果的淺層仍然存在較嚴重的噪聲(圖4c)； 值得注意的是，ALSRTM僅僅迭代5次時的成像結果(圖4b)對噪聲和采集腳印的壓制效果甚至優(yōu)于LSRTM迭代10次(圖4c)的結果。③當?shù)?0次時，常規(guī)LSRTM和ALSRTM均能有效壓制鹽上沉積層的成像噪聲，也能準確識別巖下構造(圖4e、圖4f)。

圖2 SEG/EAGE二維鹽丘模型(a)真實速度場； (b)平滑速度場； (c)反射率模型

圖3 RTM成像結果(a)拉普拉斯濾波前； (b)拉普拉斯濾波后

震源是主頻為20Hz的雷克子波，時間采樣間隔為0.5ms，時間采樣點數(shù)為4001； 總炮數(shù)為30炮，共645道檢波器接收，炮間距為215m，道間距為10m

圖4 LSRTM與ALSRTM成像結果對比(a)LSRTM迭代5次；(b)ALSRTM迭代5次；(c)LSRTM迭代10次；(d)ALSRTM迭代10次；(e)LSRTM迭代30次；(f)ALSRTM迭代30次

為了更清楚地對比常規(guī)LSRTM與ALSRTM方法的成像效果，從LSRTM和ALSRTM成像結果(圖4)中截取淺層斷塊區(qū)域(深度范圍為200～800m，水平范圍為2000～5400m)的局部放大圖進行對比(圖5)?？梢姡涸诔Ｒ?guī)LSRTM成像剖面上鹽上小斷層的接觸關系模糊，分辨率較低，當?shù)?、10次時未能有效去除鹽丘左側之上和鹽上反射層之間的低頻假象(圖5a、圖5c黑色箭頭所示)； 應用ALSRTM后(圖5b、圖5d、圖5f)，鹽上斷層更為清晰，而且僅僅在迭代5次時，就可徹底壓制鹽丘左側側翼之上和鹽上反射層之間的噪聲(圖5b)。

從LSRTM以及ALSRTM 的成像結果中分別抽取水平坐標2700m和深度400m處(圖4f中虛線位置)的地震單道振幅(圖6)對比各方法的保幅性(迭代30次)?？梢姡孩僭谒阶鴺?700m處，在模型的淺部區(qū)域由于存在低頻噪聲，常規(guī)LSRTM的振幅曲線存在明顯跳動，而ALSRTM振幅曲線抖動幅度較小，與真實反射系數(shù)振幅曲線吻合更好(圖6a的箭頭所示)；在中深部ALSRTM 與常規(guī)LSRTM的保幅性幾乎相同，兩者的振幅曲線均較好地接近真實反射系數(shù)振幅曲線(圖6a的*所示)。②截取鹽丘正上方低頻噪聲嚴重區(qū)域在深度400m處的振幅曲線(圖6b)可見，在模型的水平坐標2800、3400、4500m(被低頻噪聲嚴重干擾的反射界面)處，常規(guī)LSRTM振幅曲線與真實反射系數(shù)振幅曲線差距較大，而ALSRTM振幅曲線更接近真實反射系數(shù)振幅曲線(圖6b的箭頭所示)；在沒有反射界面的位置，ALSRTM振幅曲線較常規(guī)LSRTM振幅曲線抖動小，抗噪性更好(圖6b的*所示)。因此在淺層，ALSRTM比常規(guī)LSRTM具有更好的保幅性，尤其在低頻噪聲存在的位置，ALSRTM的優(yōu)勢更為明顯，在中深層兩者的振幅曲線均更接近于真實反射系數(shù)振幅曲線。

圖7為ALSRTM、常規(guī)LSRTM成像結果的殘差隨迭代次數(shù)變化圖。由圖可見：①ALSRTM在第3次和第6次迭代時出現(xiàn)殘差突然上升趨勢(圖7的黑色箭頭所示)；②在第7次迭代之后，隨著迭代次數(shù)增加，兩種方法的殘差都逐漸減小，在13次迭代之后，ALSRTM的殘差略微小于常規(guī)LSRTM，兩種方法的成像結果都逐漸趨近于真實模型；③隨著迭代次數(shù)的增加，兩種方法殘差的減小速度逐漸減慢。需要說明的是，第3次和第6次成像時使用式(7)，其余迭代過程均使用式(2)，式(2)是在反偏移算子和偏移算子互為共軛轉置的前提下得到的，因此數(shù)據(jù)殘差穩(wěn)定收斂，而式(7)是基于波場的動力學特征，在第3次和第6次迭代時不滿足反偏移算子和偏移算子互為共軛轉置的關系而出現(xiàn)抖動。為了測試ALSRTM和LSRTM對速度場的敏感性，在速度場中加入2%的誤差后進行成像測試。 圖8為加入2%速度誤差的成像結果。 由圖可見：隨著迭代次數(shù)增加，ALSRTM(圖8b、圖8d、圖8f)和LSRTM成像結果(圖8a、圖8c、圖8e)均能清晰展現(xiàn)淺、中、深層構造，但前者的偏移假象更少、成像精度更高。因此，ALSRTM對速度場求取不準具有更好的適應性。

圖5 圖4的局部放大(a)LSRTM迭代5次；(b)ALSRTM迭代5次；(c)LSRTM迭代10次；(d)ALSRTM迭代10次；(e)LSRTM迭代30次；(f)ALSRTM迭代30次

圖6 成像結果單道振幅對比(迭代30次)(a)水平坐標2700m處； (b)深度坐標400m處

圖7 ALSRTM、常規(guī)LSRTM成像結果的殘差隨迭代次數(shù)變化(迭代60次)

圖8 加入2%速度誤差的成像結果(a)LSRTM迭代5次；(b)ALSRTM迭代5次；(c)LSRTM迭代10次；(d)ALSRTM迭代10次；(e)LSRTM迭代30次；(f)ALSRTM迭代30次

野外地震采集數(shù)據(jù)中不可避免地含有環(huán)境及人文因素造成的隨機噪聲。為了測試ALSRTM對低信噪比數(shù)據(jù)的適應性，通過 Madagascar軟件對炮記錄添加隨機噪聲得到低信噪比單炮數(shù)據(jù)(圖9)。圖10為低信噪比數(shù)據(jù)成像結果。由圖可見：對于低信噪比數(shù)據(jù)，LSRTM(圖10a、圖10b)和ALSRTM成像結果(圖10c、圖10d)都含有類似于隨機噪聲的偏移假象，但后者更好地壓制了低頻噪聲。因此，ALSRTM對低信噪比炮記錄的適應性更好。

圖9 低信噪比單炮數(shù)據(jù)

圖10 低信噪比數(shù)據(jù)成像結果(a)LSRTM迭代5次； (b)LSRTM迭代10次； (c)ALSRTM迭代5次； (d)ALSRTM迭代10次

4 結束語

通過波場數(shù)據(jù)得到的反射角信息構建逆散射成像條件，并與最小二乘逆時偏移(LSRTM)結合，發(fā)展了一種基于角度濾波成像的最小二乘逆時偏移方法(ALSRTM)，并由波動方程的能量守恒分析了ALSRTM的可行性和保幅性。在實現(xiàn)算法的基礎上，對SEG/EAGE二維鹽丘模型的稀疏采集地震數(shù)據(jù)的成像結果表明，ALSRTM在改善常規(guī)LSRTM的成像質量及保幅性方面是有效的。為了進一步提高計算效率，可以將ALSRTM與平面波震源編碼、相位編碼、GPU加速等方法結合，使其更好地適應實際地震資料處理。

感謝中國石油大學(華東)地震波傳播與成像(SWPI)課題組的支持。感謝開源軟件Madagascar (http://www.ahay.org/RSF) 提供的技術支持。

參考文獻

[1] Youn O K,Zhou H W.Depth imaging with multiples.SEG Technical Program Expanded Abstracts,1999,18:246-255.

[2] Claerbout J F,Johnson A G.Extrapolation of time-dependent waveforms along their path of propagation.Geophysical Journal International,1971,26(3):285-293.

[3] 劉紅偉,劉洪,鄒振等.地震疊前逆時偏移中的去噪與存儲.地球物理學報,2010,53(9):2171-2180.

Liu Hongwei,Liu Hong,Zou Zhen et al.The problems of denoise and storage in seismic reverse time migration.Chinese Journal of Geophysics,2010,53(9):2171-2180.

[4] Xu S,Chen F,Tang B et al.Noise removal by migration of time-shift images.Geophysics,2014,79(3):S105-S111.

[5] Robin F F,Paul F,Phil K et al.Suppressing artifacts in prestack reverse time migration.SEG Technical Program Expanded Abstracts,2005,24:2049-2051.

[6] 李闖,黃建平,李振春等.預條件最小二乘逆時偏移方法.石油地球物理勘探,2016,51(3):513-520.

Li Chuang,Huang Jianping,Li Zhenchun et al.Preconditioned least-squares reverse time migration.OGP,2016,51(3):513-520.

[7] Guitton A,Kaelin B,Biondi B.Least-square attenuation of reverse time migration artifacts.Geophysics,2007,72(1):S19-S23.

[8] Baysal E,Dan D K,Sherwood J W C.A two-way nonreflecting wave equation.Geophysics,2012,49(2):132-141.

[9] Yoon K,Marfurt K J.Reverse-time migration using the Poynting vector.Geophysical Exploration,2006,37(1):102-107.

[10] Liu F,Zhang G,Morton S A et al.An effective imaging condition for reverse-time migration using wavefield decomposition.Geophysics,2011,76(1):S29-S39.

[11] Rocha D C J.Wavefield imaging using the energy norm.Geophysics,2015,81(4):S151-S163.

[12] Op’t Root T J P M,Stolk C C and de Hoop M V.Linearized angular filtering based on seismic reverse time migration.Journal De Mathématiques Pures Et Appliqués,2012,98(2):211-238.

[13] Whitmore N D,Crawley S.Applications of RTM inverse scattering imaging conditions.SEG Technical Program Expanded Abstracts,2012,31:1-6.

[14] Kiyashchenko D,Plessix R E,Kashtan B et al.A modified imaging principle for true-amplitude wave-equation migration.Geophysical Journal International,2007,168(3):1093-1104.

[15] Tarantola A.Linearized inversion of seismic reflection data.Geophysical Prospecting,1984,32(6):998-1015.

[16] Nemeth T,Wu C,Schuster G T.Least-squares migration of incomplete reflection data.Geophysics,1999,64(1):208-221.

[17] K?ltzsch P.Inverse problems of acoustic and elastic waves.ZAMM,1986,66(9):447-447.

[18] Li C,Huang J P,Li Z C et al.Regularized least-squares migration of simultaneous-source seismic data with adaptive singular spectrum analysis.Petroleum Science,2017,14(1):61-74.

[19] Plessix R E,Mulder W A.Frequency-domain finite-difference amplitude-preserving migration.Geophysical Journal International,2004,157(3):975-987.

[20] Dai W,Schuster J.Least-squares migration of simul-taneous sources data with a deblurring filter.SEG Technical Program Expanded Abstracts,2009,28:4338.

[21] Dai W,Schuster G T.Plane-wave least-squares reverse-time migration.Geophysics,2013,78(4):S165-S177.

[22] 黃建平,曹曉莉,李振春等.最小二乘逆時偏移在近地表高精度成像中的應用.石油地球物理勘探,2014,49(1):107-112.

Huang Jianping,Cao Xiaoli,Li Zhenchun et al.Least square reverse time migration in high resolution imaging of near surface.OGP,2014,49(1):107-112.

[23] 李振春,李闖,黃建平等.基于先驗模型約束的最小二乘逆時偏移方法.石油地球物理勘探,2016,51(4):738-744.

Li Zhenchun,Li Chuang,Huang Jianping et al.Regularized least-squares reverse time migration with prior model.OGP,2016,51(4):738-744.

[24] 李闖,黃建平,李振春等.平面波最小二乘逆時偏移編碼策略分析.石油物探,2015,54(5):592-601.

Li Chuang,Huang Jianping,Li Zhenchun et al.Analysis on the encoding strategies of plane-wave least-square reverse time migration.GPP,2015,54(5):592-601.

[25] 黃建平,孫鄖松,李振春等.一種基于分頻編碼的最小二乘裂步偏移方法.石油地球物理勘探,2014,49(4):702-707.

Huang Jianping,Sun Yunsong,Li Zhenchun et al.Least-squares split-step migration based on frequency-division encoding.OGP,2014,49(4): 702-707.

[26] Stolk C C,Symes W W.Kinematic artifacts in pre-stack depth migration.Geophysics,2004,69(2):562-575.

[27] Evans L.Partial Differential Equations(Second Edition).Wadsworth & Brooks/cole Mathematics,Berkeley,California,American,2010.

[28] 黃建平,劉培君,李慶洋等.一種棱柱波逆時偏移方法及優(yōu)化.石油物探,2016,55(5):719-727.

Huang Jianping,Liu Peijun,Li Qingyang et al.An optimized reverse time migration approach for the prism wave.GPP,2016,55(5):719-727.