一種高效的寬帶簡正波本征值計算方法

楊雪峰 ，王好忠 ，駱文于 ，胡長青

?

一種高效的寬帶簡正波本征值計算方法

楊雪峰1,2，王好忠3，駱文于4，胡長青1

(1. 中國科學院聲學研究所東海研究站，上海 201815；2. 中國科學院大學，北京 100049；3. 中國海洋大學，山東青島 266100；4. 中國科學院聲學研究所，北京 100190)

為提高寬帶簡正波本征值的計算效率，漢密爾頓方法通過公式變換抵消頻率項，將寬帶簡正波本征值的計算由頻率、本征值實部和虛部的三維尋根降低至本征值實部和虛部的二維尋根，可以一次性求解單號簡正波所有頻率對應的簡正波本征值。在已有工作的基礎上優(yōu)化和完善了漢密爾頓寬帶簡正波本征值算法，并加入并行計算方法進一步提高計算效率。以簡正波模型KRAKENC作為對比，通過若干數(shù)值算例驗證了該算法對于寬帶簡正波本征值的計算精度和計算效率。數(shù)值仿真結果顯示，在保證寬帶簡正波本征值計算精度的前提下，該方法的計算效率相對KRAKENC有著明顯的優(yōu)勢；加入并行算法后，該方法的計算效率得到大幅提高。

寬帶；簡正波方法；本征值；漢密爾頓方法；高效

0 引言

寬帶技術在地聲波形分析領域應用廣泛，但因為海洋強時變性的聲傳播環(huán)境會給傳統(tǒng)聲吶工作頻率下的聲信號施加強烈而不可預測的擾動，所以較少地應用于水聲研究領域。隨著聲吶設備工作頻率的低頻化和水聲應用領域?qū)β暡ㄐ畔⒀芯康木毣?，水聲寬帶波形預報逐漸成為熱門的研究課題[1]。

常見的寬帶技術有時域波動方程法和頻域傅里葉合成法。利用現(xiàn)有聲傳播模型進行聲場計算的頻域傅里葉合成方法實現(xiàn)簡單，但需要進行帶寬內(nèi)頻率采樣，并對每個頻率采樣點都進行一次聲場計算，效率較低。McDonald等[2]采取只計算少量頻率采樣點對應的簡正波本征值、本征函數(shù)和耦合系數(shù)的方法提高計算效率，而其他頻率對應的相關量則通過插值來獲得；Knobles等[3-4]采用類似的方法，將若干計算得到的垂直波數(shù)組成集合，通過Galerkin方法計算相鄰頻率對應的未知的垂直波數(shù)；Skarsoulis[5-6]則采用二階傅里葉合成方法，涉及計算的頻率采樣點只有一個。

上述方法都是簡化傅里葉合成方法，減少頻率采樣點的計算次數(shù)，在犧牲一定計算精度的基礎上提高效率。文獻[7-8]提出了漢密爾頓方法，將寬帶簡正波本征值的計算由頻率、本征值實部和虛部的三維尋根降低至本征值實部和虛部的二維尋根，可有效提高計算效率[7-8]。本文發(fā)展了漢密爾頓方法并對該方法寬帶簡正波本征值的計算效率進行了數(shù)值仿真驗證。

1 基本理論

漢密爾頓方法之所以可以使寬帶簡正波本征值的計算由三維尋根問題降低至二維尋根問題，是因為它通過公式變換擺脫了頻率項，只剩下由簡正波本征值實部和虛部組成的表達式。下面以等聲速剖面波導為例介紹漢密爾頓方法的基本原理。簡正波本征值的相函數(shù)依其滿足的相位條件可以定義如下[7-9]：

則

定義漢密爾頓函數(shù)

2 簡正波本征值并行計算

由式(8)中可以看到，每號簡正波對應了各自的漢密爾頓函數(shù)，因此每號簡正波本征值的求解是互相獨立的，可以采用并行算法提高計算效率。本文采用MATLAB提供的parfor進行并行計算。

parfor是MATLAB設置的針對循環(huán)并行計算的關鍵字，它將在串行計算的for循環(huán)拆解成互不相關的模塊并交給不同的worker計算，之后再將不同的計算結果匯總。其中，實際進行模塊計算的MATLAB進程叫worker，例如雙核CPU可提供2個worker同時計算。另外，parfor并行算法規(guī)定了其循環(huán)內(nèi)的變量必須有明確的定義類型，保證每次循環(huán)變量不相關。具體定義可參考相關文獻[10]。

3 數(shù)值算例

本節(jié)通過數(shù)值算例，研究漢密爾頓寬帶簡正波本征值計算方法的特性。以簡正波模型KRAKENC對比本文方法的計算精度和計算效率。

3.1 計算精度仿真驗證

圖1 Pekeris波導示意圖

首先來看簡正波本征值的計算精度，分別取聲源頻率為50、75、100 Hz和150 Hz計算簡正波本征值，KRAKENC和漢密爾頓方法的部分計算結果如表1所示(表1中，i表示虛部)。

表1 簡正波本征值計算精度對比

從表1中可以看出，漢密爾頓方法對簡正波本征值的計算與KRAKENC相比精度較高。本征值實部和虛部分別可以精確計算到小數(shù)點后六位左右和七位左右。實際上，這是由求解式(7)時采用的數(shù)值方法決定的，實部和虛部分別采用的是遍歷法和二分法，而實部步進的幅度決定了本征值計算的精度，幅度過小，計算效率就會降低；幅度過大，計算誤差就會變大。另外，漢密爾頓方法是一次性求解單號簡正波所有頻率對應的本征值，所求頻率對應的簡正波本征值需要通過插值的方法來求得，插值過程也會引入誤差。

進一步計算波導傳播損失來對比漢密爾頓方法的計算精度。全局矩陣耦合簡正波方法[11]是駱文于等人提出的精確、穩(wěn)定而高效的耦合簡正波方法，本文將漢密爾頓方法計算得到的本征值代入全局矩陣耦合簡正波方法的相關公式，計算波導傳播損失。設定圖1中波導的聲源深度為50 m，接收器深度為50 m，接收距離為4 km，分別計算聲源頻率為50 Hz和75 Hz的傳播損失，如圖2和圖3所示。圖中紅色虛線是使用了全局矩陣耦合簡正波理論的漢密爾頓方法的計算結果，藍色實線是KRAKENC計算結果。圖2中兩種方法計算結果的均方根誤差為0.10 dB，圖3中該誤差為0.55 dB?？梢钥闯?，漢密爾頓方法與KRAKENC的計算結果吻合較好，說明漢密爾頓方法在計算簡正波淺海低頻、近程聲傳播中有較好的適應性。

圖2 接收深度50 m處Pekeris波導傳播損失計算結果(聲源深度為50 m,聲源頻率為50 Hz)

圖3 接收深度50 m處Pekeris波導傳播損失計算結果，聲源深度為50 m,聲源頻率為75 Hz.

3.2 計算效率仿真驗證

本節(jié)仍采用圖1所示的波導研究漢密爾頓方法寬帶簡正波本征值的計算效率。參照模型采用修改后只計算傳播簡正波的KRAKENC模型。為討論頻率對計算效率的影響，取三種相同帶寬、不同頻率的聲源，即30～80 Hz、80～130 Hz、130～180 Hz，頻率采樣間隔都是1 Hz。幾種方法運行的PC環(huán)境均為雙核處理器，計算結果如表2所示。

表2 不同頻帶的簡正波本征值計算效率對比

從表2前3項可以看出，在相同帶寬情況下，KRAKENC計算三個頻率段所用的時間相差不大；從第1、4、5項可以看出，時間與帶寬大致成線性關系。這是因為在每個頻率采樣點，KRAKENC都需要求解一次簡正波本征值。而漢密爾頓方法可一次求解單號簡正波所有頻率采樣點對應的簡正波本征值，從表2可以看出，漢密爾頓方法的計算效率遠高于KRAKENC，且從前3項可以看出，隨著頻率的增高即簡正波號數(shù)的增多，漢密爾頓方法所用的時間也隨之增加，但仍比KRAKENC方法低幾個數(shù)量級。加入并行算法后，漢密爾頓方法的計算效率提高不大，這是因為該算例涉及的簡正波號數(shù)很少，并行算法在不同worker的數(shù)據(jù)傳遞中耗費了時間。當簡正波號數(shù)較多時，并行算法的優(yōu)勢就顯露無疑。例如，當圖1波導的水深為400 m、聲源頻帶為200～300 Hz時，涉及計算的傳播簡正波號數(shù)最高可達到56號。漢密爾頓方法和并行漢密爾頓方法的計算時間如表3所示。

表3 漢密爾頓方法并行前后的計算時間對比

可以看出，當涉及計算量較大時，并行計算在此雙核CPU環(huán)境中可以節(jié)省近一半的時間。

為討論頻率采樣間隔對計算效率的影響，取相同帶寬、不同采樣間隔的聲源作對比。設定聲源帶寬為50～90 Hz，頻率采樣間隔分別為1、0.5、0.2 Hz，PC環(huán)境不變，計算結果如表4所示。

從表4中可以看到，漢密爾頓方法不受頻率采樣間隔的影響，而KRAKENC模型則需要在每個頻率采樣點都計算一次簡正波本征值，因而所用的計算時間遠遠大于漢密爾頓方法。

表4 不同頻率采樣間隔的簡正波本征值計算時間對比

4 結論

本文進一步發(fā)展了漢密爾頓方法，用于寬帶簡正波本征值的快速計算。數(shù)值仿真顯示，漢密爾頓方法的計算效率不受頻率大小和頻率采樣間隔的影響，在保證簡正波本征值計算精度的基礎上，能夠大幅提高寬帶簡正波本征值的計算效率，且在采用并行計算方法后，計算效率提升更加明顯。漢密爾頓方法作為精確高效的寬帶聲場計算方法，為水聲寬帶波形預報提供了有力的工具。

由于漢密爾頓方法的數(shù)值求解過程采用遍歷法和二分法，計算精度和計算效率受到一定限制，因此，可在以后的工作中改進數(shù)值計算的方法。另外，本文只涉及傳播模式簡正波的計算，泄漏模式的寬帶簡正波本征值計算有待進一步研究。

[1] 唐帥, 笪良龍, 謝駿. 水聲寬帶信號波形預報技術研究[J]. 海洋科學, 2012, 36(11): 67-72. TANG Shuai, DA Lianglong, XIE Jun. Research on underwater broad-band signal prediction[J]. Marine Sciences, 2012, 36(11): 67-72.

[2] MCDONALD B E, COLLINS M D, KUPERMAN W A, et al. Comparison of data and model predictions for Heard Island acoustic transmissions[J]. J. Acoust. Soc. Am., 1994, 96(4): 2357- 2370.

[3]KNOBLES D P, KOCH R A. A time series analysis of sound propagation in a strongly multipath shallow water environment with an adiabatic normal mode approach[J]. IEEE Journal of Oceanic Engineering, 1996, 21(1): 1-13.

[4] KNOBLES D P, WESTWOOD E K, LE MOND J E. Modal time-series structure in a shallow-water environment [Hudson Canyon region][J]. IEEE Journal of Oceanic Engineering, 1998, 23(3): 188-202.

[6]SKARSOULIS E K. Fast coupled-mode approximation for broad-band pulse propagation in a range-dependent ocean[J]. IEEE Journal of Oceanic Engineering, 1999, 24(2): 172-182.

[7]王好忠. 簡正波寬帶本征波數(shù)和等效耦合矩陣計算方法[D]. 青島: 中國海洋大學, 2015. WANG Haozhong. Methods to compute wideband modal eigenwavenumber and equivalent coupled-matrix[D]. Qingdao: Ocean University of China, 2015.

[8] WANG H, WANG N, GAO D. A Hamiltonian method for finding broadband modal eigenvalues[J]. J. Acoust. Soc. Am., 2012, 131 (2): 1047-1054.

[9] TINDLE C T, CHAPMAN N R. A phase function for finding normal mode eigenvalues over a layered elastic bottom[J]. J. Acoust. Soc. Am., 1994, 96(3): 1777-1782.

[10]維. 實戰(zhàn)Matlab之并行程序設計[M]. 北京: 北京航空航天大學出版社, 2012. LIU Wei. Parallel programming of Matlab[M]. Beijing: Beihang University Press, 2012.

[11]LUO W Y, YANG C M, QIN J X, et al. A numerically stable coupled-mode formulation for acoustic propagation in range- dependent waveguides[J]. SCIENCE CHINA Physics, Mechanics & Astronomy, 2012, 55(4): 572-588.

An efficient method for broadband eigenvalue computation

YANG Xue-feng1,2, WANG Hao-zhong3, LUO Wen-yu4, HU Chang-qing1

(1. Shanghai Acoustic Laboratory, Institute of Acoustics, Chinese Academy of Sciences, Shanghai 201815, China;2. University of Chinese Academy of Sciences, Beijing 100049, China; 3. Ocean University of China, Qingdao 266100,Shandong, China;4. Institute of Acoustics, Chinese Academy of Sciences, Beijing 100190, China)

Based on existing Hamiltonian method for eigenvalue calculation of normal mode, an efficient method for broadband eigenvalue calculation is developed. Hamiltonian method gets rid of the frequency term by formula transformation to change the root finding problem from a three-dimensional problem with respect to frequency as well as real and imaginary parts of eigenvalue into a two-dimensional one. As a result, the eigenvalues of singular normal mode corresponding to all frequencies could be solved at one time. A MATLAB parallel computing method and other optimizations are included to improve the computational efficiency of Hamiltonian method. The significant advantage of Hamiltonian method in efficiency is indicated by several numerical simulation comparisons with KRAKENC, while retaining the same accuracy.

broadband; normal-mode method; eigenvalue; Hamiltonian method; efficient

P733.21

A

1000-3630(2018)-03-0201-04

10.16300/j.cnki.1000-3630.2018.03.001

2017-07-01;

2017-08-25

楊雪峰(1987－), 男, 山東煙臺人, 博士研究生, 研究方向為水聲技術。

胡長青, E-mail: hchq@mail.ioa.ac.cn