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

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

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

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

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

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

1 高斯束偏移成像公式

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

(1)

(2)

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

(3)

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

(4)

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

uGB=AGBexp(-iωTGB)。

(5)

其中

(6)

為高斯束振幅，而在全局直角坐標(biāo)系中

(7)

(8)

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

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

2 基于辛幾何算法的運(yùn)動(dòng)學(xué)射線追蹤

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

2.1 運(yùn)動(dòng)學(xué)射線追蹤方程組的三級(jí)四階辛RKN格式

運(yùn)動(dòng)學(xué)射線追蹤方程組可寫為

(9)

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

對(duì)于形如

(10)

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

(11)

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

(12)

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

(13)

更進(jìn)一步，若式還滿足i≤j時(shí)aij=0，則式為一個(gè)顯式辛RKN格式[13]。

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

(14)

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

(15)

2.2 射線追蹤方法精度與效率對(duì)比測(cè)試

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

圖1 運(yùn)動(dòng)學(xué)射線追蹤精度測(cè)試模型

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

(16)

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

表1 不同運(yùn)動(dòng)學(xué)射線追蹤方法誤差率

Note:①M(fèi)othod of kinematic ray tracing；②Error rate of fimal positionσx；③Error rate of traveling timeσt.

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

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

表2 運(yùn)動(dòng)學(xué)射線追蹤方法耗時(shí)

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

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

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

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

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

(17)

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

(18)

由其可得以下正則方程組

(19)

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

(20)

(21)

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

(22)

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

(23)

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

(24)

4 數(shù)值算例

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

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

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

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

5 結(jié)論

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

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

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