求解橢圓各向異性介質(zhì)程函方程的源點(diǎn)快速掃描算法

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

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

0 引言

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

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

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

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

1 方法原理

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

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

(1)

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

(2)

(3)

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

(4)

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

(5)

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

1.2 乘法因式分解

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

對(duì)式(3)中的旅行時(shí)T進(jìn)行乘法分解，即

T=T0τ

(6)

式中：T0為假設(shè)模型為均勻介質(zhì)時(shí)的旅行時(shí)值；τ為擾動(dòng)值，代表波場(chǎng)傳播過程中波前曲率的變化，是需要求取的參數(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為源點(diǎn)的系數(shù)。

1.3 加法因式分解

加法因式分解比乘法因式分解更為簡(jiǎn)單。將T進(jìn)行加法分解，有

T=T0+τ

(10)

則

(11)

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

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

(12)

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

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

2 源點(diǎn)快速掃描算法

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

常規(guī)快速掃描算法直接對(duì)計(jì)算域進(jìn)行全局掃描。假設(shè)震源點(diǎn)位于計(jì)算區(qū)域中心(i0，j0)(紅點(diǎn)處)，若x方向共有nx個(gè)點(diǎn)、y方向共有ny個(gè)點(diǎn)，以T(i，j)表示計(jì)算點(diǎn)(i，j)的旅行時(shí)，則常規(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)格點(diǎn)總數(shù)

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

2.2 源點(diǎn)快速掃描算法的掃描方式

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

圖2 掃描方式分析

圖3 源點(diǎn)掃描法的四種掃描方式

(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分別為源點(diǎn)的x、y方向網(wǎng)格序號(hào)

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

2.3 多個(gè)震源的源點(diǎn)快速掃描算法

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

圖4 多震源源點(diǎn)掃描法的四種掃描方式

(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)

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

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

3.1 差分方向判定因子

本文主要討論二維情況下的源點(diǎn)快速掃描法，離散化網(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)

上式去根號(hào)后，可以簡(jiǎn)化為一元二次方程形式

(18)

式中

(19)

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

(20)

(21)

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

3.3 因果性判定條件

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

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

(22)

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

(23)

離散化因果判定條件為

(24)

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

(25)

(26)

加法因式分解方法的計(jì)算公式分別為

(27)

(28)

3.4 迭代過程

3.4.1 初始化

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

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

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

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

3.4.2 旅行時(shí)計(jì)算

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

(2)若τ1*和τ2*都不符合因果條件，則利用式(25)和式(26)計(jì)算τ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|≤δ時(shí)，終止迭代。

4 數(shù)值模擬

4.1 旅行時(shí)計(jì)算過程分析

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

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

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

圖7 源點(diǎn)快速掃描算法的掃描過程分解

4.2 旅行時(shí)精度分析

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

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

計(jì)算T3與T2的平均絕對(duì)誤差

(29)

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

4.3 計(jì)算效率分析

4.3.1 單震源計(jì)算效率對(duì)比

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

圖9 兩種算法的運(yùn)行時(shí)間比較

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

圖10 不同源點(diǎn)位置對(duì)算法計(jì)算速度的影響分析

4.3.2 多個(gè)震源的計(jì)算效率

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

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

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

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

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

5 結(jié)論

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

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

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

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

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