梅勝全,鄧 飛,鐘本善,周熙襄

(成都理工大學(xué) 信息工程學(xué)院,四川成都 610059)

0 前言

三維地震射線走時(shí)計(jì)算是地震層析成像的基礎(chǔ),隨著三維高精度、高密度地震勘探的廣泛應(yīng)用,其作用也越來越突出。而走時(shí)計(jì)算的精度與效率問題[1、2],也一直是三維地震層析成像難以有效解決的難點(diǎn)。地震射線走時(shí)計(jì)算的精度與效率是一個(gè)矛盾問題,一方面較高的射線精度能更準(zhǔn)確地反映地震波場(chǎng)主能量傳播軌跡,但計(jì)算復(fù)雜度高,計(jì)算量大,計(jì)算效率會(huì)降低;另一方面高效、快速的走時(shí)計(jì)算,卻往往是以犧牲射線走時(shí)精度為代價(jià)的,這成為一個(gè)難以兩全的技術(shù)難點(diǎn)[1]。

三維地震走時(shí)計(jì)算的精度是關(guān)鍵,射線精度越高,其反映地震波場(chǎng)主能量傳播的射線路徑越合理,地震波走時(shí)擬合累積誤差越小,層析反演收斂就越快,反演的速度結(jié)構(gòu)也越接近真實(shí)情況,這在我國(guó)西部復(fù)雜近地表層析成像中尤其突出。若射線精度達(dá)不到一定的要求,其層析反演解的不適定性就嚴(yán)重,地震資料處理的成像效果就差。同時(shí),射線走時(shí)計(jì)算效率也是三維高精度三維地震勘探中的一個(gè)瓶頸問題,由于炮點(diǎn)、檢波點(diǎn)加密,道間距變小,導(dǎo)致射線走時(shí)計(jì)算數(shù)據(jù)規(guī)模呈指數(shù)級(jí)增長(zhǎng),計(jì)算效率成為三維地震層析成像實(shí)際應(yīng)用的一個(gè)巨大障礙。那么,如何有效地解決三維地震走時(shí)計(jì)算的精度與效率問題呢?

由于旅行時(shí)射線追蹤算法已經(jīng)很成熟了[5],大多以初值問題的試射法與邊值問題的彎曲法為基礎(chǔ),目前應(yīng)用廣泛的有最短路徑算法[6～8]和旅行時(shí)插值算法(LTI)[9]。眾多的文獻(xiàn)與應(yīng)用表明,雙線性旅行時(shí)插值算法能很好地模擬初至波傳播過程,計(jì)算的初至?xí)r間具有較高的精度[10、11],但其計(jì)算量大,計(jì)算效率相對(duì)不高。

作者在三維高精度地震層析成像應(yīng)用中,結(jié)合最短路徑算法和旅行時(shí)插值算法,通過對(duì)雙線性插值算法的改進(jìn),提出了一種三維射線追蹤算法。采用最短路徑算法尋找地震初至波的波前面,利用改進(jìn)的雙線性旅行時(shí)插值算法,計(jì)算網(wǎng)格單元的插值點(diǎn),進(jìn)行正演計(jì)算射線走時(shí)。這樣既保證了射線的精度,同時(shí)又簡(jiǎn)化計(jì)算,提高了計(jì)算效率。

1 三維初至波射線追蹤原理

在進(jìn)行三維射線追蹤時(shí),首先應(yīng)對(duì)模型進(jìn)行離散化處理,將地層模型轉(zhuǎn)化為以規(guī)則網(wǎng)格(一般為長(zhǎng)方體)為基本單元的三維網(wǎng)格模型,假設(shè)單元內(nèi)為常速,射線追蹤算法分為二個(gè)步驟:

(1)向前處理計(jì)算單元結(jié)點(diǎn)的時(shí)間場(chǎng)分布,從激發(fā)點(diǎn)出發(fā),基于波前面擴(kuò)展,利用初至波旅行時(shí)插值算法計(jì)算整個(gè)模型節(jié)點(diǎn)上的最小走時(shí)。

圖1 向前處理過程Fig.1 Forward processing procedure

(2)向后處理則利用互易原則,反向追蹤射線路徑,從接收點(diǎn)開始利用旅行時(shí)插值公式反向計(jì)算各單元的插值,直至激發(fā)點(diǎn),得到一條完整的初至射線路徑。

1.1 向前處理:計(jì)算模型節(jié)點(diǎn)時(shí)間場(chǎng)

向前處理過程用于計(jì)算激發(fā)點(diǎn)所覆蓋的網(wǎng)格節(jié)點(diǎn)的最小旅行時(shí)間,這是最核心也最復(fù)雜的一步,基本的算法步驟如下:

(1)在激發(fā)點(diǎn)單元,直接計(jì)算網(wǎng)格單元節(jié)點(diǎn)的旅行時(shí)間,如圖1(a)所示。

(2)以激發(fā)點(diǎn)單元為中心的,利用旅行時(shí)插值算法,逐一計(jì)算出六個(gè)相鄰單元的網(wǎng)格節(jié)點(diǎn)的旅行時(shí)間,如圖1(b)所示。

(3)使用最短路徑射線追蹤算法,求取新的次震源網(wǎng)格單元(假定次震源點(diǎn)為網(wǎng)格單元中心點(diǎn)),尋找最小走時(shí)的新的次震源點(diǎn),繼續(xù)向外插值計(jì)算或更新相鄰的非次震源單元網(wǎng)格的節(jié)點(diǎn)旅行時(shí),如圖1(c)所示。

(4)重復(fù)執(zhí)行步驟(3),直到整個(gè)模型的網(wǎng)格節(jié)點(diǎn)都作為次震源被計(jì)算完為止。

1.2 向后處理:反向追蹤射線路徑

向前處理計(jì)算了全部網(wǎng)格節(jié)點(diǎn)最小旅行時(shí)間場(chǎng)的分布。采用雙線性插值算法計(jì)算射線的入射點(diǎn),可以從長(zhǎng)方體單元中任意穿過,而不是在固定的網(wǎng)格節(jié)點(diǎn),或次震源路徑上,因此不能簡(jiǎn)單通過保留次震源點(diǎn),直接連接作為射線路徑,還需要利用旅行時(shí)插值算法,根據(jù)互換原則,反向確定滿足費(fèi)馬原理的最小旅行時(shí)節(jié)點(diǎn)位置,逐一插值計(jì)算出射線在網(wǎng)格單元的交點(diǎn),直至回到激發(fā)點(diǎn),連接而成為一條完整的射線路徑。

具體步驟如下:

(1)利用旅行時(shí)插值算法,對(duì)檢波點(diǎn)所在單元的六個(gè)面分別計(jì)算出初至射線穿出網(wǎng)格單元時(shí)與網(wǎng)格單元的交點(diǎn),找出最小旅行時(shí)的插值點(diǎn)即為射線入射點(diǎn),如圖2(a),這個(gè)最小旅行時(shí)插值就是檢波點(diǎn)的初至旅行時(shí)。

(2)根據(jù)當(dāng)前交點(diǎn)的位置,利用旅行時(shí)插值算法,繼續(xù)找出下一個(gè)網(wǎng)格單元的交點(diǎn)如圖2(b)所示,需要對(duì)交點(diǎn)的位置進(jìn)行判定。如果位于單元的一個(gè)面上,則由下一個(gè)網(wǎng)格單元決定;如果位于網(wǎng)格單元的棱線或者頂點(diǎn)上,則需要用周圍多個(gè)網(wǎng)格單元共同確定下一個(gè)插值最小的交點(diǎn)。

(3)如果射線已經(jīng)回到激發(fā)點(diǎn)所在網(wǎng)格單元,則已完成一條射線路徑的追蹤,否則繼續(xù)執(zhí)行步驟(2)。

圖2 向后處理過程Fig.2 Backward processing procedure

2 基于最短路徑的三維次震源波前面掃描算法

在向前處理過程中,作者在本文采用最短射線路徑算法不斷計(jì)算次震源點(diǎn),通過次震源點(diǎn)的最小旅行時(shí)排序,實(shí)現(xiàn)地震波前面的擴(kuò)展。并以次震源為中心,推算臨近網(wǎng)格節(jié)點(diǎn)的最小旅行時(shí),直至計(jì)算出整個(gè)模型中所有網(wǎng)格節(jié)點(diǎn)的最小旅行時(shí)。作者采用堆排序算法進(jìn)行波前面掃描,將網(wǎng)格單元分成二個(gè)集合,已經(jīng)獲得了最小旅行時(shí)的次震源集合與由其它網(wǎng)格插值計(jì)算旅行時(shí),但還未成為次震源的集合。每一次從堆中取出最小旅行時(shí)的網(wǎng)格單元,作為次震源點(diǎn),計(jì)算相鄰的六個(gè)網(wǎng)格,并加入堆中。如此反復(fù),直到所有的網(wǎng)格單元都成次震源點(diǎn),計(jì)算完所有單元最小旅行時(shí)為止。

在向后處理插值反向追射線路徑時(shí),不需要簡(jiǎn)單地插值計(jì)算六個(gè)面的最小旅行時(shí)交點(diǎn),而是先計(jì)算當(dāng)前交點(diǎn)與其相鄰次震源的最小旅行時(shí),完成次震源傳播路徑的反推,從而指導(dǎo)射線的入射方向。再應(yīng)用改進(jìn)的雙線性旅行時(shí)插值算法,快速求出下一個(gè)交點(diǎn),直至整條射線路徑追蹤至激發(fā)點(diǎn)。

這里采用Moser計(jì)算公式計(jì)算插值交點(diǎn)與鄰近網(wǎng)格次震源點(diǎn)的旅行時(shí)差[12],

其中 dij為插值交點(diǎn)與次震源點(diǎn)的歐氏距離;ˉs為相鄰網(wǎng)格計(jì)算單元的平均速度求得的平均慢度。

找出相鄰網(wǎng)格單元的次震源點(diǎn)與當(dāng)前插值交點(diǎn)的最小旅行時(shí),則相應(yīng)網(wǎng)格單元為射線下一個(gè)入射的方向。

通過這步次震源點(diǎn)掃描,可以最大程度地減少插值計(jì)算次數(shù)。

3 雙線性旅行時(shí)插值原理及改進(jìn)

旅行時(shí)插值算法的目的,是利用已知的最小旅行時(shí)間節(jié)點(diǎn),計(jì)算未知節(jié)點(diǎn)的最小旅行時(shí)間。在整個(gè)向前及向后處理過程中,都需要多次使用旅行時(shí)插值計(jì)算次震源附近的網(wǎng)格節(jié)點(diǎn)的最小旅行時(shí)間,對(duì)旅行時(shí)插值算法的優(yōu)化,可以在很大程度上提高整個(gè)算法的效率。作者在Asakawa提出的二維線性旅行時(shí)初至算法(LTI)基礎(chǔ)上[13],進(jìn)行算法推導(dǎo),得出簡(jiǎn)化的插值公式與判定條件,下面是詳細(xì)的算法說明。

在圖3中,A、B、C、D為四個(gè)已知最小旅行時(shí)的節(jié)點(diǎn),構(gòu)成一個(gè)空間矩形,在X方向的寬度為w,在Y方向的高度為h,記A點(diǎn)坐標(biāo)(0,0),最小旅行時(shí)記為t0,B點(diǎn)坐標(biāo)(1,0),最小旅行時(shí)記為t1,C點(diǎn)坐標(biāo)(1,1),最小旅行時(shí)記為t2,D點(diǎn)坐標(biāo)(0,1),最小旅行時(shí)記為t3,F點(diǎn)是待求的未知最小旅行時(shí)節(jié)點(diǎn),相對(duì)于A點(diǎn)的坐標(biāo)記為(x0,y0,z0),網(wǎng)格單元的慢度為s。

圖3 三維雙線性旅行時(shí)插值示意圖Fig.3 Three-dimensional bilinear travel-timeinterpolation scheme

在矩形ABCD上旅行時(shí)按雙線性變化,滿足

通過A、B、C、D四點(diǎn),可以求得式(3)的系數(shù):a=t0,b=t1-t0,c=t3-t0,d=t0+t2-t1-t3

在矩形ABCD內(nèi)求一點(diǎn)E,坐標(biāo)為(x,y),旅行時(shí)為Te,通過E點(diǎn)到達(dá)空間F點(diǎn)的旅行時(shí)為T。其中

要得到最小旅行時(shí)間,對(duì)式(4)求偏導(dǎo)得到方程組:

從式(5)中消去L后,可以得出:

對(duì)式(6)展開后,可得到:

從式(7)可見,最小旅行時(shí)點(diǎn)位于一個(gè)原點(diǎn)(-(cw-dx0)/2d,-(bh-dy0)/2d)的雙曲線上。圖4(見下頁)是一個(gè)可能解的示意圖,假設(shè)與矩形ABCD邊界相交有二個(gè)點(diǎn)H、K,則最小旅行時(shí)點(diǎn)應(yīng)該位于之內(nèi)。這是一簇解,一般采用迭代法求解,如牛頓～拉斐森方法,可以在有限步迭代收斂于真實(shí)解附近。

圖4 雙線性旅行時(shí)插值可行解示意圖Fig.4 The schematic diagram of bilinear travel-time interpolation feasible solution

但在三維射線追蹤,尤其在三維高密度勘探中,炮點(diǎn)數(shù)據(jù)與檢波點(diǎn)數(shù)據(jù)量大,這樣插值計(jì)算將會(huì)成為射線追蹤,耗費(fèi)大量計(jì)算時(shí)間,嚴(yán)重影響插值計(jì)算效率?；诓逯敌士紤],這里進(jìn)行插值算法改進(jìn),通過必要的近似處理,可以有效地減化運(yùn)算。

經(jīng)變換可解得:

同理可得:

這就得到了一個(gè)快速求解雙線性插值算法的公式,避免了迭代求解,可以直接快速得出插值點(diǎn),提高了計(jì)算效率。當(dāng)然與迭代法相比,在精度上有一定降低,但若控制在一個(gè)合理范圍內(nèi),則射線精度是能得到保證的。

下面給出快速插值計(jì)算的步驟:

(1)由式(7)計(jì)算與矩形ABCD的交點(diǎn),如果沒有交點(diǎn),則在邊界上應(yīng)用二維LTI得出二個(gè)解析解,比較求出最小的旅行時(shí)點(diǎn)作為插值點(diǎn)。

(2)應(yīng)用式(10)、式(11)快速求出插值點(diǎn),如果不滿足判定條件,則直接比較計(jì)算矩形ABCD的交點(diǎn)(如H、K點(diǎn))的最小旅行時(shí)點(diǎn)作為插值點(diǎn)。

改進(jìn)的雙線性插值算法的優(yōu)點(diǎn):式(10)、式(11)計(jì)算簡(jiǎn)單,通過波前面掃描算法與判定條件,能快速求出一個(gè)接近真實(shí)的解,最好情況是一個(gè)網(wǎng)格只需計(jì)算一次,最壞情況是計(jì)算邊上插值5次。相比于牛頓～拉斐森方法迭代法而言,一個(gè)網(wǎng)格的每個(gè)面平均8次～12次的插值計(jì)算,六個(gè)面的計(jì)算最好也有50次,最差70次以上。由此可以看出,計(jì)算量大幅度降低,同時(shí)射線精度也能最大限度得到保證。

圖5 一個(gè)復(fù)雜的三維山地資料Fig.5 A complex 3D mountain data example

圖6 改進(jìn)算法的射線追蹤(局部)Fig.6 The improved algorithm of ray tracing(local)

4 資料試算及誤差、效率分析

通過對(duì)一個(gè)復(fù)雜山地資料的射線追蹤實(shí)際資料處理,采用縱向速度梯度變化建立模型,地表速度為500 m/s,每10 m速度增量為5 m/s,炮點(diǎn)7 000炮,檢波點(diǎn)20 000個(gè),炮水平覆蓋范圍11 000 m,模型最大深度1 500 m,道間距25 m,最大偏移距2 500 m,采用普通微機(jī)進(jìn)行計(jì)算,平均單炮正演時(shí)間4 s左右。

為了驗(yàn)證改進(jìn)的雙線性插值算法的有效性和可靠性,通過牛頓～拉斐森方法迭代插值計(jì)算與改進(jìn)的雙線性插值算法快速計(jì)算的結(jié)果進(jìn)行對(duì)比,這里以二種算法計(jì)算的走時(shí)對(duì)比基本上能反映其走時(shí)計(jì)算的精度(如圖7及下頁圖8所示)。圖9(見下頁)是二種計(jì)算的運(yùn)算效率(每炮正演計(jì)算用時(shí))對(duì)比。

在圖7中,通過對(duì)二種算法的單炮射線走時(shí)來反映二種計(jì)算方法的差異。圖8中最大差值在20 ms,主要集中在近炮點(diǎn)附近,這說明在近炮點(diǎn)附近的計(jì)算誤差較大。當(dāng)偏移距越大,二者走時(shí)差值越小,平均在5 ms左右。同時(shí)可以看出,基于迭代計(jì)算的走時(shí)曲線比較光滑,而快速插值計(jì)算的走時(shí)曲線起伏較大,從炮點(diǎn)左右射線情況來看,炮點(diǎn)右邊射線走時(shí)差值更小。經(jīng)過分析發(fā)現(xiàn),與炮點(diǎn)左右二邊的地形起伏較大,網(wǎng)格間速度變化劇烈有較大的相關(guān)性。總的來說,二種算法計(jì)算的走時(shí)差值是在一個(gè)合理的范圍內(nèi),通過快速插值計(jì)算的走時(shí)可以代替迭代計(jì)算走時(shí),計(jì)算精度能得到保證。

圖7 單炮的各條射線走時(shí)對(duì)比(ms)Fig.7 A shot's all rays travel-time contrast(ms)

圖8 二種算法單炮射線走時(shí)差值對(duì)比(ms)Fig.8 A shot's all rays travel-time differential contrast with two algorithms(ms)

圖9是選取800炮計(jì)算用時(shí)的一個(gè)對(duì)比情況?？偟膩碚f,每炮計(jì)算用時(shí)差異不大,但二種算法的單炮計(jì)算用時(shí)則非常明顯,顯然,改進(jìn)的插值算法,計(jì)算效率得到了明顯提升。

圖9 800炮的射線追蹤運(yùn)算效率對(duì)比(s)Fig.9 The 800 shots'ray tracing operation efficiency contrast(in seconds)

5 結(jié)論

作者在本文提出的三維射線追蹤算法,在一定程度上解決了射線追蹤的精度與效率問題。在最短路徑算法的基礎(chǔ)上,提出了以次震源波前面掃描算法,代替了旅行時(shí)插值比較,從而減少了插值計(jì)算量。同時(shí)旅行時(shí)雙線性插值改進(jìn)算法原理簡(jiǎn)單,通過判定條件進(jìn)行快速插值,既保留了旅行時(shí)插值算法較高精度的優(yōu)點(diǎn),同時(shí)簡(jiǎn)化計(jì)算,效率也得到了明顯提升。通過數(shù)值計(jì)算表明,本方法可以適應(yīng)各種復(fù)雜的網(wǎng)格模型,能夠較準(zhǔn)確地計(jì)算地震波的初至?xí)r間和射線路徑,射線走時(shí)精度與效率能夠滿足層析反演需要。

在實(shí)際的復(fù)雜近地表地區(qū)資料處理中,采用此方法進(jìn)行射線追蹤正演模擬和反演成像,能獲得較精確的近地表模型。從理論模型與實(shí)際資料試算數(shù)值計(jì)算對(duì)比分析,說明了方法的計(jì)算精度和速度是可靠和有效地,也說明了在現(xiàn)有計(jì)算機(jī)硬件條件下進(jìn)行大規(guī)模的精確三維初至射線正演是可行的。

[1] JUL IAN B R,GUBBL INSD.Three-dimensional seismic ray tracing[J].J.Geophys,1977,43:95.

[2] N ISH I K.A three-d imensional robust seismic ray tracer for volcanic regions[J].Earth Planets Space,2001,53:101.

[3] 常旭,劉伊克.地震正反演與成像[M].北京:華文出版社,2001.

[4] 鄧飛,周杲,王美平,等.基于剖面重構(gòu)的三維地層建模[J].物探化探計(jì)算技術(shù),2007,29(1):19.

[5] 吳軍,賈雨.復(fù)雜介質(zhì)的三維塊狀模型快速射線追蹤[J].物探化探計(jì)算技術(shù),2008,30(6):465.

[6] 張建中,陳世軍,徐初偉.動(dòng)態(tài)網(wǎng)絡(luò)最短路徑射線追蹤[J].地球物理學(xué)報(bào),2004,47(5):146.

[7] MOSER T J.Shortest path calculation of seismic rays[J].Geophysics,1991,56(1):59.

[8] 孔選林,李錄明,羅省賢,等.各向異性介質(zhì)中地震波射線正演[J].物探化探計(jì)算技術(shù),2008,30(3):178.

[9] 張東,謝寶蓮,楊艷,等.一種改進(jìn)的線性走時(shí)插值射線追蹤算法[J].地球物理學(xué)報(bào),2009,52(1):200.

[10]高爾根,徐果明,蔣先藝,等.三維結(jié)構(gòu)下逐段迭代射線追蹤方法[J].石油地球物理勘探,2002,37(1):11.

[11]涂齊催,劉懷山.利用線性旅行時(shí)插值射線追蹤計(jì)算地表模型初至波走時(shí)[J].2006,30(2):148.

[12]MOSER T J,VANECK T,NOLET G.Hypocenter determination in strongly heterogeneous earth models using the shortest path method[J].J.Geophys.Res,1992,97:6563.

[13]ASAKAWA E,KAWANAKA T.Seismic ray tracing using linear traveltime interpolation[J].Geophysical Prospecting,1993,41(1):99.