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

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

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

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

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

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

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

1 方法原理

1.1 程函方程

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

(1)

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

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

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

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

(2)

當(dāng)|Txmin-Tymin|

(3)

式中:下標(biāo)i和j分別對應(yīng)x和y軸方向上的網(wǎng)格節(jié)點序號;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表示取括號中兩個值的較小值。

1.2 因果條件

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

(4)

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

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

(5)

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

(6)

離散形式表示為:

(7)

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

圖1 不同的波場傳播方向和波前面示意

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

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

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

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

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

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

2 方法實現(xiàn)

2.1 窄帶技術(shù)

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

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

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

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

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

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

1) 初始化

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

2) 迭代

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

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

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

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

3) 截止條件

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

3 數(shù)值模擬

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

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

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

3.1 光滑非均勻速度模型

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

(8)

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

該速度模型的解析解為:

(9)

(10)

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

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

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

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

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

E=|Tana-Tnum|

(11)

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

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

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

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

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

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

3.2 鹽丘速度模型

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

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

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

3.3 Marmousi速度模型

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

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

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

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

4 討論與分析

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

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

(12)

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

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

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

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

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

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

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

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

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

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

K=Ka+Kb+Kc

(13)

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

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

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

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

5 結(jié)論

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

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