一种高效的宽带简正波本征值计算方法

杨雪峰 ，王好忠 ，骆文于 ，胡长青



一种高效的宽带简正波本征值计算方法

杨雪峰1,2，王好忠3，骆文于4，胡长青1

(1. 中国科学院声学研究所东海研究站，上海 201815；2. 中国科学院大学，北京 100049；3. 中国海洋大学，山东青岛 266100；4. 中国科学院声学研究所，北京 100190)

为提高宽带简正波本征值的计算效率，汉密尔顿方法通过公式变换抵消频率项，将宽带简正波本征值的计算由频率、本征值实部和虚部的三维寻根降低至本征值实部和虚部的二维寻根，可以一次性求解单号简正波所有频率对应的简正波本征值。在已有工作的基础上优化和完善了汉密尔顿宽带简正波本征值算法，并加入并行计算方法进一步提高计算效率。以简正波模型KRAKENC作为对比，通过若干数值算例验证了该算法对于宽带简正波本征值的计算精度和计算效率。数值仿真结果显示，在保证宽带简正波本征值计算精度的前提下，该方法的计算效率相对KRAKENC有着明显的优势；加入并行算法后，该方法的计算效率得到大幅提高。

宽带；简正波方法；本征值；汉密尔顿方法；高效

0 引言

宽带技术在地声波形分析领域应用广泛，但因为海洋强时变性的声传播环境会给传统声呐工作频率下的声信号施加强烈而不可预测的扰动，所以较少地应用于水声研究领域。随着声呐设备工作频率的低频化和水声应用领域对声波信息研究的精细化，水声宽带波形预报逐渐成为热门的研究课题[1]。

常见的宽带技术有时域波动方程法和频域傅里叶合成法。利用现有声传播模型进行声场计算的频域傅里叶合成方法实现简单，但需要进行带宽内频率采样，并对每个频率采样点都进行一次声场计算，效率较低。McDonald等[2]采取只计算少量频率采样点对应的简正波本征值、本征函数和耦合系数的方法提高计算效率，而其他频率对应的相关量则通过插值来获得；Knobles等[3-4]采用类似的方法，将若干计算得到的垂直波数组成集合，通过Galerkin方法计算相邻频率对应的未知的垂直波数；Skarsoulis[5-6]则采用二阶傅里叶合成方法，涉及计算的频率采样点只有一个。

上述方法都是简化傅里叶合成方法，减少频率采样点的计算次数，在牺牲一定计算精度的基础上提高效率。文献[7-8]提出了汉密尔顿方法，将宽带简正波本征值的计算由频率、本征值实部和虚部的三维寻根降低至本征值实部和虚部的二维寻根，可有效提高计算效率[7-8]。本文发展了汉密尔顿方法并对该方法宽带简正波本征值的计算效率进行了数值仿真验证。

1 基本理论

汉密尔顿方法之所以可以使宽带简正波本征值的计算由三维寻根问题降低至二维寻根问题，是因为它通过公式变换摆脱了频率项，只剩下由简正波本征值实部和虚部组成的表达式。下面以等声速剖面波导为例介绍汉密尔顿方法的基本原理。简正波本征值的相函数依其满足的相位条件可以定义如下[7-9]：

则

定义汉密尔顿函数

2 简正波本征值并行计算

由式(8)中可以看到，每号简正波对应了各自的汉密尔顿函数，因此每号简正波本征值的求解是互相独立的，可以采用并行算法提高计算效率。本文采用MATLAB提供的parfor进行并行计算。

parfor是MATLAB设置的针对循环并行计算的关键字，它将在串行计算的for循环拆解成互不相关的模块并交给不同的worker计算，之后再将不同的计算结果汇总。其中，实际进行模块计算的MATLAB进程叫worker，例如双核CPU可提供2个worker同时计算。另外，parfor并行算法规定了其循环内的变量必须有明确的定义类型，保证每次循环变量不相关。具体定义可参考相关文献[10]。

3 数值算例

本节通过数值算例，研究汉密尔顿宽带简正波本征值计算方法的特性。以简正波模型KRAKENC对比本文方法的计算精度和计算效率。

3.1 计算精度仿真验证

图1 Pekeris波导示意图

首先来看简正波本征值的计算精度，分别取声源频率为50、75、100 Hz和150 Hz计算简正波本征值，KRAKENC和汉密尔顿方法的部分计算结果如表1所示(表1中，i表示虚部)。

表1 简正波本征值计算精度对比

从表1中可以看出，汉密尔顿方法对简正波本征值的计算与KRAKENC相比精度较高。本征值实部和虚部分别可以精确计算到小数点后六位左右和七位左右。实际上，这是由求解式(7)时采用的数值方法决定的，实部和虚部分别采用的是遍历法和二分法，而实部步进的幅度决定了本征值计算的精度，幅度过小，计算效率就会降低；幅度过大，计算误差就会变大。另外，汉密尔顿方法是一次性求解单号简正波所有频率对应的本征值，所求频率对应的简正波本征值需要通过插值的方法来求得，插值过程也会引入误差。

进一步计算波导传播损失来对比汉密尔顿方法的计算精度。全局矩阵耦合简正波方法[11]是骆文于等人提出的精确、稳定而高效的耦合简正波方法，本文将汉密尔顿方法计算得到的本征值代入全局矩阵耦合简正波方法的相关公式，计算波导传播损失。设定图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 计算效率仿真验证

本节仍采用图1所示的波导研究汉密尔顿方法宽带简正波本征值的计算效率。参照模型采用修改后只计算传播简正波的KRAKENC模型。为讨论频率对计算效率的影响，取三种相同带宽、不同频率的声源，即30～80 Hz、80～130 Hz、130～180 Hz，频率采样间隔都是1 Hz。几种方法运行的PC环境均为双核处理器，计算结果如表2所示。

表2 不同频带的简正波本征值计算效率对比

从表2前3项可以看出，在相同带宽情况下，KRAKENC计算三个频率段所用的时间相差不大；从第1、4、5项可以看出，时间与带宽大致成线性关系。这是因为在每个频率采样点，KRAKENC都需要求解一次简正波本征值。而汉密尔顿方法可一次求解单号简正波所有频率采样点对应的简正波本征值，从表2可以看出，汉密尔顿方法的计算效率远高于KRAKENC，且从前3项可以看出，随着频率的增高即简正波号数的增多，汉密尔顿方法所用的时间也随之增加，但仍比KRAKENC方法低几个数量级。加入并行算法后，汉密尔顿方法的计算效率提高不大，这是因为该算例涉及的简正波号数很少，并行算法在不同worker的数据传递中耗费了时间。当简正波号数较多时，并行算法的优势就显露无疑。例如，当图1波导的水深为400 m、声源频带为200～300 Hz时，涉及计算的传播简正波号数最高可达到56号。汉密尔顿方法和并行汉密尔顿方法的计算时间如表3所示。

表3 汉密尔顿方法并行前后的计算时间对比

可以看出，当涉及计算量较大时，并行计算在此双核CPU环境中可以节省近一半的时间。

为讨论频率采样间隔对计算效率的影响，取相同带宽、不同采样间隔的声源作对比。设定声源带宽为50～90 Hz，频率采样间隔分别为1、0.5、0.2 Hz，PC环境不变，计算结果如表4所示。

从表4中可以看到，汉密尔顿方法不受频率采样间隔的影响，而KRAKENC模型则需要在每个频率采样点都计算一次简正波本征值，因而所用的计算时间远远大于汉密尔顿方法。

表4 不同频率采样间隔的简正波本征值计算时间对比

4 结论

本文进一步发展了汉密尔顿方法，用于宽带简正波本征值的快速计算。数值仿真显示，汉密尔顿方法的计算效率不受频率大小和频率采样间隔的影响，在保证简正波本征值计算精度的基础上，能够大幅提高宽带简正波本征值的计算效率，且在采用并行计算方法后，计算效率提升更加明显。汉密尔顿方法作为精确高效的宽带声场计算方法，为水声宽带波形预报提供了有力的工具。

由于汉密尔顿方法的数值求解过程采用遍历法和二分法，计算精度和计算效率受到一定限制，因此，可在以后的工作中改进数值计算的方法。另外，本文只涉及传播模式简正波的计算，泄漏模式的宽带简正波本征值计算有待进一步研究。

[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]王好忠. 简正波宽带本征波数和等效耦合矩阵计算方法[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]维. 实战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