周陽， 曹俊興， 王興建， 胡江濤， 王華忠， 劉文卿

1 油氣藏地質(zhì)及開發(fā)工程國家重點(diǎn)實(shí)驗(yàn)室， 成都理工大學(xué), 成都 610059 2 同濟(jì)大學(xué)海洋與地球科學(xué)學(xué)院， 上海 200092 3 中國石油勘探開發(fā)研究院西北分院， 蘭州 730020

0 引言

通過求解雙向波方程，逆時(shí)偏移(RTM)能夠不受傳播傾角的限制對(duì)陡傾角地質(zhì)體進(jìn)行成像.同時(shí)，由于在波場反傳播過程中校正了幾何擴(kuò)散效應(yīng)，因此逆時(shí)偏移是一種經(jīng)典意義下的真振幅成像方法(Zhang and Sun, 2009).對(duì)于復(fù)雜地質(zhì)體的精確成像，除了波場傳播的幾何擴(kuò)散效應(yīng)，由于觀測采集系統(tǒng)不規(guī)則、不完備采樣以及地下復(fù)雜構(gòu)造等因素，在成像點(diǎn)處通常還會(huì)存在照明不均勻的現(xiàn)象 (Xie et al.，2006).因此，在逆時(shí)偏移中進(jìn)一步考慮照明補(bǔ)償會(huì)在較大程度上提高對(duì)地質(zhì)體的成像精度.

照明補(bǔ)償成像可以在偏移成像結(jié)果上作用Hessian算子或與其相關(guān)的點(diǎn)擴(kuò)散函數(shù)或照明算子的逆來實(shí)現(xiàn).Audebert等(2000)通過在角度域統(tǒng)計(jì)射線打擊次數(shù)近似照明算子實(shí)現(xiàn)帶照明補(bǔ)償?shù)目讼；舴蚱瞥上穹椒?，Gelius等(2002)、Lecomte(2008)通過射線計(jì)算點(diǎn)擴(kuò)散函數(shù)，并利用點(diǎn)擴(kuò)散函數(shù)對(duì)成像結(jié)果進(jìn)行反褶積來實(shí)現(xiàn)對(duì)成像結(jié)果照明不均勻的補(bǔ)償.Wu和Chen(2002)以及Xie等(2005)討論了基于單程波波動(dòng)方程的照明算子計(jì)算及對(duì)應(yīng)的照明補(bǔ)償方法.Wu等(2004)、陳生昌等(2007)、Cao和Wu(2009)進(jìn)一步將照明算子分解到入射/散射角度域進(jìn)行照明補(bǔ)償成像，取得了更優(yōu)的照明補(bǔ)償效果.Yan等(2014)利用逆時(shí)偏移在傾角域計(jì)算照明算子并提出了一套照明補(bǔ)償RTM方案.Yan和Xie(2016)利用Poynting矢量對(duì)波場進(jìn)行角度分解，并以此為基礎(chǔ)發(fā)展了一套角度域照明補(bǔ)償RTM方法.類似地，Hu等(2001)、Plessix和Mulder(2004)、Ren等(2011)討論了在最小二乘框架下基于顯式計(jì)算Hessian及其近似逆算子的照明補(bǔ)償成像方法.上述波動(dòng)方程類照明補(bǔ)償成像方法通常需要利用波動(dòng)方程顯式計(jì)算格林函數(shù)，對(duì)于大規(guī)模的三維成像問題，在計(jì)算量和存儲(chǔ)量上具有不小的挑戰(zhàn).

為了提高利用波動(dòng)方程進(jìn)行照明補(bǔ)償成像的計(jì)算效率，Rickett(2003)、Guitton(2004)提出了利用偏移-反偏移的方式估計(jì)近似Hessian算子進(jìn)行照明補(bǔ)償成像，該方式通過增加一次反偏移和一次偏移的計(jì)算代價(jià)來避免對(duì)格林函數(shù)的顯式計(jì)算和儲(chǔ)存.Gherasim等(2010)、Shen等(2011)、Li等(2012)將此種方式推廣到了角度域.相較于直接計(jì)算格林函數(shù)構(gòu)造照明算子或Hessian算子，利用偏移-反偏移的方式進(jìn)行照明補(bǔ)償成像能在一定程度上提高計(jì)算效率.對(duì)于基于RTM的照明補(bǔ)償成像，由于一次三維RTM偏移的計(jì)算量較為龐大(劉紅偉等，2010；許璐等，2019)，因此利用偏移-反偏移的方式進(jìn)行大規(guī)模數(shù)據(jù)RTM照明補(bǔ)償成像在計(jì)算量上仍然面臨較大壓力.

另一方面，在最小二乘框架下也可以利用迭代求取Hessian逆算子或其近似算子進(jìn)行成像的照明補(bǔ)償，即迭代實(shí)現(xiàn)的最小二乘類偏移成像方法技術(shù).該類方法技術(shù)的發(fā)展經(jīng)歷了射線類迭代最小二乘偏移成像方法(Schuster，1993；Nemeth et al.，1999；Duque et al.，2000；王彥飛等，2009；黃建平等，2013)、單程波算子類迭代最小二乘偏移成像方法(Kühl and Sacchi, 2003; 楊其強(qiáng)和張叔倫，2008; 沈雄君和劉能超，2012；周華敏等，2014)以及迭代實(shí)現(xiàn)的最小二乘RTM偏移成像方法(Dai et al.，2012；劉學(xué)建和劉依克，2016；李振春等，2017a，b；方修政等，2018；張攀和毛偉建，2018；周東紅等，2020).在實(shí)際資料處理中，利用迭代方式估計(jì)Hessian逆算子的精度和穩(wěn)定性易受到噪聲不滿足先驗(yàn)分布、震源子波空變且未知、初始模型不滿足近似線性化反演要求等多種因素影響，從而使得迭代最小二乘類成像方法出現(xiàn)迭代收斂慢甚至不收斂的情況，在一定程度上增加了該方法的實(shí)際應(yīng)用難度.

本文提出一種基于RTM實(shí)現(xiàn)的高效照明補(bǔ)償成像方法.在高頻近似意義下，首先從波傳播的幾何路徑出發(fā)，定義衡量成像點(diǎn)處照明均勻的雅克比矩陣，并以此為基礎(chǔ)進(jìn)一步給出照明補(bǔ)償算子的表達(dá)式.然后通過波場的邊界積分表達(dá)，將震源端照明補(bǔ)償算子利用外推波場進(jìn)行表達(dá)，避免了震源端照明算子中對(duì)格林函數(shù)顯式計(jì)算的需求.對(duì)于檢波點(diǎn)端照明補(bǔ)償算子，通過引入對(duì)參數(shù)場的合理近似，將檢波點(diǎn)端照明補(bǔ)償算子中對(duì)檢波器坐標(biāo)積分項(xiàng)進(jìn)行顯式表達(dá)，在較大程度上提高了計(jì)算效率.利用合成數(shù)據(jù)和實(shí)際數(shù)據(jù)驗(yàn)證了本文提出方法的正確性.

1 方法原理

本小節(jié)給出基于RTM的高效照明補(bǔ)償成像方法.首先引入衡量地下照明規(guī)則程度的雅克比行列式，然后基于該雅克比行列式給出照明補(bǔ)償算子的定義.進(jìn)一步，討論利用RTM計(jì)算引擎高效實(shí)現(xiàn)該照明補(bǔ)償算子的具體算法.

1.1 照明補(bǔ)償算子

為了引入照明雅克比行列式，首先對(duì)下文推導(dǎo)中涉及的符號(hào)進(jìn)行定義說明，文中涉及主要符號(hào)的幾何定義如圖1所示.圖1a為本文涉及的射線傳播參數(shù)定義.其中ps和pr分別為源端和檢波點(diǎn)端在地下成像點(diǎn)x處射線慢度矢量，其定義為震源及檢波點(diǎn)兩端波場走時(shí)場梯度.nqr為照明傾角矢量，其定義為ps與ps的矢量和(Forgues and Lambaré, 1997).ps和pr與豎直方向的夾角βs、βr分別定義為震源及檢波點(diǎn)兩端的起飛角(Operto et al., 2000).ps和pr之間的夾角θs為成像點(diǎn)處的散射角，代表當(dāng)前觀測系統(tǒng)下該成像點(diǎn)的角度照明范圍，照明傾角矢量nqr與豎直方向夾角ψx為照明傾角，代表當(dāng)前觀測系統(tǒng)下角度照明范圍的中心位置(Audebert et al., 2005).圖1b為圖1a在成像點(diǎn)x處的局部放大圖.

圖1 (a) 射線傳播參數(shù)表征及角度定義.ps、pr分別為源檢兩端出射射線在地下點(diǎn)x處的慢度矢量，nqr=ps+pr為照明傾角矢量，θs地下點(diǎn)x處的散射角且θs=βr-βs，βs、βr分別為地下點(diǎn)x處的起飛角，ψ為照明傾角且ψ= (βr+βs)/2；(b)圖(a)在x處局部放大顯示Fig.1 (a) Ray parameters for 2D asymptotic inversion imaging and illumination compensation. ps and pr are slowness vectors coming from the source and the receiver, respectively. nqr is the normalized illumination dip vector. βs and βr are the takeoff angles on source and receiver side, respectively. θs=βr-βs is the scattering angle. ψ=(βr+ βs)/2 is the illumination dip angle； (b) Enlargement of (a) at x

對(duì)于地下成像點(diǎn)處照明均勻程度，本文利用局部角度域擾動(dòng)相對(duì)于地表坐標(biāo)擾動(dòng)量大小來進(jìn)行度量.圖2展示了在直射線假設(shè)情況下利用層狀介質(zhì)和局部速度異常介質(zhì)對(duì)本文給出照明均勻程度定義的幾何解釋.圖2a 為層狀介質(zhì)直射線假設(shè)下不同深度局部角度域擾動(dòng)與地表坐標(biāo)擾動(dòng)的關(guān)系示意圖.從圖2a可以看出當(dāng)?shù)乇碜鴺?biāo)具有相同擾動(dòng)dx時(shí)，深部反射地層局部角度域擾動(dòng)dψ2要小于淺部擾動(dòng)dψ1.這意味著對(duì)于固定觀測系統(tǒng)，深部反射地層的角度覆蓋范圍要小于淺部反射地層覆蓋范圍，為了使得不同深度反射地層具有一致的照明強(qiáng)度，深部反射地層需要更多的照明補(bǔ)償.圖2b、c分別為存在及不存在局部高速異常時(shí)照明射線路徑圖.可以看出當(dāng)?shù)叵戮植拷嵌扔虼嬖谙嗤瑪_動(dòng)dψ1時(shí)，由于速度異常體對(duì)照明射線的偏折作用，存在速度異常體模型中照明射線在地表坐標(biāo)擾動(dòng)量dx2要遠(yuǎn)大于勻速介質(zhì)中照明射線在地表坐標(biāo)擾動(dòng)量dx1.這表明對(duì)于固定偏移距的觀測系統(tǒng)，速度異常體下方更加難以被地震波場照明，需要更多的照明補(bǔ)償才能夠消除照明不均衡對(duì)成像結(jié)果的影響.據(jù)以上分析，本文定義二維情況下衡量成像點(diǎn)處照明規(guī)則程度的雅克比行列式為：

圖2 (a) 在層狀介質(zhì)直射線假設(shè)下，不同深度局部角度與擾動(dòng)域地表擾動(dòng)的關(guān)系示意圖； (b) 在常速介質(zhì)中局部角度域擾動(dòng)與地表擾動(dòng)關(guān)系示意圖； (c) 在包含局部異常情景下，局部角度域擾動(dòng)與地表擾動(dòng) 關(guān)系示意圖Fig.2 (a) Illustration of relation between local angle and surface coordinate perturbation in layered media with straight ray approximation； (b) Illustration of relation between local angle and surface coordinate perturbation in constant-velocity media； (c) Illustration of relation between local angle and surface coordinate perturbation in constant-velocity media with local anomaly

(1)

其中ψx為成像點(diǎn)x處的照明傾角，|dψx/dxs|及|dψx/dxr|分別代表照明傾角相對(duì)地表炮點(diǎn)坐標(biāo)xs及檢波器坐標(biāo)xr擾動(dòng)的變化.注意到式(1)右端行列式結(jié)果為標(biāo)量，標(biāo)量Gs和Gr即為本文給出震源和檢波點(diǎn)端衡量照明規(guī)則程度的雅克比行列式.

從式(1)可以看出，給定相同照明傾角擾動(dòng)量dψx，如果地表坐標(biāo)dxs及dxr擾動(dòng)越大，則Gs和Gr越小，代表對(duì)于固定地震采集系統(tǒng)，地下照明規(guī)則程度越差，所需的照明補(bǔ)償量越大.基于以上分析，由式(1)出發(fā)可以進(jìn)一步定義成像點(diǎn)x處的照明補(bǔ)償算子Wx的形式為：

(2)

公式(2)中對(duì)檢波器的積分代表考慮共炮數(shù)據(jù)中所有檢波器接收數(shù)據(jù)的照明補(bǔ)償.對(duì)照公式(1)和(2)可以看出，當(dāng)震源和檢波點(diǎn)兩端照明越不規(guī)則，照明補(bǔ)償算子Wx越大，照明補(bǔ)償效果越強(qiáng).

公式(2)即為本文給出的照明補(bǔ)償算子抽象形式，為了能夠在RTM過程中高效地實(shí)施照明補(bǔ)償算子(2)，需要對(duì)公式(2)中元素做更具體的表征.利用圖1 所示的照明傾角與起飛角的關(guān)系ψx=(βr+βs)/2，并注意到震源端射線慢度矢量對(duì)應(yīng)的起飛角βs僅僅與炮點(diǎn)坐標(biāo)xs有關(guān)，檢波點(diǎn)端對(duì)應(yīng)的起飛角βr僅僅與檢波點(diǎn)坐標(biāo)xr有關(guān)，可以將式(2)改寫為：

(3)

據(jù)Bleistein 等(2005)，在高頻近似下，起飛角擾動(dòng)相對(duì)于地表坐標(biāo)擾動(dòng)關(guān)系為：

(4)

其中A(x,xs)為高頻近似下波場從炮點(diǎn)xs出發(fā)傳播到成像點(diǎn)x處的振幅.c0(xs)為炮點(diǎn)xs處的背景速度.

對(duì)于檢波點(diǎn)端起飛角相對(duì)地表坐標(biāo)的擾動(dòng)關(guān)系，有類似式(4)的表達(dá)形式.將式(4)代入式(3)最終可得照明補(bǔ)償算子Wx的具體表達(dá)形式為：

Wx=4IsIr,

(5)

其中：

(6)

式(6)中A(xr,x)為高頻近似下波場經(jīng)過成像點(diǎn)x反射被位于xr檢波點(diǎn)接收到的振幅，c0(xs)為檢波點(diǎn)xr處的背景速度.

1.2 照明補(bǔ)償算子的RTM實(shí)現(xiàn)

(7)

將式(7)代入式(6)并做適當(dāng)變形可得：

(8)

式(8)即為利用外推波場表達(dá)的源端照明補(bǔ)償校正算子.可以看出，由于源端波場PF(x,xs,ω)在波場外推中已知，式(8)的數(shù)值實(shí)現(xiàn)僅需額外計(jì)算成像點(diǎn)處震源端波場的起飛角βs，起飛角βs可以利用坡印廷矢量在波場外推中快速計(jì)算(Yoon and Marfurt，2006)，因此照明補(bǔ)償算子式(8)可以在震源端波場外推過程中高效地?cái)?shù)值實(shí)現(xiàn).對(duì)于檢波點(diǎn)端照明補(bǔ)償算子Gr，受到Plessix 和 Mulder(2004)工作的啟發(fā)，采用緩變介質(zhì)假設(shè)，將式(6)中對(duì)檢波器坐標(biāo)積分項(xiàng)進(jìn)行顯式表達(dá).緩變介質(zhì)的格林函數(shù)振幅項(xiàng)(Bleistein et al.，2001)為：

(9)

將式(9)代入式(6)并做適當(dāng)代數(shù)運(yùn)算，最終可以得到檢波點(diǎn)端照明補(bǔ)償算子的形式為：

(10)

2 數(shù)值實(shí)驗(yàn)

2.1 模型資料測試

首先利用一個(gè)四層速度模型來驗(yàn)證本文提出的照明補(bǔ)償方法的有效性.將利用本文給出的照明補(bǔ)償算子成像結(jié)果與Plessix 和Mulder (2004)經(jīng)典的基于對(duì)角Hessian照明補(bǔ)償算子成像結(jié)果進(jìn)行比較，特別地，對(duì)比Plessix 和 Mulder(2004)文中給出的第三類照明補(bǔ)償算子成像結(jié)果，該類照明補(bǔ)償算子相較于文中給出的其他幾類照明補(bǔ)償算子有更好的補(bǔ)償效果 (Kiyashchenko et al.，2007).圖 3a 為四層速度真實(shí)模型，四層速度分別為3000 m·s-1、 4000 m·s-1、 5000 m·s-1和 6000 m·s-1.圖3b 展示了成像算法采用的平滑模型，該平滑模型通過利用半徑為10的高斯平滑濾波器獲得.圖 3c 為在常密度假設(shè)情況下利用速度模型計(jì)算的真實(shí)零角度反射系數(shù)(Berteussen and Ursin，1983).為了更好地對(duì)比分析不同照明補(bǔ)償算子對(duì)成像結(jié)果的影響，計(jì)算了數(shù)值實(shí)驗(yàn)采用的觀測系統(tǒng)的總照明補(bǔ)償強(qiáng)度，如圖4所示.圖4a、b分別為利用Plessix和Mulder (2004)給出的照明補(bǔ)償算子及本文給出的照明補(bǔ)償算子計(jì)算得到的總照明補(bǔ)償強(qiáng)度.可以看出，本文方法計(jì)算的總照明補(bǔ)償強(qiáng)度分布和直接利用近似對(duì)角Hessian逆算子計(jì)算出的照明補(bǔ)償強(qiáng)度(Plessix and Mulder，2004)分布整體上具有較好的一致性，這在某種程度上驗(yàn)證了本文給出的照明補(bǔ)償算子確為對(duì)角Hessian逆算子在高頻近似下的一個(gè)近似.此外，從圖 4a、b的對(duì)比可以看出，相較于Plessix和Mulder(2004)給出的對(duì)角Hessian 逆算子的近似方式，利用本文方法計(jì)算得到的總照明補(bǔ)償強(qiáng)度和速度關(guān)系更加密切，不同深度照明傾角信息體現(xiàn)也更加明顯.

圖3 (a) 四層真實(shí)速度模型； (b) 四層平滑速度模型； (c) 由速度場計(jì)算得到的真實(shí)零角度反射系數(shù)模型Fig.3 (a) The true four-layer velocity model; (b) The smoothed four-layer velocity for imaging; (c) The true normal reflection coefficient for four-layer model

圖4 四層速度模型的總照明補(bǔ)償強(qiáng)度對(duì)比 (a) 利用Plessix 和 Mulder (2004)給出的照明補(bǔ)償算子計(jì)算得到的照明補(bǔ)償強(qiáng)度； (b) 利用本文給出照明補(bǔ)償算子(8)和(10)式計(jì)算得到的照明補(bǔ)償強(qiáng)度.Fig.4 Comparison of total illumination compensation intensity for four-layer model (a) The illumination compensation intensity calculated using the operators proposed in Plessix and Mulder (2004)； (b) The illumination compensation intensity calculated using the operators (8) and (10).

圖5a—c分別展示了利用不帶照明補(bǔ)償?shù)腞TM算法、利用Plessix和Mulder(2004)照明補(bǔ)償算子以及本文照明補(bǔ)償算子(8)和(10)進(jìn)行照明補(bǔ)償成像結(jié)果.在不同照明補(bǔ)償策略的RTM成像算法中，本文利用逆散射成像條件來消除RTM中的低波數(shù)噪聲(Whitmore and Crawley, 2012; Zhou and Wang, 2017).可以看出，相較于不帶照明補(bǔ)償?shù)腞TM算法，利用Plessix 和Mulder(2004)以及本文的照明補(bǔ)償算子進(jìn)行照明補(bǔ)償成像均能有效提升深層成像能量，使得成像剖面能量更加均衡.為了定量地對(duì)比不同照明補(bǔ)償算子效果，對(duì)利用不同照明補(bǔ)償策略得到的成像剖面進(jìn)行抽道并與真實(shí)零角度反射系數(shù)進(jìn)行對(duì)比，結(jié)果展示在圖5d—f中.從圖5d可以看出，深層反射層由于照明不足，不做照明補(bǔ)償直接計(jì)算得到反射系數(shù)峰值振幅與真實(shí)值誤差較大.利用Plessix 和 Mulder (2004)給出照明補(bǔ)償算子計(jì)算得到深層反射系數(shù)峰值振幅與理論值的匹配程度有較大程度的提升，如圖5e所示.利用本文給出照明補(bǔ)償算子計(jì)算得到的反射系數(shù)抽道對(duì)比結(jié)果展示在圖5f中，可以看出利用本文給出照明補(bǔ)償算子能進(jìn)一步提升深層反射系數(shù)計(jì)算精度.為了測試本文給出照明補(bǔ)償算子(8)和(10)在數(shù)據(jù)包含噪聲情況下的計(jì)算穩(wěn)定性，在模擬炮數(shù)據(jù)中加入滿足高斯分布的噪聲，原始數(shù)據(jù)與含噪數(shù)據(jù)的抽道對(duì)比如圖6a所示.圖 6c、d分別為利用本文照明補(bǔ)償算子在原始和含噪數(shù)據(jù)上單炮成像結(jié)果.可以看出相較于不包含噪聲的原始數(shù)據(jù)，含噪聲數(shù)據(jù)單炮照明補(bǔ)償成像結(jié)果包含由數(shù)據(jù)噪聲引起的額外干擾，這些額外成像干擾能量相較于真實(shí)反射層成像結(jié)果較弱.在所有炮數(shù)據(jù)成像結(jié)果疊加后，這些由數(shù)據(jù)噪聲引起額外成像干擾會(huì)被進(jìn)一步壓制.圖6b為利用本文照明補(bǔ)償算子在原始和含噪數(shù)據(jù)上所有炮數(shù)據(jù)疊加成像結(jié)果的抽道對(duì)比，可以較為清晰看出最終照明補(bǔ)償成像結(jié)果幾乎不受數(shù)據(jù)噪聲的影響，這驗(yàn)證了本文給出照明補(bǔ)償算子對(duì)含噪數(shù)據(jù)的數(shù)值計(jì)算穩(wěn)健性.

圖5 四層模型不同照明補(bǔ)償策略成像結(jié)果對(duì)比 (a) 不帶照明補(bǔ)償補(bǔ)償?shù)腞TM成像結(jié)果； (b) 利用Plessix 和 Mulder(2004)給出的照明補(bǔ)償算子得到的RTM成像結(jié)果； (c) 利用本文給出的照明補(bǔ)償算子得到的RTM成像結(jié)果； (d)—(f) 分別為(a)—(c)抽取中間道的成像結(jié)果(紅色實(shí)線)與真實(shí)零角度反射系數(shù) (藍(lán)色 實(shí)線)對(duì)比.Fig.5 Comparisons of imaging results with different illumination compensation operators for four-layer model (a), (b) and (c) Imaging results without illumination, with compensation operator of Plessix and Mulder (2004) and the compensation operator proposed in this paper, respectively；(d), (e) and (f) Comparison for selected central traces between imaging results (red solid line) and theoretical normal reflection coefficient values (blue solid line) corresponding to (a),(b) and (c), respectively.

圖6 四層速度模型包含噪聲數(shù)據(jù)測試 (a) 包含噪數(shù)據(jù)(藍(lán)色實(shí)線)與不含噪數(shù)據(jù)(紅色實(shí)線)抽道對(duì)比； (b) 包含噪數(shù)據(jù)疊加成像結(jié)果(藍(lán)色實(shí)線)與 不含音數(shù)據(jù)疊加成像結(jié)果(紅色實(shí)線)抽道對(duì)比； (c) 不含噪數(shù)據(jù)單炮成像結(jié)果； (d) 包含噪數(shù)據(jù)單炮成像結(jié)果.Fig.6 Noisy data test for four-layer model (a) Comparison of selected traces for noisy data (blue solid line) and clean data (red solid line)； (b) Comparison of selected imaged traces with noisy data (blue solid line) and clean data (red solid line)； (c) Imaging result with single shot clean data； (d) Imaging result with single shot noisy data.

為了進(jìn)一步驗(yàn)證本文照明補(bǔ)償成像方法對(duì)復(fù)雜模型的有效性，利用Marmousi模型進(jìn)行了照明補(bǔ)償成像測試.圖 7a、b分別為Marmousi真實(shí)速度以及零角度反射系數(shù)模型，利用半徑為10的高斯平滑濾波器作用在真實(shí)模型上得到平滑模型用于成像.利用不同照明補(bǔ)償算子計(jì)算得到總照明補(bǔ)償強(qiáng)度展示在圖8中.和四層速度模型測試結(jié)論類似，利用本文給出照明補(bǔ)償算子及Plessix和Mulder(2004)給出照明補(bǔ)償算子計(jì)算得到的總照明補(bǔ)償強(qiáng)度在整體分布具有較好的一致性，同時(shí)利用本文給出照明補(bǔ)償算子計(jì)算得到總照明補(bǔ)償強(qiáng)度和速度構(gòu)造有更好的一致性.圖 9a—c分別為不帶照明補(bǔ)償、利用Plessix和Mulder (2004)給出補(bǔ)償算子和利用本文給出照明補(bǔ)償算子進(jìn)行基于RTM的照明補(bǔ)償成像疊加結(jié)果.對(duì)于不做照明補(bǔ)償?shù)腞TM成像結(jié)果，可以清晰看出由于震源附近強(qiáng)波場能量相關(guān)導(dǎo)致的成像噪聲.通過在RTM過程中施加Plessix和Mulder(2004)給出的照明補(bǔ)償算子進(jìn)行成像，震源附近強(qiáng)能量噪聲得到了較好的壓制，整個(gè)成像剖面不同深度能量更加均衡，如圖9b所示.利用本文給出的照明補(bǔ)償算子成像疊加結(jié)果展示在圖9c中，可以較明顯看出震源相關(guān)噪聲得到了進(jìn)一步壓制，深層成像振幅得到了進(jìn)一步提升.圖9d—f為圖9a—c對(duì)應(yīng)的局部放大結(jié)果，可以看出利用本文給出的照明補(bǔ)償算子對(duì)復(fù)雜構(gòu)造區(qū)域由于波場未完全干涉形成的成像噪聲也具有很好的壓制作用，可以得到信噪比更好的復(fù)雜區(qū)域成像結(jié)果. Marmousi模型不同照明補(bǔ)償策略成像結(jié)果抽道對(duì)比如圖10所示.圖10a為不做照明補(bǔ)償成像結(jié)果和真實(shí)零角度反射系數(shù)對(duì)比，可以清晰地看出淺層震源處強(qiáng)相關(guān)噪聲對(duì)成像結(jié)果的影響.同時(shí)，隨著成像深度的增加，不做照明補(bǔ)償成像結(jié)果振幅峰值與理論零角度反射系數(shù)誤差逐漸增大，這與有效照明范圍隨著成像深度增加逐步減小的趨勢是一致的.利用Plessix和Mulder(2004)給出照明補(bǔ)償算子成像抽道對(duì)比結(jié)果顯示在圖10b中，可以看出利用該照明補(bǔ)償算子能明顯壓制淺層震源強(qiáng)相關(guān)噪聲，有效提升中深層成像結(jié)果峰值振幅與真實(shí)零角度反射系數(shù)匹配程度.圖10c為利用本文照明補(bǔ)償算子(8)和(10)成像抽道對(duì)比結(jié)果，相較于圖10b，淺層震源強(qiáng)相關(guān)噪聲得到了進(jìn)一步衰減，中深層成像結(jié)果峰值振幅更加接近于理論零角度反射系數(shù).

圖7 (a) Marmousi真實(shí)速度模型； (b) Marmousi模型真實(shí)反射系數(shù)Fig.7 (a) The true Marmousi velocity model； (b) The true reflection coefficient for Marmousi model

圖8 Marmousi速度模型的總照明補(bǔ)償強(qiáng)度對(duì)比 (a) 利用Plessix 和 Mulder (2004)給出的照明補(bǔ)償算子計(jì)算得到的照明補(bǔ)償強(qiáng)度； (b) 利用本文給出照明補(bǔ)償算子(8)和(10)式計(jì)算得到的照明補(bǔ)償強(qiáng)度.Fig.8 Comparison of total illumination compensation intensity for Marmousi model (a) The illumination compensation intensity calculated using the operators proposed by Plessix and Mulder (2004)； (b) The illumination compensation intensity calculated using the operators (8) and (10).

圖9 Marmousi模型不同照明補(bǔ)償策略成像結(jié)果對(duì)比 (a) 不帶照明補(bǔ)償補(bǔ)償?shù)腞TM成像結(jié)果； (b) 利用Plessix 和 Mulder(2004)給出的照明補(bǔ)償算子得到的RTM成像結(jié)果； (c) 利用本文給出的照明補(bǔ)償算子得到的RTM成像結(jié)果； (d)—(f) 分別為(a)—(c)的局部放大對(duì)比圖.Fig.9 Comparison of imaging results with different illumination compensation operators for Marmousi model (a), (b) and (c)Imaging results without illumination, with compensation operator of Plessix and Mulder (2004) and the compensation operator proposed in this paper, respectively. (d), (e) and (f) Enlarged display of (a),(b) and (c), respectively.

圖10 Marmousi 模型不同照明補(bǔ)償策略成像5250 m處抽道對(duì)比結(jié)果 (a) 不帶照明補(bǔ)償補(bǔ)償?shù)腞TM成像結(jié)果(紅色實(shí)線)與理論零角度反射系數(shù)(藍(lán)色實(shí)線)對(duì)比結(jié)果； (b) 基于Plessix 和 Mulder(2004)給出照明補(bǔ)償算子RTM成像結(jié)果(紅色實(shí)線)與理論零角度反射系數(shù)(藍(lán)色實(shí)線)對(duì)比結(jié)果； (c) 基于本文給出的照明補(bǔ)償 算子RTM成像結(jié)果(紅色實(shí)線)與理論零角度反射系數(shù)(藍(lán)色實(shí)線)對(duì)比結(jié)果.Fig.10 Comparison of selected traces between imaging results (red solid line) and theoretical normal reflection coefficient values (blue solid line) (a) RTM imaging results calculated without illumination compensation； (b) RTM imaging results calculated with illumination compensation operator of Plessix and Mulder (2004)； (c) RTM imaging results calculated with illumination compensation operator proposed in this paper.

2.2 實(shí)際資料測試

為了檢驗(yàn)本文照明補(bǔ)償成像方法在實(shí)際資料上的應(yīng)用效果，利用西部某二維實(shí)際資料進(jìn)行測試.圖11展示了實(shí)際資料的速度場.利用不同照明補(bǔ)償策略得到的RTM成像疊加剖面展示在圖12a—c中，疊加成像部分放大結(jié)果展示在圖12d—f中.可以看出，未進(jìn)行照明補(bǔ)償RTM成像疊加結(jié)果整個(gè)剖面存在由于波場能量未完全干涉造成的成像噪聲，且深層反射能量不聚焦.利用Plessix 和Mulder(2004)給出照明補(bǔ)償算子進(jìn)行照明補(bǔ)償成像能夠較好地提升中深層反射能量，減少非相干成像噪聲干擾.利用本文提出的照明補(bǔ)償算子進(jìn)行照明補(bǔ)償成像在成像噪聲壓制，成像剖面能量均衡以及深層反射層成像結(jié)果聚焦性等方面有較明顯改進(jìn)效果，相較于利用Plessix和Mulder(2004)給出照明補(bǔ)償算子計(jì)算結(jié)果在成像質(zhì)量方面有進(jìn)一步的提升.

圖11 中國西部某實(shí)際數(shù)據(jù)速度場Fig.11 Velocity field of real data from western China

圖12 (a) 不帶照明補(bǔ)償補(bǔ)償?shù)腞TM成像結(jié)果； (b) 利用Plessix 和 Mulder(2004)給出的照明補(bǔ)償算子得到的RTM 成像結(jié)果; (c) 利用本文的照明補(bǔ)償算子得到的RTM成像結(jié)果; (d)—(f) 分別為(a)—(c)對(duì)應(yīng)的局部放大顯示結(jié)果Fig.12 (a), (b) and (c) Imaging results without illumination, with illumination compensation operator given by Plessix and Mulder (2004) and the compensation operator proposed in this paper, respectively;(d), (e) and (f) Zoomed display of partial imaging result of (a), (b) and (c), respectively.

3 討論

本文給出的照明補(bǔ)償算子為對(duì)角Hessian逆算子的一個(gè)近似，在高頻近似意義下矯正了成像點(diǎn)處照明不均勻的影響，提高了成像精度.由于本文給出的照明補(bǔ)償算子并非完整Hessian算子的逆，因此在利用本文的照明補(bǔ)償算子進(jìn)行照明補(bǔ)償成像后，會(huì)有殘余的波場傳播效應(yīng)有待進(jìn)一步消除.對(duì)由于照明補(bǔ)償算子近似Hessian逆算子不精確造成的殘留波場傳播效應(yīng)問題，可將本文給出照明補(bǔ)償算子融入到迭代求解的最小二乘類成像框架中，利用其迭代最小二乘成像方法中構(gòu)造的預(yù)條件梯度項(xiàng)，在后續(xù)的迭代過程中消除殘余的波場傳播效應(yīng)，進(jìn)一步提高成像精度.Burgess和Warner(2015)利用Plessix和Mulder(2004)給出的照明補(bǔ)償算子構(gòu)造預(yù)條件梯度項(xiàng)并應(yīng)用到全波形反演(FWI)中，成功加速了FWI收斂效率.類似地，從數(shù)值實(shí)驗(yàn)對(duì)比可以看出，利用本文給出照明補(bǔ)償算子構(gòu)造預(yù)條件梯度項(xiàng)也有望提高最小二乘類成像方法的收斂速度，這是本文后續(xù)研究方向之一.

4 結(jié)論

本文給出了一種新的RTM照明補(bǔ)償方法.該照明補(bǔ)償方法利用地下角度域相對(duì)于地表坐標(biāo)的擾動(dòng)量來刻畫局部成像點(diǎn)的照明均勻程度，據(jù)此定義照明補(bǔ)償算子進(jìn)行成像過程中的照明補(bǔ)償.通過高頻近似下的波場邊界積分表達(dá)及對(duì)介質(zhì)的合理近似，給出了照明補(bǔ)償算子基于RTM的高效實(shí)現(xiàn)方式.模型和實(shí)際資料的測試對(duì)比表明，本文提出的照明補(bǔ)償算子能夠在較大程度上提升RTM成像質(zhì)量.