張朝金孫炳文1,2?

(1中國(guó)科學(xué)院水聲環(huán)境特性重點(diǎn)實(shí)驗(yàn)室 北京 100190)

(2中國(guó)科學(xué)院聲學(xué)研究所 北京 100190)

(3中國(guó)科學(xué)院大學(xué) 北京 100049)

0 引言

射線法具有成熟的理論基礎(chǔ)，適用范圍廣且其物理圖像清晰，基于射線追蹤技術(shù)的射線模型被廣泛應(yīng)用于水下聲場(chǎng)快速預(yù)報(bào)[1]、海洋聲層析[2?3]、三維聲場(chǎng)計(jì)算[4?5]、匹配場(chǎng)定位[6]等水聲領(lǐng)域。同時(shí)，將射線理論與簡(jiǎn)正波等模型相互結(jié)合，可以建立有效的適用于復(fù)雜環(huán)境中的低頻聲場(chǎng)計(jì)算模型[7]，或可用于分析復(fù)雜環(huán)境下聲場(chǎng)特性的變化規(guī)律[8]。

由于傳統(tǒng)射線方法在方程建立與求解上的限制，無(wú)法有效計(jì)算聲影區(qū)和會(huì)聚區(qū)的聲場(chǎng)。幾何衍射理論可以處理聲影區(qū)的計(jì)算問題，但該理論中聲束寬度的微小設(shè)置誤差會(huì)導(dǎo)致很大的計(jì)算問題[1]；W.K.B.近似方法可以有效計(jì)算會(huì)聚區(qū)聲場(chǎng)，但是該近似解難以推廣至二維和三維的水聲傳播情況[9]。20世紀(jì)80年代，Porter等[10]提出了高斯射線追蹤理論，通過更加寬容的高斯函數(shù)表征聲波強(qiáng)度在聲波束管中的分布，有效地解決了經(jīng)典射線方法中聲影區(qū)和會(huì)聚區(qū)聲場(chǎng)能量不準(zhǔn)確的問題和聲束寬度設(shè)定的限制。對(duì)于高頻、與距離有關(guān)(簡(jiǎn)正波、快速場(chǎng)、拋物理論不適用)的海洋聲傳播問題，該模型計(jì)算更準(zhǔn)確，在計(jì)算精度和實(shí)時(shí)性上都有很大的提高?；诟咚孤暰€束理論的Bellhop聲場(chǎng)計(jì)算模型的出現(xiàn)也加快了高斯聲線束模型的廣泛應(yīng)用。

隨著我國(guó)對(duì)海洋的研究由淺海向深海、由近場(chǎng)向遠(yuǎn)程的探索，對(duì)實(shí)時(shí)且準(zhǔn)確預(yù)估聲場(chǎng)提出了越來越高的要求。隨著計(jì)算機(jī)計(jì)算能力的飛速提高和并行計(jì)算技術(shù)的快速發(fā)展[11]，水下聲場(chǎng)并行計(jì)算成為聲場(chǎng)快速預(yù)估的有效途徑，充分利用計(jì)算機(jī)多核計(jì)算能力，將聲場(chǎng)計(jì)算并行化，是提高聲場(chǎng)計(jì)算效率最簡(jiǎn)單、直接的方法，對(duì)于該問題的研究已經(jīng)成為水下聲傳播的熱點(diǎn)問題之一[12?13]。高斯射線束(Gaussian ray bundle,GRAB)理論具有計(jì)算速度快、物理意義清晰并且適合并行化處理等優(yōu)點(diǎn)，許多學(xué)者就射線理論并行計(jì)算聲場(chǎng)展開了一系列的探究。其中較為簡(jiǎn)單的就是應(yīng)用OpenMP在共享存儲(chǔ)的模式下進(jìn)行并行計(jì)算，但OpenMP不適用于集群。因此，陳連榮等[14]利用基于消息傳遞編程模型(Message passing interface,MPI)進(jìn)行了高斯射線聲場(chǎng)模型并行算法的設(shè)計(jì)和編程，該方法可以應(yīng)用于多處理器集群系統(tǒng)。MPI適合各種機(jī)器，但是其編程模型復(fù)雜，程序可靠性差。

綜上所述，本文基于高斯射線束理論探究Bellhop模型的數(shù)值計(jì)算原理與并行計(jì)算的可行性，考慮現(xiàn)有并行方法的優(yōu)缺點(diǎn)，利用多線程技術(shù)，開發(fā)出并行化射線模型BellhopMP(Bellhop multi process)，并且對(duì)BellhopMP的并行效率進(jìn)行了測(cè)試與分析，從而快速、準(zhǔn)確地預(yù)報(bào)聲場(chǎng)。

1 并行化射線模型設(shè)計(jì)

1.1 高斯射線理論與Bellhop聲場(chǎng)計(jì)算模型

考慮二維柱坐標(biāo)系下聲場(chǎng)的Helmholtz方程：

式(1)中，c為聲速，ω為聲源處的角頻率。根據(jù)級(jí)數(shù)解形式p(r,z)=exp(iωτ)An/(iω)n，忽略高階項(xiàng)，可得程函方程和強(qiáng)度方程：

引入射線坐標(biāo)和射線束概念，可分別獲得描述射線束中心聲線的控制方程和射線束強(qiáng)度的表達(dá)式：

其中，s為射線坐標(biāo)中的聲線弧長(zhǎng)量，θ0為聲線的初射掠射角，[r(s),z(s)]描述聲線軌跡，|J(s)|為表征柱對(duì)稱情況下引入的幾何擴(kuò)展量和s處相鄰聲線間的弧長(zhǎng)變化。

根據(jù)|J(s)|所表征的物理意義，并結(jié)合射線坐標(biāo)，針對(duì)經(jīng)典的射線理論，可建立|J(s)|計(jì)算表達(dá)式和動(dòng)力學(xué)射線方程，即

其中，p和q為射線坐標(biāo)s處相鄰射線之間的弧長(zhǎng)及其相對(duì)變化。

高斯波束跟蹤方法，不僅保留了傳統(tǒng)射線的一些優(yōu)點(diǎn)，還能夠避免傳統(tǒng)射線的實(shí)際應(yīng)用困難。在水聲學(xué)研究領(lǐng)域，目前存在幾種典型的高斯波束跟

蹤方法，理論基本相似，主要區(qū)別在于波束內(nèi)部聲場(chǎng)能量的分布形式、沿波束行進(jìn)過程中能量擴(kuò)展形式等方面。高斯射線束方法，通過引入較為簡(jiǎn)單的波束寬度限制，可以較為準(zhǔn)確地處理簡(jiǎn)單的聲影區(qū)和會(huì)聚區(qū)聲場(chǎng)問題。高斯射線理論的基本思路是將聲線看作截面能量按照高斯函數(shù)分布的聲束，該聲束的中心聲線滿足標(biāo)準(zhǔn)的控制方程，如式(3a)所示，在中心聲線鄰域可以構(gòu)造出

其中，η是距離中心聲線的垂直距離；A是一個(gè)任意的常數(shù)，通過與均勻介質(zhì)聲場(chǎng)參考值比較，由聲源性質(zhì)確定；τ(s)為沿聲線的相位延遲；p和q為在上文基礎(chǔ)上由高斯波束束寬和曲率定義的復(fù)弧長(zhǎng)及相對(duì)變化。考慮點(diǎn)源問題，聲場(chǎng)的射線理論表達(dá)式為

式(6)中，?θ為聲束的夾角，θm為相應(yīng)聲線的掠射角，τ(s)在控制方程計(jì)算的基礎(chǔ)上進(jìn)行求解，p和q則在高斯波束束寬和曲率給定的初始值下根據(jù)動(dòng)力學(xué)射線方程(4b)計(jì)算。

Bellhop就是基于高斯聲線束理論開發(fā)的聲場(chǎng)計(jì)算模型，其計(jì)算準(zhǔn)確性和可靠性已經(jīng)得到重復(fù)檢驗(yàn)。該模型寬容性較好，可以應(yīng)用高斯、帽型等波束進(jìn)行聲場(chǎng)修正，且可以獲得傳播損失、本征聲線以及聲線到達(dá)時(shí)間序列等多種實(shí)用數(shù)據(jù)，同時(shí)可以處理海面、海底和聲速剖面隨距離變化的情況。Bellhop允許輸入設(shè)計(jì)好的具有指向性的聲源，可以通過輸入具體的地聲參數(shù)或者通過海底反射系數(shù)表征海底邊界進(jìn)行聲場(chǎng)計(jì)算。模型從提出至今已發(fā)展出多個(gè)不同的版本，利用不同的編程語(yǔ)言發(fā)展出Fortran版、Matlab版、Python版。本文為了方便進(jìn)行并行化模型的設(shè)計(jì)，把Fortran語(yǔ)言版本的源程序移植為C++語(yǔ)言版本：BellhopC，C++版本的Bellhop可以在不支持Fortran語(yǔ)言的平臺(tái)中使用。

1.2 BellhopMP并行化射線模型設(shè)計(jì)

高斯射線模型Bellhop的聲場(chǎng)計(jì)算中每條聲線對(duì)聲場(chǎng)的貢獻(xiàn)是獨(dú)立疊加的，即各條聲線計(jì)算之間沒有直接聯(lián)系，因而非常適宜按照不同的射線計(jì)算來分配任務(wù)。因此可以在聲線層面上進(jìn)行并行化處理，將聲線跟蹤及對(duì)聲場(chǎng)的貢獻(xiàn)分配到多個(gè)計(jì)算核心上進(jìn)行計(jì)算，然后合并聲場(chǎng)即可得到考慮所有聲線貢獻(xiàn)的聲場(chǎng)。

根據(jù)高斯射線模型的特點(diǎn)，并行化模型BellhopMP對(duì)Bellhop的并行化體現(xiàn)在聲線分組和聲場(chǎng)合并，并行過程如圖1所示。BellhopMP讀取計(jì)算參數(shù)文件，根據(jù)用戶設(shè)定的線程數(shù)量進(jìn)行多核并行計(jì)算；如果用戶不指定線程數(shù)量，則使用硬件系統(tǒng)能提供的最大并行能力。在主線程中執(zhí)行Bellhop的所有計(jì)算流程，讀取聲線初射角度信息后，根據(jù)可以使用的并行能力，將聲線角度分為多個(gè)區(qū)間，主線程只計(jì)算第一個(gè)區(qū)間的聲線及其對(duì)聲場(chǎng)的貢獻(xiàn)。其他每個(gè)線程計(jì)算某一個(gè)聲線角度區(qū)間的聲場(chǎng)結(jié)果。最后在主線程中合并所有的聲場(chǎng)計(jì)算結(jié)果，輸出到結(jié)果文件中。

圖1 聲場(chǎng)計(jì)算模型BellhopMP并行處理流程Fig.1 Flow chart of BellhopMP

2 并行計(jì)算算例分析

圖2 Pekeris波導(dǎo)情況下Bellhop、BellhopC與BellhopMP的聲場(chǎng)計(jì)算結(jié)果Fig.2 The sound f i eld computed by Bellhop,BellhopC and BellhopMP in the Pekeris waveguide

為了檢驗(yàn)并行算法的準(zhǔn)確性，實(shí)驗(yàn)過程中對(duì)Acoustic Toolbox中多個(gè)算例使用Bellhop、BellhopC和BellhopMP進(jìn)行了計(jì)算分析，得到的結(jié)果基本一致，基本驗(yàn)證了并行射線程序BellhopMP的準(zhǔn)確性和可靠性，文章中只列出了經(jīng)典情況(calibB_gb和MunkB_Coh_gb)下的結(jié)果。淺海Pekeris波導(dǎo)情況下，環(huán)境參數(shù)和幾何參數(shù)為水中聲速1500 m/s，海底為聲速1590 m/s、密度1.2 g/cm3、衰減0.5 dB/λ的液態(tài)半無(wú)限空間；海深100 m，接收距離5 km，聲源深度50 m，頻率250 Hz。Bellhop、BellhopC和BellhopMP聲場(chǎng)計(jì)算的傳播損失結(jié)果如圖2所示，Bellhop和BellhopC相對(duì)于BellhopMP的計(jì)算結(jié)果誤差分布(x軸為兩者計(jì)算得到的傳播損失的差，y軸對(duì)應(yīng)該差值的個(gè)數(shù))如圖3所示。

圖3 Pekeris波導(dǎo)情況下不同版本Bellhop模型的聲場(chǎng)計(jì)算結(jié)果誤差分布Fig.3 The error distribution of dif f erent versions of Bellhop in the Pekeris waveguide

典型深海Munk剖面情況下，聲速剖面如圖4所示，其他環(huán)境參數(shù)和幾何參數(shù)為海底為聲速1600 m/s、密度1.8 g/cm3、衰減0.8 dB/λ的液態(tài)半無(wú)限空間；海深5000 m，接收距離100 km，聲源深度1000 m，頻率50 Hz。Bellhop、BellhopC和BellhopMP聲場(chǎng)計(jì)算的傳播損失結(jié)果如圖5所示，Bellhop和BellhopC相對(duì)于BellhopMP的計(jì)算結(jié)果誤差分布如圖6所示。

圖4 Munk聲速剖面Fig.4 Munk sound speed prof i le

圖5 深海情況下Bellhop、BellhopC與BellhopMP的聲場(chǎng)計(jì)算結(jié)果Fig.5 The sound f i eld computed by Bellhop,BellhopC and BellhopMP in deep water

圖6 深海情況下不同版本Bellhop模型的聲場(chǎng)計(jì)算結(jié)果誤差分布Fig.6 The error distribution of dif f erent versions of Bellhop in deep water

數(shù)值計(jì)算的結(jié)果顯示BellhopMP與Bellhop、BellhopC計(jì)算結(jié)果一致，誤差的產(chǎn)生是由于浮點(diǎn)數(shù)的表示及累積誤差，對(duì)聲場(chǎng)計(jì)算的精度沒有影響，這驗(yàn)證了BellhopMP計(jì)算的準(zhǔn)確性。同時(shí)，由淺海和深海環(huán)境下的聲場(chǎng)計(jì)算結(jié)果可以看出，基于高斯波束的射線模型可有效實(shí)現(xiàn)水下聲場(chǎng)的計(jì)算，在淺海中聲線相互干涉從而體現(xiàn)簡(jiǎn)正波傳播特征，而在深海中可有效給出聲影區(qū)和會(huì)聚區(qū)的聲場(chǎng)值。

3 并行效率

本文使用Intel(R)Core(TM)i7-4770K CPU@3.50 GHz和Intel(R)Xeon(R)CPU E3-1505M v6@3.00 GHz兩種CPU進(jìn)行并行計(jì)算效能的測(cè)試，兩種CPU都是4個(gè)核心，在超線程技術(shù)的支持下，可以進(jìn)行8個(gè)線程的同時(shí)計(jì)算，以發(fā)揮CPU的最大計(jì)算能力。操作系統(tǒng)均為Windows 10。對(duì)于所有算例分別使用Bellhop、BellhopC、BellhopMP(依次設(shè)置線程數(shù)為1～8)、BellhopMP(默認(rèn)參數(shù)，使用全部可用的線程數(shù))計(jì)算100次，統(tǒng)計(jì)每個(gè)算例計(jì)算耗時(shí)的均值和方差。以經(jīng)典深海情況(MunkB_Coh_gb)為例，其環(huán)境參數(shù)和幾何參數(shù)見圖5的計(jì)算說明，其均值與標(biāo)準(zhǔn)差結(jié)果如圖7所示。

以上兩種CPU下的計(jì)算結(jié)果趨勢(shì)一致，如圖7所示。但是在Xeon處理器上，耗時(shí)的方差較大，這應(yīng)該與CPU的架構(gòu)差別有關(guān)。BellhopC的耗時(shí)比Bellhop的多，這與Bellhop的C++語(yǔ)言實(shí)現(xiàn)的效率有關(guān)，有進(jìn)一步優(yōu)化的可能。需要說明的是，對(duì)于計(jì)算一條聲線的聲場(chǎng)，BellhopMP并沒有使用多線程來計(jì)算，且當(dāng)算例的計(jì)算耗時(shí)很少時(shí)，由于BellhopMP需要額外的任務(wù)分配和結(jié)果合并操作，導(dǎo)致計(jì)算時(shí)間相對(duì)于BellhopC并沒有減少，有時(shí)計(jì)算時(shí)間反而會(huì)增加。

為了定量描述并行計(jì)算的能力，定義并行效率(P)為加速比與線程數(shù)目之比，反映執(zhí)行并行算法時(shí)，平均每個(gè)處理器的執(zhí)行效率。其中加速比為同一個(gè)任務(wù)在單處理器系統(tǒng)和并行處理器系統(tǒng)中運(yùn)行消耗的時(shí)間的比率：

T1是在單個(gè)處理器下的運(yùn)行時(shí)間，TN是在有N個(gè)處理器并行系統(tǒng)中的運(yùn)行時(shí)間。

根據(jù)圖7的計(jì)算結(jié)果，取T1為Bellhop的運(yùn)行時(shí)間，BellhopMP在不同線程下相對(duì)于Bellhop的并行化效率如圖8所示。BellhopMP采用大粒度的并行處理，在主線程中不僅計(jì)算聲場(chǎng)，還負(fù)責(zé)聲線區(qū)間的分割和聲場(chǎng)計(jì)算結(jié)果的合并。這導(dǎo)致了隨著線程數(shù)量的增加，時(shí)間減小得越來越不明顯，并行效率逐漸降低，如圖8中i7 CPU在使用4個(gè)線程時(shí)并行效率只有0.5，這是由于Fortran版本的Bellhop的計(jì)算時(shí)間只有12 s。

并行效率與串行的計(jì)算時(shí)間有關(guān)，串行時(shí)間越長(zhǎng)，并行效率越高。如圖9所示為使用Acousitc Toolbox中提供的測(cè)試算例(arcticB_cpp、DickinsCervenyB、MunkB_Coh_CervenyR、MunkB_Coh_gb、MunkB_Coh_SGB、MunkB_gbt)，皆使用4個(gè)線程(與CPU的核心數(shù)相同)時(shí)的并行效率。由圖中結(jié)果可以看出，隨著串行計(jì)算時(shí)間的增加，并行效率逐漸增加。如圖9(a)中可以看出，當(dāng)Bellhop計(jì)算時(shí)間為105 s時(shí)，i7 CPU的并行效率可達(dá)0.9。因而，對(duì)于遠(yuǎn)距離深海環(huán)境下的長(zhǎng)時(shí)間聲場(chǎng)計(jì)算問題，本文提出的并行模型可在保證計(jì)算準(zhǔn)確性的前提下，有效提高計(jì)算效率，減少計(jì)算時(shí)間。

圖7 各類計(jì)算方式耗時(shí)的均值和標(biāo)準(zhǔn)差。其中，F(xiàn)、C、1–8、A依次表示Bellhop、BellhopC、BellhopMP(線程數(shù)1～8)、BellhopMP(默認(rèn)參數(shù))Fig.7 The mean value and standard variance of computation time cost by dif f erent models on two platform.F,C,1–8,A represents dif f erent models:Bellhop,BellhopC,BellhopMP(thread count:1～8),BellhopMP(default)

圖8 不同線程數(shù)目時(shí)BellhopMP相對(duì)于Bellhop的并行計(jì)算效率Fig.8 The parallel efficiency of BellhopMP using dif f erent count of thread on two platforms

圖9 串行模型Bellhop的計(jì)算時(shí)間和使用4個(gè)線程的BellhopMP的并行效率之間的關(guān)系Fig.9 The relationship of serial computation time and parallel efficiency of BellhopMP using 4 threads

4 結(jié)論

本文在分析高斯射線束理論的基礎(chǔ)上，根據(jù)射線模型的特征，利用多線程技術(shù)對(duì)射線模型Bellhop進(jìn)行了并行化處理，開發(fā)了并行化射線模型BellhopMP，該模型在保證計(jì)算結(jié)果準(zhǔn)確性的同時(shí)，在多核計(jì)算機(jī)上能大幅度提高計(jì)算速度且有較好的并行效率，可以解決深海遠(yuǎn)程以及復(fù)雜情況下聲場(chǎng)計(jì)算時(shí)間過長(zhǎng)的問題。同時(shí)該模型串行時(shí)間越長(zhǎng)，并行效率越高，故對(duì)于深海超遠(yuǎn)程等長(zhǎng)時(shí)間聲場(chǎng)計(jì)算問題，BellhopMP可有效降低計(jì)算時(shí)間，而對(duì)于耗時(shí)較短的聲場(chǎng)計(jì)算，并行效率提高的程度有限。相較于其他并行方法，BellhopMP并行化處理方式簡(jiǎn)單易行，程序可靠性較高，這種并行化處理方式也適用于DSP(Digital signal processor)等多核心處理器計(jì)算平臺(tái)。