苗 賀 孫建國 王 蕤 顏鴻群 尹 暢 魏脯力(吉林大學(xué)地球探測科學(xué)與技術(shù)學(xué)院，吉林長春 130026)

1 引言

射線追蹤在地震勘探的各個(gè)領(lǐng)域得到了廣泛的應(yīng)用，同時(shí)對于油氣勘探中出現(xiàn)的復(fù)雜橫向非均勻介質(zhì)的正反演問題，提供了有力的解決手段。目前求解射線路徑的方法有很多，但大體分為兩類：一類是基于射線追蹤方程的拉格朗日法，它關(guān)注的是一點(diǎn)在地震波場中的運(yùn)動(dòng)軌跡，在計(jì)算過程中能夠同時(shí)計(jì)算出射線路徑和旅行時(shí)，但是這種方法計(jì)算速度慢，且存在盲區(qū)；另一類是對程函方程進(jìn)行數(shù)值分析的歐拉法，它關(guān)注的是固定網(wǎng)格點(diǎn)上行進(jìn)的時(shí)間，因?yàn)樵谶@種方法中唯一的控制方程是程函方程，相較于拉格朗日法，其計(jì)算速度要快得多，且不存在盲區(qū)，但是程函方程本身不能給出任何關(guān)于波傳播的方向，計(jì)算的正確性需要使用費(fèi)馬原理進(jìn)行檢查[1]，有限差分法[2,3]、最短路徑法[4-6]和旅行時(shí)線性插值法(LTI)[7-9]是三種具有代表性的方法[10]。有限差分法通過對程函方程在矩形網(wǎng)格上做差分離散近似計(jì)算給定速度模型的旅行時(shí)。最短路徑法是在計(jì)算區(qū)域內(nèi)形成一個(gè)由節(jié)點(diǎn)間的射線路徑構(gòu)成的網(wǎng)絡(luò)圖，然后根據(jù)Fermat原理選取旅行時(shí)最小的路徑作為射線路徑。Asakawa等[11]將LTI作為一種獨(dú)立的旅行時(shí)計(jì)算方法提出，并給出了相應(yīng)的射線路徑計(jì)算方法。該方法是基于網(wǎng)格單元模型計(jì)算邊界插值點(diǎn)的旅行時(shí)，并運(yùn)用Fermat原理選取具有最小旅行時(shí)的射線路徑。三維情況下LTI極小值問題為超越方程，一些學(xué)者采用網(wǎng)格界面剖分法[12]、快速插值法[13]、最速下降法[14]求解。然而，無論哪一種方法，其射線追蹤的步長都會(huì)受到網(wǎng)格限制，即次級(jí)震源點(diǎn)只能位于網(wǎng)格界面或節(jié)點(diǎn)上[10]，這會(huì)給路徑的求取帶來誤差。

旅行時(shí)場梯度法求取射線路徑可以很好地解決這一問題。旅行時(shí)場梯度法利用已知的旅行時(shí)場求得各級(jí)震源處的旅行時(shí)及其梯度。Rawlinson等[15,16]首先通過快速匹配追蹤法(Fast Marching Method，F(xiàn)MM)求取旅行時(shí)，然后在接收點(diǎn)沿旅行時(shí)梯度追蹤到震源點(diǎn)，旅行時(shí)梯度的求取通過一種類似于迎風(fēng)方法的有限差分近似得到； 張東等[10]應(yīng)用B樣條插值方法獲得連續(xù)的旅行時(shí)場，故追蹤步長不受網(wǎng)格尺寸的限制從而克服計(jì)算區(qū)域離散化所帶來的問題，但在速度突變區(qū)域計(jì)算次級(jí)源點(diǎn)的旅行時(shí)和梯度會(huì)產(chǎn)生一定的誤差[17]。通過采用函數(shù)逼近的方法能夠從根本上解決這個(gè)問題，本文應(yīng)用Chebyshev多項(xiàng)式逼近旅行時(shí)場、計(jì)算旅行時(shí)場梯度，不僅提高了射線路徑的精度，而且提高了計(jì)算效率。

2 三維旅行時(shí)場計(jì)算方法

FMM無條件穩(wěn)定且計(jì)算速度快，是一種滿足波前傳播規(guī)律的波前擴(kuò)展方式(窄帶技術(shù))[18-22]。3D情況下的程函方程為

=s2(x,y,z)

(1)

(2)

(3)

式中h為網(wǎng)格間距。因?yàn)橛?jì)算點(diǎn)的其中一個(gè)鄰居點(diǎn)一定是窄帶技術(shù)選取的局部最小旅行時(shí)點(diǎn)，而該點(diǎn)的旅行時(shí)必為已知。將三個(gè)方向上選取的差分格式代入式(2)中即可計(jì)算點(diǎn)(i,j,k)的旅行時(shí)梯度。FMM在實(shí)現(xiàn)時(shí)引入窄帶技術(shù)用于波前的近似模擬，運(yùn)用堆選排技術(shù)選取波前上(窄帶內(nèi))的最小旅行時(shí)點(diǎn)作為擴(kuò)展點(diǎn)，通過設(shè)定網(wǎng)格節(jié)點(diǎn)的屬性“凍結(jié)”被認(rèn)為已經(jīng)完成計(jì)算的網(wǎng)格節(jié)點(diǎn)，并利用窄帶內(nèi)網(wǎng)格節(jié)點(diǎn)的更新近似地模擬波前的擴(kuò)展演化[23-27]。

3 射線路徑計(jì)算方法

圖1 射線傳播路徑示意圖

3.1 一維逼近函數(shù)

函數(shù)逼近的定義為：設(shè)f(x)為[a，b]上的連續(xù)函數(shù)，由于其表達(dá)式難于求解，只能給出其在有限個(gè)點(diǎn)上的值，需要尋求一個(gè)近似函數(shù)P(x)(多項(xiàng)式)，使P(x)與f(x)在[a，b]上的誤差在某種意義下最小。逼近函數(shù)的基函數(shù)有很多種，例如：正交多項(xiàng)式、三次樣條插值、Pade逼近等。由于Chebyshev正交多項(xiàng)式能使插值余項(xiàng)達(dá)到最小，從而獲得最準(zhǔn)確的射線路徑，所以選擇Chebyshev正交多項(xiàng)式為基函數(shù)建立逼近函數(shù)。

第一類Chebyshev多項(xiàng)式Tn(x)的前四項(xiàng)為

(4)

取其前三項(xiàng)作為基函數(shù)，可得二階逼近函數(shù)為

P(x) =a0T0(x)+a1T1(x)+a2T2(x)

=a0+a1x+a2(2x2-1)

(5)

3.2 三維逼近函數(shù)

根據(jù)一維的二階逼近函數(shù)，可得三維二階逼近函數(shù)為

(6)

F(x,y,z)=a0+a1z+a2(2z2-1)+a3y+a4yz+a5y(2z2-1)+a6(2y2-1)+a7z(2y2-1)+a8(2y2-1)(2z2-1)+a9x+a10xz+

a11x(2z2-1)+a12xy+a13xyz+a14xy(2z2-1)+a15x(2y2-1)+a16(2y2-1)xz+

a17x(2y2-1)(2z2-1)+a18(2x2-1)+

a19(2x2-1)z+a20(2x2-1)(2z2-1)+

a21(2x2-1)(2y2-1)+a22(2x2-1)y+

a23(2x2-1)yz+a24(2x2-1)y(2z2-1)+

a25(2x2-1)(2y2-1)z+a26(2x2-1)×

(2y2-1)(2z2-1)

(7)

可寫為

(8)

式中Tt為正交函數(shù)簇，其系數(shù)為

(9)

(10)

為保證各項(xiàng)基函數(shù)正交，網(wǎng)格節(jié)點(diǎn)需要進(jìn)行如圖2所示局部坐標(biāo)變換

(11)

經(jīng)坐標(biāo)變換后，應(yīng)用式(9)求取多項(xiàng)式系數(shù)，再代入式(7)中可得到三維逼近函數(shù)。

與B樣條插值采用4×4×4網(wǎng)格不同，Chebyshev插值坐標(biāo)變換網(wǎng)格是以當(dāng)前追蹤點(diǎn)最臨近網(wǎng)格點(diǎn)為中心點(diǎn)的局部3×3×3網(wǎng)格。

對三維的逼近函數(shù)F(x，y，z)求導(dǎo)，即可求得旅行時(shí)場梯度

a13yz+a14y(2z2-1)+a15(2y2-1)+

a16(2y2-1)z+a17(2y2-1)(2z2-1)+

4a18x+4a19xz+4a20x(2z2-1)+

4a21x(2y2-1)+4a22xy+4a23xyz+

4a24xy(2z2-1)+4a25x(2y2-1)z+

4a26x(2y2-1)(2z2-1)

(12)

圖2 插值網(wǎng)格示意圖(a)插值網(wǎng)格； (b)插值坐標(biāo)變換網(wǎng)格

4a7yz+4a8y(2z2-1)+a12x+a13xz+

a14x(2z2-1)z+4a15xy+4a16xyz+

4a17xy(2z2-1)+4a21y(2x2-1)+

a22(2x2-1)+4a23xy+a24(2x2-1)(2z2-1)+4a25yz(2x2-1)+4a26y(2x2-1)(2z2-1)

(13)

4a14xyz+a16(2y2-1)x+4a17xz(2y2-1)+a19(2x2-1)+4a20z(2x2-1)+a23(2x2-1)y+4a24yz(2x2-1)+a25(2x2-1)(2y2-1)+4a26z(2x2-1)(2y2-1)

(14)

4 數(shù)值模擬

應(yīng)用速度連續(xù)變化模型和Marmousi模型，對Chebyshev插值射線追蹤算法與B樣條插值射線追蹤算法進(jìn)行對比分析。

4.1 速度連續(xù)變化模型

模型尺寸為500m×500m×500m，網(wǎng)格間距為10m，源點(diǎn)置于(1.0m，1.0m，1.0m)點(diǎn)處，接收點(diǎn)位于模型上表面均勻分布，間距為20m。模型的速度分布函數(shù)為(500+60z)(單位為m/s)。圖3a為Chebyshev插值法計(jì)算的射線路徑，圖3b為B樣條插值方法計(jì)算的射線路徑。圖4a和圖4b對比了兩種方法在y=25m處的xz剖面射線路徑，二者整體相差不大，只是Chebyshev方法計(jì)算的射線路徑在源點(diǎn)附近精度要高一些。圖5為兩種方法的旅行時(shí)相對誤差對比，可以看出，Chebyshev插值法精度明顯高于B樣條插值法。采用同樣的計(jì)算機(jī)硬件(酷睿i5 2450M/CPU，6G/RAM)，Chebyshev插值法用時(shí)為4.205s，B樣條插值法用時(shí)為5.191s，效率提高了約20%。

圖3 速度連續(xù)變化模型兩種方法射線路徑三維對比(a)Chebyshev插值法； (b)B樣條插值法

圖4 速度連續(xù)變化模型兩種方法射線路徑剖面對比(a)B樣條插值法； (b)Chebyshev插值法

圖5 速度連續(xù)變化模型兩種方法計(jì)算的旅行時(shí)誤差

4.2 Marmousi模型

為了檢驗(yàn)本文算法在面對復(fù)雜介質(zhì)時(shí)的有效性，由二維Marmousi模型生成的三維模型，模型尺寸為3800m×3800m×1220m，計(jì)算時(shí)采用的網(wǎng)格間距為10m，震源置于(0，0，0)處。接收點(diǎn)位于模型上表面、下底面與右側(cè)面均勻分布，間距為100m。圖6a為FMM計(jì)算的二維模型旅行時(shí)分布。圖6b和圖6c分別為B樣條插值和Chebyshev插值法計(jì)算的二維射線路徑，每條射線都垂直于等時(shí)線。 由圖5d和圖5e可見，兩種方法都能夠很好地求解復(fù)雜模型的三維射線路徑。但對比兩種方法的旅行時(shí)與原始旅行時(shí)場的相對誤差可知，Chebyshev方法計(jì)算的路徑誤差比較小(圖7)。另外，Chebyshev插值法用時(shí)為80.892s，B樣條插值法用時(shí)為102.060s，前者效率較后者提高了約20%。

圖6 Marmousi模型射線追蹤結(jié)果(a)二維模型旅行時(shí)； (b)B樣條插值法計(jì)算的二維射線路徑； (c)Chebyshev插值法計(jì)算的二維射線路徑； (d)Chebyshev插值計(jì)算的三維射線路徑； (e)B樣條插值計(jì)算的三維射線路徑

5 結(jié)束語

本文應(yīng)用Chebyshev插值法擬合旅行時(shí)場并求取旅行時(shí)場梯度，避免了對離散旅行時(shí)場的數(shù)值求導(dǎo)過程。相較于B樣條插值法，Chebyshev插值法能夠較好地得到計(jì)算區(qū)域內(nèi)任一點(diǎn)的梯度，從而獲得更準(zhǔn)確的射線路徑。數(shù)值實(shí)驗(yàn)表明該方法在構(gòu)造復(fù)雜時(shí)精度優(yōu)于B樣條插值法射線追蹤，且運(yùn)行效率高于B樣條插值法。本文射線追蹤方式為單條追蹤，如果能夠同時(shí)計(jì)算出多個(gè)接收點(diǎn)的路徑，勢必能夠大幅提升計(jì)算效率。

參考文獻(xiàn)

[1] Sun J，Sun Z，Han F.A finite difference scheme for solving the eikonal equation including surface topo-graphy.Geophysics，2011，76(4)：T53-T63.

[2] Vidale J E.Finite-difference calculation of traveltimes in 3-D.SEG Technical Program Expanded Abstracts，1989，8：1096-1098.

[3] Vidale J.Finite-difference calculation of travel times.Bulletin of the Seismological Society of America，1988，78(6)：2062-2076.

[4] Moser T J.Shortest path calculation of seismic rays.Geophysics，1991，56(1)：59-67.

[5] 唐小平，白超英，劉寬厚.分區(qū)多步最短路徑極值法多值多次反射波追蹤.地球物理學(xué)進(jìn)展，2011，26(6)：2064-2074.

Tang Xiaoping，Bai Chaoying，Liu Kuanhou.Multiva-lued and multiple reflected raytracing with extreme value based on the multistage modified shortest-path method.Progress in Geophysics，2011，26(6)：2064-2074.

[6] 趙后越，張美根.起伏地表?xiàng)l件下各向異性地震波最短路徑射線追蹤.地球物理學(xué)報(bào)，2014，57(9)：2910-2917.

Zhao Houyue，Zhang Meigen.Tracing seismic shortest path rays in anisotropic medium with roiling surface. Chinese Journal of Geophysics，2014，57(9)：2910-2917.

[7] 李培明，梅勝全，馬青坡.一種改進(jìn)的雙線性插值射線追蹤方法.石油地球物理勘探，2013，48(4)：553-558.

Li Peiming，Mei Shengquan and Ma Qingbo.An improved bilinear interpolation travel-time ray-tracing method.OGP，2013，48(4)：553-558.

[8] 韓復(fù)興，孫建國，楊昊.基于二維三次卷積插值算法的波前構(gòu)建射線追蹤.吉林大學(xué)學(xué)報(bào)(地球科學(xué)版)，2008， 38(2)：336-340.

Han Fuxing，Sun Jianguo，Yang Hao.Ray-tracing of wavefront construction by bicubic convolution interpolation.Journal of Jilin University (Earth Science Edition)，2008，38(2)：336-340.

[9] 張婷婷，張東，邱達(dá).三維層狀介質(zhì)中基于旅行時(shí)梯度的多次波射線追蹤.石油地球物理勘探，2014，49(6)：1097-1105.

Zhang Tingting，Zhang Dong and Qiu Da.Multiple ray tracing in 3D layered media using the traveltime gradient method.OGP，2014，49(6)：1097-1105.

[10] 張東，張婷婷，喬友鋒等.三維旅行時(shí)場B樣條插值射線追蹤方法.石油地球物理勘探，2013，48(4)：559-566.

Zhang Dong，Zhang Tingting，Qiao Youfeng et al.A 3-D ray tracing method based on B-spline traveltime interpolation.OGP，2013，48(4)：559-566.

[11] Asakawa E，Kanawa T.Seismic ray tracing using linear traveltime interpolation.Geophysical Prospecting，1993，41(1)：99-111.

[12] 張東，傅相如，楊艷等.基于LTI和網(wǎng)格界面剖分的三維地震射線追蹤算法.地球物理學(xué)報(bào)，2009，52(9)：2370-2376.

Zhang Dong，F(xiàn)u Xiangru，Yang Yan et al.3-D seismic ray tracing algorithm based on LTI and partition of grid interface.Chinese Journa of Geophysics，2009，52(9)：2370-2376.

[13] 劉鋒，張東，楊艷等.三維LTI射線追蹤極小值方程的快速數(shù)值解法.武漢大學(xué)學(xué)報(bào)(理學(xué)版)，2012，58(5)：395-400.

Liu Feng，Zhang Dong，Yang Yan et al.A fast nume-rical solution to minimum equation in 3-D seismic LTI ray tracing.Journal of Wuhan University (Natural Science Edition).2012，58(5)：395-400.

[14] 金昌昆，張建中.共接收面元三維折射旅行時(shí)反演方法.石油地球物理勘探，2016，51(4)：751-759.

Jin Changkun and Zhang Jianzhong.3D refraction traveltime inversion in common receiver-bins. OGP，2016，51(4)：751-759.

[15] Rawlinson N，Sambridge M.Wave front evolution in strongly heterogeneous layered media using the fast marching method.Geophysical Journal International，2004，156(3)：631-647.

[16] Rawlinson N，Sambridge M.Multiple reflection and transmission phases in complex layered media using a multistage fast marching method.Geophysics，2004，69(5)：1338-1350.

[17] 張婷婷，邱達(dá)，張東.一種改進(jìn)的三維旅行時(shí)梯度射線追蹤方法.石油地球物理勘探，2016，51(5)：916-923.

Zhang Tingting，Qiu Da and Zhang Dong.A modified 3D traveltime gradient ray tracing method.OGP，2016，51(5)：916-923.

[18] Sethian J A.Evolution，implementation and applica-tion of level set and fast marching methods for advancing fronts.Journal of Computational Physics，2001，169(2)：503-555.

[19] Sethian J A.A fast marching level set method for monotonically advancing fronts.Proceedings of the National Academy of Sciences，1996，93(4)：1591-1595.

[20] 黃興國，孫建國，孫章慶等.基于復(fù)程函方程和改進(jìn)的快速推進(jìn)法的復(fù)旅行時(shí)計(jì)算方法.石油地球物理勘探，2016，51(6)：1109-1118.

Huang Xingguo，Sun Jianguo，Sun Zhangqing et al.A perturbation approach for solving the TTI complex eikonal equation.OGP，2016，51(6)：1109-1118.

[21] 李永博，李慶春，吳瓊等.快速行進(jìn)法射線追蹤提高旅行時(shí)計(jì)算精度和效率的改進(jìn)措施.石油地球物理勘探，2016，51(3)：467-473.

Li Yongbo，Li Qingchun，Wu Qiong et al.Improved fast marching method for higher calculation accuracy and efficiency of traveltime.OGP，2016，51(3)：467-473.

[22] 王飛，曲昕馨，劉四新等.一種新的基于多模板快速推進(jìn)算法和最速下降法的射線追蹤方法.石油地球物理勘探，2014，49(6)：1106-1114.

Wang Fei，Qu Xinxin，Liu Sixin et al.A new ray tra-cing approach based on both multi-stencils fast mar-ching and the steepest descent.OGP，2014，49(6)：1106-1114.

[23] 張風(fēng)雪，吳慶舉，李永華等.FMM射線追蹤方法在地震學(xué)正演和反演中的應(yīng)用.地球物理學(xué)進(jìn)展，2010，25(4)：1197-1205.

Zhang Fengxue，Wu Qingju，Li Yonghua et al.Application of FMM ray tracing to forward and inverse problems of seismology.Progress in Geophysics，2010，25(4)：1197-1205.

[24] 鄧飛，劉超穎.三維射線快速追蹤及高斯射線束正演.石油地球物理勘探，2009，44(2)：158-165.

Deng Fei，Liu Chaoying.3-D rapid ray tracing and Gaussian ray beam forward simulation.OGP，2009，44(2)：158-165.

[25] 李曉玲，白超英，胡光義.起伏層狀TI介質(zhì)中多次波射線追蹤.石油地球物理勘探，2013，48(6)：924-931.

Li Xiaoling，Bai Chaoying and Hu Guangyi.Multiple ray tracing in an undulated layered TI media.OGP，2013，48(6)：924-931.

[26] 唐秋惠，薛霆虓，何誠.三維多步FMM射線追蹤在正態(tài)分布速度地層模型中的應(yīng)用.工程地球物理學(xué)報(bào)，2017，14(3)：346-352.

Tang Qiuhui，Xue Tingxiao，He Cheng.The application of multi-stage 3D fast marching method of ray tracing in normal distribution speed model.Chinese Journal of Engineering Geophysics，2017，14(3)：346-352.

[27] 孫章慶，孫建國，岳玉波等.復(fù)雜山地隨機(jī)介質(zhì)GMM-ULTI法射線追蹤.地球物理學(xué)報(bào)，2016，59(7)：2615-2627.

Sun Zhangqing，Sun Jianguo，Yue Yubo et al.Ray tracing using GMM-ULTI method in the random media with complex topography.Chinese Journal of Geophysics，2016，59(7)：2615-2627.