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

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

0 前言

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

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

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

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

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

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

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

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

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

1.1 向前處理:計算模型節(jié)點時間場

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

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

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

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

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

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

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

具體步驟如下:

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

同理可得:

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

下面給出快速插值計算的步驟:

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

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

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

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

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

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

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

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

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

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

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

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

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

5 結(jié)論

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

在實際的復(fù)雜近地表地區(qū)資料處理中,采用此方法進(jìn)行射線追蹤正演模擬和反演成像,能獲得較精確的近地表模型。從理論模型與實際資料試算數(shù)值計算對比分析,說明了方法的計算精度和速度是可靠和有效地,也說明了在現(xiàn)有計算機硬件條件下進(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].物探化探計算技術(shù),2007,29(1):19.

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

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

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

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

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

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

[11]涂齊催,劉懷山.利用線性旅行時插值射線追蹤計算地表模型初至波走時[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.