馬樹青，藍 強，朱勛江，邱 偉，顏 冰

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

0 引 言

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

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

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

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

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

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

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

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

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

2 水聲模型多層次并行的應用與性能

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

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

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

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

從而可以得到重要性權值更新方程：

為簡化積分可以選擇：

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

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

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

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

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

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

2.2 反演環(huán)境

圖3給出了距離不相關環(huán)境下的仿真波導環(huán)境參數(shù)。設接收陣為32個陣元的垂直均勻線陣，聲源輻射頻率為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é)果

實驗設計了兩個反演場景，一個反演動態(tài)的聲速剖面，另一個是反演靜態(tài)的聲速剖面，對于靜態(tài)聲速剖面反演實驗，EOF的前三階系數(shù)真值為20、12、7。

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

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

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

圖5 不同粒子數(shù)反演中粒子的分布及相應的EOF前三階系數(shù)的真值和粒子濾波算法的最大后驗估計結(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

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

圖6 粒子數(shù)不同情況下基于粒子濾波算法的EOF動態(tài)反演結(jié)果及真實值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ù)越多反演效果越好，但計算量也會隨之增大。本節(jié)將針對粒子濾波反演的多層次并行性能進行測試。

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

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

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

圖8 多層次并行的強擴展測試相對加速度Fig.8 Relative acceleration ratio of multi-level parallel strong expansion test

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

3 結(jié) 論

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