孫建國，李懿龍，孫章慶，苗 賀

吉林大學(xué)地球探測科學(xué)與技術(shù)學(xué)院，長春 130026

0 引言

模型參數(shù)化是對地震波傳播介質(zhì)的定量化描述[1]，影響著地震的正演和反演等各個方面。根據(jù)實際問題的需要，模型參數(shù)化的過程需要精心地設(shè)計安排。不恰當(dāng)?shù)膮?shù)化方式或是過于簡化的方式都會降低探測地球非均勻性地質(zhì)結(jié)構(gòu)的能力，使得一些有意義的信息丟失掉；而對于復(fù)雜的模型參數(shù)化設(shè)計又反而會使反演的不確定性劇烈增加，引入一些不真實的地下信息[2]。

地震波走時計算及射線路徑追蹤被廣泛應(yīng)用于地震層析成像[3]、地震波場正演、Kirchhoff偏移[4]、走時反演、模型計算等方面中，在理論研究和實際生產(chǎn)方面都具有重要意義，因此它們計算的有效性和精確性直接影響到上述各方面研究的效率和精度。計算地震波走時的方法主要分為兩大類：射線追蹤類和波前追蹤類。射線追蹤類基本理論基礎(chǔ)是費(fèi)馬原理，主要包括試射法[5]、彎曲法[6]、最短路徑射線追蹤法[7]、插值法[8]、辛幾何算法[9-10]等；波前追蹤類方法主要包括基于運(yùn)動學(xué)射線追蹤系統(tǒng)的波前構(gòu)建法[11-12]和基于惠更斯原理、解程函方程的有限差分方法[13-14]。本文走時計算采用有限差分法中的波前快速推進(jìn)(FMM)法。該方法由Sethian等[15]提出，具有靈活性強(qiáng)、計算速度快、對任何復(fù)雜速度場都具有無條件穩(wěn)定性的優(yōu)點。差分格式采用迎風(fēng)有限差分格式。目前射線路徑追蹤最具代表性的3種方法為：有限差分法、最短路徑法和走時線性插值(LTI)法[16]。本文采用的是最大旅行時梯度射線追蹤(MGT)法，其核心思想與LTI一樣，都是基于費(fèi)馬原理，但MGT更加精確高效。

本文的研究目的是在射線追蹤研究中加入模型參數(shù)化方法，進(jìn)而提高射線追蹤中走時計算及路徑追蹤的精度。模型參數(shù)化采用最小二乘參數(shù)化法，分別參數(shù)化速度模型和走時模型。以速度模型為例，先將離散的地震速度模型參數(shù)化，轉(zhuǎn)換成離散的參數(shù)模型，進(jìn)而得到連續(xù)地震速度模型，將其應(yīng)用于地震波走時計算。通過均勻梯度速度模型進(jìn)行地震波走時計算誤差分析，進(jìn)而研究模型參數(shù)化對地震波走時計算的影響。

1 FMM走時計算和MTG路徑追蹤

1.1 FMM基本理論

射線追蹤是基于對波動方程的高頻近似，其近似方程為程函方程：

|t|=s。

(1)

式中：t為地震波走時；s為傳播介質(zhì)的慢度。它們都是空間位置的函數(shù)。

FMM是基于求解程函方程來獲得地震波走時的，計算得到的是初至走時，因而可以看成是追蹤一個沿給定模型速度的法線方向擴(kuò)展演化的二維曲線或三維曲面，即追蹤波前。

構(gòu)建迎風(fēng)差分格式來求解程函方程：

(2)

(3)

式中：Δx和Δz分別表示橫向和縱向網(wǎng)格間距。

1.2 MGT基本理論

MGT法求取射線路徑的核心思想是，射線的傳播路徑是沿著旅行時最小的方向，即走時場的梯度方向傳播的。該方法在計算中是沿著走時場梯度的負(fù)方向逐步檢索次級源點，最終走到地震源點，即能夠確定出地震波在地下介質(zhì)中的射線追蹤路徑。

2 參數(shù)化方法

模型參數(shù)化過程是將介質(zhì)參數(shù)投影到某一函數(shù)空間，并基于該空間的基函數(shù)進(jìn)行展開[17]。本文采用的是最小二乘參數(shù)化法。

最小二乘法是一種數(shù)學(xué)優(yōu)化技術(shù)。它通過最小化誤差的平方和尋找數(shù)據(jù)的最佳函數(shù)匹配。這種求誤差平方和的方式可以避免正負(fù)誤差相抵，而且便于數(shù)學(xué)處理。由于我們所接觸的模型都是離散的點陣模型，因此使用最小二乘法進(jìn)行擬合，可以將其轉(zhuǎn)換成連續(xù)的函數(shù)形式。在函數(shù)的基礎(chǔ)上進(jìn)行網(wǎng)格加密和求導(dǎo)等操作更為簡單、合理。

本文采用的函數(shù)為多項式函數(shù)：

pn(x)=anxn+…+a2x2+a1x+a0。

(4)

速度模型的一般形式為v(xi,yi)(i=1,2,…,m)，其中(xi,yi)為速度模型中任意一點，m為橫縱方向模型點數(shù)。通過確定多項式的一系列系數(shù)aj(j=0,1,…,n)，使得各點函數(shù)值與實際值的偏差平方和M達(dá)到最小。數(shù)學(xué)形式為使

(5)

達(dá)到最小。由于M非負(fù)且為a0,a1,…,an的二次多項式，所以M(a0,a1,…,an)必有最小值。因此令

(6)

將地震模型v(xi,yi)帶入到式(6)中得到方程組：

sja0+sj+1a1+…+sj+nan=uj(j=0,1,…,n)。

(7)

(8)

由式(8)即可很方便地求取多項式系數(shù)aj，從而得到模型參數(shù)化函數(shù)。

3 方法實現(xiàn)及精度分析

本文采用最小二乘參數(shù)化法進(jìn)行模型參數(shù)化，對于走時計算及射線路徑追蹤分別對速度模型和走時模型進(jìn)行模型參數(shù)化，使得FMM的走時計算精度和MGT的射線路徑計算精度都得到明顯提高。下面將分別闡述模型參數(shù)化在走時計算及射線路徑追蹤中的作用。

3.1 速度模型參數(shù)化

FMM是一種理論基礎(chǔ)好、算法實現(xiàn)優(yōu)越性明顯的方法。采用差分算子來近似微分算子，在近似中會省略掉網(wǎng)格間距的n階無限小量，因此導(dǎo)致網(wǎng)格間距成為走時計算精度的主要影響因素。采用模型參數(shù)化后的模型進(jìn)行網(wǎng)格加密，縮小網(wǎng)格間距，從而提高計算精度。

在剛開始的研究中采用全列數(shù)據(jù)進(jìn)行多項式擬合，發(fā)現(xiàn)誤差極大；經(jīng)過分析研究后采取按列再按點的方式進(jìn)行細(xì)化，每五點采用最小二乘多項式擬合一次，得到速度模型所對應(yīng)的參數(shù)數(shù)組數(shù)據(jù)，化速度模型為參數(shù)模型，使得模型不再是離散點數(shù)據(jù)模型而是參數(shù)化連續(xù)速度模型。將此模型應(yīng)用到地震波走時計算上，可以極大地提高計算精度。

下面應(yīng)用均勻的速度梯度模型進(jìn)行精度分析，網(wǎng)格大小為40×40，網(wǎng)格間距為40 m，震源設(shè)置在(0, 20)。加密后模型網(wǎng)格為400×400，網(wǎng)格間距減小到4 m，震源保持不變。

圖1a為初始模型FMM地震波走時計算值與解析值的逼近程度，圖1b為模型參數(shù)化之后加密重新采樣的計算值與解析值的逼近程度，從中可以明顯看出，經(jīng)過模型參數(shù)化后的地震波走時值與解析值逼近效果比原模型效果好。圖2所示為模型參數(shù)化前后走時的相對誤差，表1為相對誤差與模型參數(shù)化程度的關(guān)系，從中可以很明顯地看出：相對誤差由原來的10-2數(shù)量級降到10-3數(shù)量級，模型參數(shù)化可以極大地提高地震波走時計算的精度；并且隨著模型參數(shù)化加密程度的提高，地震波走時精度也隨之提高。因此，我們的方法對于提高地震波走時計算精度是有用的。

圖1 模型參數(shù)化前(a)后(b)計算值與解析值的逼近程度Fig.1 Degree of approximation between the calculated and analytic values before (a) and after (b) model parameterization

圖2 模型參數(shù)化加密前(a)后(b)誤差對比Fig.2 Model parametric encryption before (a) and after (b) error contrast

Table1Comparisonoftraveltimerelativeerrorsbeforeandaftermodelparametricencryption

網(wǎng)格大小走時相對誤差40×40(原模型)0．026980×80(加密1倍)0．0188400×400(加密10倍)0．0083

3.2 走時模型參數(shù)化

對于MGT法，由于初始模型都是離散的，因此導(dǎo)致由離散走時場而獲得的旅行時梯度存在很大的不穩(wěn)定性，路徑計算的誤差較大。在此情況下，應(yīng)用本文所提出的模型參數(shù)化方法，將對其起到根本性的作用，提高計算精度。將離散走時場轉(zhuǎn)化成連續(xù)走時場，進(jìn)而可以求取任意走時點的時間梯度，使得時間梯度的求取更為合理、穩(wěn)定，提高射線路徑追蹤的精度。本文對走時場所采用Chebyshev多項式進(jìn)行最小二乘參數(shù)化。在二維走時場上采取九點擬合，獲得4個Chebyshev多項式系數(shù)，進(jìn)而將離散的走時場轉(zhuǎn)化成連續(xù)的時間場函數(shù)。依然應(yīng)用均勻速度梯度模型進(jìn)行精度分析，并將其與離散模型進(jìn)行對比，離散模型的時間場梯度由中心差分法求得。模型網(wǎng)格大小為100×100，網(wǎng)格間距為4 m，震源設(shè)置在(0, 200)處。

圖3 離散模型(a)與參數(shù)化模型(b)射線路徑對比Fig.3 Ray path contrast between discrete model (a) and parametric model (b)

圖4 參數(shù)化前后射線路徑相對誤差對比圖Fig.4 Comparison of relative errors of ray paths before and after parameterization

圖3為模型參數(shù)化前后梯度速度模型的射線追蹤路徑圖，圖4為模型參數(shù)化與離散模型的射線路徑精度對比圖，從中可以看出模型參數(shù)化對于追蹤射線路徑的精度有明顯改善，平均精度提高了14.33%。

4 計算實例

下面給出兩組經(jīng)典地震速度模型——Marmousi模型和Sigsbee 2A模型的走時計算及射線追蹤實例，以驗證模型參數(shù)化法在地震波走時計算及在射線路徑追蹤計算方面的確切可行性。

圖5為Marmousi模型參數(shù)化前后的地震波走時對比圖，該速度模型包含有極為復(fù)雜的地質(zhì)構(gòu)造，為50×50，網(wǎng)格間距為40 m，震源取在(1, 35)，參數(shù)化后重新加密采樣網(wǎng)格大小為500×500。對比參數(shù)化前后的計算結(jié)果可以看出，在模型較深的復(fù)雜介質(zhì)中，模型參數(shù)化處理之后的等時線在局部明顯被細(xì)化，地震波傳播顯得更為合理。圖6為Sigsbee 2A模型，網(wǎng)格大小為600×600，網(wǎng)格間距為4 m,參數(shù)處理之后為6 000×6 000，震源取在(1,1)。對比圖6中模型參數(shù)化前后的計算結(jié)果依然可以看出，在地下異常高速體中地震波走時的等時線在模型參數(shù)化后顯得更為細(xì)致。圖5、圖6的計算結(jié)果表明，本文提出的模型參數(shù)化法切實可行，適用于任意的復(fù)雜速度模型。

圖7為走時模型參數(shù)化在射線路徑追蹤中的應(yīng)用。圖7a為Marmousi模型，網(wǎng)格大小為500×500，網(wǎng)格間距為4 m；圖7b為Sigsbee 2A模型，網(wǎng)格大小為600×600，網(wǎng)格間距為4 m。圖7中模型的有些區(qū)域沒有射線穿過，那是由于本文所計算的是最小走時射線路徑，所以射線會盡量選擇高速層穿越，進(jìn)而達(dá)到走時最小的目的。圖7從射線追蹤計算路徑的角度證明了本文的模型參數(shù)化方法對于地下復(fù)雜介質(zhì)的有效性。

圖5 二維Marmousi模型參數(shù)化前(a)后(b)等時線對比Fig.5 Two dimensional Marmousi model parameterization before (a)and after (b) isochronous correlation

圖6 二維Sigsbee 2A模型參數(shù)化前(a)后(b)等時線對比Fig.6 Two dimensional Sigsbee2A model parameterization before (a)after (b) isochronous correlation

a. Marmousi模型；b. Sigsbee 2A模型。圖7 走時模型參數(shù)化在追蹤射線路徑中的應(yīng)用Fig.7 Application of traveltime model parameterization in tracking ray paths

5 結(jié)論

1)模型參數(shù)化是對地震波傳播介質(zhì)的定量化描述，影響著地震處理的各個過程。對速度模型進(jìn)行恰當(dāng)?shù)膮?shù)化處理，可以提高地震波走時計算的精度。通過對走時模型的參數(shù)化處理，使得求取走時場的梯度更為合理、穩(wěn)定，相對于離散模型，最終獲得的射線路徑精度更高。

2)模型參數(shù)化具有普遍適用性，且無論是速度模型還是走時模型等都可以很好地參數(shù)化。并且，本文方法對于任何復(fù)雜走時計算均適用，且具有很好的穩(wěn)定性。

[1] 李貞,郭飚,劉啟元，等. 地震層析成像中模型參數(shù)化研究進(jìn)展[J]. 地球物理學(xué)進(jìn)展,2015, 30(4): 1616-1624.

Li Zhen, Guo Biao, Liu Qiyuan, et al. Parameterization in Seismic Tomography[J]. Progress in Geophysics, 2015， 30(4): 1616-1624.

[2] Thurber C H, Aki K. Three-Dimensional Seismic Ima-ging[J]. Annual Review of Earth and Planetary Sciences, 1987，15(1): 115-139.

[3] Dziewonski A M, Woodhouse J H. Global Images of the Earth’s Interior[J]. Science, 1987, 236: 37-48.

[4] Sena A G, Toksoz M N. Kirchhoff Migration and Ve-locity Analysis Forconverted and non Converted Waves in Anisotropic media[J]. Geophysics, 1993, 58: 265-276.

[5]Julian B R, Gubbins D. Three Dimensional Seismic Ray Tracing[J]. Geophysics, 1977, 43: 95-119.

[6] Cerveny V. Seismic Ray Theory[M]. Cambridge: Cam-bridge University Press, 2001.

[7] Moser T J. Shortest Path Calculation of Seismic Rays[J]. Geophysics, 1991, 56: 59-67.

[8] Asakawa E, Kawanaka T. Seismic Ray Tracing Using Linear Traveltime Interpolation[J]. Geophysical Prospecting, 1993, 41(1): 99-111.

[9] 陳景波,秦孟兆.辛幾何算法在射線追蹤中的應(yīng)用[J].數(shù)值計算與計算機(jī)應(yīng)用,2000, 43(4): 254-265.

Chen Jingbo, Qin Mengzhao. Ray Ttacing by Symplectic Algorithm[J]. Journal on Numerical Methods and Computer Applications, 2000, 43(4): 254-265.

[10] 羅明秋，劉洪，李幼銘.地震波傳播的哈密頓表述及辛幾何算法[J].地球物理學(xué)報，2001, 44(1)： 120-127.

Luo Mingqiu, Liu Hong, Li Youming. Hamiltionian Deacription and Symplectic Method of Seismic Wave Propagtion[J]. Chinese Journal of Geophysics, 2001, 44(1): 120-127.

[11] Vinje V, Iversen E G, Joystdal H. Traveltime and Amplitude Estimation Using Wavefront Construction[J]. Geophysics, 1993, 58(8): 1157-1166.

[12] 石秀林，孫建國，孫輝，等. 基于波前構(gòu)建法的時間域深度偏移：Delta波包途徑[J].吉林大學(xué)學(xué)報(地球科學(xué)版)，2016，46(6)：1847-1854.

Shi Xiulin, Sun Jianguo, Sun Hui, et al. Depth Migration in Time Domain Using Wavefront Construction: Delta Packet Approach[J]. Journal of Jilin University (Earth Science Edition), 2016, 46(6): 1847-1854.

[13] 孫章慶. 復(fù)雜地表條件下地震波走時計算方法研究[D].長春：吉林大學(xué)，2008.

Sun Zhangqing. Study on the Computation Method of Seismic Traveltimes Including Surface Topography[D]. Changchun: Jilin University, 2008.

[14] Vidle J. Finite-Difference Calcaulation of Traveltimes[J]. Bull Seism Soc Am, 1988, 78(6): 2062-2076.

[15] Sethian J A, Popovici A M. 3-D Traveltime Com-putation Using the Fast Marching Method[J]. Geophysics, 1999, 64(2): 516-523.

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

Zhang Dong, Zhang Tingting, Qiao Youfeng, et al. A 3-D Ray Method Based on B-Spline Traveltime Interpolation[J]. Oil Geophysical Prospecting, 2013, 48(4): 559-566.

[17] Nolet G. A Breviary of Seismic Tomography: Ima-ging the Interior of the Earth and Sun[M]. Cambriage: Cambriage University Press, 2008.