崔寧城 黃光南* 李紅星 張 華 張曉峰 肖 昆

(①東華理工大學(xué)核資源與環(huán)境國家重點實驗室，江西南昌 330013；②東華理工大學(xué)地球物理與測控技術(shù)學(xué)院，江西南昌 330013)

0 引言

室內(nèi)實驗和野外實踐發(fā)現(xiàn)自然界的地層和巖石呈現(xiàn)出各向異性性質(zhì)[1-3]。常規(guī)各向異性介質(zhì)，在增加一定的假設(shè)條件后，可以近似成特定的各向異性介質(zhì)[4-6]，如：三斜各向異性介質(zhì)、單斜各向異性介質(zhì)、正交各向異性介質(zhì)、橢圓各向異性介質(zhì)和橫向各向同性介質(zhì)等[7-9]。

地震波旅行時計算結(jié)果在合成地震記錄、速度層析成像、震源定位、克希霍夫疊前深度偏移等方面均具有廣泛應(yīng)用價值[10-11]。當(dāng)介質(zhì)存在各向異性時，地震速度層析成像不僅可以運用于近地表速度建模，還可以用于裂縫探測、儲層預(yù)測等[12]。計算地震旅行時的主要方法可分為有限差分法、有限單元法和射線追蹤法，其中求解程函方程(Eikonal equation)有限差分法計算效率高、更易于編程實現(xiàn)。Vidale[13]和Van Trier等[14]提出了幾種有限差分旅行時計算方法，但目前較為流行的是快速推進法[15-19]和快速掃描法，這兩種流行算法排序理念不同：快速推進法使用堆排序算法，在運算中對旅行時進行排序，利用嚴格的因果關(guān)系強制執(zhí)行對旅行時做合理排序，這樣后求解的旅行時對先前的結(jié)果沒有影響；快速掃描法[20-27]使用Gauss-Seidel迭代，將計算后得到的值馬上代入迭代循環(huán)中，如果方程存在嚴格的因果關(guān)系，則可以保證在有限次迭代后達到收斂。這兩種方法可以得到相同精度的旅行時，后者具有較高的計算效率，但是選擇哪一種方法更優(yōu)，則需要根據(jù)求解的方程確定。另外還有一些其他形式的程函方程求解方法[28-32]。

Zhao[20]將快速掃描算法運用于各向同性介質(zhì)旅行時計算，Zhang等[21]和Qian等[22]對算法進行了一系列的改進，包括適應(yīng)各向異性介質(zhì)的快速掃描算法和高階快速掃描算法；Fomel 等[23]提出了一種基于因式分解程函方程的快速掃描算法，對各向同性的程函方程進行乘法因式分解，提高了旅行時的計算精度；Qian等[24]提出了基于靜態(tài)凸Hamilton-Jacobi方程的快速掃描算法，并推導(dǎo)了橢圓各向異性介質(zhì)程函方程。Luo等[25]針對橢圓各向異性介質(zhì)提出了加法和乘法因式分解形式的快速掃描算法，并對比了兩種分解形式的計算精度；Waheed等[26-27]推導(dǎo)了各向異性聲波近似方程的因式分解形式，并利用快速掃描算法求解。

相比各向同性介質(zhì)，各向異性介質(zhì)程函方程較為復(fù)雜，快速掃描法的效率相對較低。特別是當(dāng)模型較大時，提高算法的效率更具有實用價值。本文對橢圓各向異性程函方程進行因式分解，以壓制震源奇異性產(chǎn)生的旅行時誤差，提高旅行時計算的精度；根據(jù)地震波旅行時計算的因果關(guān)系，提出源點快速掃描方式，以去除不必要的計算過程，提高旅行時計算的效率；同時，結(jié)合迎風(fēng)差分格式，通過添加差分方向判定因子，使算法的計算更加合理和有效；最后，應(yīng)用數(shù)值模型驗證了算法可靠性和高效性。

1 方法原理

1.1 橢圓各向異性介質(zhì)程函方程

各向同性介質(zhì)的程函方程為

(1)

式中：s(x)=1/v(x)，為速度的倒數(shù)，即慢度；T為旅行時；Ω∈RN為N維的有界開集， 二維情況下N=2；x=(x,y)為計算點坐標(biāo)，x0=(x0,y0)為源點坐標(biāo)。在點震源情況下，橢圓各向異性介質(zhì)的程函方程為

(2)

(3)

式中：Tx和Ty分別為旅行時在x、y方向的導(dǎo)數(shù)；系數(shù)a(x)、b(x)、c(x)分別代表橫向、縱向和傾斜方向的各向異性強度，且滿足a(x)>0，b(x)>0，c2(x)-a(x)b(x)<0。各向異性對稱正定矩陣為

(4)

指示各向異性程度的系數(shù)為

(5)

式中λmax和λmin分別表示矩陣M最大、最小特征值。

1.2 乘法因式分解

常規(guī)求解程函方程計算旅行時的有限差分方法存在源奇異性問題[24]，求解因式分解形式的程函方程能得到更高精度的旅行時結(jié)果。

對式(3)中的旅行時T進行乘法分解，即

T=T0τ

(6)

式中：T0為假設(shè)模型為均勻介質(zhì)時的旅行時值；τ為擾動值，代表波場傳播過程中波前曲率的變化，是需要求取的參數(shù)。T的導(dǎo)數(shù)可表示為

(7)

式(3)的乘法因式分解形式為

2T0τ[aT0xτx-c(T0xτy+T0yτx)+bT0yτy]+

(8)

式中：T0x、T0y分別為T0在x和y方向上的偏導(dǎo)數(shù)；τx和τy同理。其中T0的定義為

T0(x,y)=

(9)

式中a0、b0、c0為源點的系數(shù)。

1.3 加法因式分解

加法因式分解比乘法因式分解更為簡單。將T進行加法分解，有

T=T0+τ

(10)

則

(11)

式(3)的加法分解形式為

2τx(aT0x-cT0y)+2τy(bT0y-cT0x)+

(12)

加法因式分解的T0(x，y)與乘法形式的一致，都由式(9)計算得到。

由于加法因式分解計算效果不如乘法分解，故本文應(yīng)用乘法因式分解方法。

2 源點快速掃描算法

2.1 常規(guī)快速掃描算法的掃描方式

常規(guī)快速掃描算法直接對計算域進行全局掃描。假設(shè)震源點位于計算區(qū)域中心(i0，j0)(紅點處)，若x方向共有nx個點、y方向共有ny個點，以T(i，j)表示計算點(i，j)的旅行時，則常規(guī)快速掃描算法的掃描方式如圖1所示。

圖1 常規(guī)快速掃描法的四種掃描方式

(a)i=(1，2，…，nx)，j=(1，2，…，ny)；(b)i=(nx，nx-1，…,1)，j=(1，2，…，ny)；(c)i=(1，2，…，nx)，j=(ny，ny-1，…,1)；(d)i=(nx，nx-1，…,1)，j=(ny，ny-1，…，1)。nx、ny分別為x、y方向網(wǎng)格點總數(shù)

以圖1a為例，掃描未抵達源點前，由于方程沒有源點信息作為初始條件，無法計算得到有效旅行時。當(dāng)掃描經(jīng)過源點后，才能得到有效的旅行時信息(圖2)。因此，圖中只有對灰色區(qū)域(i=(i0，i0+1，…，nx)，j=(j0，j0+1，…，ny))的掃描有效，對其余白色區(qū)域的掃描是無效的。

2.2 源點快速掃描算法的掃描方式

常規(guī)快速掃描算法的掃描方式存在大量無效掃描，為解決該問題，提出了一種源點快速掃描方法。將掃描的起始點移至震源點處，可以減少大量的無效掃描過程。源點快速掃描算法的掃描方式如圖3所示。

圖2 掃描方式分析

圖3 源點掃描法的四種掃描方式

(a)i=(i0，i0+1，…，nx)，j=(j0，j0+1，…，ny)；(b)i=(i0，i0-1，…，1)，j=(j0，j0+1，…，ny)；(c)i=(i0，i0+1，…，nx)，j=(j0，j0-1，…，1)；(d)i=(i0，i0-1，…，1)，j=(j0，j0-1，…，1)。i0和j0分別為源點的x、y方向網(wǎng)格序號

源點快速掃描算法相對于常規(guī)快速掃描算法省去了從端點處開始至源點的掃描部分，直接將源點作為起始點開始掃描，掃描過程更簡潔有效。

2.3 多個震源的源點快速掃描算法

當(dāng)存在多個震源時，選取所有震源中x和y方向上的坐標(biāo)序號最大值和最小值組成新的起始掃描點。設(shè)有O1(i1，j1)、O2(i2，j2)、…、On(in，jn)共n個源點，則x和y方向上的坐標(biāo)序號最大和最小值為imin=min(i1，i2，…，in)；imax=max(i1，i2，…，in)；jmin=min(j1，j2，…，jn)；jmax=max(j1，j2，…，jn)。當(dāng)存在多個震源點時，為覆蓋所有震源點，掃描的起始位置不再固定為某個點。圖4以兩個震源點O1(i1，j1)、O2(i2，j2)為例。

圖4 多震源源點掃描法的四種掃描方式

(a)i=(imin，imin+1，…，nx)，j=(jmin，jmin+1，…，ny)；(b)i=(imax，imax-1，…，1)，j=(jmin，jmin+1，…，ny)；(c)i=(imin，imin+1，…，nx)，j=(jmax，jmax-1，…，1)；(d)i=(imax，imax-1，…，1)，j=(jmax，jmax-1，…，1)

多個震源情況下的快速掃描算法對掃描速率的提升依賴于震源的分布情況。當(dāng)多個震源集中分布時，源點快速掃描算法可以節(jié)省更多的計算時間。

3 數(shù)值算法實現(xiàn)過程

3.1 差分方向判定因子

本文主要討論二維情況下的源點快速掃描法，離散化網(wǎng)格如圖5所示。

圖5 離散化網(wǎng)格

差分方向判定因子Sx和Sy的定義為

(13)

(14)

式(8)中的T0x、T0y和τx、τy離散化形式為

(15)

(16)

3.2 乘法因式分解形式的程函方程離散化

將T0x、T0y和τx、τy離散化形式代入式(8)，得到乘法因式分解形式程函方程的離散形式為

c(T0xSy(τi,j-τi,j+Sy)+T0ySx(τi,j-τi+Sx,j))+

2cSxSy(τi,j-τi+Sx,j)(τi,j-τi,j+Sy)+

(17)

上式去根號后，可以簡化為一元二次方程形式

(18)

式中

(19)

Syτi,j+Sy(cT0x-bT0y)]-

(20)

(21)

求解式(18)可得到擾動量τi，j， 再進一步利用式(6)可得到旅行時Ti，j。加法因式分解形式程函方程的離散化與乘法形式相似，在此不再贅述。

3.3 因果性判定條件

式(18)的解可能存在三種情況，分別為：無根、有兩個相同根和有兩個不同根。利用因果條件，可對解的正確性進行判斷。

若式(18)有根(不管兩個根相同與否)，由于程函方程屬于Hamilton-Jacobi方程，而在Hamilton-Jacobi系統(tǒng)下，方程的解應(yīng)同時滿足

(22)

式中：H表示Hamilton量；px和py為對應(yīng)x和y方向上慢度導(dǎo)數(shù)。橢圓各向異性介質(zhì)程函方程式(3)的因果條件為

(23)

離散化因果判定條件為

(24)

當(dāng)方程無根時，需要考慮波單獨沿x軸和y軸傳播的情況。

(25)

(26)

加法因式分解方法的計算公式分別為

(27)

(28)

3.4 迭代過程

3.4.1 初始化

(1)將計算區(qū)域的所有τ值都設(shè)定為一個較大值τmax(τmax大于計算區(qū)域最終計算出的最大τ值)。

(2)設(shè)定震源點處τ0值，乘法分解形式設(shè)定為τ0=1，加法分解形式設(shè)定為τ0=0。

(3)利用式(9)計算T0，迭代過程中T0保持不變。

(4)通過式(6)或式(10)計算旅行時T的初始值。

3.4.2 旅行時計算

(1)定義τ*為中間變量，求解式(18)，當(dāng)方程存在兩個解，分別為τ1*和τ2*。若τ1*和τ2*都符合因果條件式(24)，則τ*=min(τ1*,τ2*)；若僅τ1*符合因果條件，則τ*=τ1*；若僅τ2*符合因果條件，則τ*=τ2*。

(2)若τ1*和τ2*都不符合因果條件，則利用式(25)和式(26)計算τx*和τy*的值。若τx*和τy*滿足Ti,j=T0i,jτx*≥Ti+Sx,j,Ti,j=T0i,jτy*≥Ti,j+Sy則τ*=min(τx*,τy*)；若僅Ti,j=T0i,jτx*≥Ti+Sx,j，則τ*=τx*；若僅Ti,j=T0i,jτy*≥Ti,j+Sy，則τ*=τy*。

3.4.3 終止條件

設(shè)定趨近于0 的閾值δ，當(dāng)?shù)h(huán)滿足|Tnew-Told|≤δ時，終止迭代。

4 數(shù)值模擬

4.1 旅行時計算過程分析

將旅行時掃描過程分解，對比常規(guī)快速掃描算法和源點快速掃描算法在不同掃描階段的區(qū)別。測試模型選擇具有解析解的各向異性介質(zhì)模型[25]。

圖6中，利用常規(guī)快速掃描算法進行了全局掃描，每次掃描后，全局的旅行時都會向最終結(jié)果靠近一些，掃描完成后得到穩(wěn)定的旅行時結(jié)果。圖7中，源點快速掃描算法以源點作為起始點開始掃描，每次只掃描計算區(qū)域中有效的一部分，完成掃描過程后，同樣能得到穩(wěn)定的旅行時場結(jié)果。

圖6 常規(guī)快速掃描算法的掃描過程分解

圖7 源點快速掃描算法的掃描過程分解

4.2 旅行時精度分析

令T1為解析解旅行時，T2為常規(guī)快速掃描算法旅行時，T3為源點快速掃描算法旅行時，T2、T3與T1的對比如圖8所示。由圖可知，兩種方法數(shù)值解接近一致(紅色與藍色等值線幾乎重合)，且在在旅行時變化平緩區(qū)域的精度較高，在旅行時變化劇烈的區(qū)域則存在一定誤差。

圖8 均勻各向異性模型兩種方法計算的旅行時與理論旅行時等值線(單位：s)對比(a)及其局部放大顯示(b)

計算T3與T2的平均絕對誤差

(29)

由上式得到Terror=3.2458×10-5s，說明源點快速掃描算法與常規(guī)快速掃描算法的計算結(jié)果基本一致。

4.3 計算效率分析

4.3.1 單震源計算效率對比

應(yīng)用不同網(wǎng)格數(shù)各向異性模型比較常規(guī)快速掃描算法和源點快速掃描算法的計算效率。模型網(wǎng)格數(shù)分別設(shè)置為200×200、400×400、800×800、1600×1600、3200×3200。網(wǎng)格間距固定為1m，閾值δ=1×10-9，源點位置固定為(nx/2，ny/2)。兩種算法各測試5次，取均值后，記錄CPU平均計算時間如圖9所示。圖中源點快速掃描算法運行時間(藍色實線)整體上要小于常規(guī)快速掃描算法運行時間(紅色實線)，顯然源點快速掃描算法的計算效率要高于常規(guī)快速掃描算法。

圖9 兩種算法的運行時間比較

若固定模型網(wǎng)格數(shù)為1000×1000，其他參數(shù)設(shè)定不變，僅改變震源位置，測試震源位置分布對源點快速掃描算法計算效率的影響。源點位置設(shè)定為：O1(10，10)；O2(100，100)；O3(300，300)；O4(500，500)；O5(300，800)；O6(100，900)；O7(10，990)。不同震源位置兩種旅行時計算方法運行耗時如圖10所示?？梢钥闯觯袋c位置的改變對源點快速掃描算法的計算效率沒有影響，證明單個點源情況下，無論點源位置如何分布，源點快速掃描算法對計算效率的提升都是等效的。

圖10 不同源點位置對算法計算速度的影響分析

4.3.2 多個震源的計算效率

應(yīng)用各向同性介質(zhì)模型中，測試多個震源時，源點快速掃描算法的計算效率。

設(shè)均勻各向同性介質(zhì)模型速度為1km/s，網(wǎng)格數(shù)為1000×1000，網(wǎng)格間距為1m，閾值δ=1×10-9。設(shè)置四種不同的震源分布：分布Ⅰ，四個震源分別位于(600，600)、(400，600)、(600，400)、(400，400)，代表震源分布集中；分布Ⅱ，四個震源分別位于(700，700)、(300，700)、(700，300)、(300，300)，代表震源分布較集中；分布Ⅲ，四個震源分別位于(800，800)、(200，800)、(800，200)、(200，200)，代表震源分布較分散；分布Ⅳ，四個震源分別位于(900，900)、(100，900)、(900，100)、(100，100)，代表震源分布分散。

圖11為四個震源不同分布時源點快速掃描算法的旅行時計算結(jié)果，圖12為四種震源不同分布時源點快速掃描算法運行時間對比。隨著震源點的分布越來越分散，源點快速掃描算法的運行時間逐漸向常規(guī)快速掃描算法的運行時間靠近。震源分布越集中，源點快速掃描算法的計算效率提高越明顯。

圖11 均勻各向同性介質(zhì)四種震源不同分布時源點快速掃描算法的旅行時計算結(jié)果

圖12 均勻各向同性介質(zhì)四種震源不同分布時源點快速掃描算法與常規(guī)快速掃描算法的運行時間對比

5 結(jié)論

常規(guī)的快速掃描算法直接對計算區(qū)域進行全局掃描，運算過程中存在大量無效掃描。本文提出的源點快速掃描算法將掃描的起始點移至源點處，使算法的掃描過程更高效。通過求解因式分解形式的橢圓各向異性程函方程，利用源點快速掃描算法實現(xiàn)了各向異性介質(zhì)中的旅行時計算。數(shù)值模擬結(jié)果表明：

(1)計算參數(shù)相同的前提條件下，源點快速掃描算法和常規(guī)快速掃描算法的計算精度一致；

(2)單震源情況下，源點快速掃描算法的計算效率明顯高于常規(guī)快速掃描算法，且源點的位置不影響算法的計算效率；

(3)多震源情況下，震源分布越集中，源點快速掃描算法計算效率提升越大。

(4)源點快速掃描算法適用于各向同性和各向異性介質(zhì)的旅行時計算。