李沅衡 ， 王修田,2,3?? ， 宋 鵬,2,3 ， 姜秀萍,2,3 ， 趙 波,2,3

(1.中國海洋大學(xué)海洋地球科學(xué)學(xué)院， 山東 青島 266100；2. 青島海洋科學(xué)與技術(shù)試點國家實驗室 海洋礦產(chǎn)資源評價與探測技術(shù)功能實驗室， 山東 青島 266237；3.海底科學(xué)與探測技術(shù)教育部重點實驗室， 山東 青島 266100)

隨著油氣勘探開發(fā)的不斷深入，新的油氣藏愈發(fā)隱蔽，人們對油氣勘探能力提出了更加苛刻的要求。地震資料疊前深度偏移能夠?qū)?fù)雜構(gòu)造成像，是解決隱蔽油氣藏勘探的有力手段。目前疊前深度偏移方法主要分為射線類偏移與波動方程類偏移兩種，其中射線類偏移以克?；舴蚍e分法偏移為代表?？讼；舴蚍e分法疊前深度偏移具有占用計算機硬件資源小、成像精度較高且方法靈活的優(yōu)點，是目前實際資料處理中最為常用的方法，但是當?shù)叵聵?gòu)造非常復(fù)雜時，常規(guī)的克?；舴蚍e分偏移難以獲得高質(zhì)量的成像[1]。波動方程類偏移以逆時偏移為代表，其能夠?qū)?fù)雜構(gòu)造精確成像，但是運算量巨大，難以應(yīng)用于生產(chǎn)實踐[2-5]。

2001年Hill將高斯束方法應(yīng)用于疊前深度偏移，其綜合應(yīng)用高斯束運動學(xué)射線追蹤(即求取中心射線的旅行時和路徑)與復(fù)值動力學(xué)射線追蹤(即求取高斯束能量以及形態(tài))的信息，將經(jīng)傾角疊加后的局部平面地震信息映射到成像點處，不僅解決了常規(guī)射線法在焦散區(qū)失效問題，還避免了兩點射線追蹤的走時單一與射線丟失等問題，提高了在復(fù)雜構(gòu)造區(qū)域成像的精度[6]。針對Hill所提出高斯束疊前偏移方法僅僅能夠適用共偏移距道集的問題，Nowack與Gray分別將高斯束深度偏移方法應(yīng)用至炮集記錄，其中Gray所提出的方法可以看作是將Hill的高斯束深度偏移方法對炮集記錄的擴展[7-8]。為了解決傳統(tǒng)的高斯束偏移方法對振幅相位校正系數(shù)計算不精確的問題，Gray等基于廣義真振幅偏移理論，通過重新計算高斯束的振幅項，實現(xiàn)了保幅高斯束偏移[9]。目前，高斯束偏移方法可被認為是一種在計算精度與效率上均介于克希霍夫積分偏移與波動方程逆時偏移之間的成像方法，被視為是深度域克?；舴蚍e分類偏移的有效補充而受到高度重視，并在實際生產(chǎn)中得到一定的應(yīng)用[10-12]。

高斯束的運動學(xué)與動力學(xué)射線追蹤過程，從物理上講是一個典型的動力學(xué)過程。辛幾何算法是專門針對動力學(xué)過程設(shè)計的算法，可以提高動力學(xué)問題求解的精度與效率[13]。在應(yīng)用辛幾何算法進行運動學(xué)射線追蹤方面，諸多學(xué)者已進行過探索：陳景波與秦孟兆將辛幾何算法應(yīng)用于運動學(xué)射線追蹤中，并將辛幾何算法與Maslov方法結(jié)合，提出了能夠適應(yīng)焦散區(qū)的射線追蹤方法[14]。高亮應(yīng)用二階歐拉型辛差分格式進行運動學(xué)射線追蹤，并將結(jié)果與四階經(jīng)典Runge-Kutta(RK)方法比較，認為在精度差異不大的情況下，二階辛方法具有效率上的優(yōu)勢[15]。李川等將二維三次卷積模型插值算法與四階辛Partitioned-Runge-Kutta(PRK)法結(jié)合，進行了運動學(xué)射線追蹤方法研究，并認為相較于常規(guī)數(shù)值方法，辛幾何算法可以提高運動學(xué)射線追蹤的精度與效率[16]。

雖然前人應(yīng)用辛幾何算法進行了運動學(xué)射線追蹤的研究，并一致認為相較于常規(guī)算法，辛幾何算法在求解運動學(xué)射線追蹤方程組時具有一定優(yōu)勢，但是目前辛幾何算法尚未被應(yīng)用于高斯束動力學(xué)射線追蹤方程組的求解，也未實現(xiàn)基于辛幾何算法的高斯束偏移像。本文將基于辛幾何算法的運動學(xué)射線追蹤引入高斯束疊前深度偏移中，并在推導(dǎo)高斯束動力學(xué)射線追蹤方程組的三級四階辛RKN格式的基礎(chǔ)上實現(xiàn)基于辛幾何算法的高斯束疊前深度偏移。

1 高斯束偏移成像公式

本文選用炮集高斯束偏移“全波場算法”進行成像值的求取[17]。Zhang等給出了二維頻率域保幅互相關(guān)成像條件[18]

(1)

(2)

在式(2)中，psz與prz分別表示震源射線與檢波點射線在z方向上的慢度分量，G(x,xs,ω)與G(x,xr,ω)分別表示震源與檢波點波場的格林函數(shù)。將式代入式可得

(3)

在二維介質(zhì)中，用高斯束疊加積分合成格林函數(shù)的公式為

(4)

其中px′x與px′z分別表示點x′出射的高斯束水平和豎直慢度分量，uGB為高斯束表達式

uGB=AGBexp(-iωTGB)。

(5)

其中

(6)

為高斯束振幅，而在全局直角坐標系中

(7)

(8)

其中，uGB(x,xs,ω)=AGBSexp(-iωTGBS)表示震源高斯束，AGBS與TGBS為震源處高斯束的振幅以及旅行時；而檢波點處的高斯束用uGB(x,xr,ω)=AGBRexp(-iωTGBR)表示，AGBR與TGBR為檢波點高斯束的振幅與旅行時。

由高斯束偏移成像公式可知，實現(xiàn)高斯束偏移成像的關(guān)鍵是中心射線上點s坐標值、走時t(s)、慢度p、以及動力學(xué)參數(shù)Ps與Qs的求取，其中中心射線上點s坐標值、走時t(s)以及慢度p的求取過程被稱為運動學(xué)射線追蹤，而高斯束動力學(xué)參數(shù)Ps與Qs的求取過程被稱為動力學(xué)射線追蹤。

2 基于辛幾何算法的運動學(xué)射線追蹤

高斯束偏移成像需通過運動學(xué)射線追蹤獲得高斯束心射線上點s坐標值、走時t(s)和慢度p，在網(wǎng)格模型中，上述參數(shù)一般通過數(shù)值求解相應(yīng)的常微分方程組完成。射線追蹤的過程是一個動力學(xué)過程，在使用常規(guī)數(shù)值方法(例如龍格庫塔法等)求解時，難免會有人為耗散等歪曲動力學(xué)體系原特征的缺陷而影響計算精度與效率[13]。為提高運動學(xué)射線追蹤的精度與效率，本文應(yīng)用辛幾何算法求解高斯束運動學(xué)射線追蹤方程組。

2.1 運動學(xué)射線追蹤方程組的三級四階辛RKN格式

運動學(xué)射線追蹤方程組可寫為

(9)

其中：xi表示中心射線上點的坐標；pi表示中心射線的慢度；v為縱波速度；μ為任意射線參數(shù)。因為運動學(xué)射線追蹤方程組，滿足辛Runge-Kutta-Nystr?m(RKN)方法的求解條件，而四階格式精度較高，其能夠滿足運動學(xué)與動力學(xué)射線追蹤的精度要求[19]，與此同時，在相同階精度的格式中，顯式格式一般比隱式格式具有更高的效率，所以本文推導(dǎo)了針對運動學(xué)射線追蹤方程組的三級四階顯式辛RKN格式。

對于形如

(10)

的正則方程組，多級RKN方法具有如下格式

(11)

c1a11a12…a1sc2a21a22…a2s?????csas1as2…ass b1 b2 bsb1b2bs=cA bTbT。

(12)

只有當式中的f(xi)為標量函數(shù)的梯度且式中的系數(shù)滿足如下關(guān)系時，式才為辛格式。

(13)

更進一步，若式還滿足i≤j時aij=0，則式為一個顯式辛RKN格式[13]。

本文使用四階顯式算法求取運動學(xué)射線追蹤方程組。在所有四階顯式辛RKN格式中，被廣泛使用的是三級四階顯式辛RKN格式，其Butcher表為

(14)

將代入可導(dǎo)出針對運動學(xué)射線追蹤方程組的三級四階顯式辛RKN差分格式為

(15)

2.2 射線追蹤方法精度與效率對比測試

如圖1所示，運動學(xué)射線追蹤精度測試所使用的網(wǎng)格速度模型分為兩個速度層，其縱波速度分別為v1=1 500 m/s和v2=2 500 m/s，速度層分界面傾角為45°且與模型左邊界相交于30 000 m處，模型最大深度為40 000 m，網(wǎng)格間距為5 m×5 m。

圖1 運動學(xué)射線追蹤精度測試模型

測試射線以30°出射角由圖1中的點A(0,10)出射，由幾何關(guān)系可得，該射線與界面相交于點B(10 977.1,19 022.9)，由Snell定律可以計算透射射線的出射角度并最終可得射線與底邊之間的交點C的坐標為(18 383.1,40 000.0)，射線總旅行時為tsnell=23.535 s。在運動學(xué)射線追蹤試算中，所有追蹤格式步長均設(shè)為h=4 000 m2/s(注：若設(shè)Δs與v分別為空間步長與相應(yīng)速度，則h=vΔs，即在數(shù)學(xué)上相當于將空間步長Δs放大了v倍，實際計算時由于h為定值，則Δs將隨v的不同而變化)。在應(yīng)用數(shù)值方法獲得射線終點坐標的基礎(chǔ)上，求取終點橫向位置以及旅行時的誤差率

(16)

保留五位有效數(shù)字可得射線終點水平位置誤差率表(見表1)。

表1 不同運動學(xué)射線追蹤方法誤差率

Note:①Mothod of kinematic ray tracing；②Error rate of fimal positionσx；③Error rate of traveling timeσt.

分析表1可知，通過三級四階辛RKN算法獲得的射線路徑與走時誤差都是最低的，而經(jīng)典四階RK方法誤差最大，Adams預(yù)報校正法的求解精度介于二者之間。

針對運動學(xué)射線追蹤效率進行測試時，仍然選用如圖1所示的網(wǎng)格模型，震源位置仍設(shè)置為(0,10)，運動學(xué)射線追蹤步長為h=4 000 m2/s，計算每種運動學(xué)射線追蹤方法從20°～70°以0.01°為間隔的5 001根射線的總CPU耗時(運行程序的CPU平臺為主頻 2.33 GHz 的Intel Xeon E5410，測試時用C++語言編寫串行程序，編譯器版本為G++ 4.4.7，所有程序都使用O2級別優(yōu)化編譯，CPU耗時精確到毫秒)。所得運動學(xué)射線追蹤耗時見表2。

表2 運動學(xué)射線追蹤方法耗時

Note: ①Time-consuming；②No. of test.

從表2中可以看出，辛幾何算法的計算效率要高于兩種常規(guī)算法，而在常規(guī)方法中，Adams預(yù)報校正法的效率高于經(jīng)典四階RK方法。

3 基于辛幾何算法的動力學(xué)射線追蹤

高斯束偏移成像還需通過動力學(xué)射線追蹤求取高斯束動力學(xué)參數(shù)Ps與Qs。動力學(xué)射線追蹤的過程實際上是在已知中心射線的基礎(chǔ)上求解動力學(xué)射線追蹤方程組。在笛卡爾坐標系中，動力學(xué)射線追蹤方程組與運動學(xué)射線追蹤方程組具有相似的形式，因此理論上辛幾何算法同樣可以更好地針對動力學(xué)射線追蹤方程組進行數(shù)值求解。根據(jù)2.2節(jié)的數(shù)值實驗結(jié)果可以看出，三級四階辛RKN格式在精度與效率上都要優(yōu)于常規(guī)格式，因此本文選用三級四階辛RKN算法求解高斯束動力學(xué)射線追蹤方程組。

高斯束動力學(xué)參數(shù)為Ps與Qs滿足

(17)

設(shè)在笛卡爾坐標系下程函方程哈密爾頓函數(shù)形式為

(18)

由其可得以下正則方程組

(19)

對正則方程組兩邊同時對γ求取偏導(dǎo)數(shù)，并注意到?/?γ與d/dμ微分次序可調(diào)換，得

(20)

(21)

針對聲波方程，式中的各參數(shù)為

(22)

將式代入式，可得在全局直角坐標系下的動力學(xué)射線追蹤方程組

(23)

動力學(xué)射線追蹤方程組滿足辛RKN方法的應(yīng)用條件，因此可以應(yīng)用三級四階辛RKN算法求解，將式與式代入式可得動力學(xué)射線追蹤方程組的差分格式為

(24)

4 數(shù)值算例

本實驗應(yīng)用國際上使用的行業(yè)標準模型─Marmousi模型(見圖2)進行高斯束疊前深度偏移試算的數(shù)值實驗。Marmousi模型包含陡傾角斷層以及大量速度劇烈變化的復(fù)雜地層，其橫向網(wǎng)格間隔為5 m，總長度為9 200 m，縱向網(wǎng)格間隔為4 m，總深度為3 000 m；實驗中采用右邊放炮左邊接收的觀測系統(tǒng)，自2 575 m處開始放炮，炮點向右移動，炮間距為25 m，共265炮；最小偏移距為0 m，道間距為25 m，每炮104道；使用主頻為30 Hz的零相位雷克子波，應(yīng)用聲波方程有限差分法進行地震記錄正演模擬，記錄總長度為3 000 ms。成像點橫縱向間隔均為5 m。應(yīng)用全波場高斯束偏移方法，從震源和檢波點均以0.5°為間隔出射-75°～75°之間的301條高斯束，并選取高斯束初始寬度為w0=300，參考角頻率為ωr=10π進行高斯束偏移波場計算，所得疊前深度偏移剖面如圖3所示。

由圖3可以看出，基于辛幾何算法的高斯束疊前深度偏移剖面，波組特征明顯，斷層歸位精確，同相軸連續(xù)性較好。由此說明，基于辛幾何算法的高斯束疊前深度偏移具有對復(fù)雜模型精確成像的能力。

圖2 Marmousi網(wǎng)格速度模型

圖3 與Marmousi模型對應(yīng)的基于辛幾何算法的高斯束偏移剖面

5 結(jié)論

本文將基于辛幾何算法的運動學(xué)射線追蹤引入高斯束疊前深度偏移中，在推導(dǎo)了動力學(xué)射線追蹤方程組的辛差分格式基礎(chǔ)上，實現(xiàn)了基于辛幾何算法的高斯束疊前深度偏移。模型實驗表明：

(1)基于辛幾何算法的運動學(xué)射線追蹤，其效率與精度相比常規(guī)算法都具有一定優(yōu)勢，能夠滿足高斯束成像的運動學(xué)射線追蹤要求。

(2)基于辛幾何運動學(xué)射線追蹤與動力學(xué)射線追蹤算法的高斯束疊前深度偏移方法可適合于復(fù)雜構(gòu)造模型的精確成像。