基于Chebyshev多項式的三維射線追蹤

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

1 引言

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

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

2 三維旅行時場計算方法

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

=s2(x,y,z)

(1)

(2)

(3)

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

3 射線路徑計算方法

圖1 射線傳播路徑示意圖

3.1 一維逼近函數(shù)

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

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

(4)

取其前三項作為基函數(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)

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

(11)

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

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

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

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)插值坐標變換網(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ù)值模擬

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

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

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

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

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

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

4.2 Marmousi模型

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

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

5 結(jié)束語

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

參考文獻

[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ū)多步最短路徑極值法多值多次反射波追蹤.地球物理學進展，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] 趙后越，張美根.起伏地表條件下各向異性地震波最短路徑射線追蹤.地球物理學報，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] 李培明，梅勝全，馬青坡.一種改進的雙線性插值射線追蹤方法.石油地球物理勘探，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] 韓復興，孫建國，楊昊.基于二維三次卷積插值算法的波前構(gòu)建射線追蹤.吉林大學學報(地球科學版)，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] 張婷婷，張東，邱達.三維層狀介質(zhì)中基于旅行時梯度的多次波射線追蹤.石油地球物理勘探，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] 張東，張婷婷，喬友鋒等.三維旅行時場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)格界面剖分的三維地震射線追蹤算法.地球物理學報，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ù)值解法.武漢大學學報(理學版)，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] 金昌昆，張建中.共接收面元三維折射旅行時反演方法.石油地球物理勘探，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] 張婷婷，邱達，張東.一種改進的三維旅行時梯度射線追蹤方法.石油地球物理勘探，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] 黃興國，孫建國，孫章慶等.基于復程函方程和改進的快速推進法的復旅行時計算方法.石油地球物理勘探，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] 李永博，李慶春，吳瓊等.快速行進法射線追蹤提高旅行時計算精度和效率的改進措施.石油地球物理勘探，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] 王飛，曲昕馨，劉四新等.一種新的基于多模板快速推進算法和最速下降法的射線追蹤方法.石油地球物理勘探，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] 張風雪，吳慶舉，李永華等.FMM射線追蹤方法在地震學正演和反演中的應用.地球物理學進展，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)分布速度地層模型中的應用.工程地球物理學報，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] 孫章慶，孫建國，岳玉波等.復雜山地隨機介質(zhì)GMM-ULTI法射線追蹤.地球物理學報，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.