声场波数积分截断波数自适应选取方法.pdf

第 42 卷第 4 期2023 年 8 月Vol.42,No.4Aug.,2023声学技术Technical Acoustics声场波数积分截断波数自适应选取方法李承帮1,2，江鹏飞1，孙军平1，衣雪娟1，林建恒1(1.中国科学院声学研究所北海研究站，山东青岛 266114；2.中国科学院大学，北京 100049)摘要：在采用波数积分法进行声场计算的过程中，需要选取合适的积分截断波数，文章提出一种应用于流体介质中宽带声场波数积分计算的截断波数自适应选取方法。首先根据波数域格林函数的衰减特性构造一个数学模型，然后利用卡尔曼滤波器对该数学模型的拟合参数进行跟踪和预测，最后根据预测的模型参数计算截断波数。仿真试验结果表明，该方法实现了给定精确度下的积分截断波数自适应选取，能够克服现有方法不能兼顾低频段精确度和高频段计算量的问题，并且不会引入太多额外的计算量。关键词：声场计算；波数积分；自适应；截断波数；海洋声学中图分类号：P733.21 文献标志码：A 文章编号：1000-3630(2023)-04-0533-08Adaptive selection of truncated wavenumber in sound field calculation with wavenumber integration methodLI Chengbang1,2,JIANG Pengfei1,SUN Junping1,YI Xuejuan1,LIN Jianheng1(1.Qingdao Branch,Institute of Acoustics,Chinese Academy of Sciences,Qingdao 266114,Shandong,China;2.University of Chinese Academy of Sciences,Beijing 100049,China)Abstract:It is necessary to select an appropriate integral truncated wavenumber in sound field calculation with wavenumber integration method.In this paper,an adaptive truncated wavenumber selection method for broadband sound filed calculation in fluid waveguide is proposed,which means that the wavenumber integral can be truncated at a given accuracy.A mathematical model is established to fit the attenuation characteristics of Greens function in wavenumber domain.The fitting parameters in the mathematical model are tracked and predicted by Kalman filters.The truncated wavenumber is estimated according to the mathematical model and its predicted parameters.The experimental results show that the truncated wavenumber selected by the proposed adaptive method can make the integral truncated at a given accuracy,which solves the problem that the existing methods cannot take into account the calculation accuracy at low frequency and the calculation amount at high frequency,and no much additional computation is introduced.Key words:sound field calculation;wavenumber integration;adaptive;truncated wavenumber;ocean acoustics0引 言波数积分法是一种常用的水下声场计算方法，此方法对波数域格林函数进行直接的数值积分，因而被认为是一种精确的声场计算方法。波数积分法由Pekeris引入到水下声场的计算中1，此后很多学者对其具体实现方法进行了研究。DiNapoli等2基于波数积分理论提出了声场计算效率非常高的快速场程序(Fast Field Program)，Schmidt3-4提出了直接全局矩阵法实现声场的波数积分计算，骆文于等5-6提出了一种可稳定计算Pekeris波导中声场的波数积分方法。波数积分法的一个重要步骤是对波数域格林函数进行逆汉克尔变换，此处会面临半无限区间的积分问题。一般来说，当水平波数大于某一个临界值时，波数域格林函数的幅值会随着水平波数的增大而衰减，并趋向于0。因此通常的做法是，当波数域格林函数幅值衰减到一定程度时，对波数积分区间进行截断，忽略大于此波数积分区间的声场。若截断波数取值太大，会增加不必要的运算量；若截断波数取值太小，则会给声场的计算带来误差。积分截断相当于给波数域格林函数加矩形窗，在截断处引入了不连续点，导致出现截断效应7。然而频率、介质类型、波导参数、声源与接收点的相对深度等参数会对波数域格林函数的衰减特性产生很大引用格式：李承帮,江鹏飞,孙军平,等.声场波数积分截断波数自适应选取方法J.声学技术,2023,42(4):533-540.LI Chengbang,JIANG Pengfei,SUN Junping,et al.Adaptive selection of truncated wavenumber in sound field calculation with wavenumber integration methodJ.Technical Acoustics,2023,42(4):533-540.DOI:10.16300/ki.1000-3630.2023.04.018收稿日期:2022-02-18;修回日期:2022-03-20基金项目:中 国 科 学 院 声 学 研 究 所 前 沿 探 索 基 金 项 目(QYTS202008),中国科学院海洋信息技术创新研究院前沿基础研究项目(QYJC201911)。作者简介:李承帮(1997),男,广东肇庆人,硕士研究生,研究方向为水下噪声。通信作者:林建恒,E-mail:2023 年声学技术的影响，在计算出波数域格林函数之前很难对截断波数进行合适的选取。在不同介质类型的海洋波导中，波数域格林函数会有不同的衰减特性。对于不含冰盖层且海底为液态的海洋波导，波数域格林函数的衰减特性比较简单，当水平波数大于临界值时，波数域格林函数不存在极点；而对于存在弹性海底或冰盖层的海洋波导，在波数域格林函数衰减较慢，在水平波数大于临界值的较大范围内，很远的位置仍包含声场的重要成分7。本文主要关注的是各层均为流体的海洋波导中，宽带声场波数积分计算的截断波数选取问题。当海洋波导中各层介质都是流体时，目前最常用的截断波数选取方法是把波导内最小声速对应的波数k0=cmin作为一个参考波数，截断波数取值为kt=k0，其中取一个固定的值。文献8认为可取1.1；文献9认为大多数情况下取1.11.2时即可包含波数域格林函数中绝大部分能量。在宽带声场计算过程中，若取固定值，则不能兼顾低频段的积分精确度和高频段的计算量。文献10提出一种消除截断效应的方法，在波数域格林函数截断点后，添加埃尔米特(Hermite)多项式函数，消除截断处的不连续性并使其迅速衰减。该方法虽然能避免出现截断效应，但是会给声场计算带来一定的误差。文献11提出了基于预估-校正思想的最大截止波数自动选取算法。但此方法需要在波数域格林函数计算过程中判断是否进行积分截断，这不利于实现波数域格林函数的并行计算。针对上述问题，本文提出一种应用于宽带点声源声场计算的积分截断波数自适应选取方法。首先构造一个数学模型拟合水平波数大于参考波数时的波数域格林函数，然后利用卡尔曼滤波器对模型的参数进行跟踪和预测，再根据已计算频率的模型参数估计出下一频率的模型参数，最后将估计的模型参数代入数学模型求解出下一个频率的截断波数，从而实现给定精确度下的积分截断。1截断波数自适应选取方法在利用波数积分法计算频率为f的点声源产生的声场时，需要对波数域格林函数G(kr;f)进行逆汉克尔变换：g(r;f)=0G(kr;f)J0(krr)krdkr12re-i40G(kr;f)kreikrrdkr(1)式中：J0(krr)为零阶贝塞尔函数，r为水平距离，kr为水平波数。当kr大于某个临界值时，G(kr;f)随kr的增大而衰减，在G(kr;f)明显衰减后进行积分截断不会带来很大误差。本文将积分截断波数记为kt，kt的取值取决于波数域格林函数的衰减特性，即G(kr;f)开始衰减的位置以及衰减的速率。考虑频率为fn=f0+nDf的宽带点声源声场的计算问题，其中n=1N，Df是频率间隔。由1.1节的算例可看出，随频率fn的变化，G(kr;fn)呈现出 一 定 的 变 化 规 律。因 而 可 根 据 已 计 算 的G(kr;f0)G(kr;f1)G(kr;fn)的衰减特性，估计G(kr;fn+1)的衰减特性。但由于G(kr;f)是数值计算的结果，为了便于跟踪和预测它的衰减特性，需要构造一个合适的数学模型对其进行拟合，并利用自适应滤波器对此数学模型的拟合参数进行跟踪和预测。根据模型参数的预测结果，可计算出下一个频率满足精度要求的截断波数。1.1波数域格林函数的衰减特性本节分析波数域格林函数的衰减特性，以选取合适的数学模型对其进行拟合。在Pekeris波导中，水层的波数域格林函数G1(kr;f)可表示为声源项与上行波、下行波的叠加：G1(kr;f)=A+1(kr)eikz1z+A-1(kr)e-ikz1z+Seikz1|z-zs4ikz1(2)式中：A+1()kr、A-1()kr分别表示上行波、下行波的波谱数；kz1=k21-k2r，k1是水层的波数；z是接收深度，zs是声源深度；S是声源强度。文献12中指出当水平波数kr比较大时，在包含声源的介质层中，式(2)中的声源项具有最慢的衰减速率，因而G1(kr;f)的衰减速率主要由声源项决定。以下讨论G1(kr;f)中声源项的衰减特性。用s表示式(2)中声源项的对数幅值，表达式为s(kr;f)=ln|Seikz1|z-zs4ikz1|(3)令=krk1，则kr=k1。当krk1时，1，且kz1=ik2r-k21。把上述kz1和kr代入式(3)中并化简可得：s(;f)=-k1|z-zs|2-1-ln2-1+ln|S4|(4)534第 4 期李承帮等：声场波数积分截断波数自适应选取方法当1时，-12-1，将其代入式(4)，可得：-k1|z-zs|-ln+ln|S4|s(;f)1时，声源项的对数幅值s(;f)随的变化关系近似满足一次线性项与对数项的叠加。在海底为液态且没有冰盖层的海洋波导中，波数域格林函数G(kr;f)也会呈现与Pekeris波导类似的衰减特性。令参考波数k0(f)为水平分层介质中最小声速对应的波数：k0(f)=2fminc1cmcM(7)其中：cm表示第m层介质中的声速。当krk0(f)时，任意一层介质的垂直波数kzm=ik2r-k2m都为纯虚数。与 Pekeris 波导中的情况类似，当krk0(f)时，G(kr;f)会随kr的增大而减小，并趋向于0。对G(kr;f)的幅值取对数，并且把水平波数kr与参考波数k0(f)的比作为自变量，令：(;f)=lg|G(kr;f)|(8)其中，水平波数比的表达式为=krk0()f(9)(;f)能较好地反映G(kr;f)的衰减特性。以图1中的波导环境为例，首先计算出不同频率的波数域格林函数G(kr;f)，然后计算出对应的(;f)，计算结果如图 2 所示。可以看出当1，即krk0(f)时，G(kr;f)开始衰减。频率越高，G(kr;f)衰减得越快。设期望的积分截断精确度参数为a，表示期望在(;f)=-a处 进 行 积 分 截 断。当a分 别 取3456时，计算在期望截断处t的取值随频率变化，结果如图3所示。结合图2的分析结果，可以看出在较低频段，截断波数比t需要取远大于1的值才能获得足够的积分精确度。其原因是在低频段(;f)衰减得比较慢，而在较高的频段(;f)衰减得比较快，截断波数比t取稍大于1的值就能获得足够的积分精确度。目前常用的截断波数选取方法是令t等于某个固定的值。在宽带声场计算中，若t取一个较大的值，可保证在低频声场计算时获得足够的积分精确度，但会给高频声场计算增加不必要运算量；若t取定一个较小的值，可减少高频声场的计算量，但会增大低频声场计算误差。图3 截断波数与参考波数的比随频率变化Fig.3 The ratio of truncated wavenumber to reference wave-number versus frequency图1 海洋波导环境特例Fig.1 A special case of marine waveguide environment50 Hz100 Hz200 Hz300 Hz400 Hz500 Hz600 Hz700 Hz800 Hz900 Hz1 000 Hz图2 不同频率波数域格林函数对数幅值Fig.2 Logarithmic magnitude of Greens function in wave-number domain at different frequencies5352023 年声学技术1.2波数域格林函数的模型拟合由1.1节的分析，在海底为液态且没有冰盖层的海洋波导中，当krk0(f)时，G(kr;f)开始衰减并趋向于0。在保证精确度的条件下，为了减少计算量，截断波数的取值应尽可能小，所以应当更关注当kr不太大时G(kr;f)的衰减特性。对于格林函数中的声源项，根据本文1.1节中式(4)(6)，当1时，s(;f)随的变化近似满足对数项与线性项叠加的关系，且当1时，s(;f)；对于格林函数中的齐次项，当kr的值不太大时，其影响也需要考虑，但由于齐次项包含数值计算结果A+1(kr)和A-1(kr)，难以对其变化特性进行理论推导分析。图2中格林函数对数幅值(;f)的数值计算结果显示，在12的范围内，当接近于 1 时(;f)衰减速率较快；随着的增大，(;f)的衰减速率逐渐减慢。(;f)的这种变化特征与声源项s(;f)所满足的对数和线性变化特征比较相似。另外，注意到式(2)中格林函数每一项都包含kz1，因而(;f)中每一项都包含2-1，且=1是(;f)的一个极点。而在数值计算中，=1处格林函数的计算结果不是无穷大，且极点的位置可能会存在偏差。因此，构造数学模型时，不应完全按照格林函数表达式把=1确定为极点，而应该把该极点的位置作为一个可调整的参数，否则数学模型可能与格林函数的数值计算结果失配。基于上述对格林函数衰减特性理论和数值计算结果的分析，并参照式(5)和式(6)中s(;f)的近似结果，构造了以下包含对数项和一次线性项的数学模型：h1();=1lg()2+3+4+5(10)或h();=1lg()+2+3+4(11)其中：=1 2 3 4 5T和=1 2 3 4T为模型参数，t为截断处的水平波数比，1t。根据对数函数的性质，可证明式(10)和式(11)是等效的，而式(11)的参数更少，因此采用式(11)作为格林函数的拟合模型。在该模型中可以通过调整参数2去调整极点的位置，避免了因=1附近极点位置的偏差而导致的模型失配。利用此模型通过最小二乘法对格林函数对数幅值进行拟合：=argmin=1th();-();f2(12)图 4 是利用式(11)中的数学模型对不同频率(;f)的拟合结果，图中粗虚线对应图2中各频率的波数域格林函数对数幅值(;f)，实线是拟合结果。所有频率的(;f)与对应的拟合曲线几乎重合，这表明利用式(11)对(;f)进行拟合取得了较好的效果。本文构造的式(11)中的数学模型符合格林函数的变化趋势，具有一定的物理意义，可调整极点的位置，避免模型在极点处失配，并且复杂程度较低。因此，本文将采用式(11)中的数学模型对格林函数对数幅值(;f)进行拟合。1.3模型参数自适应滤波卡尔曼滤波器是一种最优线性状态估计方法，被广泛应用于工程领域，本文利用卡尔曼滤波器对波数域格林函数的拟合参数进行跟踪和预测。1.2节构造的数学模型h(;)中包含4个参数j(n)，j=1234，n是频率序号。这4个参数之间是相互独立的，用4个卡尔曼滤波器分别对其进行跟踪和预测。式(11)中第j个参数的状态方程为sj(n)=Asj(n-1)+uj(n)(13)其中：状态sj(n)=j(n)j(n)T，j(n)表示j(n)对频率的偏导数；uj(n)为系统噪声，设其协方差矩阵是Qj；状态转移矩阵A=1Df01，Df为频率间隔。设状态误差协方差矩阵为Mj=E(sj-sj)(sj-sj)T，其中sj表示状态矢量的真实值，E表示对变量求期望。第j个参数的观测方程为zj(n)=Hsj(n)+vj(n)(14)(f)h(f)图4 波数域格林函数对数幅值拟合结果Fig.4 Fitting results of the logarithmic magnitude of Greens function in wavenumber domain536第 4 期李承帮等：声场波数积分截断波数自适应选取方法其中：zj(n)是状态的观测值；观测矩阵H=1 0T；vj(n)为观测噪声，设观测噪声的协方差矩阵是Rj。设第j个参数、第n步的卡尔曼增益矩阵为Kj(n)。卡尔曼滤波的每个迭代过程，都包含了预测与修正两个环节，在这两个环节中都会对状态矢量sj和状态误差协方差矩阵Mj进行更新。为了对两个环节的变量加以区别，用sj(n|n-1)表示根据第n-1步状态矢量所预测的第n步状态矢量，Mj(n|n-1)表示根据第n-1步状态误差协方差矩阵所预测的第n步状态预测误差协方差矩阵；用sj(n|n)表示第n步修正后的状态矢量，Mj(n|n)表示第n步修正后的状态误差协方差矩阵。对每个模型参数j(n)分别进行卡尔曼滤波，其具体计算为13(1)初始条件 sj()0|0=Esj()0Mj()0|0=Esj()0 sTj()0(15)(2)预测sj(n|n-1)=Asj(n-1|n-1)(16)(3)状态预测误差协方差矩阵Mj(n|n-1)=AMj(n-1|n-1)AT+Qj(17)(4)卡尔曼增益Kj(n)=Mj(n|n-1)HTRj+HMj()n|n-1 HT-1(18)(5)修正sj(n|n)=sj(n|n-1)+Kj(n)zj()n-Hsj()n|n-1(19)(6)状态误差协方差矩阵Mj(n|n)=I-Kj()n H Mj(n|n-1)(20)其中：I为单位矩阵。1.4截断波数的计算设期望的积分截断精确度参数为a，表示期望在波数域格林函数的幅值衰减到10-a时进行积分截断，即：|G(kt;f)|=10-a(21)其中：kt为截断波数。把自适应滤波器输出的第n+1个频率的参数估计结果(n+1)代入数学模型式(11)中，可估计出第n+1个频率的波数域格林函数：h;()n+1();fn+1=lg|G()kr;fn+1(22)把式(21)代入式(22)中可得第n+1个频率的截断波数t(fn+1)即为式(23)的解。ht;()n+1+a=0(23)由式(11)可知方程(23)是超越方程，需要通过数值方法求解。1.5截断波数自适应选取流程图5是本文提出的截断波数自适应选取方法的流程图。本文提出的截断波数自适应选取方法包含以下步骤：(1)初始截断波数t(f0)，t(f0)的取值可以较大，以保证声场计算精度。(2)利用波数积分法计算声场和(;fn)。(3)利用式(11)中的数学模型对进行(;fn)拟合，得到参数(n)。把上一次迭代自适应滤波器估计的参数(n)作为求解此最优化问题的初值，以加快拟合速度。(4)把参数(n)输入卡尔曼滤波器，估计下一个频率的模型参数(n+1)。(5)把(n+1)代入方程(23)并求解此超越方程，得到下一频率的截断波数t(fn+1)。求解此超越方程时，把上一次迭代的求出的截断波数t(fn)作为初值，以加快求解速度。(6)返回步骤(2)计算下一个频率的声场。图5 截断波数自适应选取流程Fig.5 The flow chart of adaptive selcetion process of truncat-ed wavenumber5372023 年声学技术2仿真验证计算图1所示海洋波导的波数域格林函数，并利用本文提出的方法选取积分截断波数，以验证本文提出方法的有效性。初始截断波数取值是t(f0)=5，计算频率范围是501 000 Hz。卡尔曼滤波器的参数如表1所示。对精确度参数a取不同值、频率间隔Df取不同值的情况进行试验。图6是Df分别为1、5和10 Hz 时，a分 别 为 3，4，5，6 的 情 况 下，(t(f);f)随频率变化，其中t(f)为本文方法选取的截断波数。图6中的结果表明在上述试验条件下，本文提出的方法都能自适应地选取积分截断波数，实现在期望的精确度附近进行积分截断。在较低频段，实际截断处幅值t()f;f会出现波动，随着频率的增大，波动逐渐减小。这是因为确定截断精确度a后，当频率较低时，(;f)随频率f的变化比较剧烈，因而模型参数的预测误差比较大；当频率较高时，(;f)随频率f的变化比较缓慢，因而模型参数的预测误差较小，截断精确度误差变小。图6说明了本文提出的自适应截断波数选取方法能在宽带声场波数积分计算过程中，把截断处的格林函数幅值稳定地保持在某一个精度附近。作为对比，图7给出了当t取固定值时，声场计算过程中截断处波数域格林函数的对数幅值随频率变化。可以看出当t取固定值时，声场计算过程中截断精确度随频率的增大而不断提高。为了评估本文方法对声场计算速度的影响，在其他条件都相同的情况下，分别利用常用的t取固图6采用本文方法计算声场时截断处波数域格林函数的对数幅值随频率变化图Fig.6Logarithmic magnitude of Greens function in wave-number domain at the truncation versus frequency when the method proposed in this paper is used for sound field calculation表1卡尔曼滤波器参数的取值Table 1Values of Kalman filter parameterssj(0|0)j(0)0TMj(0|0)10-50010-4Qj 10-50010-4Rj10-8t=1.1t=1.2t=2.0t=3.0图7 采用t取固定值法计算声场时截断处波数域格林函数的对数幅值随频率变化图Fig.7 Logarithmic magnitude of Greens function in wavenumber domain at the truncation versus frequency when the fixed t method is used for sound field calculation538第 4 期李承帮等：声场波数积分截断波数自适应选取方法定值的方法和本文提出的自适应截断波数选取方法，计算501 000 Hz宽带声场100次，对比两种方法的平均用时。为了量化对比两种方法的积分截断精确度，计算宽带平均精确度amean：amean=-lg1Nn=1N10t()fn;fn(24)其中：n为频率序号，N为频点个数，t()fn;fn为截断处格林函数对数幅值。表2和表3分别为采用t取固定值的方法和本文提出方法进行声场计算的平均用时和带宽平均精确度。在平均用时比较接近的情况下(如表2中t=1.1 1.2，表3中a=3,4,5,6的情况)，采用本文方法选取截断波数能获得更高的宽带平均精确度；在宽带平均精确度比较接近的情况下(如表2中t=2.0和表3中a=3的情况)，采用本文方法选取截断波数所需平均时间更少。这是因为t取固定值的方法无法兼顾低频精确度和高频计算量。由图7可知，当t取固定值1.1或1.2时，虽然计算量不大，但低频声场计算的积分截断精确度太低，导致宽带平均精确度较低；而当t取固定值2.0或3.0时，虽然声场计算整体的积分截断精度足够高，但高频声场计算量较大，导致计算的平均用时较大。本文提出的方法在宽带声场计算过程中，积分截断精确度可稳定地维持在某一个值上，在保证了低频积分截断精确度的同时不会给高频声场计算增加不必要的计算量。3结 论由于频率、介质类型、波导参数、声源与接收点的相对深度等参数会对波数域格林函数的衰减特性产生很大的影响，在采用波数积分法计算不同海洋波导的声场时需要选取不同的积分截断波数。针对波数积分宽带声场计算过程中，现有积分截断波数选取方法所存在的问题，本文提出一种应用于流体介质声场计算的积分截断波数的自适应选取方法。利用数学模型拟合波数域格林函数的衰减特性，根据已计算的波数域格林函数预测下一频率的波数域格林函数的衰减特性，然后估计出积分截断波数。仿真试验结果表明，本文提出的方法实现了在给定的精确度附近进行积分截断，在保证积分精确度的同时不会引入太多额外的计算量。本文提出的波数积分法截断波数自适应选取方法除了可以用于宽带多频点的声场计算外，还可以拓展到其他参数渐变的声场计算场景。本文关注的是各层都是流体介质的海洋波导的情况，当海底是弹性介质或波导中包含冰盖层时，波数域格林函数的衰减特性会有所变化。能否把本文提出的截断波数自适应选取方法，应用到包含弹性介质的海洋波导的波数积分声场计算，有待进一步研究。构造函数模型拟合图3中的曲线，可以得到给定精确度时截断波数随频率变化的经验公式，这可能是解决此问题的另一种思路。但图3中的曲线会随波导参数的改变而改变，如何建立曲线参数与波导参数之间的关系有待深入研究。参考文献1PEKERIS C L.Theory of propagation of explosive sound in shallow waterJ.Geological Society of America Memoirs,1948:1-116.2DINAPOLI F R,DEAVENPORT R L.Theoretical and nu-merical Greens function field solution in a plane multilayered mediumJ.The Journal of the Acoustical Society of Ameri-ca,1980,67(1):92-105.3SCHMIDT H,JENSEN F B.A full wave solution for propaga-tion in multilayered viscoelastic media with application to Gaussian beam reflection at fluid solid interfacesJ.The Journal of the Acoustical Society of America,1985,77(3):813-825.4SCHMIDT H,TANGO G.Efficient global matrix approach to the computation of synthetic seismogramsJ.Geophysical Journal International,1986,84(2):331-359.5骆文于,于晓林,张仁和.一种可稳定计算Pekeris波导中声场的波数积分方法J.声学学报,2016,41(3):321-329.LUO Wenyu,YU Xiaolin,ZHANG Renhe.A wave number 表3本文提出的方法声场计算平均用时和宽带平均精确度Table 3Average time and broadband average accuracy of sound field calculation with the method proposed in this paper期望精确度a3456平均用时/s140.0141.4142.9144.5宽带平均精确度2.996 83.994 94.993 15.991 4表2水平波数比取固定值方法声场计算平均用时和宽带平均精确度Table 2Average time and broadband average accuracy of sound field calculation with the fixed horizontal wavenumber ratio method水平波数比t1.11.22.03.0平均用时/s141.7145.5175.4216.8宽带平均精确度0.910 51.356 13.053 94.446 25392023 年声学技术integration method for stable calculation of sound field in Pe-keris waveguideJ.Acta Acustica,2016,41(3):321-329.6LUO W Y,YU X L,YANG X F,et al.Analytical solution based on the wavenumber integration method for the acoustic field in a Pekeris waveguideJ.Chinese Physics B,2016,25(4):198-209.7FINN B J,WILLIAM A K,MICHAEL B P,et al.Computa-tional Ocean AcousticsM.2nd ed.New York,USA:Spinger,2011:126,262.8FINN B J,WILLIAM A K,MICHAEL B P,et al.Computa-tional Ocean AcousticsM.2nd ed.New York,USA:Spinger,2011:323.9SCHMIDT H.SAFARI:seismo-acoustic fast field algorithm for range-independent environments.users guideR.SACLANT Undersea Research Centre,1988:31.10 SCHMIDT H.SAFARI:seismo-acoustic fast field algorithm for range-independent environments.users guideR.SACLANT Undersea Research Centre,1988:36.11 刘巍,肖汶斌,程兴华,等.声场波数积分最大截止波数自动选取算法J.国防科技大学学报,2019,41(4):177-181.LIU Wei,XIAO Wenbin,CHENG Xinghua,et al.Automatic selection algorithm for maximum truncate wavenumber of acoustic wavenumber integration methodJ.Journal of Na-tional University of Defense Technology,2019,41(4):177-181.12 FINN B J,WILLIAM A K,MICHAEL B P,et al.Computa-tional Ocean AcousticsM.2nd ed.New York,USA:Spinger,2011:261.13 STEVENM.KAY.统计信号处理基础:估计与检测理论(卷、卷合集)M.罗鹏飞,张文明,刘忠,等,译.北京:电子工业出版社,2014.540