王華忠,胡江濤,2,郭 頌

(1.波現(xiàn)象與反演成像研究組(WPI),同濟(jì)大學(xué)海洋與地球科學(xué)學(xué)院,上海200092;2.油氣藏地質(zhì)及開發(fā)工程國(guó)家重點(diǎn)實(shí)驗(yàn)室,成都理工大學(xué),四川成都610059)

?

最小二乘疊前深度偏移成像理論與方法

王華忠1,胡江濤1,2,郭 頌1

(1.波現(xiàn)象與反演成像研究組(WPI),同濟(jì)大學(xué)海洋與地球科學(xué)學(xué)院,上海200092;2.油氣藏地質(zhì)及開發(fā)工程國(guó)家重點(diǎn)實(shí)驗(yàn)室,成都理工大學(xué),四川成都610059)

地震波成像的目的首先是定位反射界面或散射點(diǎn)的空間位置,描述地下地質(zhì)體的幾何結(jié)構(gòu);然后是估計(jì)地下介質(zhì)的物性參數(shù),主要是與速度和密度相關(guān)的參數(shù);最終是與巖石物理結(jié)合描述含油氣儲(chǔ)層。以地震波傳播的散射場(chǎng)表達(dá)以及逆散射成像為切入點(diǎn),首先給出散射波的表達(dá)形式,討論了逆散射成像的本質(zhì)與逆散射成像所需要的疊前地震觀測(cè)數(shù)據(jù)及對(duì)應(yīng)的觀測(cè)方式;然后給出一般的線性化最小二乘疊前深度偏移(Least-squares Prestack Depth Migration,LS_PSDM)基本理論框架,只要能給出用Green函數(shù)計(jì)算的波傳播算子,就可以實(shí)現(xiàn)線性化的LS_PSDM;接著討論了當(dāng)前廣受關(guān)注的最小二乘逆時(shí)深度偏移(Least-squares Reverse Time Migration,LS_RTM)方法,給出了最小二乘偏移成像的迭代實(shí)現(xiàn)方法和保反射界面結(jié)構(gòu)的總變差正則化方法以及數(shù)值結(jié)果;最后分析了基于反射波估計(jì)反射系數(shù)的成像方法和基于散射波估計(jì)散射強(qiáng)度的成像方法的差異,指出針對(duì)勘探地震介質(zhì)特征和波場(chǎng)特征的成像方法是最具實(shí)用性的方法,逆散射成像應(yīng)該在此基礎(chǔ)上進(jìn)行?；谏⑸洳ǖ某上穹椒芊駶M足儲(chǔ)層描述的要求,有待進(jìn)一步研究工作的檢驗(yàn)。

散射場(chǎng)表達(dá);逆散射成像;反射波成像;最小二乘疊前深度偏移;全波形反演

油氣勘探地震學(xué)的最終目標(biāo)是儲(chǔ)層描述。實(shí)現(xiàn)儲(chǔ)層描述與刻畫的途徑有兩條，一是常規(guī)的、基于反射波場(chǎng)的保幅或保真角度反射系數(shù)的成像。該方法基于保真的角度成像道集,利用1D波阻抗反演和/或AVA疊前反演得到彈性參數(shù),與巖石物理、測(cè)井結(jié)合,進(jìn)行儲(chǔ)層識(shí)別和描述,是目前常用的從地震成像走向儲(chǔ)層刻畫的方法技術(shù)體系[1]。該技術(shù)體系的核心包括三方面的關(guān)鍵因素:①滿足速度估計(jì)和保真成像要求的數(shù)據(jù)體;②較準(zhǔn)確的背景速度模型建立技術(shù)和保真角度反射系數(shù)成像技術(shù);③基于保真成像道集的儲(chǔ)層參數(shù)估計(jì)及儲(chǔ)層描述方法。實(shí)現(xiàn)儲(chǔ)層描述與刻畫的另一條途徑是當(dāng)前石油工業(yè)界十分關(guān)注的全波形反演(Full Waveform Inversion,FWI)成像——基于散射和逆散射的反演成像技術(shù),它的本質(zhì)目標(biāo)是估計(jì)寬(全)波數(shù)帶的彈性參數(shù)場(chǎng),基于此反演結(jié)果直接進(jìn)行儲(chǔ)層描述。但該方法僅僅利用地表觀測(cè)的有限帶寬、有限孔徑和無假頻數(shù)據(jù),其實(shí)用化并不成功,僅僅在低頻和長(zhǎng)偏移距地表觀測(cè)數(shù)據(jù)情況下能得到背景速度的估計(jì),彈性參數(shù)躍變(高波數(shù))成分的估計(jì)仍然采用最小二乘RTM。然而,依賴散(反)射波振幅的LS_RTM也很難有效地收斂并給出滿足儲(chǔ)層描述要求的反演成像結(jié)果。

在一些數(shù)學(xué)家看來,石油地震勘探問題很簡(jiǎn)單,就是用波動(dòng)方程作為實(shí)測(cè)數(shù)據(jù)預(yù)測(cè)器,若基于一個(gè)假設(shè)的參數(shù)模型求解波動(dòng)方程得到的合成數(shù)據(jù)與實(shí)測(cè)疊前地震數(shù)據(jù)的總誤差最小,則該假設(shè)模型就是所求的地下巖石彈性參數(shù)估計(jì)結(jié)果。事實(shí)上,這完全低估了勘探地震中彈性參數(shù)反演問題的復(fù)雜性[2-3]。從勘探地震學(xué)出現(xiàn)伊始,勘探地震學(xué)家走的就是一條實(shí)際問題驅(qū)動(dòng)的路線:針對(duì)成像問題,先由CMP道集估計(jì)背景初始速度,進(jìn)行疊后偏移成像;然后進(jìn)入疊前偏移成像,由成像道集進(jìn)行偏移速度分析,或進(jìn)行層析成像建立速度模型。迭代進(jìn)行疊前偏移成像和速度模型的更新是石油工業(yè)界發(fā)展出的實(shí)用化地震波成像路線,經(jīng)實(shí)踐檢驗(yàn)是行之有效的。更細(xì)致的能用于儲(chǔ)層描述的參數(shù)估計(jì)是基于1D波阻抗反演和/或AVA角度道集反演[4]。事實(shí)上,一直有很多學(xué)者在探討是否能從更數(shù)學(xué)化的角度來解決成像問題,并且分成兩條路線,一條是以BLEISTEIN為代表的基于反射波的線性反演路線[5-6];另一條是以TARANTOLA為代表的基于散射波的全波形擬合非線性反演路線[7]。前者是當(dāng)前勘探地震反射波成像的完整理論框架,一直在指導(dǎo)反射波成像的具體實(shí)踐,取得了很好的應(yīng)用效果；后者是最近幾年勘探地震中的熱點(diǎn)研究問題,利用低頻透射波估計(jì)背景速度在合適的探區(qū)獲得了比較顯著的應(yīng)用效果,但以LS_RTM為代表的彈性參數(shù)高波數(shù)擾動(dòng)量的估計(jì)還沒有看到具體的實(shí)用效果。

基于此,本文以地震波傳播的散射場(chǎng)表達(dá)與逆散射成像的本質(zhì)為切入點(diǎn),討論了一般意義下疊前偏移成像與逆散射成像之間的區(qū)別與聯(lián)系,重點(diǎn)分析了最小二乘疊前深度偏移成像方法存在的問題,認(rèn)為勘探地震中的地震波成像要充分考慮地下介質(zhì)以分層沉積為主的特點(diǎn)以及波在這樣的介質(zhì)中傳播以反射波為主的特點(diǎn)構(gòu)建合適的成像方法技術(shù),而不是完全基于地下介質(zhì)是背景+散射的模式構(gòu)建逆散射成像技術(shù)。在方位角度反射系數(shù)估計(jì)結(jié)果比較準(zhǔn)確的基礎(chǔ)上再進(jìn)行逆散射成像應(yīng)該是LS_RTM方法走向?qū)嵱没暮侠砺肪€。

1 地震波傳播的散射場(chǎng)表達(dá)

FWI和波動(dòng)理論地震波層析成像本質(zhì)上是針對(duì)地下介質(zhì)為背景+散射體情形的,因此,散射體引起的散射場(chǎng)的表達(dá)是其中最基本的問題。

1.1 基本的聲波方程和聲散射勢(shì)的引入

1.1.1 基本的聲波方程

一般地,我們假設(shè)波在地下介質(zhì)中的傳播可以由如下聲波方程控制:

(1)

式中:u(x,t)為聲場(chǎng)勢(shì),它是時(shí)間t∈R和空間坐標(biāo)x∈Ω?R的函數(shù);v(x)代表速度場(chǎng);ρ(x)代表密度場(chǎng);S(x,t)代表體積場(chǎng)源。

在頻率-空間域,上式變?yōu)?

(2a)

其中,

(2b)

為聲波傳播算子。

定義背景介質(zhì)中的聲波傳播滿足如下方程:

(3a)

其中,

(3b)

為背景介質(zhì)中的聲波傳播算子。

1.1.2 聲散射勢(shì)的引入

根據(jù)(2)式和(3)式,定義聲散射勢(shì)為:

Va(x,ω)=La(x,ω)-La0(x,ω)

(4)

其中,

(5a)

(5b)

1.2 標(biāo)量波方程和標(biāo)量波散射勢(shì)的引入

勘探地震學(xué)很多情況下關(guān)注標(biāo)量波方程(或密度不變情況下的聲波方程)。頻率-空間域標(biāo)量波方程形式為:

(6a)

其中,

(6b)

為標(biāo)量波傳播算子。

定義背景介質(zhì)中的聲波傳播滿足如下方程:

(7a)

其中,

(7b)

為背景介質(zhì)中的標(biāo)量波傳播算子。

根據(jù)(6)式和(7)式,定義標(biāo)量波散射勢(shì)為:

(8a)

其中,

(8b)

稱為散射強(qiáng)度。(8a)式定義的散射勢(shì)在速度擾動(dòng)不大時(shí)還可以寫為如下形式:

(9)

上述形式引入δv(x)=v(x)-v0(x),更清楚地反映出速度的局部變化是散射勢(shì)的本質(zhì)。從散射勢(shì)的定義也可以看出,估計(jì)散射強(qiáng)度的繞(散)射波動(dòng)理論層析成像與估計(jì)(或定位)反射系數(shù)的疊前偏移成像是有嚴(yán)格區(qū)別的。根據(jù)(9)式定義的散射勢(shì),可以寫出散射勢(shì)引起的場(chǎng)的傳播控制方程:

(10)

求解該方程可以得到散射勢(shì)引起的場(chǎng)us(x,ω),進(jìn)而進(jìn)行波動(dòng)理論的層析成像。這就是在一定的反演理論下的線性化散射波FWI(LS_RTM)的正問題基礎(chǔ)。

1.3Lippmann-Schwinger方程的導(dǎo)出及散射場(chǎng)的具體表達(dá)

L=L0+L0VL

(11)

假如波傳播算子L用Green函數(shù)G表示,(11)式可以重寫為:

G=G0+G0VG

(12)

(12)式是非線性的,它能被進(jìn)一步重寫為:

(13)

對(duì)(13)式右端項(xiàng)進(jìn)行如下的級(jí)數(shù)展開,得:

G0VG0VG0+G0VG0VG0VG0+…

(14)

(12)式和(14)式被稱為L(zhǎng)ippmann-Schwinger方程[8-9]。

很顯然,假如j≤2,(14)式就僅描述一個(gè)忽略了二階以上散射波的波傳播過程。即線性化的波傳播算子僅僅描述波傳播過程中的一階散射，也就是:

G≈G0+G0VG0

(15)

(15)式就是Born近似下的波傳播算子。其中方程右端的第一項(xiàng)描述了背景場(chǎng)中的波傳播現(xiàn)象,第二項(xiàng)描述了一階散射波的傳播現(xiàn)象。(12)至(15)式清楚地說明了Born近似的物理含義,只有當(dāng)散射勢(shì)比較弱時(shí),(14)式中高于二階的項(xiàng)才可以被忽略不計(jì)。從(14)式可以清楚地看出,高階散射場(chǎng)能否被忽略,除了散射勢(shì)比較弱的條件外,散射體的體積也很重要,因?yàn)樯⑸鋱?chǎng)是散射勢(shì)的體積積分結(jié)果。小散射體和弱散射勢(shì)是Born近似成立的條件。

從(15)式可以看出,Born近似線性化算子模擬出的總波場(chǎng)包括兩部分:第一部分是由背景Green函數(shù)所描述的背景波場(chǎng);第二部分是由散射勢(shì)引起的、方程(14)中的第二項(xiàng)所描述的散射場(chǎng)。很顯然,這里存在一個(gè)背景速度場(chǎng)如何獲取才能得到符合Born近似要求的背景場(chǎng)的問題，原則上的要求是入射場(chǎng)中不能存在散射場(chǎng)成分。事實(shí)上,在特別復(fù)雜、變化特別劇烈的速度場(chǎng)中,一階Born近似的應(yīng)用效果是很值得考慮的。這也是人們逐步在考慮引入更高階的Born近似、降低一階Born近似弱散射要求的原因。但這樣做也會(huì)引入高階散射場(chǎng)計(jì)算上的復(fù)雜性。

假定Born近似假設(shè)成立,根據(jù)(15)式,線性化的總波場(chǎng)可以寫成

(16)

其中,Born近似后的散射波場(chǎng)記為:

(17)

Born近似的物理實(shí)質(zhì)是:在背景場(chǎng)u0(ξ,η,ω)中沒有散射場(chǎng);在擾動(dòng)場(chǎng)中,只存在一次散射場(chǎng)。盡管用這種方式描述的波傳播與實(shí)際的波傳播不同,但當(dāng)二階及二階以上的散射場(chǎng)能量比較弱時(shí),用Born近似散射場(chǎng)進(jìn)行最小二乘疊前深度偏移、估計(jì)散射強(qiáng)度還是比較有效的。這是當(dāng)前最小二乘疊前深度偏移成像的基礎(chǔ),也預(yù)示著未來的最小二乘疊前深度偏移成像需要考慮更高階次的散射波才能得到更好的結(jié)果[10]。

關(guān)于子波問題,可以認(rèn)為已經(jīng)包含在背景場(chǎng)u0(x,η,ω)中,也可以用G0(x,η,ω)S(ω)來表示,其中S(ω)代表子波的譜。但是,在最小二乘疊前深度偏移成像中,假設(shè)子波已知,僅需估計(jì)散射強(qiáng)度。事實(shí)上,子波的振幅譜和相位譜不正確都會(huì)影響散射強(qiáng)度估計(jì)的精度。因此,子波研究是決定最小二乘疊前深度偏移反演成像結(jié)果好壞的重要因素。

2 逆散射問題分析

2.1 常速背景介質(zhì)與平面波入射情形下散射場(chǎng)與散射勢(shì)間的關(guān)系

首先給出常速介質(zhì)中自由空間Green函數(shù):

(18)

式中:x代表場(chǎng)點(diǎn)位置，x′代表散射源點(diǎn)位置，|x-x′|代表源點(diǎn)與場(chǎng)點(diǎn)之間的距離;k代表場(chǎng)的波數(shù)。將(18)式代入(17)式定義的散射場(chǎng)中,得:

(19)

一般地,在遠(yuǎn)場(chǎng)近似假設(shè)下,入射場(chǎng)在常速介質(zhì)中傳播可以用平面波近似:

(20)

式中:s0是平面波入射方向的單位矢量。將(20)式代入(19)式得:

(21)

(21)式描述的波的傳播關(guān)系如圖1所示。由(16)式知,此時(shí)總場(chǎng)可表示為:

(22)

圖1 空間局部散射體引起的散射波傳播關(guān)系

由圖1可知,當(dāng)場(chǎng)點(diǎn)P距離散射體比較遠(yuǎn),散射體積源分布空間V有限時(shí),用Q代表散射體積源內(nèi)的點(diǎn),O代表坐標(biāo)原點(diǎn),則可以做如下近似推導(dǎo):

(23)

(24)

式中:s是散射波傳播的局部平面波方向的單位矢量。將(24)式代入(22)式,得:

(25)

其中,

(26)

稱為散射波方向譜。定義散射勢(shì)的Fourier變換為:

(27)

據(jù)此,可以總結(jié)出如下定理:在Born近似下,相對(duì)于散射體的遠(yuǎn)場(chǎng)區(qū)域,散射勢(shì)的Fourier變換等于散射波方向譜,且有關(guān)系式

(28)

成立[11]。

由此可知,在常速介質(zhì)情況下,所謂逆散射問題就是根據(jù)散射波方向譜確定散射勢(shì)的Fourier分量。散射勢(shì)的所有Fourier分量都確定后,通過Fourier逆變換就可以完全精確地恢復(fù)散射勢(shì)。

上述定理盡管是在平面波入射和常速背景情況下得到的,但其理論指導(dǎo)意義十分明顯。從k′的表達(dá)式中可以看出,確定k′的低波數(shù)部分需要大角度散射波和低頻數(shù)據(jù);確定k′的高波數(shù)部分需要小角度散射波和高頻數(shù)據(jù)?？傊?震源子波的頻帶范圍、觀測(cè)系統(tǒng)和背景速度場(chǎng)的分布決定了散射強(qiáng)度估計(jì)結(jié)果的分辨率。因此,寬帶、寬方位觀測(cè)數(shù)據(jù)是高精度反演成像所必需的。盡管(28)式中沒有明確體現(xiàn),無假頻的高密度觀測(cè)也是必要的。

2.2 高頻近似下的Born近似的物理含義

Born近似下的散射場(chǎng)(17)式也可以表示為:

(29)

用漸進(jìn)Green函數(shù)(即Green函數(shù)的WKBJ近似解)代入上式,得:

us(ξ,η,ω)=ω2∫dxf(x)A(x,η)A(ξ,x)·

exp{-iω[T(x,η)+T(ξ,x)]}

(30)

在散射點(diǎn)周圍,對(duì)旅行時(shí)進(jìn)行線性化,即圍繞中心射線進(jìn)行Taylor展開:

(31)

式中:p(x0,x′)表示旅行時(shí)場(chǎng)T(x,x′)在點(diǎn)x0處的一階導(dǎo)數(shù),其方向?yàn)橹行纳渚€的切線方向。因此,(30)式中的指數(shù)部分可以化為:

(32)

其中,

(33)

代表波場(chǎng)的一種波數(shù)關(guān)系[12]。其幾何意義如圖2所示。

圖2 地震波傳播過程中波矢量在反射點(diǎn)處的幾何關(guān)系

圖2中p=ω[p(x0,η)+p(x0,ξ)]=2ncosθ/v0(x0)。其中n為界面的外法線方向單位向量;θ為射線的入射角;v0(x0)為x0處的背景速度。根據(jù)圖2 所示幾何關(guān)系有:

(34)

注意(34)式與(28)式的一致性,但這是從反射的角度來解釋。從k的表達(dá)式中可以看出,確定k的低波數(shù)部分需要大角度散射波和低頻數(shù)據(jù),確定k的高波數(shù)部分要小角度散射波和高頻數(shù)據(jù)。同樣地,此處希望滿足k規(guī)定的譜的完整性。這進(jìn)一步證明了高精度的逆散射反演成像需要寬頻帶和寬方位的地震數(shù)據(jù)觀測(cè)。

3 Bayes框架下基于Born近似散射場(chǎng)的散射勢(shì)估計(jì)方法

(17)式定義了一個(gè)Fredholm第一類積分方程,它描述的是一個(gè)線性問題。此處,我們將該Fredholm第一類積分方程的求解轉(zhuǎn)化為一個(gè)Bayes估計(jì)問題,就是將線性問題的求解轉(zhuǎn)化為利用泛函變分求極值的問題。基本內(nèi)容表述如下。

假設(shè)(17)式定義的積分方程寫成算子表達(dá)形式:

Km=d

(35)

式中:K是Fredholm第一類積分方程的積分核。由(17)式知,積分核的定義為:

K[v0(x),x,ω]=ω2G0(ξ,x,ω)u0(x,η,ω)

(36)

式中:m代表由f(x)作為元素形成的矢量;d代表每個(gè)炮檢點(diǎn)對(duì)應(yīng)的地震道的頻譜值。K:H1→H2被認(rèn)為是一個(gè)有界的線性算子,定義如下Tikhonov泛函:

(37)

對(duì)(37)式求解的首要問題是計(jì)算泛函關(guān)于參數(shù)變化的導(dǎo)數(shù),即泛函的梯度。為此,利用如下的變分方法。

對(duì)于任意的介質(zhì)擾動(dòng)δm∈H1和擾動(dòng)步長(zhǎng)τ∈R,總存在如下的泛函變分:

(38)

對(duì)(38)式定義的變分求偏導(dǎo)數(shù),即:

=[(〈Km-d,Kδm〉H2+α〈m,δm〉H1)+

=(〈Km-d,Kδm〉H2+α〈m,δm〉H1)

=[〈K*(Km-d),δm〉H1+α〈m,δm〉H1]

=〈K*(Km-d)+αm,δm〉H1

=〈(K*K+αI)m-K*d,δm〉H1

=〈GradI(m),δm〉H1

(39)

其中,

(40)

是誤差泛函相對(duì)于參數(shù)模型的梯度。K*是波場(chǎng)正傳播算子的共軛算子,代表波場(chǎng)的反傳播。經(jīng)典地,我們利用該梯度構(gòu)造最速下降迭代方法求解Born近似后的線性問題的解。迭代格式為:

(41)

其中,μk為第k個(gè)迭代步長(zhǎng)。當(dāng)然,更一般性的解法是令(40)式定義的梯度等于0,得到如下線性系統(tǒng)的法方程:

(42)

法方程的求解一般用共軛梯度方法。但是,由于勘探地震中模型參數(shù)描述采用網(wǎng)格化方式,而模型參數(shù)個(gè)數(shù)非常多,(42)式右端的Hessian矩陣異常龐大。若用NX×NY×NZ表示要估計(jì)的模型參數(shù)的個(gè)數(shù),那么Hessian矩陣的規(guī)模為(NX×NY×NZ)×(NX×NY×NZ),當(dāng)前的計(jì)算機(jī)技術(shù)還不能滿足如此巨大規(guī)模的線性方程組的求解。再者,形成Hessian矩陣本身也是一個(gè)巨大的挑戰(zhàn)。所以,當(dāng)前解法方程、進(jìn)行最小二乘偏移成像的基本做法是先進(jìn)行常規(guī)的疊前深度偏移成像,然后用對(duì)角矩陣元素(地震波照明能量)對(duì)成像疊加后的結(jié)果進(jìn)行照明補(bǔ)償,就是用Hessian矩陣的對(duì)角元素去除成像結(jié)果。也可以用角度Hessian矩陣的對(duì)角元素去除每個(gè)共角度成像剖面,然后進(jìn)行疊加,生成最終剖面。當(dāng)然,也可以不疊加,產(chǎn)生照明補(bǔ)償后的角度成像道集。(42)式代表的問題可以認(rèn)為是圖像反褶積問題。右端項(xiàng)是常規(guī)偏移成像得到的結(jié)果,它是模糊的圖像。用共軛梯度法解(42)式相當(dāng)于圖像去模糊運(yùn)算,可以顯著地提高最終圖像的質(zhì)量。求解(42)式最關(guān)鍵的問題是難以形成Hessian矩陣,根據(jù)Hessian矩陣元素的物理含義(相鄰空間點(diǎn)Green函數(shù)的互相關(guān)),學(xué)者們提出了很多簡(jiǎn)化的Hessian矩陣生成方法。

實(shí)質(zhì)上,(40)式計(jì)算的梯度不會(huì)等于0,因?yàn)檎輸?shù)據(jù)和實(shí)測(cè)數(shù)據(jù)的殘差永遠(yuǎn)不會(huì)等于0。這是最小二乘偏移成像中迭代解法占主要地位的根本原因。迭代法涉及殘差計(jì)算rk+1=Kmk-dobs,這需要一步反偏移;(41)式中的K*rk+1=K*(Kmk-dobs)代表一迭代步中殘差的反向傳播,需要一步殘差的反向傳播。同時(shí),迭代法要考慮步長(zhǎng)的計(jì)算并用Hessian矩陣的逆對(duì)梯度項(xiàng)進(jìn)行預(yù)條件處理。以上只是一步迭代需要的計(jì)算量,最小二乘偏移成像一般需要n步迭代計(jì)算,因此,迭代法的計(jì)算量也很大。

迭代的線性最小二乘反演方法中并不更新核函數(shù),因?yàn)槲覀兗俣ê撕瘮?shù)K[v0(x),x,ω]僅僅與背景速度v0(x)有關(guān)。但這僅僅是從計(jì)算效率方面考慮的結(jié)果,本質(zhì)上還是應(yīng)該更新核函數(shù),因?yàn)楸尘八俣扰c擾動(dòng)速度是相互耦合的。存在較強(qiáng)多次散射波(反射波)時(shí),要么事先消除多次波,要么構(gòu)建非線性的最小二乘反演成像方法。

(41)式或(42)式中計(jì)算出的m(x)代表速度的一種相對(duì)變化,只有在速度躍變的界面或繞射點(diǎn)處才會(huì)出現(xiàn),因此,它反映了速度場(chǎng)中的高頻變化部分。這部分高波數(shù)速度擾動(dòng)可以用來代表界面位置和繞射位置的圖像,與一般的疊前深度偏移類似。也可以精確估計(jì)這部分高波數(shù)速度擾動(dòng),用來刻畫儲(chǔ)層的性質(zhì)。但是,我們認(rèn)為,只有提供滿足(28)式或(34)式要求的疊前數(shù)據(jù)體時(shí),反演成像的結(jié)果才比較可靠,才可以用來進(jìn)行有效的儲(chǔ)層解釋。目前,僅僅利用地表觀測(cè)的反射(散射)波場(chǎng),加上噪聲的存在,m(x)的估計(jì)結(jié)果很難滿足精確描述儲(chǔ)層的要求。線性的最小二乘反演成像可以認(rèn)為是背景速度正確時(shí)的反(散)射波FWI。因此,目前的FWI要得到適用于儲(chǔ)層解釋的反演結(jié)果也是不現(xiàn)實(shí)的。速度(更應(yīng)該是波阻抗)躍變的精確估計(jì)是一個(gè)比背景速度估計(jì)更難的問題。

需要指出的是,本節(jié)討論的內(nèi)容適用于利用各種波傳播算子(包括Kirchhoff積分算子、Gauss-Beam算子、單向波傳播算子以及雙向波傳播算子)計(jì)算Green函數(shù),進(jìn)行線性化的最小二乘疊前深度偏移成像。

4 線性與非線性最小二乘逆時(shí)深度偏移(LS_RTM)

FWI的本質(zhì)目的是估計(jì)寬或全波數(shù)帶的速度(最好是波阻抗)。但是當(dāng)前勘探地震數(shù)據(jù)采集并沒有按(28)式或(34)式的要求進(jìn)行,而是僅僅在地表進(jìn)行有限帶寬和有限孔徑的觀測(cè)。這是當(dāng)前FWI分成估計(jì)低波數(shù)背景速度和高波數(shù)速度擾動(dòng)兩部分的主要原因。FWI背景速度估計(jì)主要是利用低頻長(zhǎng)偏移距的透射波;FWI速度擾動(dòng)的估計(jì)是在假設(shè)背景速度正確的基礎(chǔ)上利用反(散)射波進(jìn)行的。速度(或波阻抗)擾動(dòng)的估計(jì)必須利用反(散)射波的振幅信息,基于目前的數(shù)據(jù)采集,它們的估計(jì)精度更難以保障?；谡穹母卟〝?shù)彈性參數(shù)的估計(jì)是更困難的問題。目前進(jìn)行速度(或波阻抗)擾動(dòng)估計(jì)的主要技術(shù)是線性或非線性的最小二乘逆時(shí)深度偏移方法。

4.1 非線性的LS_RTM

我們首先列出非線性的最小二乘逆時(shí)深度偏移成像(LS_RTM)公式。正向波場(chǎng)外推方程為:

(43a)

式中:fb(Ω,t)代表邊界條件,一般用PML邊界條件,Ω代表計(jì)算波場(chǎng)在空間上的邊界。

反向波場(chǎng)外推方程為:

(43b)

式中:η代表觀測(cè)點(diǎn)坐標(biāo)值,uobs(η,t)代表觀測(cè)點(diǎn)處的實(shí)測(cè)波場(chǎng),ucal(η,t)代表觀測(cè)點(diǎn)處的正演合成波場(chǎng)。FWI進(jìn)行波場(chǎng)殘差外推時(shí),僅有觀測(cè)點(diǎn)處有波場(chǎng)殘差,因此有(43b)這樣的邊界條件的定義。實(shí)際上,波場(chǎng)逆時(shí)外推并不限于這樣的邊界條件。

由此可見,(43a)式是一個(gè)初值問題,沿時(shí)間正向外推計(jì)算;(43b)式是一個(gè)邊值問題,逆時(shí)間反向外推計(jì)算。

誤差泛函梯度公式為:

(44)

式中:ξ代表不同的炮點(diǎn)位置;η代表不同的檢波點(diǎn)位置。假定初始背景速度mB是已知的，δm代表要估計(jì)的高波數(shù)速度擾動(dòng)量;δmB代表要估計(jì)的背景速度的擾動(dòng)量。對(duì)應(yīng)的迭代公式為:

(45a)

(45b)

與前述基于一階Born近似公式導(dǎo)出的LS_PSDM相比，非線性LS_RTM(或散/反射波FWI)的梯度項(xiàng)中包含了對(duì)所有波現(xiàn)象的預(yù)測(cè)。當(dāng)背景速度不準(zhǔn)時(shí),不能預(yù)測(cè)的波現(xiàn)象會(huì)干擾梯度的計(jì)算,導(dǎo)致收斂慢,甚至不收斂的情況發(fā)生。比較而言,一階Born近似僅僅考慮一次散射波的預(yù)測(cè),梯度計(jì)算相對(duì)要更準(zhǔn)確,因此LS_PSDM的收斂性更好。但LS_PSDM要求ucal(η,t)中僅僅包含一次散射(反射)波場(chǎng),必須對(duì)實(shí)測(cè)數(shù)據(jù)uobs(η,t)做很好的預(yù)處理,盡可能使其僅保留一次散射波場(chǎng)。

非線性LS_RTM的計(jì)算步驟如下:

1) 給定較準(zhǔn)確的背景偏移速度mB;

2) 計(jì)算波場(chǎng)ucal(η,t)和波場(chǎng)殘差r(η,t),并產(chǎn)生照明矩陣(即Hessian矩陣的對(duì)角陣);

3) 計(jì)算梯度g[δm(x),δmB(x),mB(x)],并給出RTM成像結(jié)果,當(dāng)成像道集拉平或數(shù)據(jù)擬合差達(dá)到要求時(shí),迭代停止;

4) 分離梯度g[δm(x),δmB(x),mB(x)]為gh[δm(x),δmB(x),mB(x)]和gl[δm(x),δmB(x),mB(x)];

5) 考慮梯度預(yù)條件和步長(zhǎng),利用(45)式得到δm和δmB當(dāng)前步的修正量;

6) 修改模型m(x)=mB(x)+δmB(x)+δm(x);

7) 當(dāng)修正量的范數(shù)‖δm‖和‖δmB‖小于預(yù)設(shè)值時(shí)計(jì)算結(jié)束,或返回步驟2)。

從上述步驟可以看出,非線性LS_RTM包含了速度擾動(dòng)量的反演和背景擾動(dòng)量的反演兩種成分,這兩種成分的聯(lián)合反演是目前地震波反演的導(dǎo)向性思想。

4.2 基于Born線性近似的LS_RTM

4.2.1 基于Born線性近似的LS_RTM方法

目前關(guān)于LS_RTM的討論主要基于(16)式,我們稱之為線性化的LS_RTM。為估計(jì)散射強(qiáng)度,定義如下誤差泛函:

(46)

目標(biāo)泛函的梯度是構(gòu)造最速下降法的核心,其定義為:

(47)

利用(47)式定義的梯度,可以進(jìn)行最速下降法的最小二乘偏移成像計(jì)算:

(48)

為了加速收斂,可以用Hessian矩陣的逆對(duì)梯度方向進(jìn)行預(yù)條件處理。Hessian矩陣定義為:

(49)

該Hessian矩陣是由背景介質(zhì)中的Green函數(shù)規(guī)定的,其中并不包括高階反射/繞射效應(yīng),即不包括非線性部分。雙向波計(jì)算Green函數(shù)時(shí)可以包含這部分的作用。由Hessian矩陣的逆對(duì)梯度方向進(jìn)行預(yù)條件處理后,最速下降法的迭代公式變?yōu)?

(50)

迭代步長(zhǎng)μ(k)可以采用試探的方式確定,一般采用試探加拋物擬合的方法來選取最優(yōu)下降步長(zhǎng),保證誤差的下降。

迭代的每一步需要在更新后對(duì)模型進(jìn)行一次反偏移,然后將反偏移得到的正演數(shù)據(jù)從原始數(shù)據(jù)中減去,判斷誤差是否達(dá)到要求,迭代是否終止。實(shí)質(zhì)上,用反偏移來產(chǎn)生與實(shí)測(cè)波場(chǎng)逼近的波場(chǎng)這個(gè)過程存在太多的問題。

具體的計(jì)算步驟為:

1) 給定偏移速度場(chǎng)vmig(x);

4) 用殘差波場(chǎng)進(jìn)行RTM,同時(shí)產(chǎn)生照明矩陣;

5) 用(50)式計(jì)算迭代成像結(jié)果;

6) 進(jìn)行RTM;

7) 轉(zhuǎn)步驟3)。

步驟3)也可以是解法方程。但是,形成法方程本身就需要巨大的計(jì)算量,解法方程也需要大量的內(nèi)存空間和計(jì)算量。解法方程的優(yōu)點(diǎn)是正則化和預(yù)條件的方法比較多?？梢钥闯?此處的梯度計(jì)算和波場(chǎng)殘差計(jì)算在理論上是不夠嚴(yán)謹(jǐn)?shù)摹?/p>

4.2.2 數(shù)值結(jié)果

圖3展示了用公式(50)計(jì)算的聲介質(zhì)最小二乘反射波疊前偏移成像結(jié)果?？梢钥闯?最小二乘偏移補(bǔ)償了照明不足,提高了成像的分辨率,使成像結(jié)果的橫向一致性更好,更有利于后續(xù)的儲(chǔ)層解釋。

圖4對(duì)比了聲介質(zhì)和粘聲介質(zhì)下最小二乘疊前深度偏移結(jié)果,以及Tikhonov正則化和加權(quán)正則化情形下迭代最小二乘疊前深度偏移結(jié)果?？梢钥闯?最后的粘聲介質(zhì)迭代20次的加權(quán)正則化疊前深度偏移結(jié)果與真實(shí)模型的反射系數(shù)很接近[13]。理論上講,如此高分辨率和高保真的疊前深度成像結(jié)果完全滿足巖性油氣藏的精細(xì)刻畫要求。但我們稱上述數(shù)值試驗(yàn)及結(jié)果是一個(gè)反演騙局。對(duì)實(shí)際數(shù)據(jù)而言,成像結(jié)果遠(yuǎn)不會(huì)如此理想,具體原因包括:①地震波正演不能很好地模擬實(shí)際觀測(cè)數(shù)據(jù);②預(yù)測(cè)誤差不是高斯分布的;③子波是未知的;④背景速度沒有滿足線性反演的要求;⑤觀測(cè)數(shù)據(jù)不完整(有限孔徑、不規(guī)則等);⑥影響振幅的因素復(fù)雜等。本質(zhì)上,線性化的LS_RTM是Born近似條件成立時(shí)的FWI。在當(dāng)前實(shí)際地震數(shù)據(jù)狀況下,用這樣的方式估計(jì)高波數(shù)參數(shù)擾動(dòng)量的真實(shí)值,其估計(jì)精度很難評(píng)價(jià)。

圖3 聲介質(zhì)迭代法最小二乘疊前深度偏移結(jié)果對(duì)比a 理論反射系數(shù); b 背景速度模型; c 常規(guī)疊前深度偏移結(jié)果; d 最小二乘疊前深度偏移結(jié)果

圖4 粘聲介質(zhì)最小二乘疊前深度偏移結(jié)果對(duì)比a 背景速度模型; b 品質(zhì)因子模型; c 聲介質(zhì)中模擬的單炮道集; d 粘聲介質(zhì)中模擬的單炮道集; e 粘聲介質(zhì)數(shù)據(jù)進(jìn)行聲介質(zhì)最小二乘疊前深度偏移的結(jié)果; f 粘聲介質(zhì)數(shù)據(jù)進(jìn)行粘聲介質(zhì)最小二乘疊前深度偏移的結(jié)果; g 一般Tikhonov正則化情況下的粘聲介質(zhì)最小二乘疊前深度偏移結(jié)果; h 加權(quán)正則化最小二乘疊前深度偏移結(jié)果

4.2.3 總變差(TV)正則化約束下的LS_RTM

對(duì)于實(shí)際數(shù)據(jù)而言,多重因素的制約使得最小二乘反演成像成為一個(gè)不適定的反問題。譬如:輸入的疊前數(shù)據(jù)中包含噪聲(非一次反射波)、不準(zhǔn)確的背景速度、正演(偏移)算子不能很好地模擬波現(xiàn)象、子波未知等因素會(huì)影響最小二乘反演成像迭代過程中的數(shù)據(jù)匹配,使不準(zhǔn)確的數(shù)據(jù)殘差反投影并混疊在梯度中,進(jìn)而影響最小二乘反演成像的收斂率。正則化是使不適定(非線性較強(qiáng)或凸性差)的地球物理反問題的解更穩(wěn)定、更有地質(zhì)意義的一類方法。常用的Tikhonov正則化可以提升反問題解的穩(wěn)定性,但這種正則化方法基于L2模,得到的是一光滑解,而我們所期望得到的是“棱角”更為分明的成像剖面,能突出地刻畫地下反射結(jié)構(gòu)。所以,具備保留圖像邊界性質(zhì)的總變差(TV)正則化是更為合適的選擇。

TV正則化形式如下:

R(m)=‖m‖TV

(51)

可利用如下迭代流程將TV正則化加入LSM迭代步驟中[14]。

初始化:u0=m0,k=0

k=k+1

上述迭代流程中,對(duì)每一步迭代中的像進(jìn)行TV去噪,并將去噪結(jié)果作為下一步迭代的先驗(yàn)信息,可以有效地去除像中的偏移假象和隨機(jī)噪聲,提高成像質(zhì)量。

圖5為實(shí)際數(shù)據(jù)常規(guī)LS-PSDM結(jié)果與TV正則化LS-PSDM結(jié)果對(duì)比,可以看出,TV正則化去除了成像結(jié)果中的隨機(jī)噪聲干擾,提高了同相軸的連續(xù)性,從而提高了成像質(zhì)量。

圖5 實(shí)際數(shù)據(jù)LS-PSDM結(jié)果對(duì)比a 常規(guī)LS-PSDM結(jié)果; b TV正則化LS-PSDM結(jié)果

5 針對(duì)反射系數(shù)的成像與針對(duì)散射強(qiáng)度的成像差異分析

到目前為止,本文討論的內(nèi)容都建立在地下介質(zhì)分布是背景+散射的基礎(chǔ)上。實(shí)質(zhì)上,勘探地震所面對(duì)的介質(zhì)可以描述為在空間上廣泛分布的層狀沉積層,加上火山活動(dòng)、構(gòu)造運(yùn)動(dòng)等引起的大、小尺度速度異常體(或波阻抗變化的異常體)。針對(duì)這樣的介質(zhì)情況,用背景+反射界面的方式來描述更為合適,體現(xiàn)了地下介質(zhì)的特征和在其中傳播的波場(chǎng)的特征。我們認(rèn)為,基于地下介質(zhì)特征和波場(chǎng)特征構(gòu)建的成像方法和技術(shù)才是最實(shí)用的。

前已述及,針對(duì)散射的成像需要滿足(28)式和(34)式規(guī)定的野外采集方式得到的數(shù)據(jù)。當(dāng)前,僅僅在地表進(jìn)行有限帶寬和有限孔徑的觀測(cè)顯然無法滿足逆散射成像的要求。事實(shí)上,到目前為止，勘探地震成像技術(shù)都是基于反射波的,成像的目標(biāo)是定位和估計(jì)反射界面的反射系數(shù)。勘探地震中,當(dāng)前反射波成像理論可以總結(jié)為:在最小二乘成像理論框架下,假設(shè)地下介質(zhì)是層狀的(每層都終結(jié)于無窮遠(yuǎn)邊界上)、觀測(cè)系統(tǒng)是無窮寬方位和長(zhǎng)偏移距及無假頻采樣的、背景速度正確(滿足Born近似)且子波已知,則方位角度反射系數(shù)可以準(zhǔn)確地估計(jì)出來并有一個(gè)半解析的表達(dá)形式[15]。COHEN等[16]和BLEISTEIN[17]等給出了上述成像方法完整的理論描述,這是勘探地震學(xué)家給出的反射波地震成像的最完整表達(dá)。

FWI(包括上述LS_RTM)都是針對(duì)散射體的成像,試圖估計(jì)散射強(qiáng)度f(x)。僅有地表觀測(cè)數(shù)據(jù),根據(jù)當(dāng)前FWI(LS_RTM)反演理論,很難給出反演解的精度評(píng)價(jià)。這是當(dāng)前基于散射和逆散射成像理論的主要困境。

研究人員已經(jīng)意識(shí)到了勘探地震的介質(zhì)特點(diǎn)和波傳播特點(diǎn),提出了改進(jìn)的LS_RTM方法:在LS_RTM的初始迭代中先進(jìn)行反射系數(shù)的估計(jì),即進(jìn)行反射波最小二乘疊前深度偏移成像,然后再進(jìn)行散射強(qiáng)度的估計(jì)。這一改進(jìn)可以有效地改善LS_RTM的收斂效率,提高最終成像質(zhì)量。

6 總結(jié)與討論

油氣地震勘探的目標(biāo)是含油氣儲(chǔ)層的識(shí)別與描述。最核心的技術(shù)是地震波成像,得到地下介質(zhì)的幾何特征圖像和彈性參數(shù)(主要是反射系數(shù)和波阻抗)。疊前深度偏移給出介質(zhì)的幾何特征圖像,FWI(LS_RTM)試圖給出彈性參數(shù)的估計(jì)。

基于散射波理論的FWI根本目標(biāo)是給出寬(全)波數(shù)帶的彈性參數(shù)估計(jì)。但是,它需要的數(shù)據(jù)是目前地表觀測(cè)得到的有限帶寬和有限孔徑的無假頻數(shù)據(jù)所無法滿足的。依據(jù)這樣不完整的數(shù)據(jù)得到的反演結(jié)果精度也很難評(píng)價(jià)。目前的FWI僅僅可以在低頻長(zhǎng)偏移距數(shù)據(jù)存在的情況下給出背景速度的估計(jì);LS_RTM是在假設(shè)背景速度正確情況下的散(反)射波FWI。數(shù)據(jù)觀測(cè)的不完整、疊前數(shù)據(jù)中的噪聲、不準(zhǔn)確的背景速度、正演(偏移)算子不能很好地模擬波現(xiàn)象、子波未知等因素導(dǎo)致LS_RTM很難收斂,解的精度難以評(píng)價(jià)。這是LS_RTM的困境,估計(jì)會(huì)一直伴隨LS_RTM方法技術(shù)的發(fā)展過程,無法得到根本解決。

方位角度反射系數(shù)的估計(jì)是勘探地震學(xué)所獨(dú)有的、針對(duì)勘探地震所面對(duì)的層狀沉積地質(zhì)特點(diǎn)以及波在層狀介質(zhì)中傳播的特點(diǎn)而發(fā)展起來的成像方法,勘探地震學(xué)發(fā)展到現(xiàn)在主要是基于這樣的反射波成像理論。

基于散射及逆散射成像的方法理論,到目前為止還沒有得到實(shí)用性的有效證明。但是,首先進(jìn)行基于反射系數(shù)的成像,在此基礎(chǔ)上進(jìn)行散射強(qiáng)度的成像無疑是可行的方法,可以有效地提高成像精度。針對(duì)散射波成像的正則化方法也有待進(jìn)一步深入研究。

總之,盡管FWI方法具有看起來完美的反演理論,但基于散(反)射波振幅的高波數(shù)參數(shù)擾動(dòng)量的估計(jì)目前缺乏有效的技術(shù)甚至缺乏明確的技術(shù)發(fā)展方向。這正是地震波成像方法理論研究需要著力的地方。

[1] WAPENAAR C P A.Inversion versus migration:a new perspective to an old discussion[J].Geophysics,1996,61(3):804-814

[2] CLAYTON R,STOLT R.A Born-WKBJ inversion method for acoustic reflection data[J].Geophysics,1981,46(11):1559-1567

[3] DEVANEY A J.Mathematical foundations of imaging,tomography and wavefield inversion[M].Cambridge:Cambridge University Press,2012:229-429

[4] SHUEY R.A simplification of the Zoeppritz equation[J].Geophysics,1985,50(4):609-614

[5] BLEISTEIN N.On the imaging of reflectors in the earth[J].Geophysics,1987,52(7):931-942

[6] BLEISTEIN N,COHEN J K,STOCKWELL J W.Mathematics of multidimensional seismic imaging,migration,and inversion[J].Applied Mechanics Reviews,2001,54(5):B94

[7] TARANTOLA A.Inversion of seismic reflection data in the acoustic approximation[J].Geophysics,1984,49(8):1259-1266

[8] STOLT R H,WEGLEIN A B.Migration and inversion of seismic data[J].Geophysics,1985,50(12):2458-2472

[9] STOLT R H,WEGLEIN A B.Seismic imaging and inversion[M].Cambridge:Cambridge University Press,2012:122-134

[10] WEGLEIN A B,VIOLETTE P B,KEHO T H.Using multiparameter Born theory to obtain certain exact multiparameter inversion goals[J].Geophysics,1986,51(5):1069-1074

[11] BORN M,WOLF E.Principles of optics[M].7th ed.Cambridge:Cambridge University Press,1999:695-732

[12] VAN LEEUWEN T,MULDER W A.A correlation-based misfit criterion for wave-equation traveltime tomography[J].Geophysical Journal International,2010,182(3):1383-1394

[13] 胡江濤.最小二乘逆時(shí)偏移及角度道集提取方法研究[D].上海:同濟(jì)大學(xué),2015 HU J T.Research on least-squares reverse time migration and angle gathers generation method[D].Shanghai:Tongji Universtiy,2015

[14] LIN Y,HUANG L.Least-squares reverse-time migration with modified total-variation regularization[J].Expanded Abstracts of 85thAnnual Internat SEG Mtg,2015:4264-4269

[15] LAMBARE G,OPERTO S,PODVIN P,et al.3D ray+Born migration/inversion,part I:theory[J].Geophysics,2003,68(4):1348-1356

[16] COHEN J K,HAGIN F,BLEISTEIN N.Three dimensional Born inversion with an arbitrary reference[J].Geophysics,1986,51(8):1552-1558

[17] BLEISTEIN N,COHEN J K,HAGIN F.Two and one-half dimensional Born inversion with an arbitrary reference[J].Geophysics,1987,52(1):26-36

(編輯：戴春秋)

Theory and method of least square prestack depth migration and imaging

WANG Huazhong1,HU Jiangtao1,2,GUO Song1

(1.Wavephenomenaandinversionimaginggroup(WPI),TongjiUniversity,Shanghai200092,China;2.StateKeyLaboratoryofOilandGasReservoirGeologyandExploitation,ChengduUniversityofTechnology,Chengdu610059,China)

The purpose of seismic imaging is firstly to determine the spatial location of the reflectivity or diffraction targets.After the structure of the reflection is settled,the rock parameters,related to velocity and density,can be estimated.Combined with the estimated elastic parameters and rock physics,the oil and gas reservoir can be delineated.We start from the seismic wave propagation representation of the scattered wave field and inverse scattering imaging.the formula of the scattered wave field is firstly given.The essence of the inverse scattering imaging is discussed,including the prestack seismic observed data and the corresponding observation mode required by inverse scattering imaging.Then,the theoretical framework of general linearized least-squares prestack depth migration (LS-PSDM) is given.The linearized LS-PSDM can be realized as long as the wave propagation operator calculated by Green functions can be reached.The least-squares reverse time migration (LS-RTM) is also discussed and the iterative realization method and the structure preserved total variation regularization method are given with numerical results.Finally,the differences between imaging methods based on reflectivity and scattering intensity estimation are analyzed.The most practical imaging methods aiming at the characteristics of medium and wave fields are stated.The inverse scattering imaging should be based on the methods.Further studies are required to testify whether the imaging methods based on scattered wave can satisfy the requirement of reservoir description.

scattered wave field representation,inverse scattering imaging,reflected wave imaging,least-squares prestack depth migration (LS_PSDM),full waveform inversion(FWI)

2016-10-09;改回日期：2016-12-06。

王華忠(1964—),男,教授,主要從事地震波傳播、地震數(shù)據(jù)分析和地震波反演成像方面的研究工作。

胡江濤(1987—),男,博士,主要從事地震信號(hào)分析、地震波成像和反演方面的研究工作。

國(guó)家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃(973計(jì)劃)(2011CB201002)、國(guó)家自然科學(xué)基金項(xiàng)目(41374117,41604100)和國(guó)家科技重大專項(xiàng)(2011ZX05005-005-008HZ,2011ZX05006-002,2011ZX05023)共同資助。

P631

A

1000-1441(2017)02-0159-12

10.3969/j.issn.1000-1441.2017.02.001

This research is financially supported by the National Key Basic Research and Development Program of China (973 Program) (Grant No.2011CB201002),the National Natural Science Foundation of China (Grant Nos.41374117,41604100),the National Science and Technology Major Project (Grant Nos.2011ZX05005-005-008HZ,2011ZX05006-002,2011ZX05023).

王華忠,胡江濤,郭頌.最小二乘疊前深度偏移成像理論與方法[J].石油物探,2017,56(2):-170

WANG Huazhong,HU Jiangtao,GUO Song.Theory and method of least square prestack depth migration and imaging[J].Geophysical Prospecting for Petroleum,2017,56(2):-170