基于高頻近似波傳播的實(shí)用化保真RTM成像方法

周 陽,王華忠

(1.成都理工大學(xué)地球物理學(xué)院,四川成都610059;2.波現(xiàn)象與智能反演成像研究組(WPI),同濟(jì)大學(xué)海洋與地球科學(xué)學(xué)院,上海200092)

地震波反演成像是勘探地球物理學(xué)的核心環(huán)節(jié)[1],其終極目標(biāo)為獲得地下全波數(shù)帶的參數(shù)場。全波形反演(FWI)是理論完備的地震波反演成像方法,理想情況下能得到地下全頻帶參數(shù)場的估計(jì)。經(jīng)典的FWI是強(qiáng)非線性反演問題,易受到各種客觀實(shí)際因素制約而出現(xiàn)收斂到局部極值或不收斂的情況。因此,在實(shí)際地震數(shù)據(jù)處理中,通常將FWI分解為一系列更為線性化的子問題求解[2-3]。估計(jì)背景參數(shù)場的層析成像和估計(jì)高波數(shù)參數(shù)場的保真成像是最為常見的分解方法。本文主要關(guān)注估計(jì)高波數(shù)參數(shù)擾動(dòng)的保真成像方法。

最小二乘偏移是一種具有代表性的保真成像技術(shù),其試圖通過匹配模擬數(shù)據(jù)和觀測數(shù)據(jù)的振幅值來估計(jì)參數(shù)的高波數(shù)擾動(dòng)。然而,最小二乘偏移的理想實(shí)現(xiàn)需將Hessian逆算子作用到常規(guī)成像結(jié)果上。在三維實(shí)際資料處理過程中,由于Hessian矩陣規(guī)模巨大,通常采用迭代方式求取Hessian的逆算子,因而使得最小二乘偏移成像成為一種計(jì)算昂貴的保真成像方法。另外,Hessian逆算子的估計(jì)精度很容易受到疊前數(shù)據(jù)觀測不完備、正算子對物理波現(xiàn)象的描述不夠準(zhǔn)確以及地震子波未知等多種因素影響,從而出現(xiàn)迭代收斂慢甚至不收斂情況,這進(jìn)一步增加了最小二乘偏移成像在大規(guī)模實(shí)際資料保真成像中應(yīng)用的難度。RICKETT[4]和GUITTON[5]提出的利用偏移、反偏移和再偏移的方式估計(jì)近似Hessian逆矩陣的做法,提高了估計(jì)近似Hessian逆的計(jì)算效率。利用此種方式能估計(jì)出近似對角Hessian逆矩陣中對應(yīng)的低秩的部分[6],其代價(jià)是增加了兩次偏移的計(jì)算量。

真振幅成像是另一種較為常見的保真成像方法。真振幅成像是一種線性化的反演成像方法,其通過在高頻近似下校正波傳播過程中的幾何擴(kuò)散期望得到參數(shù)擾動(dòng)的帶限估計(jì)。根據(jù)不同模型參數(shù)化方式,真振幅成像可以依賴于ray+Born或是ray+Kirchhoff兩種近似下的線性化正問題加以實(shí)現(xiàn)[7]。其中Born和Kirchhoff近似分別對應(yīng)于體散射和面散射的情況。從線性化后的正問題出發(fā),可以從積分算子理論或是優(yōu)化理論構(gòu)造對應(yīng)反問題的解,得到所謂真振幅成像公式[8-10]。真振幅成像公式可以基于射線實(shí)現(xiàn)[11]、Gaussian-Beam實(shí)現(xiàn)[12]、單程波實(shí)現(xiàn)[13]以及雙程波實(shí)現(xiàn)[14]。在真振幅成像的原始推導(dǎo)中,由于無需考慮中低波數(shù)參數(shù)成分估計(jì)問題,隱含假設(shè)源檢兩端波場是單向傳播的,即只存在源端的下行波場和檢端的上行波場,因此,在基于雙向波傳播算子實(shí)現(xiàn)真振幅成像方法過程中,會出現(xiàn)源檢兩端波場相關(guān)能量中對應(yīng)大散射角部分的低波數(shù)噪聲。高通濾波或是角度域衰減是消除這些低波數(shù)噪聲的常規(guī)方法[14-18],然而這些方法需要在成像過程外增加額外計(jì)算步驟,另外這些方法沒有考慮低波數(shù)噪聲去除過程中對目標(biāo)反射層波形的畸變效應(yīng)。STOLK等[19]在ray+Born近似下給出了能自動(dòng)消除低波數(shù)噪聲的基于逆時(shí)偏移(RTM)實(shí)現(xiàn)的真振幅成像公式。該方法的優(yōu)點(diǎn)在于隱式考慮波場角度域?yàn)V波,無需顯式進(jìn)行波場角度域分解,在成像過程中能自動(dòng)去除低波數(shù)噪聲,另外,由于該方法基于嚴(yán)格的反問題框架導(dǎo)出,自動(dòng)考慮了角度域?yàn)V波算子非單位化的問題,輸出結(jié)果具有明確物理意義。WHITMORE等[20]進(jìn)一步從波傳播角度解釋了STOLK等[19]的研究成果中角度濾波因子的物理意義。

在常密度聲學(xué)介質(zhì)中和高頻近似意義下,如果忽略透射損失,影響成像保真度主要有兩大因素。第一個(gè)因素為波場傳播過程中的幾何擴(kuò)散,這是高頻近似下波場傳播所要滿足的物理規(guī)律之一,它會影響射線路徑及其領(lǐng)域內(nèi)振幅分布,從而影響成像結(jié)果的振幅值。第二個(gè)因素為觀測系統(tǒng)對地下成像點(diǎn)照明炮檢射線對空間分布的均勻程度。對地面觀測系統(tǒng)來說,影響成像點(diǎn)處照明射線分布的主要因素為地震觀測系統(tǒng)的排布以及地下介質(zhì)參數(shù)場的分布。傳統(tǒng)真振幅成像方法僅僅校正波傳播過程的幾何擴(kuò)散效應(yīng),對于觀測系統(tǒng),總假設(shè)為連續(xù)的無限孔徑觀測,即對于地下任意成像點(diǎn),對其有照明的炮檢射線對分布總是一致的,這顯然與實(shí)際情況不符。傳統(tǒng)真振幅成像對觀測系統(tǒng)過于理想化的假設(shè)會在較大程度上影響成像結(jié)果的保真度[21]。

本文首先回顧了Bayes反演框架下的ray+Born和ray+Kirchhoff兩種近似情況下的反演成像公式及其基于波動(dòng)方程的實(shí)現(xiàn),然后針對傳統(tǒng)基于波動(dòng)方程實(shí)現(xiàn)的高頻漸進(jìn)反演成像公式中存在的不足,給出了基于波場隱式角度分解的低波數(shù)噪聲消除方法以及高效穩(wěn)健的角度域照明補(bǔ)償方法。同時(shí),本文還簡單討論了直接輸出角度域反射系數(shù)和阻抗擾動(dòng)的方法。本文還簡單討論了直接輸出角度域反射系數(shù)和阻抗擾動(dòng)的方法。本文給出的是二維情況下的具體計(jì)算公式,但本文方法可以很方便地推廣到三維情況。

1 高頻近似下的反演成像及其波動(dòng)方程實(shí)現(xiàn)

1.1 散射場的表達(dá)

在勘探地震中,介質(zhì)參數(shù)的高波數(shù)擾動(dòng)量((方位角度)反射系數(shù)、速度擾動(dòng)和波阻抗擾動(dòng))的估計(jì),理論上都?xì)w結(jié)為最小二乘線性反演理論。更確切地,都?xì)w結(jié)為利用線性化散射場正演模型預(yù)測擾動(dòng)波場與實(shí)測擾動(dòng)波場之間的預(yù)測誤差滿足Gauss分布時(shí)的Bayes估計(jì)理論框架下的線性反演問題。

正問題是求解反問題的基礎(chǔ),地震波線性化反演成像對應(yīng)的正問題為地震波線性化散射場的計(jì)算。散射場計(jì)算問題可以由表示定理導(dǎo)出。考慮由邊界?D圍成的計(jì)算區(qū)域D中的散射場計(jì)算問題,且假設(shè)波場在計(jì)算區(qū)域二次連續(xù)可微,假設(shè)子波的頻譜F(ω)為常數(shù)。對于常密度聲學(xué)介質(zhì),在點(diǎn)源情況下表示定理形式可以寫為[22]:

(1)

式中:uS(xr,xs,ω)為炮點(diǎn)xs激發(fā),檢波點(diǎn)xr接收到的頻率域散射場;g(xr,x,ω)為地下任意點(diǎn)x到檢波點(diǎn)xr的格林函數(shù);uI(x,xs,ω)為炮點(diǎn)xs激發(fā),傳播經(jīng)過x處的入射波場;n代表邊界?D的外法線方向;c0(x)為x處背景速度;α(x)為慢度平方擾動(dòng);dS為?D上的面積微元;dV為D中的體積微元。公式(1)即為常密度聲學(xué)情況下的表示定理。從(1)式可以看出,在最廣義的情況下,散射場可以由一個(gè)面積分加一個(gè)體積分來表達(dá)。在不同的散射模式假設(shè)下,(1)式的體積分或面積分會起到不同的主導(dǎo)作用,從而可以導(dǎo)出不同的散射模式下散射波的計(jì)算公式。下面分別討論公式(1)中面積分和體積分起主導(dǎo)作用時(shí)散射波場的計(jì)算。

在體散射假設(shè)下,邊界?D不產(chǎn)生散射場,散射場全部由?D內(nèi)部的非均勻體產(chǎn)生,則(1)式退化為:

(2)

類似的,假設(shè)在?D圍成的計(jì)算區(qū)域D中,?D內(nèi)部沒有散射源,所有散射源僅來自于邊界?D上,則(1)式退化為:

(3)

為了導(dǎo)出線性化反演成像中使用的線性化的散射場計(jì)算公式,需要進(jìn)一步做線性化近似。對于體散射積分公式(2),引入Born近似[23],即假設(shè)散射場較弱,入射場與散射場的和約等于入射場有:

(4)

將公式(4)代入公式(2),有:

(5)

公式(5)描述了散射場和速度擾動(dòng)強(qiáng)度之間的線性化關(guān)系,我們將其簡稱為Born近似下散射場計(jì)算公式。對于面散射公式(3),引入Kirchhoff近似[23],即在?S1處有如下關(guān)系式成立:

(6)

式中:R為反射系數(shù)。(6)式的負(fù)號表示相對于?S1的外法線方向入射波與反射波的傳播方向相反。將公式(6)代入公式(3),有:

(7)

公式(7)即為線性化的Kirchhoff積分法散射波計(jì)算公式[21]。對比(5)式和(7)式可以看出,Born近似下散射場(體散射假設(shè))計(jì)算公式和Kirchhoff近似下散射場(面散射假設(shè))計(jì)算公式的模型空間假設(shè)是不同的。在Born近似下,體狀的散射勢產(chǎn)生散射場,數(shù)學(xué)上接近于Heaviside階躍函數(shù)。在Kirchhoff近似下,層狀分布的角度反射系數(shù)產(chǎn)生散射場,數(shù)學(xué)上接近單位脈沖函數(shù)。對于實(shí)際的帶限數(shù)據(jù)反演,Born近似下和Kirchhoff近似下反演的模型參數(shù)可以認(rèn)為是帶限的階躍函數(shù)和帶限的脈沖函數(shù)[24]。

基于線性化的散射場計(jì)算公式(5)或者公式(7),可以構(gòu)造相應(yīng)的反演成像公式。然而,由前面討論可知,直接估計(jì)完整Hessian的逆在實(shí)際情況下較為困難。Hessian矩陣是一個(gè)主對角占優(yōu)的帶狀矩陣,其特征值主要由對角元素貢獻(xiàn),利用線性Hessian矩陣的主對角元素代替整個(gè)矩陣是一種較為可行的反演策略[25]。本文在高頻近似框架下實(shí)現(xiàn)Hessian矩陣的對角化,因此,需要對線性化的體散射和面散射公式(5)以及公式(7)做進(jìn)一步的高頻近似展開。

對線性化體散射積分公式(5),在二維情況下其高頻近似的形式(2D ray+Born近似)[26]:

(8)

(9)

在(9)式中,Anl(xr,x,xs)表示從炮點(diǎn)xs出發(fā)到地下任一點(diǎn)x然后被xr處接收到射線振幅的變化。相位項(xiàng)Tnl(xr,x,xs)的表達(dá)式為:

(10)

式中:τ(xr,x,xs)為從炮點(diǎn)xs出發(fā)到地下任一點(diǎn)x然后被xr處接收到射線的總旅行時(shí)。K(xr,x,xs)為KMAH標(biāo)號,其描述射線從炮點(diǎn)xs出射,被檢波點(diǎn)xr接收到所經(jīng)歷的一階焦散點(diǎn)的總個(gè)數(shù),在每一個(gè)一階焦散點(diǎn)處,高頻近似下的波場發(fā)生π/2的相移[27]。

對線性化面散射積分公式,在二維情況下其高頻近似的形式為(2D ray+Kirchhoff近似)[7]:

(11)

(12)

這里|qnl(xr,x,xs)|為偏移傾角矢量的模[26],其定義為:

(13)

1.2 高頻近似下線性化的反演成像框架

考慮如下Bayes框架下的反問題:

① 線性問題d=Gm或可在初始模型mprior處線性展開的非線性問題g(mprior);

② 不考慮正算子的不確定性,通常也不考慮模型的先驗(yàn)信息;

據(jù)TARANTOLA[28],基于以上四點(diǎn)考慮,可以得到最小二乘意義下法方程的解形式為:

(14)

其中梯度項(xiàng)為:

(15)

(16)

注意:梯度和擬Hessian記號中的上標(biāo)代表其為對偶空間的元素[27]。進(jìn)一步地,可以構(gòu)建法方程解的迭代格式(擬Newton算法):

(17)

式中:mn,mn+1分別為當(dāng)前模型和下一次更新模型;μn為步長。本文僅考慮非迭代格式的求解。

1.3 高頻近似下的反演成像公式

將ray+Born近似下散射場的計(jì)算公式(8)代入公式(14)～公式(16),得到ray+Born近似在最小二乘意義下法方程解的具體形式為[26]:

(18)

(19)

式中:J(xs,xr,ω)→(kx,kz,θs)為笛卡爾坐標(biāo)向相空間坐標(biāo)轉(zhuǎn)化的行列式,其表達(dá)式為[29-30]:

(20)

式中:βxs,βxr分別為源檢兩端的起飛角。公式(18)～公式(20)即為ray+Born近似下的反演成像公式。類似地,基于公式(11)和公式(14)～公式(16)可以得到ray+Kirchhoff近似最小二乘意義下法方程解的具體形式為[7]:

(21)

(22)

1.4 高頻近似反演的RTM實(shí)現(xiàn)-保真的RTM偏移成像方法

本文主要討論基于反射波的保真成像,因此本小節(jié)僅討論ray+Kirchhoff真振幅成像公式的波動(dòng)方程實(shí)現(xiàn)方法,ray+Born真振幅成像公式的波動(dòng)方程實(shí)現(xiàn)方法可以利用類似的方式得出。據(jù)BLEISTEIN等[31],公式(21)可以利用如下波場相關(guān)運(yùn)算實(shí)現(xiàn):

(23)

式中:PF(x,xs,ω),PB(x,xs,ω)分別為頻率域的正傳播及反傳播波場。我們將(23)式反變換回時(shí)間域有:

(24)

這里H代表Hilbert變換。注意到ZHANG等[13]給出的輸出為R′(x,θr)=R(x,θr)/|q|,即有:

(25)

公式(25)即為ZHANG等[13]給出的二維情況下利用波動(dòng)方程實(shí)現(xiàn)的ray+Kirchhoff真振幅成像公式。注意到|q|的作用是為了校正在水平層狀介質(zhì)假設(shè)下子波峰值振幅隨角度的變化[32],對保真角度反射系數(shù)的獲取是有必要的[33]。因此在下文中,采用公式(24)作為二維ray+Kirchhoff近似下RTM真振幅成像條件。

2 基于波場隱式分解的RTM真振幅成像條件

高頻近似下成像點(diǎn)處波矢量的表達(dá)式可以表示為[34]:

(26)

其中:

(27)

(27)式中的n為單位化的偏移傾角。從公式(26)和公式(27)可以看出,當(dāng)震源端和檢波點(diǎn)端波場對應(yīng)的射線慢度矢量之間散射角較大時(shí),利用成像公式(24)進(jìn)行基于波場相關(guān)成像會產(chǎn)生低波數(shù)能量。對于估計(jì)參數(shù)場中的高波數(shù)部分,這些低波數(shù)成分是成像噪聲應(yīng)當(dāng)去除[18]。

為了消除這些低波數(shù)噪聲,STOLK等[19]在ray+Born近似意義下給出了能自動(dòng)消除低波數(shù)噪聲的基于逆時(shí)偏移(RTM)實(shí)現(xiàn)的真振幅成像公式,在其它文獻(xiàn)中也被稱之為逆散射成像條件[20]。該成像條件能夠在不顯式的進(jìn)行波場角度分解情況下完成角度域?yàn)V波,是一種高效穩(wěn)健的低波數(shù)噪聲消除方法。對于ray+Kirchhoff近似下的RTM,我們可以構(gòu)造類似能自動(dòng)消除低波數(shù)噪聲的反演成像條件。

3 高效穩(wěn)健的角度域照明補(bǔ)償

經(jīng)典的高頻近似下基于hit-count方式進(jìn)行角度域照明補(bǔ)償成像可以表示如下[33]:

(28)

式中:ψ為對應(yīng)于照明傾角的積分變量;C(x,θr,ψ)為照明圓單位面積上落入的單位偏移傾角矢量的數(shù)量;R(x,θr,ψ)為未校正前成像點(diǎn)x對應(yīng)于反射角θr的共照明傾角道集;Rc(x,θr)為成像點(diǎn)x處對應(yīng)于反射角θr校正后的成像結(jié)果。在三維情況下,有類似公式:

(29)

式中:φr為反射方位角;λ為照明方位角。關(guān)于公式(28)及公式(29)中角度的定義可以參見文獻(xiàn)[35]。由公式(28)和公式(29)可知,傳統(tǒng)基于照明傾角道集的照明校正方法物理上直觀清晰,并且由于校正過程僅需顯式計(jì)算偏移傾角矢量的空間分布密度,因此計(jì)算較為穩(wěn)健。然而,基于照明傾角道集的校正方法需要在高維空間進(jìn)行,計(jì)算量和存儲量相較于傳統(tǒng)真振幅成像方法大大增加。即使將照明校正投影到地表坐標(biāo)上完成,也僅能減少1維的數(shù)據(jù)處理量,對于3維觀測仍然需要處理6維數(shù)據(jù),同時(shí),將照明校正投影到地表坐標(biāo)上完成也會引入新的計(jì)算量[36],計(jì)算和存儲量依然較大。另外,如果考慮用波動(dòng)方程算子去實(shí)現(xiàn)傳統(tǒng)照明校正方法,在目前的計(jì)算機(jī)硬件條件下對3維實(shí)際問題幾乎不可行。基于與公式(28)和公式(29)相同的物理意義,可以在波動(dòng)方程類成像過程中高效穩(wěn)健實(shí)施角度域照明補(bǔ)償方法,這些方法將在另文單獨(dú)討論。

4 角度反射系數(shù)的輸出及阻抗成像

從公式(24)可以看出,直接利用基于ray+Kirchhoff反演類的波動(dòng)方程成像條件得到結(jié)果是角度域的反射系數(shù),因此可以直接輸出該角度域反射系數(shù),即所謂的角度道集。同時(shí),為了在不輸出角度反射系數(shù)的情況下得到有物理意義的疊加結(jié)果,也可以利用常密度聲學(xué)近似下的AVA近似關(guān)系式[37]將利用公式(24)得到的不同角度反射系數(shù)結(jié)果校正到零角度進(jìn)行疊加成像[38],最終輸出零角度反射系數(shù)成像剖面。

對于勘探地震學(xué)中的儲層參數(shù)建模,波阻抗是比反射系數(shù)更加貼近儲層描述的參數(shù)。在已知零角度反射系數(shù)的情況下,可以通過道積分來得到波阻抗分布[39]。但在零角度反射系數(shù)存在一定的誤差時(shí),直接利用道積分計(jì)算波阻抗容易產(chǎn)生數(shù)值噪聲。為了更加穩(wěn)健地利用零角度反射系數(shù)估計(jì)波阻抗,可以利用稀疏反演方法由零角度反射系數(shù)計(jì)算波阻抗[40-41]。

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

首先利用層狀模型驗(yàn)證本文提出的成像條件的有效性。模型第一層速度為3000m/s,第二層速度為4000m/s。模型的x和z方向范圍均為0～4000m,兩個(gè)方向上網(wǎng)格大小均為10m。利用真實(shí)模型正演得到的數(shù)據(jù)作為輸入炮記錄,利用半徑為10的Gussian函數(shù)平滑真實(shí)速度后的模型作為成像算法采用的背景速度。圖1a為利用不帶角度濾波的ray+Kirchhoff近似下真振幅成像條件(公式(24))得到的成像結(jié)果,從成像結(jié)果上可以看出,目標(biāo)反射層的成像為帶限的脈沖函數(shù),這與ray+Kirchhoff真振幅成像的期望輸出一致[24]。同時(shí),在目標(biāo)反射層上方存在可以預(yù)期的低波數(shù)噪聲。圖1b為利用本文提出的成像條件得到的成像結(jié)果,可以看出,低波數(shù)噪聲被消除,同時(shí)成像結(jié)果的振幅量級與圖1a保持一致。圖1c為利用本文提出的成像條件濾出的低波數(shù)噪聲。圖1d為利用帶角度濾波不帶振幅校正成像條件得到的成像結(jié)果,在這種情況下,雖然低波數(shù)噪聲被消除了,但是由于角度域?yàn)V波算子的非單位化效應(yīng)未被消除,成像結(jié)果的振幅量級相較于原始成像結(jié)果(圖1a)發(fā)生了明顯變化。圖2為單層模型不同成像條件成像結(jié)果抽道的對比。圖2a為本文提出成像條件成像結(jié)果與不帶角度濾波的ray+Kirchhoff近似下真振幅成像條件得到成像結(jié)果抽道結(jié)果對比,圖2b為本文提出成像條件濾出的低波數(shù)噪聲結(jié)果與不帶角度濾波的ray+Kirchhoff近似下真振幅成像條件得到成像結(jié)果抽道結(jié)果對比。我們可以更加清晰地看出利用不帶角度濾波的ray+Kirchhoff近似下真振幅成像條件得到成像結(jié)果淺部存在的低波數(shù)噪聲,利用本文提出的成像條件得到的反射層成像結(jié)果以及濾出的低波數(shù)噪聲在振幅和相位上與不帶角度濾波成像結(jié)果中對應(yīng)成分具有很好的對應(yīng)關(guān)系。這說明了本文提出的成像條件在自動(dòng)消除低波數(shù)成像噪聲的同時(shí)較好保持了目標(biāo)反射層的振幅和相位信息。

圖1 單層模型不同成像條件成像結(jié)果的對比a 利用不帶角度濾波的ray+Kirchhoff近似下真振幅成像條件((24)式)成像結(jié)果; b 利用本文提出的成像條件成像結(jié)果; c 圖1a被濾出的低波數(shù)噪聲; d 利用帶角度濾波不帶振幅校正的成像條件的成像結(jié)果

對于角度反射系數(shù)輸出,我們利用Marmousi模型進(jìn)行數(shù)值實(shí)驗(yàn),并對輸出角度反射系數(shù)結(jié)果進(jìn)行對比。圖3和圖4展示了數(shù)值實(shí)驗(yàn)結(jié)果。從中可以看出利用本文方法輸出角度反射系數(shù)與理論值具有較好的匹配關(guān)系。Marmousi模型數(shù)值實(shí)驗(yàn)結(jié)果表明了本文提出成像條件輸出角度反射系數(shù)對復(fù)雜模型的適應(yīng)性。對于角度域照明補(bǔ)償及阻抗成像,首先利用層狀模型進(jìn)行測試。層狀模型參數(shù)見圖5和圖6。圖7展示了在真振幅RTM基礎(chǔ)上不進(jìn)行和進(jìn)行角度域照明補(bǔ)償?shù)膸薹瓷湎禂?shù)估計(jì)結(jié)果。可以看出,在正確校正幾何擴(kuò)散的基礎(chǔ)上進(jìn)一步校正成像點(diǎn)處照明補(bǔ)償不均勻效應(yīng),能夠大大提高帶限反射系數(shù)的估計(jì)精度。在得到的帶限反射基礎(chǔ)上,進(jìn)一步得到了阻抗擾動(dòng)的估計(jì),結(jié)果見圖8和圖9。很明顯,在考慮了照明補(bǔ)償后,阻抗擾動(dòng)估計(jì)結(jié)果的精度也有較大程度的提高。為了驗(yàn)證本文討論的角度域照明補(bǔ)償及阻抗成像在復(fù)雜模型上的測試結(jié)果,我們進(jìn)一步利用部分Sigsbee模型進(jìn)行測試。圖10為部分Sigsbee模型參數(shù)。帶限反射系數(shù)以及阻抗擾動(dòng)成像結(jié)果分別見圖11、圖12以及圖13。由圖可見,與前面層狀模型測試結(jié)果類似,在真振幅成像基礎(chǔ)上進(jìn)一步考慮角度域照明補(bǔ)償能夠大大提高反射系數(shù)以及阻抗擾動(dòng)反演成像精度。

圖2 單層模型不同成像條件成像結(jié)果抽道的對比a 利用不帶角度濾波的真振幅成像條件((24)式)成像結(jié)果(藍(lán)色實(shí)線)與利用本文提出的成像條件成像結(jié)果(紅色實(shí)線)對比; b 利用不帶角度濾波的真振幅成像條件((24)式)成像結(jié)果與利用本文提出成像條件濾出的低波數(shù)噪聲(紅色實(shí)線)對比

圖3 Marmousi模型所有CDP的理論角度道集(a)和計(jì)算得到的角度道集(b)

圖4 Marmousi模型角度反射系數(shù)輸出對比a CDP119～CDP121理論角度道集; b CDP119～CDP121計(jì)算得到的角度道集; c CDP120、深度824m處理論AVA曲線(藍(lán)色實(shí)線)與計(jì)算AVA曲線(紅色曲線)對比; d CDP120、深度1572m處理論AVA曲線與計(jì)算AVA曲線對比

圖5 層狀模型模型參數(shù)場a 真實(shí)速度模型; b 真實(shí)速度擾動(dòng)

圖6 成像采用的平滑模型(a)及真實(shí)零角度反射系數(shù)(b)

圖7 層狀模型帶限反射系數(shù)成像對比a 不考慮照明校正的零角度反射系數(shù)成像結(jié)果; b 利用本文提出的照明補(bǔ)償算子進(jìn)行照明校正的零角度反射系數(shù)成像結(jié)果; c,d分別為圖7a、圖7b抽取中間一道零角度反射系數(shù)和理論零角度系數(shù)的對比

圖8 層狀模型利用零角度反射系數(shù)反演速度擾動(dòng)結(jié)果(Ⅰ)a 基于不考慮照明校正的零角度反射系數(shù)估計(jì)結(jié)果反演得到的速度擾動(dòng)場; b 抽取中間道和理論速度擾動(dòng)的對比

圖9 層狀模型利用零角度反射系數(shù)反演速度擾動(dòng)結(jié)果(Ⅱ)a 基于本文提出照明補(bǔ)償算子進(jìn)行照明校正的零角度反射系數(shù)估計(jì)結(jié)果反演得到的速度擾動(dòng)場; b 抽取中間道和理論速度擾動(dòng)的對比

圖10 部分Sigsbee模型參數(shù)場a 真實(shí)速度模型; b 成像采用的平滑模型; c 真實(shí)速度擾動(dòng); d 真實(shí)零角度反射系數(shù)

圖11 部分Sigsbee模型帶限反射系數(shù)成像對比(Ⅰ)a 不考慮照明校正的零角度反射系數(shù)成像結(jié)果; b 抽取中間一道的零角度反射系數(shù)和理論零角度反射系數(shù)的對比結(jié)果

圖12 部分Sigsbee模型帶限反射系數(shù)成像對比(Ⅱ)a 利用本文提出的照明補(bǔ)償算子進(jìn)行照明校正的零角度反射系數(shù)成像結(jié)果; b 抽取中間一道的零角度反射和理論零角度反射系數(shù)的對比結(jié)果

6 結(jié)論與討論

本文系統(tǒng)地討論了體散射以及面散射模式下的散射場計(jì)算。在高頻近似框架下,本文進(jìn)一步討論了不同散射模式下的反演成像公式及其波動(dòng)方程實(shí)現(xiàn)。在反射模式下,我們給出了改進(jìn)的真振幅RTM成像條件,用以在成像過程中自動(dòng)且保真地消除低波數(shù)噪聲。在此基礎(chǔ)上,我們還討論了角度反射系數(shù)的輸出、角度域的照明補(bǔ)償以及阻抗成像方法。模型數(shù)據(jù)的數(shù)值試驗(yàn)結(jié)果證明了本文方法相較于傳統(tǒng)基于反射模式真振幅成像條件的優(yōu)勢。本文結(jié)果可以推廣到三維情況,相應(yīng)的分析測試工作是本文下一步的工作重點(diǎn)之一。

致謝:感謝中石油勘探開發(fā)研究院及西北分院、中海油研究院和湛江分公司、中國石油化工股份有限公司石油物探技術(shù)研究院和勝利油田分公司對波現(xiàn)象與智能反演成像研究組(WPI)研究工作的資助與支持。