基于Poynting矢量行波分離的最小二乘逆時(shí)偏移

王曉毅， 張江杰， 許宏橋， 田寶卿

1 中國科學(xué)院地質(zhì)與地球物理研究所， 中國科學(xué)院油氣資源研究重點(diǎn)實(shí)驗(yàn)室， 北京 100029 2 中國科學(xué)院地球科學(xué)研究院， 北京 100029 3 中國科學(xué)院大學(xué)， 北京 100049 4 中國科學(xué)院地質(zhì)與地球物理研究所， 中國科學(xué)院頁巖氣與地質(zhì)工程重點(diǎn)實(shí)驗(yàn)室， 北京 100029

0 引言

作為一種基于雙程波動(dòng)方程進(jìn)行波場延拓的成像方法，逆時(shí)偏移具有諸多優(yōu)點(diǎn)：可以利用反射波、多次波、回轉(zhuǎn)波和棱柱波等各種波型，不受地層傾角的限制，并且對起伏地表和復(fù)雜速度分布有很好的適應(yīng)性.但是相較于單程波偏移算子，逆時(shí)偏移面臨著由背向反射所引起的低頻噪聲發(fā)育的問題.所謂背向反射，是指在震源波場或者檢波器波場延拓過程中，當(dāng)波遇到波阻抗界面時(shí)所發(fā)生的次生反射.背向反射波與正常傳播的波之間的互相關(guān)在地震剖面上所表現(xiàn)出的結(jié)果即低頻噪聲(劉紅偉等，2010；杜啟振等，2013).一般而言，界面波阻抗差異越大，背向反射越發(fā)育，低頻噪聲越強(qiáng).

常規(guī)的偏移方法受低覆蓋的觀測系統(tǒng)，波的幾何擴(kuò)散、吸收衰減等效應(yīng)引起的弱振幅，有限偏移孔徑造成的低分辨率，以及震源子波的震蕩等諸多因素的影響，串?dāng)_噪聲和偏移假象比較發(fā)育(Schuster，2017).相比之下，最小二乘偏移(LSM)通過觀測數(shù)據(jù)和理論數(shù)據(jù)的不斷擬合和校正，取得的成像結(jié)果振幅保真性好，分辨率和信噪比更高(馬方正等，2016；李振春等，2017).最初，LSM主要用于射線類和單程波類偏移算法，效果顯著.后來，最小二乘逆時(shí)偏移將最小二乘思想和逆時(shí)偏移方法結(jié)合，既可以利用各種波攜帶的信息，又可以得到能量更加均衡且消除了采集腳印的地下結(jié)構(gòu)剖面.

在逆時(shí)偏移(RTM)中，關(guān)于低頻噪聲的壓制有著諸多辦法，例如Laplace濾波、成像角度限制等(Zhang et al., 2007; Yoon and Marfurt, 2006).相應(yīng)地，這些技術(shù)也被引入LSRTM中，并取得了良好的效果(楊凱和張劍鋒, 2017; Yang and Zhang, 2018).以Laplace濾波為例，我們可以在每次迭代求取梯度之后對其進(jìn)行濾波消除噪聲，但是這改變了振幅和相位信息，因此需要預(yù)先對震源子波進(jìn)行二次積分并對求取的梯度進(jìn)行振幅補(bǔ)償.近些年來，一種基于行波分離的成像方法被提出(Liu et al., 2011; Fei et al., 2015).根據(jù)入射波和出射波之間的幾何關(guān)系，我們可以發(fā)現(xiàn)兩者是分別沿著下行和上行方向的.Liu等(2011)提出只保留震源波場和檢波器波場沿不同垂直方向的組分之間的互相關(guān)作為新的成像條件，從而有效地壓制地震剖面中的低頻噪音.Fei等(2015)進(jìn)一步提出只保留震源波場的下行波和檢波器波場的上行波之間的互相關(guān)運(yùn)算，以避免回轉(zhuǎn)波產(chǎn)生的偏移假象.王一博等(2016)分別將震源波場和檢波器波場分離為左上、左下、右上和右下四個(gè)組分，并將震源波場和檢波器波場的這些組分兩兩互相關(guān)運(yùn)算，得到了16種成像結(jié)果，有效地分析和分離了低頻噪聲和串?dāng)_噪聲.值得注意的是，以上的波場分離主要是在波數(shù)域?qū)崿F(xiàn)的，涉及到Hilbert變換和Fourier變換，因而計(jì)算量較大，存儲成本較高.

Poynting矢量提供了波場傳播的方向信息，因此無論是在聲波還是彈性波地震成像都有著廣泛而充分的應(yīng)用.一方面，Poynting矢量被用于壓制低頻噪聲(Yoon and Marfurt, 2006)；借助震源波場和檢波器波場的方向信息，我們可以得到入射波與出射波之間的張角，并在成像條件中添加張角(或者入射角)加權(quán)項(xiàng)對大角度進(jìn)行衰減或截?cái)?，以達(dá)到壓制噪聲的目的.另一方面，Poynting矢量還可以被用于提取角度域共成像點(diǎn)道集(ADCIG)，包括張角道集和傾角道集(王保利等，2013; Liu and Zhang, 2018).這些道集為AVA、MVA等分析提供了依據(jù).在本文中，借助Poynting矢量所包含的方向信息，我們以較低的代價(jià)獲取波場的上行波組分和下行波組分(黃杰等，2018).特別是為了使Poynting矢量求取的魯棒性更強(qiáng)，我們采用了一種局部最小二乘意義下的求取結(jié)果(蘆永明等，2017).

本文結(jié)構(gòu)如下：首先，我們回顧了LSRTM的基本原理，對Born近似和迭代反演進(jìn)行了闡述.其次，通過對雙程波動(dòng)方程類偏移方法成像條件的分析，解釋了低頻噪聲的產(chǎn)生和壓制方法，并提出了對LSRTM梯度算子的改造.然后，我們引入Poynting矢量并利用其方向指示性實(shí)現(xiàn)了行波分離，提高了算法實(shí)現(xiàn)的效率.最后，通過數(shù)值算例驗(yàn)證了我們方法的有效性.

1 基本理論

最小二乘逆時(shí)偏移通過迭代反演得到的最終剖面具有更高的分辨率.為取得更快的收斂速率，我們通過在梯度構(gòu)建中對震源和檢波器波場分別選擇合適的組分進(jìn)行互相關(guān)運(yùn)算以壓制低頻噪聲和偏移假象.特別地，Poynting矢量的引入使得我們方案的實(shí)現(xiàn)變得方便快捷.

1.1 最小二乘逆時(shí)偏移

按照擾動(dòng)理論，真實(shí)速度v可以分為偏移速度v0和擾動(dòng)速度δv兩部分，對應(yīng)的波場分別為背景波場p0和擾動(dòng)波場δp.根據(jù)Born近似，擾動(dòng)波場由背景波場和速度擾動(dòng)共同決定，即：

(1)

其中，xs表示震源所在的位置，x是任一空間坐標(biāo)點(diǎn)，t表示時(shí)間，f是震源子波.m(x)=2δv(x)/v0(x)表示擾動(dòng)模型，也叫真實(shí)反射率，刻畫了偏移速度與真實(shí)速度之間的偏離程度.在檢波器位置xr上記錄擾動(dòng)波場，我們可以得到模擬數(shù)據(jù)d，即：

d(t,xr)=δp(xr,t;xs).

(2)

由以上敘述可知，借助于Born近似，我們可以建立擾動(dòng)模型與觀測記錄之間的線性關(guān)系.以矩陣的形式對此加以表達(dá)，我們有：

d=Lm,

(3)

在這里，L為正演算子，也被稱為反偏移算子.

在已有觀測系統(tǒng)和采集記錄dobs的情況下，求取擾動(dòng)模型m需要用到正演算子的逆L-1.然而，對于地震勘探而言，遇到的反問題本身就是病態(tài)的,L-1往往并不存在.另一方面勘探范圍大，處理的數(shù)據(jù)也是海量的，因此直接求逆從運(yùn)算量和存儲量上也不可行.所謂最小二乘方法，就是通過迭代的方法來逼近反問題的最優(yōu)解，從而避免直接的取逆.目標(biāo)函數(shù)選用預(yù)測數(shù)據(jù)d與觀測數(shù)據(jù)dobs之間殘差的L2范數(shù)，殘差反傳使得模型不斷被校正，直至預(yù)測數(shù)據(jù)和觀測數(shù)據(jù)足夠接近.在成像過程中，梯度算子LT也被稱為偏移算子，按照伴隨狀態(tài)法可以表示成如下形式(Plessix, 2006)：

(4)

其中，q是伴隨波場，滿足：

(5)

對于常規(guī)的逆時(shí)偏移，通過一次震源波場正向延拓和檢波器波場反向延拓且以(4)式為成像條件可以得到地下結(jié)構(gòu)的成像剖面mmig.在最小二乘逆時(shí)偏移中，以(4)式作為梯度對反演結(jié)果不斷進(jìn)行校正，可以消除成像剖面中由采集腳印和非均勻照明引起的偏移假象.相較于常規(guī)偏移，最小二乘偏移的運(yùn)算量需求更加龐大.在每次迭代過程中，僅梯度計(jì)算就需要一次震源波場(包括背景波場和擾動(dòng)波場)正向延拓和檢波器(伴隨)波場反向延拓.在這里需要指出的是，我們每次迭代更新的是反射率模型m，而非速度模型v.

求解最小二乘問題的最優(yōu)解涉及到的迭代方法包括最速下降法(SD)，共軛梯度法(CG)以及擬牛頓法(QN)等.在本文中我們采用帶有預(yù)條件的共軛梯度算法來對反問題進(jìn)行求解.這一方法的優(yōu)點(diǎn)在于其既避免了對Hessian陣顯式求逆，又克服了最速下降法收斂慢的缺點(diǎn)，并且所需存儲量少，穩(wěn)定性高(Schuster, 2017).

1.2 低頻噪聲的產(chǎn)生與去除

波動(dòng)方程是數(shù)值模擬地震波傳播過程的理論依據(jù)，一方面它可以生成包括一次反射波和多次波在內(nèi)的多種類型的波，另一方面它具有同時(shí)刻畫波沿各個(gè)方向傳播的能力.因此，相較于一般的單程波偏移算法，逆時(shí)偏移具有如下優(yōu)勢：處理回轉(zhuǎn)波、棱柱波等特殊波型十分方便，不受地層傾角限制且具有更大的照明區(qū)域.然而，強(qiáng)振幅的低頻噪聲往往會污染偏移結(jié)果，特別是在淺層和強(qiáng)反射界面的上部.當(dāng)背向反射發(fā)育時(shí)，對震源波場和檢波器波場采用傳統(tǒng)的互相關(guān)條件會沿著射線路徑成像，從而在偏移剖面上產(chǎn)生強(qiáng)烈的干擾.

對于雙程波偏移方法，基于全波的互相關(guān)成像條件可以寫為

(6)

其中，us和ur分別表示震源波場和檢波器波場，上標(biāo)u和d則分別表示上行波場和下行波場.由式(6)可知，除了震源波場的下行波和檢波器波場的上行波可以成像外，其他組分也參與成像.其中，沿相同方向傳播的組分互相關(guān)運(yùn)算會產(chǎn)生較大幅值的噪聲，這種噪聲在頻率域表現(xiàn)出低頻特征(杜啟振等, 2013)，因此可以用Laplace濾波加以消除.但這也會改變地震剖面同相軸的振幅和相位特征，因此需要對震源的預(yù)處理和額外的矯正.為了提高成像質(zhì)量，可以改造成像條件，即省略(6)式的后兩項(xiàng)來壓制低頻噪聲.從而，我們得到了新的成像條件(Liu et al., 2011)：

(7)

對于一次反射或者散射，震源波場的入射波下行而檢波器波場的出射波上行，如圖1a中的A點(diǎn)和B點(diǎn)所示.特別地，對于多次波的一些反射點(diǎn)也滿足這一特征，如圖1b中的C點(diǎn).式(7)中右側(cè)的第一項(xiàng)可以很好地將這部分波加以利用，使其在正確位置成像.對于如圖1b中的D點(diǎn)，式(7)中右側(cè)第二項(xiàng)滿足對其成像的要求.但是當(dāng)偏移速度變化特別劇烈時(shí)，回轉(zhuǎn)波會在剖面上產(chǎn)生假的同相軸(Fei et al.，2015)，如圖2所示.這會對后續(xù)的處理和解釋造成困擾.因此，式(7)又可近一步簡化為

(8)

最小二乘逆時(shí)偏移與常規(guī)的逆時(shí)偏移一樣，都面臨著雙程波成像所帶來的低頻噪聲和串?dāng)_噪聲問題.因此，我們可以仿照上述成像條件對梯度算子進(jìn)行改造.與式(7)和(8)相對應(yīng)，我們分別有：

(9)

圖1 波的傳播路徑(a) 反射波和繞射波的傳播；其中A點(diǎn)為反射點(diǎn)，B點(diǎn)為繞射點(diǎn)； (b) 多次波的傳播路徑.C點(diǎn)和D點(diǎn)均為多次反射點(diǎn)，但是適用于不同的成像條件.Fig.1 Wave propagation paths(a) Propagation of reflected and diffracted waves. Point A is the reflection point and point B is the diffraction point； (b) The propagation path of multiples. Points C and D are multiple reflection points, but they are suitable for different imaging conditions.

(10)

1.3 Poynting矢量與上下行波分離

從物理意義上來說，Poynting矢量是方向性的能流密度矢量，其大小為單位時(shí)間穿過單位面積的能量.這一矢量在電磁勘探和地震勘探領(lǐng)域有著重要的意義和作用，特別是其良好的方向指示性對波傳播的刻畫和反演成像有著很大的幫助(王保利等，2013).

圖2 偏移假象產(chǎn)生示意圖當(dāng)偏移速度存在劇烈變化區(qū)域時(shí)，能夠在正確的位置成像，而則會生成虛假的軸.Fig.2 Schematic of the generation of migration artifactsWhen there is a sharp change in the migration velocity, energy is generated in the correct position by using but the false image is also produced by

對于聲波介質(zhì)，Poynting矢量可以按照式(11)進(jìn)行計(jì)算(Yoon and Marfurt, 2006)：

(11)

基于(11)式的Poynting矢量求取面臨計(jì)算的穩(wěn)定性問題，也就是當(dāng)聲壓場達(dá)到峰值時(shí)，s的幅值為0，從而喪失了其方向性.為使算法更加穩(wěn)健，局部最小二乘平均意義下的Poynting矢量被提出(蘆永明等，2017),即：

(12)

其中，i表示代求點(diǎn)，Ω表示i點(diǎn)的鄰域.

圖3a是聲波在速度為1000 m·s-1的均勻模型中傳播的波場快照,其傳播時(shí)間為0.2 s.紅色箭頭即為按照(11)式和(12)式所求的Poynting矢量，其指向聲波的傳播方向，也就是波前面的外法線方向.按照Poynting矢量的指向，我們可以非常方便快捷地完成波場的上下行波分離.如圖3b、c所示，波場被分為上下行波兩部分.對于復(fù)雜波場我們也可以按照這種方式對波場進(jìn)行分離.圖4a是一復(fù)雜模型，將震源置于模型中間位置激發(fā)，在傳播時(shí)間為0.6 s記錄下的波場快照如圖4b所示.與此同時(shí)，我們求取了Poynting矢量，并在圖4c、d中分別展示其x分量和z分量.按照Poynting矢量的指向，我們對復(fù)雜的波場快照進(jìn)行了分離(圖4e、f).值得注意的是，對于波場重疊的區(qū)域，(11)式只能給出單一方向，這意味著這種簡單的分離對于重疊區(qū)域是存在誤差的.當(dāng)偏移速度足夠光滑時(shí)，波的混疊并不嚴(yán)重，因而對成像造成的影響并不大.

1.4 方法流程

盡管通過改造梯度算子，可以避免低頻噪聲和假象的產(chǎn)生，但是這也有可能會造成部分有效信號的丟失.特別地，Poynting矢量在重疊區(qū)域只能給出一個(gè)方向，導(dǎo)致部分能量被錯(cuò)誤地歸位.為了降低這部分的影響，我們將改造后的梯度算子只用在早期的迭代中，并根據(jù)迭代次數(shù)或者數(shù)據(jù)誤差的閾值在后邊調(diào)整回全波的梯度算子，從而充分利用所有波所攜帶的地下信息.

結(jié)合上述理論，我們將LSRTM的整體工作流程總結(jié)如下：

(1)根據(jù)偏移速度v0正演得到模擬數(shù)據(jù)d0，并讓原始數(shù)據(jù)dobs減去其得到dres；

(2)設(shè)定初始反射率模型m0,并基于Born近似計(jì)算模擬數(shù)據(jù)Lm0，從而得到殘差r=Lm0-dres；

(3)利用偏移算子(9)或者(10),將殘差反傳得到gk+1=LTr;

(4)根據(jù)帶有預(yù)條件的共軛梯度法，利用gk、gk+1和zk分別計(jì)算模型更新方向zk+1和迭代步長α，即：

(13)

其中，P算子表示震源照明補(bǔ)償.隨后，完成模型迭代mk+1=mk-αzk+1；

(5)計(jì)算數(shù)據(jù)殘差r=Lmk+1-dres，重復(fù)第3步和第4步,直至‖r‖2小于某個(gè)閾值或者迭代次數(shù)達(dá)到預(yù)定次數(shù);

(6)利用偏移算子(4),將殘差反傳得到gk+1=LTr,重復(fù)第4步，直至誤差低于某個(gè)閾值或者迭代次數(shù)達(dá)到設(shè)定值.

圖3 均勻介質(zhì)中聲波Poynting矢量的求取及上下行波分離(a) 0.2 s時(shí)刻的波場快照；紅色箭頭表示歸一化后的Poynting矢量，沿波的外法線方向； (b) 上行波場； (c) 下行波場.Fig.3 Calculation of Poynting vector for acoustic waves and the up/down wavefield decomposition in the homogeneous medium(a) Snapshot of acoustic wave at 0.2 s. The red arrow represents the normalized Poynting vector, which is along the outer normal direction of the wave； (b) Upgoing wavefield； (c) Downgoing wavefield.

圖4 復(fù)雜介質(zhì)中聲波Poynting矢量的求取及上下行波分離(a) 復(fù)雜介質(zhì)速度模型； (b) 0.6 s時(shí)刻的波場快照； (c) Poynting矢量的x分量； (d) Poynting矢量的z分量； (e) 上行波場； (f) 下行波場.對于波場發(fā)生重疊的區(qū)域，Poynting矢量只能給出一個(gè)方向.Fig.4 Calculation of Poynting vector for acoustic waves and the up/down wavefield decomposition in a complex medium(a) Velocity of the complex model； (b) Snapshot of acoustic wave at 0.6 s； (c) x-component of the Poynting vector； (d) z-component of the Poynting vector； (e) Upgoing wavefield； (f) Downgoing wavefield. The Poynting vector gives only one direction for the overlapping area of the wavefield.

2 數(shù)值算例

為了解釋方法的有效性，我們用兩個(gè)數(shù)值算例進(jìn)行了驗(yàn)證.其中既包括簡單的雙層模型，又包括相對復(fù)雜的Marmousi模型.偏移所需要的地震記錄數(shù)據(jù)集由時(shí)間二階，空間八階的有限差分算法模擬得到.

2.1 雙層介質(zhì)

模型界面位于500 m處，上層介質(zhì)速度為1000 m·s-1，而下層介質(zhì)速度為1500 m·s-1.震源子波設(shè)置為主頻25 Hz的雷克子波，將其在地表激發(fā)且間隔設(shè)為250 m.所以，總炮數(shù)為5，從左至右依次排列.檢波器置于地表，空間間隔5 m,時(shí)間采樣間隔為0.5 ms，總記錄時(shí)長為2 s.

偏移所用的速度模型是由真實(shí)模型經(jīng)高斯光滑所得，并由此計(jì)算出真實(shí)的反射率模型.在進(jìn)行LSRTM之前，我們選取第三炮的數(shù)據(jù)對正演算子和偏移算子進(jìn)行了點(diǎn)積測試.表1詳細(xì)列出了我們點(diǎn)積測試的結(jié)果，可見無論是采用公式(4)還是(9)或者(10)作為偏移算子，都近似滿足〈Lm,Lm〉=〈m,LTLm〉以及〈LTLm,LTLm〉=〈Lm,LLTLm〉.

表1 雙層介質(zhì)點(diǎn)積測試結(jié)果Table 1 Dot product test result for a two-layer model

此外，因?yàn)楣?9)、(10)不能保證像公式(4)一樣與正演算子滿足充分的伴隨關(guān)系，其測試結(jié)果的相等性要差于公式(4)測試結(jié)果的相等性.圖6a、b、c分別是三種成像條件對應(yīng)的偏移結(jié)果，可以看出相比于全波梯度算子(4)，成像條件(9)和(10)產(chǎn)生的單炮剖面更為干凈，沒有低頻噪聲發(fā)育.

對于簡單的雙層模型，我們使用(9)式作為梯度算子.圖7a是真實(shí)反射率模型，圖7b是常規(guī)全波RTM偏移成像的結(jié)果，而圖7c、d分別是使用常規(guī)LSRTM方法和我們的方案迭代10次得到的反射率模型.顯然，與RTM相比，LSRTM的偏移剖面振幅更加均衡，偏移假象也比較少.特別地，我們將常規(guī)LSRTM和基于Poynting矢量行波分離的LSRTM在收斂速率上進(jìn)行了比較，如圖8所示.可見，在相同的迭代次數(shù)條件下我們的方法模型預(yù)測誤差更小.

2.2 Marmousi模型

本模型深度為1.5 km，寬度為4.0 km，具體速度參數(shù)如圖9a所示.本實(shí)驗(yàn)采用主頻為15 Hz，時(shí)間采樣間隔為1 ms的雷克子波.炮點(diǎn)置于地表0 km至4.0 km處，炮間距設(shè)為200 m.檢波點(diǎn)在地表均勻排列，間隔設(shè)為20 m，總記錄時(shí)長為4 s.反演所用背景速度為圖9b，模型真實(shí)反射率如圖9c所示.

圖5 雙層介質(zhì)中重構(gòu)的震源波場(a)和經(jīng)AGC處理后的檢波器波場(b)其中藍(lán)色箭頭所指為背向反射波，紅色線條表示反射界面.Fig.5 Reconstructed source wavefield (a) and receiver wavefield processed by AGC (b) in a two-layer modelBlue arrows point to the back-reflected waves, and red lines represent the reflective interfaces.

圖6 不同成像條件單炮數(shù)據(jù)的成像結(jié)果(a)、(b)和(c)對應(yīng)的偏移算子分別為(4)式、(9)式和(10)式.相較于全波成像條件，后兩者沒有低頻噪聲發(fā)育，因此剖面質(zhì)量更高.Fig.6 Imaging results of single-shot data under different imaging conditionsThe migration operators corresponding to (a), (b) and (c) are equations (4), (9) and (10), respectively. Compared with full-wave imaging conditions, the latter two have higher-quality profiles without low-frequency noises.

圖7 (a) 真實(shí)反射率模型； (b) 基于全波成像條件的RTM偏移結(jié)果； (c) 常規(guī)LSRTM迭代10次的反演結(jié)果； (d) 基于Poynting矢量行波分離的LSRTM方法迭代10次的反演結(jié)果Fig.7 (a) True reflectivity model； (b) Result of RTM with the full-wave imaging condition； (c) Inversion result after 10 iterations using the conventional LSRTM； (d) Inversion result after 10 iterations using LSRTM with the wavefield decomposition based on Poynting vector

圖8 對于雙層介質(zhì)，常規(guī)LSRTM和基于Poynting矢量行波分離的LSRTM收斂性的對比(顯然，我們方法的收斂速率更高，殘差也更小)Fig.8 Comparison of convergence between conventional LSRTM and LSRTM with the wavefield decomposition based on the Poynting vector for the two-layer model(Note that our scheme has a higher convergence rate and a minor residual)

圖9 Marmousi模型(a) 真實(shí)速度； (b) 偏移速度； (c) 真實(shí)反射率.Fig.9 Marmousi model(a) True velocity； (b) Migration velocity； (c) True reflectivity.

圖10 Marmousi模型中的波場快照(a) 震源波場； (b) 經(jīng)AGC處理后顯示的震源波場； (c) 檢波器波場.Fig.10 Snapshots of wavefield in the Marmousi model(a) The source wavefield; (b) The displayed source wavefield after the AGC process; (c) The receiver wavefield.

圖11 Marmousi模型的偏移成像結(jié)果(a) 基于全波成像條件的RTM； (b) 常規(guī)LSRTM； (c) 基于Poynting矢量行波分離的LSRTM.其中，LSRTM的迭代次數(shù)為20.Fig.11 Migration results of the Mamousi model(a) RTM with the full-wave imaging condition; (b) The conventional LSRTM; (c) LSRTM with the wavefield decomposition based on Poynting vectors. Here, the iteration numbers of two LSRTM schemes are both 20.

圖12 對于Marmousi模型，常規(guī)LSRTM和基于Poynting矢量行波分離的LSRTM的收斂曲線Fig.12 Comparison of convergence between conventional LSRTM and LSRTM with the wavefield decomposition based on the Poynting vector for the Marmousi model

圖13 常規(guī)LSRTM和基于Poynting矢量行波分離的LSRTM達(dá)到指定誤差閾值所需時(shí)間Fig.13 Time required for conventional LSRTM and LSRTM based on the Poynting-vector wavefield separation to reach the specified error threshold

為了展示在光滑的背景速度下背向反射的發(fā)育程度和波場的混疊程度，我們分別記錄了震源位于地表中央處時(shí)1.15 s時(shí)刻的震源波場和檢波器波場.如圖10a、c所示，因?yàn)榻缑娓浇淖杩共町惪s小，導(dǎo)致波場背向反射的能量弱化.這也是常規(guī)RTM得到的剖面上低頻噪聲并沒有十分嚴(yán)重的原因(如圖11a所示).為了展示背向反射，我們對震源波場快照進(jìn)行了自動(dòng)增益控制以增加較弱的反射能量，得到如圖10b所示的結(jié)果.從中，我們可以發(fā)現(xiàn)在部分區(qū)域，波場的混疊比較嚴(yán)重，如圖中所標(biāo)注的A和B區(qū)域.以A區(qū)域?yàn)槔?，兩種波發(fā)生了混疊但是都沿相同垂直向上方向，在此求出來的Poynting矢量仍能保持向上，因此其分離結(jié)果是準(zhǔn)確的.而在界面處(B區(qū)域)，入射波和反射波同時(shí)存在且垂直方向不同，但是因?yàn)槿肷洳芰棵黠@優(yōu)于反射波能量，其Poynting矢量方向與入射波傳播方向相近，因此采用(9)式或(10)式仍可以在此成像.

在此，我們分別使用RTM，常規(guī)LSRTM和基于Poynting矢量行波分離的LSRTM對地震數(shù)據(jù)進(jìn)行了反演成像，并將結(jié)果展示在圖11中.其中，LSRTM的迭代次數(shù)為20，我們的方法在前五次迭代中采用了(10)式作為梯度算子.相較于RTM方法，LSRTM的結(jié)果分辨率和信噪比得到了極大的改善.值得注意的是，與常規(guī)LSRTM結(jié)果相比，使用我們的LSRTM方法得到的剖面淺部噪聲更少，成像軸清晰.圖12展示了常規(guī)LSRTM和基于Poynting矢量的LSRTM在收斂性上的比較，結(jié)果發(fā)現(xiàn)我們的方案對于復(fù)雜的模型在收斂速率上依然優(yōu)于常規(guī)的LSRTM.

除此之外，我們還對兩種方案的運(yùn)行時(shí)間進(jìn)行了比較，結(jié)果發(fā)現(xiàn)：采用(9)式或者(10)式作為梯度算子，其運(yùn)行時(shí)間大致為采用全波成像條件(4)式的1.375倍.相較于傳統(tǒng)的全波LSRTM，我們的方案只在早期的迭代中增加運(yùn)算量.假設(shè)最大迭代次數(shù)是n, 我們在前m次迭代中以(9)式或者(10)式作為梯度算子，常規(guī)方法中每次迭代需要時(shí)間為t，所以常規(guī)LSRTM所需總時(shí)長為nt,而我們的方法所需要的總時(shí)長則為nt+0.375mt.當(dāng)n?m時(shí)，程序所增加的時(shí)間可忽略不計(jì).另外，將標(biāo)準(zhǔn)誤差閾值分別設(shè)置為0.5、0.4、0.3和0.2，我們計(jì)算了傳統(tǒng)LSRTM和基于Poynting矢量的LSRTM程序所需的時(shí)間長度，并將其展示如圖13.可見我們的方法在加快收斂速率上具有明顯效果.

3 結(jié)論

本文提出并實(shí)現(xiàn)了基于Poynting矢量行波分離的最小二乘逆時(shí)偏移方案.作為一種基于雙程波的地震成像方法，LSRTM面臨著背向反射所引發(fā)的強(qiáng)振幅噪聲和虛假軸的干擾.低頻噪聲來自于背向反射波與正常傳播的波之間的互相關(guān)，因此在梯度構(gòu)建時(shí)選擇不同垂直方向的波場組分進(jìn)行互相關(guān)可以有效壓制低頻噪聲.進(jìn)一步地，在成像條件中只保留震源波場的下行組分和檢波器波場的上行組分之間的互相關(guān)，雖然會使逆時(shí)偏移退化為類似于單程波偏移的方法，但是可以進(jìn)一步地去除偏移假象.

借助Poynting矢量的方向指示性，對震源波場和檢波器波場進(jìn)行上下行波分離具有簡單直接，方便快捷的優(yōu)點(diǎn).并且考慮到Poynting矢量對于重疊區(qū)域的效果不佳且不完全的成像條件也會損失波的部分有效信息，我們可以只在迭代的早期使用改造后的梯度算子.這樣我們既充分利用了地震波所攜帶的有效信息，又降低了方法所需的運(yùn)算成本.數(shù)值算例證明通過本文的方法，在減少迭代次數(shù)(運(yùn)算量)的基礎(chǔ)上，我們可以得到振幅均衡性好且信噪比更高的成像結(jié)果.