崔寧城,黃光南,李紅星,肖 昆

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

地震旅行時(shí)層析成像技術(shù)可用于預(yù)測(cè)地下地質(zhì)體的分布情況,該技術(shù)已廣泛應(yīng)用于石油勘探、礦產(chǎn)勘查和工程勘探等領(lǐng)域[1-6]。旅行時(shí)正演模擬是地震旅行時(shí)層析成像技術(shù)的重要組成部分,正演模擬的精度和計(jì)算效率影響著最終的反演效果[7]。射線追蹤類方法、程函方程的有限元解法和有限差分解法是常見(jiàn)的3類地震旅行時(shí)正演模擬方法[8-13],其中有限差分法相對(duì)有限元法更易于實(shí)現(xiàn),且能計(jì)算出目標(biāo)區(qū)域內(nèi)所有網(wǎng)格節(jié)點(diǎn)的旅行時(shí),不存在射線追蹤類方法無(wú)法到達(dá)的陰影區(qū)域,因此該方法在地震波旅行時(shí)計(jì)算中得到廣泛應(yīng)用。

1988年VIDALE[14]提出了盒式擴(kuò)展法,利用有限差分的盒式擴(kuò)展形式計(jì)算程函方程的數(shù)值解,首次快速高效地計(jì)算了復(fù)雜速度模型的地震波旅行時(shí)場(chǎng),但這種方法的擴(kuò)展規(guī)律不符合程函方程的因果關(guān)系條件,因此穩(wěn)定性較差。1991年VAN等[15]改進(jìn)了盒式擴(kuò)展法的有限差分格式,并結(jié)合迎風(fēng)差分格式提高了地震波旅行時(shí)計(jì)算方法的穩(wěn)定性。1992年QIN等[16]以波前面的波前點(diǎn)作為擴(kuò)展要素,提出了波前擴(kuò)展方法,該方法的求解方式存在缺陷,計(jì)算復(fù)雜速度模型時(shí)穩(wěn)定性較差。后經(jīng)不斷改進(jìn),有限差分計(jì)算方法求解程函方程的結(jié)果穩(wěn)定性有了明顯改善。1996年SETHIAN[17]結(jié)合迎風(fēng)有限差分方法和窄帶技術(shù)提出了快速推進(jìn)算法,2004年ZHAO[18]綜合利用程函方程的因果關(guān)系條件和Gauss-Seidel迭代法提出了快速掃描算法。上述兩種有限差分旅行時(shí)計(jì)算方法都具有計(jì)算精度高、計(jì)算速度快且無(wú)條件穩(wěn)定的特點(diǎn),是現(xiàn)今較為成熟的旅行時(shí)計(jì)算方法[19-21]。

快速推進(jìn)算法在考慮因果關(guān)系條件的前提下,采用了窄帶技術(shù)[22],該技術(shù)的思路是將計(jì)算區(qū)域內(nèi)的網(wǎng)格節(jié)點(diǎn)分為接受點(diǎn)(旅行時(shí)已知的網(wǎng)格點(diǎn))、窄帶點(diǎn)(旅行時(shí)等待計(jì)算的網(wǎng)格點(diǎn))和遠(yuǎn)離點(diǎn)(旅行時(shí)尚未計(jì)算的網(wǎng)格點(diǎn))。窄帶技術(shù)先將遠(yuǎn)離點(diǎn)逐步轉(zhuǎn)化為窄帶點(diǎn),再進(jìn)一步轉(zhuǎn)換為接受點(diǎn),最終使未知點(diǎn)全部轉(zhuǎn)換為已知點(diǎn)。在轉(zhuǎn)換過(guò)程中,需要頻繁查找窄帶點(diǎn)中的最小旅行時(shí)點(diǎn),并將其代入后續(xù)的迭代運(yùn)算。

高效的排序方法能極大地提高快速推進(jìn)算法的計(jì)算效率,但相關(guān)的研究成果并不多見(jiàn)。常規(guī)的快速推進(jìn)算法一般采用計(jì)算機(jī)技術(shù)中排序能力強(qiáng)的方法對(duì)窄帶點(diǎn)進(jìn)行排序。堆排序方法如二叉樹和三叉樹排序方法是快速推進(jìn)算法中常用的排序方法,能有效提高算法的計(jì)算效率[23-24]。但堆排序方法仍然存在一定的不足:堆排序法屬于非穩(wěn)定類排序方法,計(jì)算精度不如穩(wěn)定類排序方法;堆排序法的實(shí)現(xiàn)步驟較為繁瑣,需要額外的運(yùn)算成本。針對(duì)上述問(wèn)題,YATZIV等[25]提出了一種基于非嚴(yán)格優(yōu)先列隊(duì)的快速推進(jìn)算法,以提高旅行時(shí)算法的穩(wěn)定性。在此基礎(chǔ)上,楊昊[26]對(duì)該算法進(jìn)行了改進(jìn),提出了分段最小排序的思路,進(jìn)一步提升了排序方法的穩(wěn)定性。雖然這兩種改進(jìn)后的排序方法都屬于穩(wěn)定類排序法,但二者只在排序?qū)用嬗懻摿怂惴ǖ目尚行?未深入發(fā)掘旅行時(shí)場(chǎng)本身所隱含的有序性。本文通過(guò)分析程函方程的因果關(guān)系條件,從原理上論證了旅行時(shí)場(chǎng)的有序性,并以此為基礎(chǔ),設(shè)計(jì)了適用于快速推進(jìn)算法的插入排序方法。

本文首先介紹了利用有限差分格式求解程函方程的基本過(guò)程;然后從程函方程的因果關(guān)系條件出發(fā),論證了旅行時(shí)場(chǎng)的有序性;接著介紹了窄帶技術(shù)的具體實(shí)現(xiàn)過(guò)程,設(shè)計(jì)了相應(yīng)的插入排序方法對(duì)窄帶點(diǎn)進(jìn)行排序,進(jìn)而實(shí)現(xiàn)了插入排序快速推進(jìn)算法的應(yīng)用;然后將本文提出的插入排序快速推進(jìn)算法與常規(guī)的三叉樹堆排序快速推進(jìn)算法以及作為穩(wěn)定對(duì)照組的常規(guī)快速掃描方法分別進(jìn)行測(cè)試并比較,展示了插入排序快速推進(jìn)算法的優(yōu)勢(shì)和不足;最后,詳細(xì)分析和討論了影響插入排序快速推進(jìn)算法計(jì)算效率的因素,以及算法計(jì)算效率與旅行時(shí)場(chǎng)的有序程度之間的相關(guān)性,拓展了插入排序方法的適用領(lǐng)域。

1 方法原理

1.1 程函方程

各向同性介質(zhì)的二維程函方程為[18-21]:

(1)

式中:S為慢度;(x,y)為計(jì)算區(qū)域內(nèi)任意一點(diǎn)的空間坐標(biāo);T為地震波旅行時(shí)。

利用迎風(fēng)有限差分格式將程函方程離散,可以得到網(wǎng)格節(jié)點(diǎn)(i,j)的旅行時(shí)計(jì)算公式:

當(dāng)|Txmin-Tymin|≥Sh時(shí):

Ti,j=min(Txmin,Tymin)+Sh

(2)

當(dāng)|Txmin-Tymin|

(3)

式中:下標(biāo)i和j分別對(duì)應(yīng)x和y軸方向上的網(wǎng)格節(jié)點(diǎn)序號(hào);h代表網(wǎng)格間距(假設(shè)網(wǎng)格為等間距網(wǎng)格);Txmin=min(Ti+1,j,Ti-1,j),Tymin=min(Ti,j+1,Ti,j-1),min表示取括號(hào)中兩個(gè)值的較小值。

1.2 因果條件

程函方程屬于哈密頓-雅克比方程(Hamilton-Jacobi equation),二維哈密頓-雅克比方程因果關(guān)系條件為[27]:

(4)

式中:H為哈密頓算子;px,py分別為x軸和y軸的廣義坐標(biāo)。

二維地震波程函方程對(duì)應(yīng)的哈密頓算子可以表示為:

(5)

其對(duì)應(yīng)的因果條件為:

(6)

離散形式表示為:

(7)

式中:Ti,j由(2)式、(3)式計(jì)算得到。(7)式的物理意義如圖1所示,圖中波前面(虛線)與波場(chǎng)傳播方向垂直,波前面到達(dá)每個(gè)網(wǎng)格節(jié)點(diǎn)的先后順序決定了該點(diǎn)的旅行時(shí)。

圖1 不同的波場(chǎng)傳播方向和波前面示意

圖1展示了旅行時(shí)計(jì)算的局部規(guī)律,由程函方程因果關(guān)系條件可知,隨著計(jì)算的進(jìn)行,旅行時(shí)有逐漸增大的趨勢(shì)。均勻速度模型(速度為1km/s)的旅行時(shí)場(chǎng)展示了旅行時(shí)場(chǎng)的有序性。其整體規(guī)律表現(xiàn)為:網(wǎng)格節(jié)點(diǎn)上的旅行時(shí)值隨網(wǎng)格節(jié)點(diǎn)位置與源點(diǎn)間的距離增大而增大(圖2)。

選取快速推進(jìn)算法的排序方法時(shí),應(yīng)當(dāng)考慮以下3個(gè)方面的因素:①排序方法的排序能力強(qiáng)弱;②實(shí)現(xiàn)排序方法所需的運(yùn)算成本大小;③排序目標(biāo)是否存在特殊規(guī)律,可否簡(jiǎn)化排序過(guò)程。

常規(guī)快速推進(jìn)算法主要考慮排序方法的排序能力強(qiáng)弱,因此認(rèn)為“堆排序方法”的排序效果較好[23-24]。實(shí)際上,雖然堆排序方法的排序能力強(qiáng),但其運(yùn)算成本高。因果關(guān)系條件表明,作為排序?qū)ο蟮穆眯袝r(shí)場(chǎng)存在由小到大的分布規(guī)律。根據(jù)這一規(guī)律,本文采用運(yùn)算成本較小的插入排序方法替換堆排序方法,以取得更好的計(jì)算效果。

圖2 均勻速度模型的旅行時(shí)場(chǎng)(速度為1km/s)

對(duì)比插入排序方法和堆排序方法的時(shí)間和空間復(fù)雜度[24,27],結(jié)果見(jiàn)表1,O為不同數(shù)值方法的復(fù)雜度函數(shù),n代表模型網(wǎng)格節(jié)點(diǎn)的數(shù)量。插入排序方法在時(shí)間復(fù)雜度(最佳)上優(yōu)于堆排序方法。這說(shuō)明在適當(dāng)?shù)臈l件下,插入排序方法的效率高于堆排序方法。

表1 排序算法的時(shí)間和空間復(fù)雜度

2 方法實(shí)現(xiàn)

2.1 窄帶技術(shù)

快速推進(jìn)算法采用窄帶技術(shù)將計(jì)算區(qū)域內(nèi)的所有網(wǎng)格節(jié)點(diǎn)分為接受點(diǎn)、窄帶點(diǎn)和遠(yuǎn)離點(diǎn)(圖3)。接受點(diǎn)代表計(jì)算完成的已知旅行時(shí)點(diǎn),在后續(xù)計(jì)算過(guò)程中不再變化。遠(yuǎn)離點(diǎn)代表尚未計(jì)算的網(wǎng)格節(jié)點(diǎn),等待轉(zhuǎn)換為窄帶點(diǎn)接受計(jì)算。窄帶點(diǎn)處于接受點(diǎn)與遠(yuǎn)離點(diǎn)之間,是準(zhǔn)備接受計(jì)算的點(diǎn),類似于地震波傳播過(guò)程中的波前面。地震波場(chǎng)的傳播方向由上風(fēng)區(qū)向下風(fēng)區(qū)擴(kuò)散傳播。在這一傳播過(guò)程中,首先將窄帶點(diǎn)轉(zhuǎn)化為接受點(diǎn),再將與窄帶點(diǎn)相鄰的遠(yuǎn)離點(diǎn)轉(zhuǎn)化為窄帶點(diǎn),使地震波場(chǎng)逐漸向外擴(kuò)展,從而達(dá)到將遠(yuǎn)離點(diǎn)全部轉(zhuǎn)換為接受點(diǎn)的目的。

2.2 插入排序快速推進(jìn)算法

在快速推進(jìn)算法的迭代過(guò)程中,需要頻繁地查找窄帶點(diǎn)中的最小旅行時(shí)點(diǎn)并將其代入迭代計(jì)算中。常規(guī)快速推進(jìn)算法通常采用堆排序的方法,將窄帶點(diǎn)中的數(shù)據(jù)按二叉樹或三叉樹形式進(jìn)行排序。當(dāng)排序的數(shù)據(jù)量較大時(shí),該方法能有效提升排序效率。

圖3 窄帶技術(shù)示意

堆排序方法屬于非穩(wěn)定排序方法。非穩(wěn)定排序方法定義如下:若序列中存在數(shù)據(jù)點(diǎn)a=b,且原序列中a排列在b之前,那么排序之后a可能會(huì)出現(xiàn)在b之后[27]。這種非穩(wěn)定性對(duì)排序結(jié)果造成了一定程度的不利影響。插入排序法屬于穩(wěn)定的排序方法,作為線性類型的排序方法,其實(shí)現(xiàn)過(guò)程簡(jiǎn)單穩(wěn)定?？紤]到旅行時(shí)場(chǎng)具有有序性,理論上需要排序的窄帶點(diǎn)數(shù)組的數(shù)據(jù)量不會(huì)太大,本文設(shè)計(jì)并實(shí)現(xiàn)了插入排序快速推進(jìn)算法。

插入排序快速推進(jìn)算法的具體流程如下。

1) 初始化

首先令源點(diǎn)處的旅行時(shí)T0=0,并將其設(shè)置為接受點(diǎn);然后計(jì)算出源點(diǎn)的相鄰點(diǎn)旅行時(shí),并將相鄰點(diǎn)設(shè)為窄帶點(diǎn);再將窄帶點(diǎn)數(shù)組中的元素按旅行時(shí)值由小到大排序;最后將其它網(wǎng)格節(jié)點(diǎn)設(shè)置為某個(gè)大于全局可能旅行時(shí)值的較大值,將其屬性設(shè)為遠(yuǎn)離點(diǎn)。

2) 迭代

在迭代過(guò)程中,若窄帶點(diǎn)數(shù)組(T1,T2,…,Tk,…,Tn)中共有n個(gè)元素,則以k=1,2,3,…,n的順序從窄帶點(diǎn)數(shù)組中選取Tk點(diǎn)代入循環(huán)。此時(shí)窄帶點(diǎn)數(shù)組中的元素由小到大排列,因此每次選取的Tk點(diǎn)必定為窄帶點(diǎn)數(shù)組中的最小值點(diǎn)。

選定Tk點(diǎn)后,利用(2)式和(3)式更新Tk點(diǎn)的旅行時(shí)值,并將Tk點(diǎn)屬性從窄帶點(diǎn)轉(zhuǎn)換為接受點(diǎn),更新與Tk點(diǎn)相鄰的遠(yuǎn)離點(diǎn)的旅行時(shí)值得到Tnew,將這些點(diǎn)的屬性轉(zhuǎn)換為窄帶點(diǎn),并將其插入窄帶點(diǎn)數(shù)組中。其它與Tk點(diǎn)相鄰的非遠(yuǎn)離點(diǎn)不作處理。

將更新得到的旅行時(shí)Tnew插入至窄帶點(diǎn)數(shù)組的方法為:首先從數(shù)組末尾向前,查找數(shù)組中剛好小于旅行時(shí)Tnew的點(diǎn),在數(shù)組的第m點(diǎn)處,若Tm>Tnew,則將窄帶點(diǎn)數(shù)組中的Tm點(diǎn)向后移動(dòng)一位,即Tm+1=Tm;接著令m=m-1,對(duì)比下一個(gè)點(diǎn)的Tm值,直到m點(diǎn)的旅行時(shí)值滿足Tm

當(dāng)所有與Tk相鄰的遠(yuǎn)離點(diǎn)都完成了插入操作后,令k=k+1,選取窄帶點(diǎn)數(shù)組中的下一個(gè)最小值點(diǎn)Tk代入循環(huán),直到滿足迭代終止條件。

3) 截止條件

隨著波前面的擴(kuò)展和推進(jìn),當(dāng)計(jì)算區(qū)域內(nèi)的所有網(wǎng)格節(jié)點(diǎn)都變?yōu)榻邮茳c(diǎn)時(shí),算法終止迭代。

3 數(shù)值模擬

本文選用3種有限差分旅行時(shí)算法(插入排序快速推進(jìn)算法(insert-sorting fast marching method,簡(jiǎn)稱FMM 1)、三叉樹堆排序快速推進(jìn)算法(ternary tree fast marching method,簡(jiǎn)稱FMM 2)和常規(guī)快速掃描算法(fast sweeping method,FSM))參與測(cè)試。因常規(guī)快速掃描算法經(jīng)多次迭代掃描可以保證結(jié)果達(dá)到穩(wěn)定收斂狀態(tài),故將其作為一個(gè)穩(wěn)定的參考解。將快速掃描算法的迭代終止的閾值設(shè)置為δ=10-9,即相鄰兩次迭代間的旅行時(shí)場(chǎng)誤差小于10-9時(shí),終止迭代。

測(cè)試的3個(gè)速度模型分別為:光滑非均勻速度模型、鹽丘速度模型和Marmousi速度模型。其中光滑非均勻速度模型具有解析解,可直接用于測(cè)試算法的精度;其它兩個(gè)復(fù)雜速度模型可用于測(cè)試算法的穩(wěn)定性。

為對(duì)比3種有限差分算法的計(jì)算效率,我們將每個(gè)速度模型都采樣成6種網(wǎng)格離散模型,網(wǎng)格節(jié)點(diǎn)數(shù)分別為:201×201,401×401,801×801,1601×1601,3201×3201,6401×6401。對(duì)不同網(wǎng)格的模型測(cè)試3種算法,每種算法測(cè)試5次,并統(tǒng)計(jì)CPU時(shí)間,取CPU時(shí)間的均值作為最終的運(yùn)行時(shí)間。

3.1 光滑非均勻速度模型

參考等梯度速度模型,令模型速度值隨計(jì)算點(diǎn)與源點(diǎn)之間的距離呈等梯度變化,得到光滑非均勻速度模型[28]。該速度模型S(x)可表示為:

(8)

式中:S0為初始源點(diǎn)的慢度;G0為梯度向量;x表示計(jì)算點(diǎn)的空間位置;x0表示震源點(diǎn)的空間位置。

該速度模型的解析解為:

(9)

(10)

式中:a代表任意輸入變量。

令S0=1,G0=(0.2,0.2),計(jì)算區(qū)域?yàn)?0km×10km,震源點(diǎn)為(5km,0),圖4為網(wǎng)格節(jié)點(diǎn)數(shù)為401×401時(shí)光滑非均勻速度模型的旅行時(shí)場(chǎng)。

圖4中紅色、藍(lán)色和黑色等值線分別代表插入排序快速推進(jìn)算法(FMM 1)、三叉樹堆排序快速推進(jìn)算法(FMM 2)和常規(guī)快速掃描算法(FSM)計(jì)算得到的旅行時(shí),紫色等值線代表(9)式的解析解。圖4中4種顏色的等值線幾乎完全重合,難以分辨采用3種不同算法計(jì)算得到的旅行時(shí)場(chǎng)之間的差別。

圖4 光滑非均勻速度模型的旅行時(shí)場(chǎng)(網(wǎng)格節(jié)點(diǎn)數(shù)為401×401)

為了更詳細(xì)地分析3種算法的計(jì)算精度,根據(jù)如下公式:

E=|Tana-Tnum|

(11)

計(jì)算解析解Tana和數(shù)值解Tnum之間的誤差E。

為方便比較,將圖5中的色標(biāo)(數(shù)值)范圍固定為0～0.04s,經(jīng)逐點(diǎn)驗(yàn)證可知,采用FMM 1和FSM得到的數(shù)值解與解析解的旅行時(shí)誤差計(jì)算結(jié)果完全一致,因此圖5a可同時(shí)表示FMM 1、FSM的數(shù)值解與解析解之間的旅行時(shí)誤差。由于直角坐標(biāo)系下的常規(guī)有限差分旅行時(shí)算法存在源點(diǎn)奇異性問(wèn)題,因此算法的計(jì)算誤差較大,且誤差區(qū)域主要分布在假設(shè)波前與實(shí)際波前不匹配的區(qū)域(相對(duì)源點(diǎn)的傾斜方向上)[29]。對(duì)比圖5a和圖5b可以看出,圖5a的旅行時(shí)誤差明顯更小,且分布更為規(guī)律和平滑;圖5b的旅行時(shí)誤差分布一定程度上呈現(xiàn)鋸齒狀且較為散亂,這證明采用FMM 2得到的旅行時(shí)計(jì)算結(jié)果收斂程度不夠。

圖5 光滑非均勻速度模型中采用不同算法得到的數(shù)值解與解析解的旅行時(shí)誤差a FMM 1、FSM的數(shù)值解與解析解的旅行時(shí)誤差; b FMM 2的數(shù)值解與解析解的旅行時(shí)誤差

固定計(jì)算區(qū)域不變,測(cè)試不同網(wǎng)格大小的光滑非均勻速度模型下各個(gè)算法的L2誤差和CPU時(shí)間。圖6a 展示了FMM 1和FMM 2的L2誤差(均方差)統(tǒng)計(jì)結(jié)果[28]?？梢钥闯稣w上采用FMM 1得到的L2誤差明顯低于采用FMM 2得到的L2誤差,并且隨著網(wǎng)格節(jié)點(diǎn)數(shù)的增加,采用FMM 1得到的L2誤差曲線下降更快。

比較圖6b中3種算法的CPU時(shí)間可知,在網(wǎng)格節(jié)點(diǎn)數(shù)小于801×801的情況下,采用FMM 1所需的CPU時(shí)間約為其它兩種方法耗時(shí)的一半,說(shuō)明在相同條件下FMM 1的運(yùn)算效率明顯更高。隨著網(wǎng)格節(jié)點(diǎn)數(shù)的增加,FMM 1的運(yùn)算效率逐漸降低。這是由于常規(guī)旅行時(shí)計(jì)算方法存在源點(diǎn)奇異性問(wèn)題,網(wǎng)格點(diǎn)數(shù)增加導(dǎo)致波前窄帶點(diǎn)的數(shù)量增多且無(wú)序性增強(qiáng)。

圖6 不同網(wǎng)格大小的光滑非均勻速度模型下各個(gè)算法的L2誤差和CPU時(shí)間a FMM 1和FMM 2的L2誤差; b 3種算法的CPU時(shí)間

3.2 鹽丘速度模型

通過(guò)增加鹽丘速度模型的復(fù)雜度,來(lái)測(cè)試復(fù)雜模型下插入排序快速推進(jìn)算法的穩(wěn)定性和計(jì)算效率。鹽丘速度模型的計(jì)算區(qū)域?yàn)?.50km×4.99km,震源點(diǎn)位置為(3.25km,0),圖7a為鹽丘速度模型網(wǎng)格節(jié)點(diǎn)數(shù)為651×500時(shí)的測(cè)試結(jié)果。

鹽丘速度模型的旅行時(shí)場(chǎng)如圖7a所示,采用FMM 1、FMM 2和FSM計(jì)算得到的旅行時(shí)一致性良好(三者的等值線幾乎完全重合),在一定程度上證明插入排序快速推進(jìn)算法可在復(fù)雜模型中保持良好的穩(wěn)定性。

由于模型的復(fù)雜度增加,采用FSM時(shí),需要進(jìn)行6次迭代才能滿足迭代終止條件,這一方面保證了旅行時(shí)計(jì)算結(jié)果足夠穩(wěn)定和精確,另一方面也導(dǎo)致在該模型下采用FSM所需的CPU時(shí)間比采用FMM 2所需的CPU時(shí)間更長(zhǎng)(圖7b)。此外,相較于簡(jiǎn)單模型,復(fù)雜模型中FMM 1的運(yùn)算時(shí)間隨模型網(wǎng)格節(jié)點(diǎn)數(shù)的增加上升得更快。

3.3 Marmousi速度模型

Marmousi速度模型的計(jì)算區(qū)域?yàn)?.36km×7.49km,源點(diǎn)位置為(3.68km,0),圖8a展示了該模型下網(wǎng)格節(jié)點(diǎn)數(shù)為737×750時(shí)的測(cè)試結(jié)果。

如圖8a所示,Marmousi速度模型中,采用3種算法所計(jì)算的旅行時(shí)結(jié)果整體匹配良好,這進(jìn)一步證明了FMM 1算法的穩(wěn)定性。由于Marmousi速度模型的復(fù)雜度高于鹽丘模型,因此采用FSM需要9次迭代才能滿足迭代終止的條件。如圖8b所示,采用FSM所需的CPU時(shí)間明顯大于采用FMM 2所需的CPU時(shí)間。結(jié)合圖6b、圖7b和圖8b,進(jìn)一步證明了模型的復(fù)雜度影響了FMM 1的運(yùn)算效率。隨著模型復(fù)雜度的增加,FMM 1的CPU時(shí)間隨網(wǎng)格節(jié)點(diǎn)數(shù)增加而增長(zhǎng)的趨勢(shì)愈發(fā)強(qiáng)烈。

圖7 鹽丘速度模型的測(cè)試結(jié)果a 3種算法的旅行時(shí)場(chǎng); b 3種算法的CPU時(shí)間

圖8 Marmousi速度模型的測(cè)試結(jié)果a 3種算法的旅行時(shí)場(chǎng); b 3種算法的CPU時(shí)間

4 討論與分析

從上述的數(shù)值模擬結(jié)果不難發(fā)現(xiàn),模型網(wǎng)格節(jié)點(diǎn)數(shù)增多后,插入排序快速推進(jìn)算法的計(jì)算效率有所下降。為了進(jìn)一步分析影響插入排序方法計(jì)算效率的因素,我們統(tǒng)計(jì)了模型網(wǎng)格節(jié)點(diǎn)數(shù)均為401×401時(shí),各個(gè)模型中每個(gè)點(diǎn)在插入排序過(guò)程中向前移動(dòng)的次數(shù),統(tǒng)計(jì)結(jié)果如圖9所示。圖中某一點(diǎn)的移動(dòng)次數(shù)越多,則該點(diǎn)在運(yùn)算過(guò)程中所需的排序時(shí)間越長(zhǎng)。

為了更好地判斷旅行時(shí)場(chǎng)的有序程度與計(jì)算效率之間的關(guān)系,我們對(duì)不同模型中旅行時(shí)的有序程度進(jìn)行量化,并定義其計(jì)算公式為:

(12)

式中:μ為有序程度量,μ越小則表明旅行時(shí)場(chǎng)越有序,反之則越無(wú)序;τi為i點(diǎn)在排序過(guò)程中向前移動(dòng)的次數(shù);n為模型的總網(wǎng)格節(jié)點(diǎn)數(shù)。

圖9a的色標(biāo)范圍為0～140次,圖9b,圖9c和圖9d 的色標(biāo)范圍為0～700次。對(duì)比圖9a和圖5,可以發(fā)現(xiàn)向前移動(dòng)次數(shù)較多的點(diǎn)集中在因源點(diǎn)奇異性導(dǎo)致的誤差較大的區(qū)域(圖9a中相對(duì)源點(diǎn)呈45°方向的明亮區(qū)域)。圖9c和圖9d分別與圖7a和圖8a存在相似的分布規(guī)律,這證明源點(diǎn)奇異性不僅會(huì)影響插入排序快速推進(jìn)算法的計(jì)算精度,還會(huì)影響其計(jì)算效率。圖9中3個(gè)模型的有序程度量分別為18.902,122.310,133.302。根據(jù)有序程度量的定義可知,圖9c 和圖9d的兩個(gè)速度模型的旅行時(shí)場(chǎng)比光滑非均勻速度模型的旅行時(shí)場(chǎng)更無(wú)序。回顧之前的CPU時(shí)間統(tǒng)計(jì)結(jié)果可知,相較于簡(jiǎn)單模型(圖6b),復(fù)雜模型(圖7b和圖8b)中采用FMM 1所需的CPU時(shí)間曲線隨模型網(wǎng)格節(jié)點(diǎn)數(shù)的增加上升得更快,證明了旅行時(shí)場(chǎng)越無(wú)序,插入排序快速推進(jìn)算法的效率越低。

圖9中排序次數(shù)較多的區(qū)域明顯集中在由源點(diǎn)奇異性導(dǎo)致的誤差較大的區(qū)域,證明了源點(diǎn)奇異性與排序次數(shù)之間存在相關(guān)性。由于極坐標(biāo)網(wǎng)格形態(tài)與點(diǎn)源地震波場(chǎng)的形態(tài)相似,因此當(dāng)震源點(diǎn)位于極坐標(biāo)原點(diǎn)時(shí),源點(diǎn)奇異性問(wèn)題可以被很好地壓制[28]。為了排除源點(diǎn)奇異性對(duì)計(jì)算結(jié)果的影響,從而更好地分析排序方法與計(jì)算效率之間的關(guān)系,我們?cè)跇O坐標(biāo)系下對(duì)插入排序快速推進(jìn)算法展開了測(cè)試。除了極坐標(biāo)插入排序快速推進(jìn)算法之外,其它能壓制源點(diǎn)奇異性問(wèn)題的方法還包括在震源點(diǎn)周圍局部網(wǎng)格加密和因式分解方法等[22,29-31]。

圖9 各模型中各計(jì)算點(diǎn)在插入排序過(guò)程中向前移動(dòng)的次數(shù)a 光滑非均勻速度模型(未固定色標(biāo)范圍); b 光滑非均勻速度模型(固定色標(biāo)范圍); c 鹽丘速度模型(固定色標(biāo)范圍); d Marmousi速度模型(固定色標(biāo)范圍)

利用光滑非均勻速度模型測(cè)試極坐標(biāo)插入排序快速推進(jìn)算法。計(jì)算區(qū)域設(shè)置為半徑7.1km的半圓,本文只展示10km×10km的矩形計(jì)算區(qū)域,模型網(wǎng)格節(jié)點(diǎn)數(shù)為401×401,源點(diǎn)位于極坐標(biāo)原點(diǎn)(5.0km,0)處,模型的速度分布情況與圖4一致,測(cè)試結(jié)果如圖10所示。

圖10a中的色標(biāo)范圍為0～0.04s(與圖5一致),對(duì)比圖5和圖10a可知,采用極坐標(biāo)插入排序快速推進(jìn)算法得到的誤差遠(yuǎn)小于直角坐標(biāo)插入排序快速推進(jìn)算法得到的誤差。前者得到的數(shù)值解與解析解之間的誤差主要由模型網(wǎng)格化過(guò)程中導(dǎo)致的截?cái)嗾`差組成,誤差分布光滑,極大地壓制了由源點(diǎn)奇異性引起的誤差。因此,采用極坐標(biāo)插入排序快速推進(jìn)算法得到的旅行時(shí)計(jì)算結(jié)果更接近理論結(jié)果。圖10b中,窄帶點(diǎn)數(shù)組中的全部計(jì)算點(diǎn)向前移動(dòng)的次數(shù)全部為0,即有序程度量為0,理論上旅行時(shí)場(chǎng)處于高度有序狀態(tài)。

我們又統(tǒng)計(jì)了光滑非均勻速度模型在不同網(wǎng)格點(diǎn)數(shù)條件下,采用極坐標(biāo)插入排序快速推進(jìn)算法(PFMM 1)和極坐標(biāo)三叉樹堆排序快速推進(jìn)算法(PFMM 2)所需的CPU時(shí)間,并將統(tǒng)計(jì)結(jié)果與之前的統(tǒng)計(jì)結(jié)果進(jìn)行了比較,結(jié)果分別如圖11a和圖11b所示。

根據(jù)圖11的統(tǒng)計(jì)結(jié)果進(jìn)一步分析快速推進(jìn)方法的適用范圍。該算法的運(yùn)算成本可以粗略地拆分為:

圖10 光滑非均勻速度模型下極坐標(biāo)插入排序快速推進(jìn)算法的測(cè)試結(jié)果a 解析解與數(shù)值解之間的誤差; b 插入排序過(guò)程中全部計(jì)算點(diǎn)向前移動(dòng)的次數(shù)

圖11 光滑非均勻速度模型下各個(gè)快速推進(jìn)算法的CPU時(shí)間統(tǒng)計(jì)結(jié)果a PFMM 1和PFMM 2的CPU時(shí)間; b 5種算法的CPU時(shí)間

K=Ka+Kb+Kc

(13)

式中:K表示實(shí)現(xiàn)快速推進(jìn)算法所需的總運(yùn)算成本,Ka表示實(shí)現(xiàn)排序方法所需的運(yùn)算成本;Kb表示不同排序方法對(duì)窄帶點(diǎn)排序所需的運(yùn)算成本;Kc表示計(jì)算旅行時(shí)所需的運(yùn)算成本。

設(shè)插入排序快速推進(jìn)算法的總運(yùn)算成本為K1,且分別由Ka1(實(shí)現(xiàn)插入排序方法所需的運(yùn)算成本),Kb1(插入排序方法對(duì)窄帶點(diǎn)排序所需的運(yùn)算成本)和Kc1(插入排序快速推進(jìn)算法計(jì)算旅行時(shí)所需的運(yùn)算成本)組成。設(shè)堆排序快速推進(jìn)算法的總運(yùn)算成本為K2,分別由Ka2(實(shí)現(xiàn)堆排序方法所需的運(yùn)算成本),Kb2(堆排序方法對(duì)窄帶點(diǎn)排序所需的運(yùn)算成本)和Kc2(堆排序快速推進(jìn)算法計(jì)算旅行時(shí)所需的運(yùn)算成本)組成。由于插入排序方法的實(shí)現(xiàn)過(guò)程比堆排序方法的更簡(jiǎn)單(方法實(shí)現(xiàn)的運(yùn)算成本更少),因此Ka1Kb2;由于旅行時(shí)計(jì)算方法完全一樣,這兩種算法計(jì)算旅行時(shí)所需的時(shí)間相差無(wú)幾,因此Kc1≈Kc2。

綜上所述,分析各測(cè)試模型的CPU時(shí)間統(tǒng)計(jì)結(jié)果,可以發(fā)現(xiàn)影響算法計(jì)算效率的關(guān)鍵因素是Ka和Kb在總運(yùn)算成本K中所占比重。對(duì)于常規(guī)的直角坐標(biāo)快速推進(jìn)算法,由于源點(diǎn)奇異性問(wèn)題的存在,造成窄帶點(diǎn)數(shù)組的有序程度低,因此增加模型網(wǎng)格節(jié)點(diǎn)數(shù)量會(huì)增加Kb在運(yùn)算成本中所占的比重,且由于Kb1>Kb2,最終導(dǎo)致隨著網(wǎng)格節(jié)點(diǎn)數(shù)增加插入排序方法所需的運(yùn)算成本顯著增加。在極坐標(biāo)系下,快速推進(jìn)算法克服了源點(diǎn)奇異性問(wèn)題,Kb在總運(yùn)算成本中所占的比重小,故受網(wǎng)格點(diǎn)數(shù)變化的影響也小,且由于Ka1

對(duì)比圖11b中5種算法的CPU時(shí)間可知,PFMM 1的計(jì)算效率明顯優(yōu)于其它方法,說(shuō)明在克服了源點(diǎn)奇異性問(wèn)題的情況下,快速推進(jìn)算法無(wú)需采用復(fù)雜的排序方法對(duì)窄帶點(diǎn)數(shù)組進(jìn)行排序,即可通過(guò)簡(jiǎn)單的插入排序方法可取得最佳的計(jì)算效果。

5 結(jié)論

本文對(duì)程函方程的因果條件進(jìn)行分析,探討了旅行時(shí)場(chǎng)隱含的有序性。在此基礎(chǔ)上利用插入排序方法傳統(tǒng)的堆排序方法對(duì)窄帶點(diǎn)排序,實(shí)現(xiàn)了基于插入排序方法的快速推進(jìn)算法。采用數(shù)值模擬方法對(duì)插入排序快速推進(jìn)算法、三叉樹堆排序快速推進(jìn)算法和快速掃描算法進(jìn)行測(cè)試,統(tǒng)計(jì)和比較3種算法的精度和效率,結(jié)果表明:插入排序快速推進(jìn)算法的計(jì)算精度要高于常規(guī)的三叉樹堆排序快速推進(jìn)算法,計(jì)算結(jié)果與作為穩(wěn)定參考組的快速掃描算法計(jì)算結(jié)果完全相同,達(dá)到高度收斂狀態(tài),三叉樹堆排序快速推進(jìn)算法存在計(jì)算結(jié)果不夠收斂的問(wèn)題。在計(jì)算效率的研究中,分析插入排序快速推進(jìn)算法的運(yùn)算效率與旅行時(shí)場(chǎng)的有序性之間的關(guān)系可知,對(duì)于一般(未解決源點(diǎn)奇異性問(wèn)題)的直角坐標(biāo)快速推進(jìn)算法,源點(diǎn)奇異性問(wèn)題會(huì)導(dǎo)致旅行時(shí)場(chǎng)無(wú)序,從而導(dǎo)致插入排序快速推進(jìn)算法的計(jì)算效率不能達(dá)到理想狀態(tài),此時(shí)插入排序快速推進(jìn)算法的運(yùn)算效率只在模型網(wǎng)格點(diǎn)數(shù)較少的情況下具有優(yōu)勢(shì)。但在壓制了源點(diǎn)奇異性問(wèn)題后,旅行時(shí)場(chǎng)的有序性接近理論狀態(tài),此時(shí)極坐標(biāo)插入排序快速推進(jìn)算法的計(jì)算效率全面優(yōu)于常規(guī)的堆排序快速推進(jìn)算法。

同樣的技術(shù)應(yīng)用于不同領(lǐng)域時(shí),應(yīng)當(dāng)考慮不同領(lǐng)域自身的特殊性。在選取快速推進(jìn)算法的排序方法時(shí),應(yīng)當(dāng)考慮到旅行時(shí)場(chǎng)潛在的規(guī)律性,而不應(yīng)僅考慮排序方法的排序能力強(qiáng)弱。雖然對(duì)于存在源點(diǎn)奇異性問(wèn)題的快速推進(jìn)算法而言,隨著模型網(wǎng)格點(diǎn)數(shù)的增加,堆排序方法的計(jì)算效率優(yōu)于插入排序方法,但造成這一現(xiàn)象的原因并不只是堆排序方法的排序能力優(yōu)于插入排序方法,其基本前提還包括了旅行時(shí)場(chǎng)受源點(diǎn)奇異性影響導(dǎo)致其有序性被破壞。隨著旅行時(shí)算法的不斷發(fā)展,源點(diǎn)奇異性問(wèn)題被解決,因此采用插入排序方法替換原算法中的堆排序方法將是更好的選擇。