一種基于逆時偏移的高效照明補償成像方法

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

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

0 引言

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

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

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

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

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

1 方法原理

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

1.1 照明補償算子

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

圖1 (a) 射線傳播參數(shù)表征及角度定義.ps、pr分別為源檢兩端出射射線在地下點x處的慢度矢量，nqr=ps+pr為照明傾角矢量，θs地下點x處的散射角且θs=βr-βs，βs、βr分別為地下點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

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

圖2 (a) 在層狀介質直射線假設下，不同深度局部角度與擾動域地表擾動的關系示意圖； (b) 在常速介質中局部角度域擾動與地表擾動關系示意圖； (c) 在包含局部異常情景下，局部角度域擾動與地表擾動 關系示意圖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為成像點x處的照明傾角，|dψx/dxs|及|dψx/dxr|分別代表照明傾角相對地表炮點坐標xs及檢波器坐標xr擾動的變化.注意到式(1)右端行列式結果為標量，標量Gs和Gr即為本文給出震源和檢波點端衡量照明規(guī)則程度的雅克比行列式.

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

(2)

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

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

(3)

據(jù)Bleistein 等(2005)，在高頻近似下，起飛角擾動相對于地表坐標擾動關系為：

(4)

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

對于檢波點端起飛角相對地表坐標的擾動關系，有類似式(4)的表達形式.將式(4)代入式(3)最終可得照明補償算子Wx的具體表達形式為：

Wx=4IsIr,

(5)

其中：

(6)

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

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

(7)

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

(8)

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

(9)

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

(10)

2 數(shù)值實驗

2.1 模型資料測試

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

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

圖5 四層模型不同照明補償策略成像結果對比 (a) 不帶照明補償補償?shù)腞TM成像結果； (b) 利用Plessix 和 Mulder(2004)給出的照明補償算子得到的RTM成像結果； (c) 利用本文給出的照明補償算子得到的RTM成像結果； (d)—(f) 分別為(a)—(c)抽取中間道的成像結果(紅色實線)與真實零角度反射系數(shù) (藍色 實線)對比.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ù)(藍色實線)與不含噪數(shù)據(jù)(紅色實線)抽道對比； (b) 包含噪數(shù)據(jù)疊加成像結果(藍色實線)與 不含音數(shù)據(jù)疊加成像結果(紅色實線)抽道對比； (c) 不含噪數(shù)據(jù)單炮成像結果； (d) 包含噪數(shù)據(jù)單炮成像結果.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.

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

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

圖8 Marmousi速度模型的總照明補償強度對比 (a) 利用Plessix 和 Mulder (2004)給出的照明補償算子計算得到的照明補償強度； (b) 利用本文給出照明補償算子(8)和(10)式計算得到的照明補償強度.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模型不同照明補償策略成像結果對比 (a) 不帶照明補償補償?shù)腞TM成像結果； (b) 利用Plessix 和 Mulder(2004)給出的照明補償算子得到的RTM成像結果； (c) 利用本文給出的照明補償算子得到的RTM成像結果； (d)—(f) 分別為(a)—(c)的局部放大對比圖.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 模型不同照明補償策略成像5250 m處抽道對比結果 (a) 不帶照明補償補償?shù)腞TM成像結果(紅色實線)與理論零角度反射系數(shù)(藍色實線)對比結果； (b) 基于Plessix 和 Mulder(2004)給出照明補償算子RTM成像結果(紅色實線)與理論零角度反射系數(shù)(藍色實線)對比結果； (c) 基于本文給出的照明補償 算子RTM成像結果(紅色實線)與理論零角度反射系數(shù)(藍色實線)對比結果.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 實際資料測試

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

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

圖12 (a) 不帶照明補償補償?shù)腞TM成像結果； (b) 利用Plessix 和 Mulder(2004)給出的照明補償算子得到的RTM 成像結果; (c) 利用本文的照明補償算子得到的RTM成像結果; (d)—(f) 分別為(a)—(c)對應的局部放大顯示結果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 討論

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

4 結論

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