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

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

?

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

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

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

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

散射場表達;逆散射成像;反射波成像;最小二乘疊前深度偏移;全波形反演

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

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

基于此,本文以地震波傳播的散射場表達與逆散射成像的本質為切入點,討論了一般意義下疊前偏移成像與逆散射成像之間的區(qū)別與聯(lián)系,重點分析了最小二乘疊前深度偏移成像方法存在的問題,認為勘探地震中的地震波成像要充分考慮地下介質以分層沉積為主的特點以及波在這樣的介質中傳播以反射波為主的特點構建合適的成像方法技術,而不是完全基于地下介質是背景+散射的模式構建逆散射成像技術。在方位角度反射系數(shù)估計結果比較準確的基礎上再進行逆散射成像應該是LS_RTM方法走向實用化的合理路線。

1 地震波傳播的散射場表達

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

1.1 基本的聲波方程和聲散射勢的引入

1.1.1 基本的聲波方程

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

(1)

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

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

(2a)

其中,

(2b)

為聲波傳播算子。

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

(3a)

其中,

(3b)

為背景介質中的聲波傳播算子。

1.1.2 聲散射勢的引入

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

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

(4)

其中,

(5a)

(5b)

1.2 標量波方程和標量波散射勢的引入

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

(6a)

其中,

(6b)

為標量波傳播算子。

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

(7a)

其中,

(7b)

為背景介質中的標量波傳播算子。

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

(8a)

其中,

(8b)

稱為散射強度。(8a)式定義的散射勢在速度擾動不大時還可以寫為如下形式:

(9)

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

(10)

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

1.3Lippmann-Schwinger方程的導出及散射場的具體表達

L=L0+L0VL

(11)

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

G=G0+G0VG

(12)

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

(13)

對(13)式右端項進行如下的級數(shù)展開,得:

G0VG0VG0+G0VG0VG0VG0+…

(14)

(12)式和(14)式被稱為Lippmann-Schwinger方程[8-9]。

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

G≈G0+G0VG0

(15)

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

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

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

(16)

其中,Born近似后的散射波場記為:

(17)

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

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

2 逆散射問題分析

2.1 常速背景介質與平面波入射情形下散射場與散射勢間的關系

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

(18)

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

(19)

一般地,在遠場近似假設下,入射場在常速介質中傳播可以用平面波近似:

(20)

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

(21)

(21)式描述的波的傳播關系如圖1所示。由(16)式知,此時總場可表示為:

(22)

圖1 空間局部散射體引起的散射波傳播關系

由圖1可知,當場點P距離散射體比較遠,散射體積源分布空間V有限時,用Q代表散射體積源內的點,O代表坐標原點,則可以做如下近似推導:

(23)

(24)

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

(25)

其中,

(26)

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

(27)

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

(28)

成立[11]。

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

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

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

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

(29)

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

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

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

(30)

在散射點周圍,對旅行時進行線性化,即圍繞中心射線進行Taylor展開:

(31)

式中:p(x0,x′)表示旅行時場T(x,x′)在點x0處的一階導數(shù),其方向為中心射線的切線方向。因此,(30)式中的指數(shù)部分可以化為:

(32)

其中,

(33)

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

圖2 地震波傳播過程中波矢量在反射點處的幾何關系

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

(34)

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

3 Bayes框架下基于Born近似散射場的散射勢估計方法

(17)式定義了一個Fredholm第一類積分方程,它描述的是一個線性問題。此處,我們將該Fredholm第一類積分方程的求解轉化為一個Bayes估計問題,就是將線性問題的求解轉化為利用泛函變分求極值的問題?；緝热荼硎鋈缦?。

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

Km=d

(35)

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

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

(36)

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

(37)

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

對于任意的介質擾動δm∈H1和擾動步長τ∈R,總存在如下的泛函變分:

(38)

對(38)式定義的變分求偏導數(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)

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

(41)

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

(42)

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

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

迭代的線性最小二乘反演方法中并不更新核函數(shù),因為我們假定核函數(shù)K[v0(x),x,ω]僅僅與背景速度v0(x)有關。但這僅僅是從計算效率方面考慮的結果,本質上還是應該更新核函數(shù),因為背景速度與擾動速度是相互耦合的。存在較強多次散射波(反射波)時,要么事先消除多次波,要么構建非線性的最小二乘反演成像方法。

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

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

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

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

4.1 非線性的LS_RTM

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

(43a)

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

反向波場外推方程為:

(43b)

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

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

誤差泛函梯度公式為:

(44)

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

(45a)

(45b)

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

非線性LS_RTM的計算步驟如下:

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

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

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

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

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

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

7) 當修正量的范數(shù)‖δm‖和‖δmB‖小于預設值時計算結束,或返回步驟2)。

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

4.2 基于Born線性近似的LS_RTM

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

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

(46)

目標泛函的梯度是構造最速下降法的核心,其定義為:

(47)

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

(48)

為了加速收斂,可以用Hessian矩陣的逆對梯度方向進行預條件處理。Hessian矩陣定義為:

(49)

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

(50)

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

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

具體的計算步驟為:

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

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

5) 用(50)式計算迭代成像結果;

6) 進行RTM;

7) 轉步驟3)。

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

4.2.2 數(shù)值結果

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

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

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

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

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

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

TV正則化形式如下:

R(m)=‖m‖TV

(51)

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

初始化:u0=m0,k=0

k=k+1

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

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

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

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

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

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

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

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

6 總結與討論

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

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

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

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

總之,盡管FWI方法具有看起來完美的反演理論,但基于散(反)射波振幅的高波數(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] 胡江濤.最小二乘逆時偏移及角度道集提取方法研究[D].上海:同濟大學,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—),男,博士,主要從事地震信號分析、地震波成像和反演方面的研究工作。

國家重點基礎研究發(fā)展計劃(973計劃)(2011CB201002)、國家自然科學基金項目(41374117,41604100)和國家科技重大專項(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