馬樹青，藍(lán) 強(qiáng)，朱勛江，邱 偉，顏 冰

(1.國(guó)防科技大學(xué)氣象海洋學(xué)院，湖南長(zhǎng)沙 410073；2.中國(guó)人民解放軍96863部隊(duì)，河南洛陽 471000)

0 引 言

匹配場(chǎng)反演[1]與匹配場(chǎng)定位[2]一直以來是水聲領(lǐng)域的研究重點(diǎn)，其基本思想是將通過聲場(chǎng)計(jì)算得到的仿真計(jì)算結(jié)果與實(shí)際測(cè)量結(jié)果進(jìn)行匹配，并運(yùn)用反問題方法得到反演結(jié)果。在水聲匹配反演中，聲場(chǎng)計(jì)算一直是制約其運(yùn)算效率的重要原因：對(duì)于寬帶反演，由于需要計(jì)算不同頻率聲源激發(fā)的聲場(chǎng)結(jié)果，計(jì)算通常需要上百次以上；對(duì)于時(shí)序反演，如卡爾曼濾波[3]、粒子濾波[4-8]等方法，通常需要多次迭代計(jì)算聲場(chǎng)。由于每次聲場(chǎng)計(jì)算時(shí)間在秒或分鐘量級(jí)，整體反演耗時(shí)較長(zhǎng)，從而限制了人們的深入研究。

與此同時(shí)，水聲模型種類較多，通常包括了射線模型[9]、簡(jiǎn)正波模型[10]、拋物方程模型[11]等，不同模型存在不同特點(diǎn)與適用范圍。以往研究人員常針對(duì)某類模型進(jìn)行并行優(yōu)化，以提高其計(jì)算效率，如針對(duì)RAM模型的MPIRAM研究[12]、針對(duì)射線模型的并行優(yōu)化[13]等，針對(duì)特定種類模型的優(yōu)化可利用OPENMP[14]或MPI[15]技術(shù)，實(shí)現(xiàn)單機(jī)多核或多節(jié)點(diǎn)并行，以提高計(jì)算效率。但是在實(shí)際應(yīng)用過程中，常需要多類模型進(jìn)行比較性研究，以選擇最為合適的模型提高精度。因此本文提出了一種跨模型的多級(jí)并行框架，在水聲模型計(jì)算的上層采用多級(jí)并行技術(shù)，實(shí)現(xiàn)匹配反演的上層并行計(jì)算與底層水聲模型計(jì)算的解耦合，可應(yīng)用至不同類型的水聲模型中，從而簡(jiǎn)化并行開發(fā)的復(fù)雜度。本文將其應(yīng)用至海洋聲速剖面反演的典型應(yīng)用中，評(píng)估了其運(yùn)行效率與并行優(yōu)化效率。

1 水聲模型多層次并行思想

在水聲領(lǐng)域，基于匹配反演應(yīng)用的典型計(jì)算模式是對(duì)海洋環(huán)境參數(shù)以及聲源參數(shù)進(jìn)行大規(guī)模采樣，根據(jù)采樣結(jié)果生成水聲模型輸入文件，調(diào)用水聲模型進(jìn)行計(jì)算。由于需要對(duì)反演參數(shù)進(jìn)行采樣，會(huì)造成水聲模型的大規(guī)模循環(huán)計(jì)算。本文提出了多層次并行框架，在水聲模型計(jì)算的循環(huán)層次實(shí)現(xiàn)多核與多節(jié)點(diǎn)的混合并行，具體包含了配置文件生成模塊、多層次并行控制模塊以及水聲傳播模型計(jì)算模塊，水聲模型多層次并行計(jì)算框架如圖1所示。

圖1 水聲模型多層次并行計(jì)算框架Fig.1 Block diagram of multi-level parallel computation of underwater acoustic model

當(dāng)前開源的水聲傳播模型種類多樣，但這些模型編譯后通常只可以使用單線程進(jìn)行計(jì)算，當(dāng)遇到大規(guī)模計(jì)算問題時(shí)耗時(shí)較長(zhǎng)，但是如果針對(duì)各類模型進(jìn)行并行優(yōu)化，則不利于模型比較等研究，且對(duì)硬件平臺(tái)的適用性不強(qiáng)。本文采用了一種多層次并行框架來控制不同的聲傳播模型，該框架僅需較少地修改模型輸入輸出接口，即將模型的輸出文件名與輸入文件名設(shè)置為相同，并編譯為可執(zhí)行文件即可。該框架無需改變水聲模型的使用方式、環(huán)境輸入文件格式以及輸出文件格式等，保證了其適用性。

框架中配置文件生成模塊將根據(jù)匹配反演算法的需要，生成不同名稱的環(huán)境輸入配置文件，并生成環(huán)境輸入文件的描述配置文件ref.txt，描述配置文件包含了水聲模型程序路徑、環(huán)境文件總數(shù)目以及每個(gè)環(huán)境文件名稱。

多層次并行控制模塊，主要面向多核集群系統(tǒng)體系結(jié)構(gòu)平臺(tái)，以環(huán)境輸入文件的描述配置文件作為輸入，根據(jù)環(huán)境文件的數(shù)目即計(jì)算任務(wù)數(shù)，采用MPI+OpenMP的混合并行方法，實(shí)現(xiàn)計(jì)算任務(wù)的多核以及多節(jié)點(diǎn)兩級(jí)并行映射。節(jié)點(diǎn)間的并行是通過MPI編程方法實(shí)現(xiàn)的，主要完成計(jì)算任務(wù)在節(jié)點(diǎn)間的任務(wù)劃分，節(jié)點(diǎn)內(nèi)多核的并行主要通過OpenMP編程方法實(shí)現(xiàn)，是針對(duì)每個(gè)節(jié)點(diǎn)負(fù)責(zé)的計(jì)算任務(wù)在多核間進(jìn)行任務(wù)劃分。

水聲傳播模型模塊負(fù)責(zé)處理各種具體方法的聲場(chǎng)計(jì)算，可以根據(jù)需要使用不同的水聲傳播模型，將多層次并行控制模塊與水聲模型計(jì)算模塊進(jìn)行解耦合，提高了多層次并行計(jì)算框架的靈活性和適應(yīng)性。

2 水聲模型多層次并行的應(yīng)用與性能

2.1 粒子濾波反演聲速剖面理論

一般的粒子濾波步驟包括了初始化粒子、重要性權(quán)值計(jì)算、重采樣保留權(quán)值較大的粒子、計(jì)算估計(jì)結(jié)果等步驟。在k時(shí)刻傳感器量測(cè)為yk，聲源狀態(tài)與環(huán)境參數(shù)的聯(lián)合待估計(jì)量表示為θk。粒子濾波的步驟為：

(1) 初始化，在k=0時(shí)刻，需初始化采樣粒子，可以從先驗(yàn)分布p(θ0)中產(chǎn)生粒子群，各粒子重要性權(quán)值ωi0=1/N，其中下標(biāo)表示時(shí)刻，上標(biāo)表示粒子。

(2) 重要性權(quán)值計(jì)算，根據(jù)設(shè)定每個(gè)粒子的權(quán)重大小ωik，實(shí)現(xiàn)“優(yōu)質(zhì)”粒子的復(fù)制以及“劣質(zhì)”粒子的淘汰，并實(shí)現(xiàn)粒子的狀態(tài)轉(zhuǎn)移與量測(cè)更新。

從而可以得到重要性權(quán)值更新方程：

為簡(jiǎn)化積分可以選擇：

(3) 重采樣：重采樣是抑制粒子退化的重要手段，可以采用隨機(jī)重采樣、系綜重采樣、多項(xiàng)式重采樣、殘差重采樣等多種方法。重采樣決策可以按照Geweke等提出的相對(duì)效率(Relative Numerical Efficiency, RNE)來進(jìn)行門限判斷，估計(jì)有效抽樣尺度為

越小退化越嚴(yán)重，可通過手動(dòng)設(shè)定門限Nth與其進(jìn)行比較，判斷是否進(jìn)行重采樣步驟，自動(dòng)化的設(shè)定則需要嘗試采用其他方法，在本文后續(xù)中，為更好觀察參數(shù)后驗(yàn)分布情況，使用門限為75%的確定值，小于該值后進(jìn)行系綜重采樣。

(4) 計(jì)算結(jié)果輸出：可以按照后驗(yàn)均值估計(jì)、最大后驗(yàn)(Maximum A Posteriori, MAP)估計(jì)或后驗(yàn)誤差標(biāo)準(zhǔn)差估計(jì)等方法得到估計(jì)的后驗(yàn)概率密度，MAP表示為

粒子濾波算法應(yīng)用于環(huán)境參數(shù)反演等反問題求解時(shí)，常用的算法流程圖如圖2所示，將待反演的參數(shù)設(shè)置為粒子集合，通過初始化粒子值后開始狀態(tài)預(yù)測(cè)與更新的步驟。本文中對(duì)經(jīng)驗(yàn)正交函數(shù)(Empirical Orthogonal Function, EOF)的前三階系數(shù)進(jìn)行反演，根據(jù)先驗(yàn)觀測(cè)數(shù)據(jù)，前三階系數(shù)的測(cè)量噪聲的方差選取分別為22，1.22，0.72，測(cè)量值為聲壓的幅度，粒子濾波器初始值根據(jù)先驗(yàn)信息在設(shè)定范圍區(qū)間進(jìn)行隨機(jī)選取。

圖2 粒子濾波算法流程圖Fig.2 Flow chart of the particle filter algorithm

其中主要的計(jì)算部分為似然概率的計(jì)算。似然概率的求解需要對(duì)環(huán)境參數(shù)進(jìn)行采樣，然后循環(huán)地調(diào)用水聲傳播模型進(jìn)行多次聲場(chǎng)計(jì)算，根據(jù)傳感器接收到的真實(shí)聲場(chǎng)大小計(jì)算出似然概率。因此，本文采用了多層次并行框架，實(shí)現(xiàn)大量聲場(chǎng)計(jì)算的多級(jí)并行。

2.2 反演環(huán)境

圖3給出了距離不相關(guān)環(huán)境下的仿真波導(dǎo)環(huán)境參數(shù)。設(shè)接收陣為32個(gè)陣元的垂直均勻線陣，聲源輻射頻率為120 Hz，深度位于50 m，距離接收陣5 000 m。該區(qū)域水深為200 m，海底的密度，吸收系數(shù)及聲速分別為ρb=1.6 g·cm-3，αb=0.05 dB·λ-1，cb=1 700 m·s-1。本文使用的聲傳播模型為拋物方程模型RAM。

圖3 海洋環(huán)境仿真參數(shù)Fig.3 Parameters of the simulated marine environment

2.3 仿真結(jié)果

實(shí)驗(yàn)設(shè)計(jì)了兩個(gè)反演場(chǎng)景，一個(gè)反演動(dòng)態(tài)的聲速剖面，另一個(gè)是反演靜態(tài)的聲速剖面，對(duì)于靜態(tài)聲速剖面反演實(shí)驗(yàn)，EOF的前三階系數(shù)真值為20、12、7。

在聲速剖面的靜態(tài)反演實(shí)驗(yàn)中，進(jìn)行了4次測(cè)試，粒子濾波算法的粒子數(shù)量分別設(shè)為24、120、480和960。圖4為粒子數(shù)分別取24、120、480、960時(shí)，對(duì)應(yīng)EOF的前3陣系數(shù)a1、a2、a3的反演結(jié)果隨著時(shí)間步長(zhǎng)的變化情況，當(dāng)粒子數(shù)較少時(shí)，粒子的前三階系數(shù)并沒有出現(xiàn)收斂的現(xiàn)象，而當(dāng)粒子數(shù)達(dá)到480或者960時(shí)，可以看到反演粒子很快收斂到了EOF前三階系數(shù)的真值附近。

圖4 粒子數(shù)不同時(shí)各粒子的三個(gè)分量隨著時(shí)間步長(zhǎng)的變化Fig.4 Variations of three components of each particle with time step for different numbers of particles

圖5為不同粒子數(shù)下，代表EOF前三階系數(shù)的粒子分布情況以及EOF前三階系數(shù)的真值和粒子濾波算法的最大后驗(yàn)估計(jì)結(jié)果。圖5中紅色的豎直線代表的是EOF系數(shù)真值，黑色豎直線代表的是后驗(yàn)概率估計(jì)。從圖5中可以看出，隨著粒子數(shù)的增多，粒子采樣越充分，最后得到EOF前三階系數(shù)的后驗(yàn)概率估計(jì)與EOF系數(shù)的真值越接近。

圖5 不同粒子數(shù)反演中粒子的分布及相應(yīng)的EOF前三階系數(shù)的真值和粒子濾波算法的最大后驗(yàn)估計(jì)結(jié)果Fig.5 The particle distributions and the corresponding true values and MAP results of the three-order coefficients EOF in the particle inversion simulation with different numbers of particles

此外，本文還進(jìn)行了聲速剖面的動(dòng)態(tài)反演測(cè)試，設(shè)置3階EOF系數(shù)隨時(shí)間變化規(guī)律滿足正弦的變化規(guī)律，其中值為20、12、7。圖6為粒子集合分別為24、120、480、960的粒子濾波算法對(duì)聲速剖面前三階EOF系數(shù)的反演結(jié)果。圖6中時(shí)長(zhǎng)均為200 s，實(shí)線表示真實(shí)的EOF的系數(shù)變化情況，虛線表示濾波反演得到的結(jié)果。從圖6中可以看到，粒子濾波算法在真實(shí)值劇烈變化的情況下，增大初始隨機(jī)采樣的粒子集合，能夠較好地跟蹤上系數(shù)的時(shí)間變化。當(dāng)粒子數(shù)較少時(shí)，從圖6(a)可以看到，反演效果較差，與真實(shí)值相差很大，隨著粒子集合的增大，第二階和第三階系數(shù)已經(jīng)可以較好地跟蹤上，但第一階系數(shù)在時(shí)間步長(zhǎng)100 s以前，都不能很好地收斂，只有在粒子集合數(shù)達(dá)到480時(shí)，也就是圖6(d)才能很好地收斂。這主要是因?yàn)榈谝浑A系數(shù)所占能量較多，因此第一階系數(shù)的真值變化與方差大小會(huì)對(duì)聲速剖面產(chǎn)生劇烈擾動(dòng)，而后兩階系數(shù)逐級(jí)減小，所占的能量較少，進(jìn)而對(duì)接收信號(hào)的聲壓影響較小。

圖6 粒子數(shù)不同情況下基于粒子濾波算法的EOF動(dòng)態(tài)反演結(jié)果及真實(shí)值Fig.6 EOF inversion results based on the particle filter algorithm for different numbers of particles and their true values

2.4 多層次并行性能

從2.3節(jié)中可以看出基于粒子濾波的反演方法中，粒子數(shù)越多反演效果越好，但計(jì)算量也會(huì)隨之增大。本節(jié)將針對(duì)粒子濾波反演的多層次并行性能進(jìn)行測(cè)試。

本文首先在天河2號(hào)超級(jí)計(jì)算平臺(tái)的單個(gè)節(jié)點(diǎn)上進(jìn)行了多核性能測(cè)試。天河2號(hào)每個(gè)計(jì)算節(jié)點(diǎn)中有2個(gè)CPU、共24個(gè)計(jì)算內(nèi)核，本文采用的線程數(shù)分別為1、2、4、8、16、24，當(dāng)固定問題的規(guī)模保持粒子總數(shù)為1 584時(shí)，圖7為不同線程下的加速比。從圖7可以看到當(dāng)節(jié)點(diǎn)中24個(gè)處理器內(nèi)核都用上時(shí)，加速比達(dá)到了20.9，并行效率達(dá)到了87.5%，取得較好的加速性能。

圖7 不同線程下粒子濾波算法的多核并行運(yùn)行時(shí)間Fig.7 Comparison of the multi-core parallel computation runtime of particle filter algorithm under different threads

為了測(cè)試更大規(guī)模的粒子數(shù)，在天河2號(hào)上進(jìn)行了強(qiáng)擴(kuò)展實(shí)驗(yàn)，即保持粒子數(shù)為12 288，使用的節(jié)點(diǎn)數(shù)分別為1、2、4、8、16、32、64、128、256，圖8為不同節(jié)點(diǎn)數(shù)下相對(duì)于1個(gè)節(jié)點(diǎn)的加速比，可以看出，當(dāng)使用的節(jié)點(diǎn)數(shù)在64個(gè)以下時(shí)，基本上取得線性的加速比，然而當(dāng)超過64個(gè)節(jié)點(diǎn)時(shí)，強(qiáng)擴(kuò)展性有所下降，主要是由于代碼中的串行部分占的時(shí)間占主導(dǎo)，雖然最后的加速效果有所下降，但最多取得了109倍的加速性能，這將有效提高大規(guī)模粒子濾波反演的效率。

圖8 多層次并行的強(qiáng)擴(kuò)展測(cè)試相對(duì)加速度Fig.8 Relative acceleration ratio of multi-level parallel strong expansion test

當(dāng)結(jié)合OpenMP與MPI并行后，為了更好對(duì)比串行與并行間的差異，假設(shè)粒子濾波算法中設(shè)置3 072個(gè)粒子，時(shí)長(zhǎng)為400 s，串行運(yùn)行時(shí)間需要5天以上；若使用64個(gè)節(jié)點(diǎn)，每節(jié)點(diǎn)使用24線程，則并行運(yùn)行時(shí)間僅需400 s左右。

3 結(jié) 論

粒子濾波方法在聲速剖面的EOF系數(shù)反演中取得了很好的效果，然而粒子濾波算法隨著粒子數(shù)的逐漸增大，對(duì)計(jì)算的需求也逐漸增大。本文提出的多層次并行方法，采用了MPI+OpenMP的混合并行技術(shù)，實(shí)現(xiàn)了粒子濾波算法向多核集群體系結(jié)構(gòu)的多級(jí)并行映射。通過在天河2號(hào)超級(jí)計(jì)算平臺(tái)的并行測(cè)試，驗(yàn)證了多層次并行方法的有效性。本文研究的多層次并行計(jì)算方法能夠在未來基于粒子濾波方法的多參數(shù)反演應(yīng)用中發(fā)揮更大的作用。