陳生昌,魯方正,劉亞楠,周華敏

(1.浙江大學地球科學學院,浙江杭州310027;2.長江水利委員會長江科學院,湖北武漢430010)

在波形成像論文的第一篇中,對于給定的地震波運動學光滑模型,將擾動法應用于波動方程,得到了描述非均勻體擾動波場的波動方程。依據(jù)地震波長與非均勻體大小或非均勻體物理性質(zhì)空間變化特征尺度之間的相對大小關系,將非均勻體劃分為在空間上具有有限延續(xù)度可產(chǎn)生散射波的局部散射體或空間上具有一定延續(xù)度可產(chǎn)生反射波的反射體,并將相應的擾動波場的波動方程轉(zhuǎn)化為散射波動方程和反射波動方程。對于散射波,利用小擾動假定和Born近似,得到基于散射體物性參數(shù)相對擾動的散射數(shù)據(jù)線性正演表達式;對于反射波,利用波阻抗相對擾動的小值假定、一次波近似和高頻近似(地震波長相對于非均勻體大小或其物理性質(zhì)空間變化特征尺度為較小量),得到基于反射體波阻抗相對擾動的反射數(shù)據(jù)線性正演表達式。為了描述反射體邊界對反射波的作用和反射體邊界上物理量在反射數(shù)據(jù)偏移成像中的廣泛應用,在高頻近似條件(光滑模型的空間變化相對于地震波長尺度可視為緩慢變化甚至不變化)下定義反射體波阻抗相對擾動沿入射波傳播方向的方向?qū)?shù)為反射體邊界局部反射系數(shù),推導出了基于反射體局部反射系數(shù)的反射數(shù)據(jù)線性正演表達式,并且給出了標量波、聲波和彈性波方程下散射數(shù)據(jù)和反射數(shù)據(jù)的線性正演表達式。

本文是波形成像論文的第二篇,將第一篇推導出的地震一次散射數(shù)據(jù)和一次反射數(shù)據(jù)線性正演表達式作為地震數(shù)據(jù)線性反演的正演方程,利用線性反演方法提出包括散射數(shù)據(jù)波形成像和反射數(shù)據(jù)波形成像的地震數(shù)據(jù)波形成像方法理論。首先根據(jù)第一部分論述中提出的地震數(shù)據(jù)線性正演表達式,提出用于散射數(shù)據(jù)波形成像和反射數(shù)據(jù)波形成像的地震波傳播通用方程、散射通用方程和反射通用方程,再利用線性反演理論建立地震數(shù)據(jù)波形成像的基本方法理論框架,最后在波形成像的線性反演中利用地震波傳播算子的伴隨算子和最小二乘逆算子分別提出波形成像的波形偏移方法和最小二乘波形偏移方法。對于地震反射數(shù)據(jù)在兩種不同線性正演表達下的波形成像,首先給出面向地層波阻抗相對擾動的巖性成像和反射面局部反射系數(shù)的構造成像,然后針對具體的標量波散射數(shù)據(jù)、聲波反射數(shù)據(jù)和彈性波反射數(shù)據(jù)分別給出它們的波形偏移和最小二乘波形偏移具體計算公式,最后進行合成數(shù)據(jù)的波形成像數(shù)值試驗,試驗結果表明波形偏移結果相對于逆時偏移結果具有更加準確的相位和更高的分辨率。

1 地震數(shù)據(jù)波形成像的方法理論框架

地震數(shù)據(jù)的正演表達式是地震數(shù)據(jù)反演成像的基礎,不同的地震數(shù)據(jù)正演表達式會導致不同的反演目標和反演方法。當前地震數(shù)據(jù)波動方程偏移方法缺少與之對應的嚴謹?shù)牡卣饠?shù)據(jù)正演方程,使得偏移方法中描述波場傳播的方程與描述波場成像的成像公式不匹配,皆因波場傳播方程和成像公式是由不同條件下的波動方程導出的。偏移成像公式應該和產(chǎn)生散射波或反射波的虛源有關,即應該和散射波或反射波與入射波之間所滿足的數(shù)學物理方程有關。本文將描述散射波或反射波與入射波之間關系的數(shù)學物理方程稱之為散射(波)方程或反射(波)方程。在第一部分論述中關于地震數(shù)據(jù)線性正演表達的基礎上,首先建立與地震數(shù)據(jù)波形成像有關的地震波傳播、散射、反射的通用方程,然后在此基礎上利用反演理論構建地震數(shù)據(jù)波形成像的方法理論框架。

1.1 用于散射數(shù)據(jù)波形成像的地震波傳播與散射通用方程

當前的地震數(shù)據(jù)波動方程偏移成像不區(qū)分散射和反射,將無窮大平邊界平面波的鏡面反射方程(第一部分論述中的(6)式)作為描述散射波與入射波之間關系的散射方程,用于散射數(shù)據(jù)偏移成像。第一部分論述中將相對于地震波長尺度可近似為點源的虛源產(chǎn)生的地震波定義為散射波,因此我們認為在散射波偏移成像中采用第一部分論述中的(6)式描述散射不合適,它不能準確地描述局部的散射現(xiàn)象和散射波與入射波之間的相位變化。

根據(jù)第一部分論述中在光滑模型下建立的散射數(shù)據(jù)線性正演表達公式,我們建立用于描述地震波傳播與一次散射的通用波動方程,即:

L(mb)ui(x,t;xs)=f(x,t)

(1)

L(mb)us(x,t;xs)=-δLs(mb,δm)ui(x,t;xs)

(2)

式中:mb為已知的光滑模型;L(mb)為光滑模型下的通用波動算子;ui(x,t;xs)為光滑模型下的背景波場或入射波場;f(x,t)為震源函數(shù);xs為震源位置;us(x,t;xs)為一次散射波場;δLs(mb,δm)為散射算子,它的具體表達式與δm/mb和波動算子L(mb)的具體形式有關;δm為模型擾動,其空間尺寸大小或其物理性質(zhì)空間變化特征尺度相對于地震波長為較小量,滿足散射體條件。波動方程(1)描述了光滑模型下的入射波場及其傳播,它所對應的Green函數(shù)也是光滑模型下散射波的波場傳播算子。波動方程(2)不僅描述了散射波場與入射波場之間應滿足的數(shù)學物理關系,也描述了散射波的產(chǎn)生(即產(chǎn)生散射波的虛源)與散射波的傳播。我們將(1)式(或其對應的齊次方程)和(2)式統(tǒng)稱為用于散射數(shù)據(jù)波形成像的地震波傳播與散射通用方程,它們是由相同條件下的波動方程推導而來的,因此二者相互匹配,由此可以進一步推導出一次散射波數(shù)據(jù)的線性正演表達,該方程也是利用一次散射波對模型相對擾動δm/mb或模型擾動δm進行反演成像的基礎方程。

1.2 用于反射數(shù)據(jù)波形成像的地震波傳播與反射通用方程

在當前的地震數(shù)據(jù)波動方程逆時偏移方法中,將無窮大平邊界平面波的鏡面反射方程(第一部分論述中的(6)式)作為一般情況下反射波的反射方程,并用于反射數(shù)據(jù)偏移成像。這樣的鏡面反射方程不能準確地描述非平邊界、非平面波情況下的反射以及反射波與入射波之間的相位差異。根據(jù)第一部分論述中在光滑模型下利用高頻近似和局部平面波建立的反射數(shù)據(jù)線性正演表達公式,可以得到地震波傳播與反射的通用方程。對于基于反射體波阻抗相對擾動產(chǎn)生的一次反射,可以得到:

L(mb)ui(x,t;xs)=f(x,t)

(3)

L(mb)ur(x,t;xs)=-δLr(mb,Ir,σ)ui(x,t;xs)

(4)

式中:ur(x,t;xs)表示一次反射波場;δLr(mb,Ir,σ)表示隨角度σ變化的反射算子,它的表達式與波阻抗相對擾動Ir和波動算子L(mb)的關系式有關;σ表示入射方向與反射方向之間的開角。(3)式描述了光滑模型下的入射波場及其傳播規(guī)律,所對應的Green函數(shù)也是光滑模型下反射波的波場傳播算子。(4)式不僅反映了反射波場與入射波場之間應滿足的數(shù)學物理關系,也描述了基于相對擾動Ir的反射波產(chǎn)生(即產(chǎn)生反射波的虛源)與傳播規(guī)律。我們將(3)式(或其對應的齊次方程)和(4)式稱為基于波阻抗相對擾動Ir的反射數(shù)據(jù)波形成像地震波傳播與反射通用方程,它們均由相同條件下的波動方程推導而來,因此二者是相互匹配的,由此可以推導出基于波阻抗相對擾動Ir的一次反射波數(shù)據(jù)的線性正演表達公式,該方程也是利用一次反射波數(shù)據(jù)對波阻抗相對擾動Ir進行反演成像的基礎方程。

考慮基于反射體邊界局部反射系數(shù)產(chǎn)生的反射,推導得到如下的地震波傳播與反射的通用方程:

L(mb)ui(x,t;xs)=f(x,t)

(5)

L(mb)ur(x,t;xs)=-δLr(mb,r,σ)ui(x,t;xs)

(6)

式中:r表示反射體邊界局部反射系數(shù);δLr(mb,r,σ)表示基于r隨角度σ變化的反射算子,它的具體表達式與r和波動算子L(mb)的具體形式有關。(6)式不僅反映反射波場與入射波場之間應滿足的數(shù)學物理關系,也描述了基于局部反射系數(shù)的反射波的產(chǎn)生與傳播規(guī)律。我們將(5)式(或其對應的齊次方程)和(6)式稱為利用基于局部反射系數(shù)的反射數(shù)據(jù)波形成像的地震波傳播與反射通用方程,它們也是由相同條件下的波動方程推導而來,因此二者是相互匹配的,它們的結合可作為基于反射體邊界局部反射系數(shù)產(chǎn)生的一次反射波數(shù)據(jù)的線性正演表達,由此可以推導出利用一次反射波對反射體邊界局部反射系數(shù)r進行反演成像的基礎方程。

由上述分析可知,對于反射波和散射波而言,它們具有相同的入射方程和傳播方程(或傳播算子),卻具有不同的散射方程和反射方程。

1.3 地震數(shù)據(jù)波形成像的方法理論框架

根據(jù)上文中的地震波傳播通用方程、散射通用方程和反射通用方程,再利用(1)式所對應的Green函數(shù),可得到下述通用的散射數(shù)據(jù)或反射數(shù)據(jù)線性積分表達式:

(7)

式中:u(xg,t;xs)代表散射數(shù)據(jù)或反射數(shù)據(jù);xg代表檢波點坐標;Ω表示散射波或反射波虛源的空間分布范圍;g(xg,x,t)為Green函數(shù),代表從Ω內(nèi)x到xg的波場傳播算子;*t代表時間褶積;s(x,t;xs)為(2)式、(4)式或(6)式右端的散射波或反射波虛源項,該虛源項包含了光滑模型下的入射波場ui(x,t;xs)、產(chǎn)生散射或反射的模型參數(shù)和散射算子或反射算子。

利用(7)式左端項的地震數(shù)據(jù)u(xg,t;xs),估計右端積分中的虛源項s(x,t;xs),即為求解第一類線性積分方程,也稱為反源問題(inverse source problem)[1]。將(7)式代表的積分方程離散為矩陣方程形式,如下:

d=Lm

(8)

式中:d代表共炮點道集地震數(shù)據(jù)(即一次散射數(shù)據(jù)或一次反射數(shù)據(jù))矩陣;L代表由波場傳播算子組成的波場傳播矩陣;m表示待求的地下虛源矩陣。

對于矩陣方程(8),有下述的最小二乘解:

m=(L*L)-1L*d=L-gd

(9)

式中:L-g代表波場傳播矩陣L的最小二乘逆矩陣;L*代表矩陣L的伴隨矩陣。如果(L*L)為酉矩陣或忽略(L*L)-1的作用,即用L*作為L的逆矩陣,則可得到矩陣方程(8)的近似解ma,即:

ma=L*d

(10)

對于大規(guī)模反演成像而言,利用最小二乘方法求解(9)式存在逆矩陣計算困難和不穩(wěn)定性問題,但得到的解保真性和分辨率高。(10)式的求解方法稱為反投影法或反傳播法,具有計算量小、穩(wěn)定的特點,所得到的解保真性和分辨率不及最小二乘法所得結果。對于求解(9)式中存在的逆矩陣計算困難和解不穩(wěn)定問題,均可由迭代方法解決[2]。

求解(9)式或(10)式均可以得到產(chǎn)生散射波或反射波的虛源s(x,t;xs)的解估計。為了得到產(chǎn)生散射波的散射體物性參數(shù)擾動或產(chǎn)生反射波的反射體波阻抗相對擾動或反射體邊界局部反射系數(shù)的成像結果,我們還需要利用波動方程(1)計算虛源s(x,t;xs)中的入射波場ui(x,t;xs),然后再從s(x,t;xs)的解估計中消除ui(x,t;xs),同時還要考慮散射算子或反射算子的具體形式,并將其從s(x,t;xs)的解估計中消除。這些操作類似于常規(guī)逆時偏移中的波場成像運算,但在常規(guī)逆時偏移方法中不區(qū)分散射算子和反射算子之間的差異,并利用無窮大平界面平面波的鏡面反射算子R(σ)I近似本文中的散射算子和反射算子,R(σ)為隨角度σ變化的反射系數(shù)[3],I為單位算子。因此,偏移成像公式應該與散射或反射虛源中的散射算子或反射算子(或散射或反射虛源的表達式)有關。

根據(jù)上述內(nèi)容可得到如下的地震數(shù)據(jù)波形成像方法理論框架:

1) 利用觀測到的地震數(shù)據(jù)(散射數(shù)據(jù)或反射數(shù)據(jù))求解反源問題,得到對地下散射波或反射波虛源的解估計;

2) 利用光滑模型下的波動方程(1)求解地下入射波場;

3) 利用地下入射波場、散射算子或反射算子的具體形式,根據(jù)散射虛源或反射虛源的解估計得到散射體物性參數(shù)擾動或反射體波阻抗相對擾動或反射體邊界局部反射系數(shù)的成像。

在波形成像方法理論框架的反源問題求解中,利用(10)式得到虛源的解估計,對應的地震數(shù)據(jù)波形成像方法稱為地震數(shù)據(jù)波形偏移方法。在波形成像方法理論框架的反源問題求解中,利用(9)式得到虛源的解估計,則對應的地震數(shù)據(jù)波形成像方法稱為地震數(shù)據(jù)最小二乘波形偏移方法。

利用上述最小二乘波形偏移方法,在理論上可以得到散射體物性參數(shù)擾動、反射體波阻抗相對擾動、反射體邊界局部反射系數(shù)的高分辨率和高保真成像結果。但由于實際地震數(shù)據(jù)的頻帶局限性(特別是缺少低頻成分)和觀測范圍的有限性(缺失大偏移距觀測),導致得到的最小二乘波形偏移結果通常缺少低波數(shù)成分,僅僅是一種對散射體物性參數(shù)擾動或反射體波阻抗相對擾動的視反演,因此需要在最小二乘波形偏移中加入有關散射體物性參數(shù)擾動或反射體波阻抗相對擾動低波數(shù)成分的先驗信息,才能實現(xiàn)對散射體物性參數(shù)擾動值或反射體波阻抗相對擾動值的高分辨率和高保真反演。

2 標量波散射數(shù)據(jù)的波形成像

根據(jù)上節(jié)給出的地震數(shù)據(jù)波形成像方法理論框架,對于不同波動方程所描述的地震數(shù)據(jù)線性正演表達式,給出不同的地震波傳播方程和散射或反射方程,進而可以得到不同的地震數(shù)據(jù)波形成像公式。在給定的光滑速度模型下,對于標量波動方程所描述的散射線性數(shù)據(jù)正演表達式,存在下述具體的地震波傳播方程和散射方程,即:

(11)

(12)

(13)

式中:g(xg,x,t)表示光滑速度模型下方程(11)的Green函數(shù)。

將上節(jié)的地震數(shù)據(jù)波形成像方法理論框架及其波形偏移的伴隨求解方法以及(10)式應用于求解線性積分方程(13),再經(jīng)過一些簡單的推導可得到下述求取速度相對擾動av(x)近似估計的標量波散射數(shù)據(jù)波形偏移計算公式。

1) 地下入射波場的計算:

(14)

2) 基于邊值條件和波場逆時傳播的地下散射虛源近似反演:

(15)

3) 對av(x)的近似估計:

(16)

上述公式中:us(xg,t;xs)表示共炮點道集散射數(shù)據(jù);T表示地震數(shù)據(jù)的最大記錄時間。由于我們將產(chǎn)生散射波的散射體相對于地震波長近似為一個點源體,因此利用(16)式得到的av(x)的近似估計可以實現(xiàn)對散射體空間位置的成像[4]。

(16)式相較于常規(guī)逆時偏移成像公式的主要區(qū)別是多了對震源波場的時間二階導數(shù),因此(16)式得到的成像結果相對于常規(guī)逆時偏移結果存在180°的相位差異(即反極性),而且分辨率也會更高。

為了得到具有高保真性和高分辨率的av(x)反演結果,根據(jù)上節(jié)的地震數(shù)據(jù)波形偏移方法理論框架及最小二乘波形偏移的最小二乘公式,再結合上述的散射數(shù)據(jù)波形偏移計算公式,可通過迭代方式求解av(x)的散射數(shù)據(jù)最小二乘波形偏移計算公式。

① 初始化。

入射波場計算:利用公式(14)計算入射波場。

② 迭代過程

計算數(shù)據(jù)殘差:

(17)

數(shù)據(jù)殘差的逆時外推:

(18)

(19)

k=k+1,迭代結束。

修正公式(19)中的αk為第k次迭代步長因子(下同),可采用簡單的線性搜索法獲得[5],也可通過震源方向和檢波點方向的地震波雙向照明強度近似求解[6]得到。

3 聲波反射數(shù)據(jù)的波形成像

根據(jù)前文所述的反射數(shù)據(jù)波形成像的地震波傳播與反射通用方程,針對速度、密度變化的聲波動方程所描述的地震反射數(shù)據(jù),在給定光滑速度vb(x)、密度ρb(x)的模型下,對于基于波阻抗相對擾動產(chǎn)生的反射,有下述具體的地震波傳播方程和反射方程:

(20)

(21)

(22)

式中:g(xg,x,t)表示光滑速度、密度模型下聲波方程(20)的Green函數(shù)。

如果不考慮密度的變化,即δρ(x)=0,則與角度有關的聲波阻抗相對擾動簡化為與角度無關的速度相對擾動,即Ir(x,σ)=av(x),方程(20)至方程(22)簡化為下述描述基于速度相對擾動產(chǎn)生的反射所對應的方程:

(20-1)

(21-1)

(22-1)

對于基于反射體邊界局部反射系數(shù)產(chǎn)生的反射,可以得到具體的地震波傳播方程、反射方程和反射數(shù)據(jù)積分表達式,如下:

(23)

(24)

(25)

式中:r(x,σ)表示反射體邊界的局部反射系數(shù),是聲波阻抗相對擾動Ir(x,σ)沿入射波傳播方向的方向?qū)?shù),r(x,σ)=(iω)/[vb(x)]Ir(x,σ)=(iω)/[vb(x)][av+aρ(1+cosσ)]。方程(24)的右端項是產(chǎn)生反射波的虛源,包含局部反射系數(shù)、反射算子和入射波場,其中反射算子由光滑速度、密度模型和時間一階導數(shù)組成,即-1/[ρb(x)vb(x)]?/?t。

對于方程(23)至方程(25),如果不考慮密度的變化,則δρ(x)=0,可簡化為如下方程:

(23-1)

(24-1)

(25-1)

式中:r(x)代表反射體邊界的局部反射系數(shù),是速度相對擾動沿入射波傳播方向的方向?qū)?shù),r(x)=iω/vb(x)av(x)。

將上文的地震數(shù)據(jù)波形成像方法理論框架及波形偏移的伴隨求解方法公式(10)分別應用于線性積分方程(22)和方程(25),再經(jīng)過一些簡單的推導可以得到求取聲波阻抗相對擾動Ir(x,σ)和反射體邊界局部反射系數(shù)r(x,σ)近似估計的聲波反射數(shù)據(jù)波形偏移計算公式,如下。

1) 地下入射波場的計算:

(26)

2) 基于邊值條件和波場逆時傳播的地下反射波虛源近似反演:

(27)

3)對波阻抗相對擾動Ir(x,σ)和局部反射系數(shù)r(x,σ)的近似估計:

(28)

(29)

式中:ur(xg,t;xs)為共炮點道集反射數(shù)據(jù)。(28)式和(29)式相較于常規(guī)逆時偏移成像公式的主要區(qū)別在于多了時間二階導數(shù)和時間一階導數(shù),因此它們得到的成像結果相對于常規(guī)逆時偏移結果分別存在180°和90°的相位差異,此外分辨率有所提高。

為了得到σ角度域的Ir(x,σ)和r(x,σ),在進行(28)式和(29)式右邊的運算前需要將波場ur(x,t;xs)和ui(x,t;xs)轉(zhuǎn)換到相應的角度域[7-8]。由(28)式可知,利用角度域的共成像點道集Ir(x,σ)數(shù)據(jù)和Ir(x,σ)的表達式可以進行類似于基于Zoeppritz方程的振幅隨入射角變化(AVA)速度擾動、密度擾動和波阻抗擾動分析和反演。但是波阻抗相對擾動Ir(x,σ)不同于Zoeppritz方程中隨入射角變化的反射系數(shù),在地震波長相對于反射體的大小或其物理性質(zhì)空間變化的特征尺度滿足高頻近似條件下,其表達式是由波動方程導出的,而Zoeppritz方程是由波動方程在平邊界和平面波條件下導出的;再者Ir(x,σ)中的角度σ是根據(jù)局部切平面和局部平面波的概念定義的,是x點處的局部角度,不同于Zoeppritz方程中假定無窮大平邊界和平面波條件下定義的角度。由于本文定義的局部反射系數(shù)與入射波的局部波數(shù)有關(即與頻率有關),因此我們目前還不能推導出r(x,σ)與速度變化、密度變化、波阻抗變化和角度之間的具體函數(shù)關系,因此利用波形偏移得到的角度域共成像點道集r(x,σ)難以進行AVA速度分析、密度分析和反演。

如果不考慮密度的變化,即δρ(x)=0,將上文中的地震數(shù)據(jù)波形成像方法理論框架及其波形偏移的伴隨求解方法公式(10)應用于線性積分方程(22-1)和方程(25-1),再經(jīng)過一些簡單的推導可分別得到常密度聲波反射數(shù)據(jù)波形偏移計算公式、速度相對擾動av(x)和反射體邊界局部反射系數(shù)r(x)近似估計公式。

1) 地下入射波場的計算:

(26-1)

2) 基于邊值條件和波場逆時傳播的地下反射波虛源近似反演:

(27-1)

3) 對速度相對擾動av(x)和局部反射系數(shù)r(x)的近似估計:

(28-1)

(29-1)

由于標量波方程(或常密度聲波方程)所描述的反射波與入射角無關[9],所以標量波散射數(shù)據(jù)的波形成像公式(14)式至(16)式和常密度聲波反射數(shù)據(jù)的波形成像速度擾動公式(26-1)式至(28-1)式在形式上完全一致。

1) 初始化。

入射波場計算:利用(26)式計算入射波場。

2) 迭代過程。

計算數(shù)據(jù)殘差:

(30)

數(shù)據(jù)殘差的逆時外推:

(31)

(32)

k=k+1,迭代結束。

如果不考慮密度的變化,即δρ(x)=0,可以得到如下的速度相對擾動av(x)的反射數(shù)據(jù)最小二乘波形偏移。

1) 初始化。

入射波場計算:利用(26-1)式計算入射波場。

2) 迭代過程。

計算數(shù)據(jù)殘差:

(33)

數(shù)據(jù)殘差的逆時外推:

(34)

(35)

k=k+1,迭代結束。

1) 初始化。

入射波場計算:利用(26)式計算入射波場。

2) 迭代過程。

計算數(shù)據(jù)殘差:

(36)

數(shù)據(jù)殘差的逆時外推:

(37)

(38)

k=k+1,迭代結束。

如果不考慮密度的變化,即δρ(x)=0,得到下述局部反射系數(shù)下r(x)反射數(shù)據(jù)最小二乘波形偏移成像結果,r(x)是速度相對擾動沿入射波傳播方向的方向?qū)?shù)。

1) 初始化。

令k=0,rk(x)=r0(x)=0。

入射波場計算:利用(26-1)式計算入射波場。

2) 迭代過程。

計算數(shù)據(jù)殘差:

(39)

數(shù)據(jù)殘差的逆時外推:

(40)

對rk(x)的修正:

(41)

k=k+1,迭代結束。

4 彈性波反射數(shù)據(jù)的波形成像

對于縱、橫波速度和密度變化的各向同性彈性波動方程所描述的地震反射數(shù)據(jù),假定地震波在給定的縱、橫波速度和密度光滑模型αb(x),βb(x),ρb(x)下不發(fā)生散射、反射和波型轉(zhuǎn)換,利用反射數(shù)據(jù)波形成像的地震波傳播與反射通用方程,可以得到如下的彈性入射波方程和一次反射波方程:

Lbui(x,t;xs)=f(t,xs)

(42)

Lbur(x,t;xs)=ΔLrui(x,t;xs)

(43)

(44)

式中:g(xg,x,t)為光滑模型下彈性波方程(42)所對應的矢量Green函數(shù)。

由于笛卡爾坐標系下的彈性波場中P波和S波是耦合的,在彈性波反射數(shù)據(jù)的逆時偏移成像中,為了避免彈性波場的分量對分量成像產(chǎn)生的交叉噪聲以及得到易于物理解釋的成像結果,需要在成像前將彈性波場解耦為P/S波型波場[10-15]。本文將笛卡爾坐標系下的彈性波場轉(zhuǎn)換成以P、SH和SV波極化方向為坐標系下的P/S波型彈性波場,實現(xiàn)彈性波場的P波和S波解耦。

將光滑模型下的彈性波傳播算子Lb和反射算子ΔLr轉(zhuǎn)換成與P/S波型對應的算子,根據(jù)波形成像論文的第一篇,本文(43)式可表示為:

(45)

(46)

(47)

(48)

(49)

(50)

(51)

(52)

(53)

P波入射、P波反射的波場局部反射系數(shù)為:

(54)

SH波入射、SH波反射的波場局部反射系數(shù)為:

(55)

SV波入射、SV波反射的波場局部反射系數(shù)為:

(56)

P波入射、SV波反射的波場局部反射系數(shù)為:

(57)

SV波入射、P波反射的波場局部反射系數(shù)為:

(58)

由上述的rmn(m,n=p,sh,sv)表達式可知,rsvp,rpsv是關于σ的奇函數(shù),所以rsvp,rpsv在σ=0的兩側也會出現(xiàn)極性反轉(zhuǎn)。

利用上述的P波、SH波和SV波極化方向單位矢量和P波、SH波和SV波標量波場(本文的P波、SH波和SV波分量),還可以得到P、SH和SV波的矢量波場[12-13,17]。對入射波場ui(x,t;xs)而言,P波的矢量波場為:

(59)

SH波的矢量波場為:

(60)

SV波的矢量波場為:

(61)

對于反射波場ur(x,t;xs),P波的矢量波場為:

(62)

SH波的矢量波場為:

(63)

SV波的矢量波場為:

(64)

1) 地下入射彈性波場的計算:

Lbui(x,t;xs)=f(t,xs)

(65)

2) 基于邊值條件和波場逆時傳播的地下反射波虛源近似反演:

(66)

式中:ur(xg,t;xs)代表彈性反射波的炮集數(shù)據(jù)。

(67)

(68)

(69)

(70)

1) 初始化。

入射波場計算:利用公式(65)計算入射波場。

反射波場逆時外推:

(71)

2) 迭代過程。

計算分量位波場殘差:

(72)

(73)

k=k+1,迭代結束。

1) 初始化。

入射波場計算:利用公式(65)計算入射波場。

反射波場逆時外推:利用公式(71)進行反射波場逆時外推。

2) 迭代過程。

計算分量位波場殘差:

(74)

(75)

k=k+1,迭代結束。

(76)

(77)

(78)

(79)

由上述的標量波、聲波和彈性波數(shù)據(jù)的波形成像計算公式可知,波形成像方法中的波形偏移、最小二乘波形偏移的計算過程與逆時偏移、最小二乘逆時偏移的計算過程相同,計算量相當。對于地震數(shù)據(jù)的波形成像處理而言,如果波形偏移成像處理的目標是散射體,則將地震數(shù)據(jù)視為散射數(shù)據(jù)進行本文提出的散射數(shù)據(jù)波形偏移,可以得到散射體的波形偏移成像結果;如果偏移成像處理的目標是反射體(或反射面),則將地震數(shù)據(jù)視為反射數(shù)據(jù)進行本文提出的反射數(shù)據(jù)波形偏移,可以得到反射體(或反射面)的波形偏移成像結果。

上述提出的反射數(shù)據(jù)波形成像不僅能實現(xiàn)反射界面的準確成像,還能得到地層波阻抗相對變化(巖性)的成像,綜合應用波阻抗相對擾動的角度域共成像點道集數(shù)據(jù)和波阻抗相對擾動的表達式還可以進行類似于基于Zoeppritz方程的AVA分析和反演。

5 數(shù)值試驗

為了驗證上述波形成像方法理論,我們首先采用不考慮密度變化的標量波動方程合成地震數(shù)據(jù)進行標量波地震數(shù)據(jù)的波形成像試驗,然后再采用考慮密度變化的聲波動方程合成地震數(shù)據(jù)進行聲波地震數(shù)據(jù)的波形成像試驗。

5.1 標量波地震數(shù)據(jù)的波形成像試驗

5.1.1 Sigsbee2A模型試驗

采用SEG提供的Sigsbee2A模型數(shù)據(jù)驗證本文提出的散射數(shù)據(jù)波形偏移方法和反射數(shù)據(jù)局部反射系數(shù)波形偏移方法效果。圖1a為Sigsbee2A的地層模型,該模型中不僅包含高速鹽體、精細的地層,還包含許多局部高速散射體和底部的水平高速反射面,十分適合用于驗證散射數(shù)據(jù)波形偏移和反射數(shù)據(jù)波形偏移方法的效果。圖1b為用于波形偏移試驗的光滑速度模型,模型大小為1067×1201個網(wǎng)格點,水平和垂直方向的網(wǎng)格間距分別為22.86m和7.62m。該模擬數(shù)據(jù)的觀測系統(tǒng)為海上拖纜觀測系統(tǒng),共有500炮,炮間距為45.72m,拖纜長度為7932.42m,檢波點間距為22.86m,各炮的道數(shù)不均,最少為57道,最多為348道,時間采樣間隔為8ms,時間采樣點數(shù)為1500。

圖1 Sigsbee2A模型

分別應用常規(guī)逆時偏移方法(波形成像第一部分論述中的(1)式至(3)式)、本文提出的散射數(shù)據(jù)波形偏移方法計算公式(本文(14)式至(16)式)和反射數(shù)據(jù)局部反射系數(shù)波形偏移方法計算公式(本文(26-1)式至(29-1)式)對Sigbee2A模型數(shù)據(jù)進行試驗。圖2 為應用常規(guī)逆時偏移方法、反射數(shù)據(jù)波形偏移方法和散射數(shù)據(jù)波形偏移方法得到的Sigsbee2A模型偏移結果。

圖2 Sigbees2A模型偏移結果

由圖2中的偏移結果可以看出,3種偏移方法都取得了不錯的成像結果。對比圖2中的3種偏移結果可知,對于高速散射體的成像,散射數(shù)據(jù)波形偏移方法得到的散射體偏移結果要優(yōu)于常規(guī)偏移方法和反射波形偏移方法的偏移成像結果;對于反射面的成像,反射數(shù)據(jù)波形偏移方法得到的反射面的偏移成像結果要優(yōu)于常規(guī)偏移方法和散射波形偏移方法的偏移成像結果。

為了仔細對比常規(guī)逆時偏移、反射數(shù)據(jù)波形偏移和散射數(shù)據(jù)波形偏移對散射體和平反射面的成像結果的差異,我們分別對圖1所示的紅色方框中的高速散射體和藍色方框中的高速平反射面的成像結果進行對比。圖3為圖1紅色方框中的高速散射體模型和分別利用逆時偏移方法、反射數(shù)據(jù)波形偏移方法、散射數(shù)據(jù)波形偏移方法得到的成像結果。由圖3可以清楚地看出,散射數(shù)據(jù)波形偏移方法(其偏移結果見圖3d)對高速散射體成像(成像位置和零相位)準確;逆時偏移方法(其偏移結果見圖3b)對高速散射體成像位置準確,但出現(xiàn)了180°的相位差(反極性),將高速散射體成像為低速散射體;反射數(shù)據(jù)波形偏移方法(其偏移結果見圖3c)對高速散射體成像出現(xiàn)了90°的相位差,致使成像位置變淺,有1/4波長的深度誤差,高速散射體模型位于成像結果的波峰至波谷過渡帶的中間位置。由圖3還可以看出,散射數(shù)據(jù)波形偏移結果的分辨率最高,而逆時偏移結果的分辨率最低。圖3的試驗結果與本文的理論(見波形成像第一部分論述中逆時偏移方法中的(3)式和本文的(16)式(29-1)式)相符,波形成像第一部分論述中的(3)式與本文的(16)式相差一個時間二階導數(shù),即產(chǎn)生180°的相位差,本文的(29-1)式與(16)式相差一個負號和時間一階導數(shù),即產(chǎn)生-90°的相位差。本文的(16)式和(29-1)式中入射波場的時間導數(shù)有助于提高偏移結果的分辨率。

圖4為圖1藍色方框中的高速平反射面模型和分別利用逆時偏移方法、反射數(shù)據(jù)波形偏移方法、散射數(shù)據(jù)波形偏移方法得到的成像結果。由圖4可以清楚地看出,采用反射數(shù)據(jù)波形偏移方法(其偏移結果見圖4c)得到的高速平反射面成像結果(成像位置和零相位)準確;采用逆時偏移方法(其偏移結果見)圖4b 得到的高速平反射面成像結果存在90°的相位差,致使成像位置變淺,并存在1/4波長的深度誤差,高速平反射面模型位于成像結果的波峰至波谷過渡帶的中間位置;采用散射數(shù)據(jù)波形偏移方法(其偏移結果見圖4d)對高速平反射面成像出現(xiàn)了-90°的相位差,致使成像位置變深,并存在1/4波長的深度誤差,高速平反射面模型位于成像結果的波谷至波峰過渡帶的中間位置。圖4中的試驗結果與本文的理論相符,波形成像第一部分論述中的(3)式與本文的(28-1)式相差一個一階時間導數(shù),即產(chǎn)生90°的相位差,本文的(16)式與(29-1)式相差一個負號和時間一階導數(shù),即產(chǎn)生-90°的相位差。

圖4 高速平反射面模型(a)、逆時偏移(RTM)方法(b)、反射數(shù)據(jù)波形偏移(RWFM)方法(c)和散射數(shù)據(jù)波形偏移(SWFM)方法(d)成像結果

對比上述偏移成像結果可知,將本文提出的散射數(shù)據(jù)波形偏移方法應用于散射體成像、反射數(shù)據(jù)波形偏移方法應用于反射面成像,均可以取得理想的效果,它們的成像結果相對于逆時偏移結果不僅分辨率高,而且位置準確、相位(極性)正確。散射數(shù)據(jù)波形偏移方法有利于散射體的準確成像,反射數(shù)據(jù)波形偏移方法有利于反射面的準確成像。

5.1.2 反射數(shù)據(jù)的局部反射系數(shù)最小二乘波形偏移試驗

為了驗證本文提出的反射數(shù)據(jù)局部反射系數(shù)最小二乘波形偏移方法,我們截取Sigsbee2A速度模型(圖5a)的一部分進行試驗。對于圖5a的速度模型采用標量波動方程的有限差分方法合成地震數(shù)據(jù)。本次試驗合成了300炮地震數(shù)據(jù),炮點范圍是0～6km,炮間距20m,為中間放炮兩邊接收,每炮由601道接收,道間距10m,時間采樣間隔是0.5ms,時間采樣點數(shù)為3000,震源為主頻18Hz的Ricker子波。由于圖5a的速度模型為層狀結構,因此合成地震數(shù)據(jù)主要為反射數(shù)據(jù)。圖5b為用于最小二乘波形偏移試驗的光滑速度模型,它是對圖5a的速度模型進行光滑處理得到的。

在試驗中,我們首先應用反射數(shù)據(jù)局部反射系數(shù)波形偏移方法對合成的炮道集數(shù)據(jù)進行局部反射系數(shù)波形偏移,圖5c為得到的波形偏移結果。由圖5c可以看出,將局部反射系數(shù)波形偏移方法應用于該模型可獲得很好的成像結果。圖5d為應用本文提出的反射數(shù)據(jù)局部反射系數(shù)最小二乘波形偏移方法得到的偏移結果,共進行了35次迭代。對比圖5c和圖5d可以看出,圖5d的最小二乘波形偏移結果的分辨率相較于圖5c的波形偏移結果明顯得到了提高,并且層位面與斷層面的成像更清晰。由圖5e所示的模型真實垂向反射率可以看出,圖5d的最小二乘波形偏移結果的振幅保真性相較于圖5c的波形偏移結果也存在明顯提升。

圖5 部分Sigsbee2A模型及其波形偏移結果、最小二乘波形偏移結果和模型的理論垂向反射率

5.2 聲波地震數(shù)據(jù)的波形成像

為了驗證聲波地震數(shù)據(jù)的波形成像方法,采用考慮密度變化的聲波動方程合成地震數(shù)據(jù)進行聲波阻抗相對擾動的聲波反射數(shù)據(jù)波形偏移試驗。

圖6a為用于合成地震數(shù)據(jù)的速度模型,它是在常速為2000m/s的背景模型上疊加一些理論速度值在(-100m/s,100m/s)變化的隨機層。圖6b為用于合成地震數(shù)據(jù)的密度模型,它是在常密度為2000kg/m3的背景模型上疊加一些理論密度值在(-50kg/m3,50kg/m3)變化的隨機層。圖6c為利用圖6a和圖6b的速度、密度擾動、背景速度和密度得到的聲波阻抗相對擾動模型。利用聲波方程的有限差分方法合成地震數(shù)據(jù),本試驗合成了50炮地震數(shù)據(jù),炮點范圍為0～2km,炮間距為40m,為中間放炮兩邊接收,每炮由201道接收,道間距為10m,時間采樣間隔為1ms,時間采樣點數(shù)為3000,震源為主頻20Hz的Ricker子波。由于圖6的模型為層狀結構,因此合成的地震數(shù)據(jù)為反射數(shù)據(jù)。用于波形成像的速度和密度模型分別為速度為2000m/s的常速模型和密度為2000kg/m3的常密度模型。

利用(26)式至(28)式所表示的波阻抗相對擾動波形偏移計算公式對合成地震反射數(shù)據(jù)進行波阻抗相對擾動的波形偏移。圖7為波形偏移得到的結果,其中圖7a為波形偏移的角度域疊加結果,圖7b為x=1000m處的角度域共成像點道集,圖7c為x=1000m處角度域共成像點道集零度角對應的波阻抗相對擾動波形偏移與波阻抗相對擾動模型值對比結果(右上)和角度域共成像點道集中60°角對應的速度相對擾動波形偏移與速度相對擾動模型值對比結果(右下)。對比圖6和圖7中的波形偏移結果可以看出,本文提出的聲波阻抗相對擾動波形偏移方法可以很好地實現(xiàn)對聲波阻抗相對擾動的成像。由圖6c 和圖7a驗證了在(28)式中不進行波場角度分解而直接得到角度域的平均波阻抗相對擾動是一種快速的波阻抗相對擾動成像方法,可用來對地層波阻抗變化進行成像。由圖7c中的對比曲線可以看出,波阻抗相對擾動的角度域共成像點道集中的小角度道和大角度道可分別作為波阻抗相對擾動和速度相對擾動的近似估計。圖7b的波阻抗相對擾動角度域共成像點道集為進一步的角度域疊前巖性分析和反演提供了基礎數(shù)據(jù)。

圖6 用于波形成像試驗的速度模型(a)、密度模型(b)和相對聲波阻抗相對擾動模型(c)

圖7 波阻抗相對擾動的波形偏移結果

上述合成數(shù)據(jù)的數(shù)值試驗結果驗證了本文提出的標量波和聲波的散射數(shù)據(jù)波形偏移方法、反射數(shù)據(jù)波形偏移方法和最小二乘波形偏移方法的正確性和有效性。根據(jù)波形成像方法的計算公式,我們認為本文提出的波形成像方法具有很強的實用性。彈性波數(shù)據(jù)的波形成像數(shù)值試驗結果將在后續(xù)的文章中展示。

6 結論

對于波形成像的光滑背景模型,在地震波不發(fā)生散射、反射和波型轉(zhuǎn)換的假定下,基于地震數(shù)據(jù)的線性正演表達,本文首先給出了地震數(shù)據(jù)波形成像的地震波通用傳播方程、散射數(shù)據(jù)波形成像的通用散射方程和反射數(shù)據(jù)波形成像的通用反射方程,然后利用線性反演理論建立了面向散射數(shù)據(jù)和反射數(shù)據(jù)的地震數(shù)據(jù)波形成像方法理論框架。散射數(shù)據(jù)波形成像是對散射體物性參數(shù)擾動的線性反演,適合相對于地震波長較小尺寸的散射體和繞射體的成像。反射數(shù)據(jù)波形成像是對反射體波阻抗相對擾動或反射體邊界局部反射系數(shù)的線性反演,適合相對于地震波長較大尺寸的反射體物性和地層反射界面的成像。反射數(shù)據(jù)波形成像不僅能實現(xiàn)對反射界面位置的準確成像,還能得到對地層波阻抗相對變化的巖性成像。在波形成像中,如果利用波傳播算子的伴隨算子進行地震波的反傳播,則波形成像退化為波形偏移,如果利用波傳播算子的最小二乘逆算子進行地震波的反傳播,則波形成像退化為最小二乘波形偏移。波形偏移計算量小穩(wěn)定性好,但其偏移成像結果會因數(shù)據(jù)觀測系統(tǒng)和地下復雜構造的影響在保真性和分辨率方面存在不足。最小二乘波形偏移計算量大,但能有效校正數(shù)據(jù)觀測系統(tǒng)和地下復雜構造的影響,因此其偏移成像結果保真性好、分辨率高。我們給出了標量波散射數(shù)據(jù)、聲波反射數(shù)據(jù)和彈性波反射數(shù)據(jù)波形成像的計算式,并進行了標量波地震數(shù)據(jù)和聲波地震數(shù)據(jù)的波形成像數(shù)值試驗,取得了與理論相一致的試驗結果。我們提出的波形成像方法彌補了當前地震數(shù)據(jù)偏移成像方法理論上存在的不足,可以準確地利用地震數(shù)據(jù)的波形信息,修正當前逆時偏移方法存在的相位誤差,得到的偏移結果沒有相位誤差,分辨率也有所提高。與常規(guī)的逆時偏移和最小二乘逆時偏移的計算相比,波形偏移和最小二乘波形偏移僅增加了對入射波場的時間導數(shù)計算。反射體波阻抗相對擾動隨反射開角、巖性變化之間的定量表達式和反射數(shù)據(jù)波形成像得到的波阻抗相對擾動角度域共成像點道集為角度域成像數(shù)據(jù)的巖性分析和反演提供了方法理論和數(shù)據(jù)基礎,這是我們下一步研究工作的方向。