一种水平变化波导中声传播问题的耦合模态法

刘娟 李琪

1) (哈尔滨工程大学， 水声技术重点实验室， 哈尔滨 150001)

2) (海洋信息获取与安全工业和信息化部重点实验室(哈尔滨工程大学)， 哈尔滨 150001)

3) (哈尔滨工程大学水声工程学院， 哈尔滨 150001)

1 引 言

水平变化波导中的声传播问题是水下声传播研究中的热点问题.受季节和日照的影响， 以及海洋中的中小尺度现象和内波的作用， 海洋环境参数空间分布不均匀.作为海洋波导的上下界面， 海面在风力驱动作用下随机起伏， 海底地形粗糙不平.声波在水平变化波导中传播会发生随机散射， 严重影响能量和信息的传输.本文针对介质参数和海底边界两种水平变化因素， 研究水平变化波导中的声传播问题.

耦合模态法是求解水平变化波导中声传播问题的常用方法之一， 由Pierce[1]和Milder[2]引入到水声学中.Pierce-Milder 耦合模态理论将水平变化波导中任意位置处的声压表示为本地模态的叠加， 模态系数与水平位置有关.通过对声压满足的Helmholtz 方程在本地模态上作投影， 将Helmholtz 方程转化为关于模态系数的二阶耦合模态方程组.波导环境的水平变化对声场的作用由两个耦合矩阵描述.本地模态的计算精度直接影响耦合模态法的求解精度.对于介质参数均匀的波导， 本地模态可解析计算.对于介质参数复杂的波导， 需借助数值手段求解， 如谱方法[3]、多模态法[4]、有限元法[5]、有限差分法[6]及标准简正波程序KRAKEN[7]等.Rutherford 和Hawker[8]指出， Pierce-Milder耦合模态理论在处理边界倾斜波导的声传播问题时能量不守恒.这是由于本地模态并不符合严格的边界条件.Fawcett[9]将正确的边界条件引入到耦合模态方程组的推导中， 提出一种能量守恒的耦合模态理论.但该理论不仅要计算两个耦合矩阵， 还需计算两个界面矩阵.这些矩阵与本地模态在水平方向上的导数有关， 只能近似求解.另外， 直接积分求解二阶耦合模态方程组可能会遇到数值发散的问题[4].本地模态和其水平导数的计算复杂度及二阶耦合模态方程的数值求解问题使得Fawcett的方法在实际中难以广泛应用.Abawi[10]通过忽略高阶耦合项和后向散射场的能量， 结合抛物方程方法推导得到耦合简正波-抛物方程理论， 有效解决了耦合模态方程的求解问题.在此基础上， 彭朝晖和张仁和[11]采用广义相积分理论和波束位移射线简正波理论计算本地模态， 实现了声场的快速计算.莫亚枭等[12]通过忽略高阶耦合项并只取前向场， 提出了一种大步长格式的二维耦合模态方法.这些单向耦合模态理论对于环境特性水平不变或缓变波导是正确的， 但不适用于后向场能量不可忽视的情况[13].针对求解二阶耦合模态方程中可能出现的数值发散问题， Pagneux 等[4，14]提出一种数值收敛的多模态导纳法， 分析了变截面声管中的声传播问题.多模态导纳法通过引入导纳矩阵将Helmholtz 方程的边值问题等价成两个一阶演化方程的初值问题， 采用Runge-Kutta 法[4，15]或Magnus 法[14，16，17]数值计算获得稳定声场解.目前研究海底地形变化的模态耦合理论已相对成熟， 而介质参数水平变化波导中的声传播问题仍有待研究.对于介质参数水平缓变(相比声波波长)的海洋环境， 绝热近似耦合模态理论[18]通过忽略模态间的耦合作用， 得到形式简单的耦合模态方程.然而， 对于垂直于内波波峰方向上的声传播问题， 环境特性会迅速变化， 绝热近似理论通常无效[19].因此， 有必要推导一种双向耦合模态方法， 考虑介质参数水平变化对模态耦合及声场的影响.

前文所述均为连续耦合模态理论(continuous coupled mode method)， Evans[20]将水平变化区域分为若干段水平不变波导， 提出阶梯耦合模态理论(stepwise coupled mode method).骆文于等[7，21]引入全局矩阵算法， 推导出一种数值稳定的阶梯耦合模态方法.阶梯耦合模态方法数值实现相对简单， 在水声领域有广泛应用.相比阶梯耦合模态理论， 连续耦合模态方法有两个优势: 1)给出了耦合矩阵的具体表达式， 直观地体现水平变化因素对模态间能量交换及声场的影响; 2)方法可以延伸应用于更实际的海洋波导中的声传播问题， 如: 三维模型[22]、含散射体[23，24]等复杂波导环境.

鉴于上述问题， 本文将结合Fawcett 耦合模态理论和多模态导纳法， 提出一种能量守恒且数值收敛的耦合模态方法， 研究介质参数和边界水平变化波导中的声传播问题.全文结构如下: 首先， 针对介质参数及边界水平变化波导， 推导基于多模态导纳法的耦合模态方法， 给出耦合模态方程及求解耦合模态方程的Magnus 数值积分方法; 其次， 使用该方法计算水平变化波导中的声场， 与COMSOL的计算结果比较， 验证方法的准确性; 最后， 给出讨论和总结.

2 理论方法

针对水平变化波导， 本节给出基于多模态导纳法的耦合模态方法(coupled mode method)， 下文中简称CMM.

2.1 双向耦合模态理论

考虑水平方向半无穷、垂直方向受限的二维平面波导， 如图1 所示.上边界为水平不变的空气-海水绝对软边界， 下边界为水平变化的海水-岩石绝对硬边界， 以连续函数 H (x) 表 示.声速 c (x，y) 和密度 ρ (x，y) 均 是空 间 坐标 的连 续 函 数.假 设 区 域Ω1:x ∈[0，xr]为水平变化区域， 声波在该区域中发生散射.区域 Ω2:x ＞xr为水平不变区域， 声波在该区域中全部透射.入射波从 x =0 处向右入射.

图1 水平变化波导示意图Fig.1.Configuration of a waveguide with range-dependent environments.

省略时间因子 e-iωt， 声压满足Helmholtz 方程

和边界条件

其中p 是声压， ω 是角频率， n 表示外法线方向.

根据耦合模态理论， 波导中任意位置处的声压可用本地模态叠加求和表示， 这里本地模态是指上下边界分别满足绝对软和绝对硬边界条件， 声速和密度分别等于本地声速 c (y;x) 和本地密度 ρ (y;x) 的水平不变波导中的简正波.而依据多模态理论， 本地模态可用一组正交完备的本地本征函数叠加求和表示， 因此， 波导中的声压可表示为

其中: Pn(x) 为 模态系数; φn(y;x) 为本地本征函数，代表上下边界分别满足绝对软和绝对硬边界条件.波导深度为本地深度 H (x) 的均匀(声速和密度均为常数)波导中的简正波.用本地本征函数φn(y;x)作基函数的优势在于 φn(y;x) 及其导数的解析表达式容易计算， 而真正的本地模态及其导数往往不易解析求解.φn(y;x) 的表达式为φn(y;x)=满足正交性φn(y;x)dy =δmn.注意(3)式中的求和号表示在数值实现中要对级数求和作截断处理其中N 表示截断数.对于水平不变波导， 求和项通常只需包含全部传播模态.对于水平变化波导， 由于衰逝模态对水平变化区域中近场声压的贡献不可忽略， 求和项中需包含全部的传播模态和部分衰逝模态[14].

代入声压展开(3)式， 得到二阶耦合模态方程组

写成矩阵形式

其中， 列向量P 中元素为 Pn(x) ， 矩阵A， B， C，D， E， F， G， J 和K 中的元素分别为

(6)式中的各个矩阵描述了水平变化因素对声场的影响.矩阵A， C 和K 是由于基函数不是真正的本地模态产生的系数矩阵.B， D， J 和G 是耦合矩阵， 描述模态间的耦合作用.E 和F 是边界矩阵， 来自于计算的过程中加入的严格下边界条件， 即 ∂yp(H)=H′∂xp(H) ， 目的是修正本地本征函数的边界条件， 使二阶耦合模态方程(6)满足能量守恒， 详细过程可参考文献[9].

至此， Helmholtz 方程引导的声传播问题被转化为关于模态系数P 的二阶耦合模态方程的边值问题.由于边界矩阵E 和F 修正了边界条件， (6)式已满足边界条件， 对(6)式的求解只需考虑声源条件和辐射条件.但根据辐射条件从右向左积分求解(6)式时， 由于衰逝模态在这个方向上是指数发散的， 导致数值计算结果不收敛， 故而(6)式并不是理想形式.根据多模态导纳法的思想， 将二阶微分方程形式下的耦合模态方程重构为两个一阶演化方程.定义列向量s 为

(7)式对x 作全微分， 有s′=AP′′+(A′+D)P′+D′P ， 矩阵 A′中的元素为

即A′=E+B+D+DT.矩阵 D′中的元素为

即 D′=F +G+J +L ， 其中L 中的元素为Lmn=因此， 有

根据(6)—(8)式， 有

其中用到了等式 L =DTA-1D.下面给出该等式的证明， 利用本地本征函数 φn的完备性， 将函数代 入 矩 阵D 中 有写成矩阵形式，得到 D =AWT， 其中 W 中的元素为 wnn′， 同理可得即L=W AWT， 因 此 DTA-1D =W ATA-1AWT，由于矩阵A 为对称阵， 有DTA-1D =L.

最终， 得到两个一阶耦合模态方程

(10)式可通过引入导纳矩阵获得稳定的声场解， 并且无需求解边界矩阵E 和F.导纳矩阵Y 代表模态域的Dirichlet to Neumann (DtN)算子[17，25]， 定义为 s =Y P.将导纳矩阵代入(10)式， 得到导纳矩阵满足的Riccati 方程

和声压的模态系数满足的一阶微分方程

(11)式的初值条件为区域 Ω2中的辐射条件Y(xr)=事实上， 矩阵A-1(xr)Y(xr)的本征值等于 i kxn， 其中 kxn对应波导在 Ω2区域中水平方向上的波数.(12)式的初值条件为 x =0 处的声压对应的模态系数.工程应用及数值实现中，声源通常为入射波， 而 x =0 处的声压包含入射波成分及反射波成分.为了计算反射波， 定义模态域的反射矩阵R 为 P-(0)=RP+(0) ， 其中 P+和P-代表前向传播和后向传播的模态系数分量， 有

根据(12)—(14)式， 得到

其中 P+(0) 为入射波的模态系数.本文中的导纳矩阵和声压模态系数均采用Magnus 数值积分算法求解， 详细步骤见2.3 节.

另外， 散射矩阵是描述散射区域性质的主要参量， 是连接入射信息和出射信息的重要桥梁.针对图1 所示模型， 散射矩阵包含反射矩阵R 和透射矩阵T.透射矩阵定义为 P+(xr)=T P+(0).引入传播矩阵Q， 满足 P (xr)=Q(x)P(x) ， 根据(10)式，有 Q′=-Q(-A-1D+A-1Y) ， 初值条件为Q(xr)等于单位阵.透射矩阵则可写为 T =Q(0)(I +R).

总而言之， 对于一个水平变化波导， 双向耦合模态方法(two-way CMM)通过辐射条件确定任意水平位置处的导纳矩阵， 进而计算散射区域对应的反射矩阵， 再根据入射条件推出声压的模态展开系数， 最后代入声压的展开式中得到波导中的稳态声场.

2.2 单向近似与绝热近似耦合模态理论

当波导环境参数水平连续缓变(相比波长)时，单向传播模型及绝热耦合模态理论能够近似求解声场[13].本节将给出基于多模态导纳法的单向近似(one-way CMM)和绝热近似耦合模态(adiabatic CMM)方法.

单向近似理论假设后向散射场的能量相比前向场的能量是可忽略不计的高阶小量， 等价于反射矩阵近似为零 R ≈0， 即 Y0≈Y (0) ， 表示任意位置处的导纳矩阵无需迭代计算， 近似为Y(x)=根据(10)式， 有

初值条件为入射波对应的模态系数 P (0)=P+(0).该单向耦合模态方法避免了从无穷远辐射条件迭代计算导纳矩阵的过程， 无需计算反射矩阵， 相比双向耦合模态方法提高了计算速度.

绝热近似理论假设模态间绝热耦合， 不发生能量交换.因此， 忽略耦合矩阵项， 得到水平变化波导中的绝热近似耦合模态方程

同样地， 初值条件为入射波对应的模态系数P (0)=P+(0).

2.3 Magnus 数值积分方法

本文通过Magnus 数值积分方法求解导纳矩阵及声压模态系数.Magnus 法是一种高效的数值积分算法， 即使是剧烈振荡的解函数也只需稀疏的离散点， 这种大步长的计算特性使得Magnus 法在高频条件下能够实现快速计算[14].Magnus 法计算导纳矩阵与声压模态系数的具体步骤如下.

首先将散射区域 Ω1中的横坐标从右向左离散为， 其 中，=0 ， 离散间隔为＜0.将(10)式重写为

(19)式的Magnus 数值解为

其中Magnus 矩阵Γj的形式取决于不同阶的Magnus 积分方法， 二阶Magnus 矩阵为

四阶Magnus 矩阵为

将 eΓj写为分块矩阵， 即

最终得到导纳矩阵和声压的模态系数

3 数值结果

本节将利用第2 节提出的耦合模态方法， 求解水平变化波导中分布源及点源声传播问题， 并与COMSOL 参考解比较， 验证耦合模态方法的正确性.针对水平缓变波导的声传播问题， 讨论单向近似及绝热近似耦合模态方法.数值验证双向耦合模态方法的能量守恒特性.最后针对浅海孤立子内波中的低频声传播问题， 分析模态间的耦合作用.

3.1 节至3.3 节中考虑如图1 所示的含声速和密度 剖 面 的 水平 变化 波 导， 密度 ρ (x，y) (单 位 为kg/m3)和声速 c (x，y) (单位为 m /s )均是空间坐标的函数， 表达式分别为

其中， xr为水平变化距离， ρ1和 ρ2分别代表密度垂直和水平变化率， c1和 c2分别代表声速垂直和水平变化率.上边界水平不变， 下边界的几何参数为

其中 h1代表海底倾斜率.

3.1 分布源和点源声传播问题

首先计算分布源激发声场.环境参数选取为ρ1=0.2 ， ρ2=0.2 ， c1=0.2 ， c2=0.4 ， h1=-0.1.H =200 m， xr=1600 m.声源为从 x =0 处向右入射的分布源 pi=εφ1(y;0)=sin(1.5πy/200) ， 其中ε=10 为入射波归一化系数.频率 f =20 Hz.在该频率下， 波导在 x =0 m 处有五阶可传播模态， 在水平不变区域中有三阶可传播模态.图2(a)为利用双向CMM((10)式)计算得到的声压幅值(Pa)分布， 截断数选取为 N =10 ， 图2(b)为使用COMSOL 计算得到的声压幅值(Pa)分布， 图2(c)为深度 y =20 m 处声压幅值沿x 轴的分布， 图2(d)为深度 y =50 m 处声压幅值沿x 轴的分布.可以看出水平变化区域中存在明显的后向散射效应， 模态间耦合作用剧烈， 两种方法的计算结果一致.

图2 水平变化波导中的声场(声源为分布源， 频率为20 Hz) (a) 利用CMM 计算得到的声压幅值分布， 截断数 N =10 ; (b) 利用COMSOL 计算得到的声压幅值分布; (c) 深度为20 m 处， 声压幅值的水平分布; (d) 深度为50 m 处， 声压幅值的水平分布Fig.2.Sound fields in a range-dependent waveguide (the source is a distributed source at 20 Hz): (a) Sound field computed by CMM where the truncation number is N =10 ; (b) sound field computed with COMSOL; (c) sound field distribution along x direction at depth 20 m; (d) sound field distribution along x direction at depth 50 m.

CMM 方法处理点源(相当于对称三维问题中的线源)声传播问题的关键在于用本地本征函数展开表示点源传播函数—Green 函数， 即

对于水平变化波导， 相应的模态系数g 可写为[24，26]

图3 所示为点源激发声场.水平变化波导的环境参数为 ρ1=0.2 ， ρ2=0.2 ， c1=0.2 ， c2=0.4 ，h1=-0.05.H =200 m， xr=1600 m.声源位于(0，10) m 处.频率为 f =20 H z.图3(a)为利用双向CMM((10)式)计算得到的声压幅值(Pa)， 截断数选取为 N =50 ， 图3(b)为使用COMSOL 计算得到的声压幅值(Pa)， 图3(c)为深度 y =71 m 处声压幅值沿x 轴的分布， 图3(d)为深度 y =101 m处声压幅值沿x 轴的分布.可以看出明显的后向散射作用.两种方法的计算结果一致， 说明双向CMM能够准确计算水平变化波导中的声场.两种方法的计算偏差来源于两方面: 1) CMM 在数值实现中对级数求和(3)式作截断处理， 导致的部分和与真值之间的误差; 2)本文中采用Clenshaw-Curtis 数值积分方法[26]计算系数矩阵及耦合矩阵， 离散点采用Chebyshev 插值点， 而COMSOL 计算采用三角形网格， CMM 与COMSOL 离散方式不同导致两种结果出现偏差.

图3 水平变化波导中的声场(声源为位于 ( 0，10) m 处的点源， 频率为20 Hz) (a) 利用双向CMM 计算得到的声压幅值分布，截断数 N =50 ; (b) 利用COMSOL 计算得到的声压幅值分布; (c) 深度为71 m 处， 声压幅值的水平分布; (d) 深度为101 m 处， 声压幅值的水平分布Fig.3.Sound fields in a range-dependent waveguide generated by a point source at ( 0，10) m (the frequency is 20 Hz): (a) Sound field computed by CMM where the truncation number is N =50 ; (b) sound field computed with COMSOL; (c) sound field distribution along x direction at depth 71 m; (d) sound field distribution along x direction at depth 101 m.

3.2 水平缓变波导: 单向近似与绝热近似

水平缓变是指环境参数在一个波长内的水平变化量远小于波长， 对于本节的计算模型， 水平缓变代表 ρ1≪1 ， c1≪1 及 h1≪1.考虑水平缓变波导算例， 假设 ρ1=0.01 ， ρ2=0.2 ， c1=0.01 ， c2=0.4 ， h1=-0.001 ， H =200 m， xr=1600 m.声源为从 x =0 处向右入射的分布源pi=εφ1(y;0)=sin(1.5πy/200).频率 f =20 Hz.图4(a)—(c)分别为利用双向CMM 理论((10)式)、单向近似CMM理论((17)式)和绝热近似CMM 理论((18)式)计算得到的声压幅值(Pa)分布， 图4(d)和图4(e)分别为接收点在60 m 和120 m 深度处， 声压幅值的水平分布.三种方法的数值实现采用相同的离散点和截断数N = 10.可以看出， 单向近似和绝热近似耦合模态理论均可较为准确地计算水平缓变波导中的近距离声场， 但计算误差随水平距离的递推逐渐累积.

3.3 能量守恒验证

本节讨论数值声场解的能流守恒， 理论验证将在第4 节讨论部分中给出.根据能流定义

代入声压展开式(3)， 写成矩阵形式， 得到模态(本地本征函数)域的能流表达式

其中上标H 代表共轭转置.

图4 水平缓变波导中的声场(声源为分布源， 频率为20 Hz) (a) 利用双CMM 计算得到的声压幅值分布; (b) 利用单向近似CMM 计算得到的声压幅值分布; (c) 利用绝热近似CMM 计算得到的声压幅值分布; (d) 深度为60 m 处， 声压幅值的水平分布;(e) 深度为120 m 处， 声压幅值的水平分布Fig.4.Sound fields in a waveguide with weak range dependence generated by a distributed source (the frequency is 20 Hz): (a) Sound field computed by two-way CMM; (b) sound field computed with one-way CMM; (c) sound field computed with adiabatic CMM;(d) sound field distribution along x direction at depth 60 m; (e) sound field distribution along x direction at depth 120 m.

为直观起见， 考虑边界水平不变， 密度均匀，声速水平变化波导中的声传播.假设 ρ1=0 ，ρ2=0 ， c1=0.1 ， c2=0 ， h1=0 ， H =200 m.频 率f =10 Hz.在 x =0 m 处波导中有三阶可传播模态.随着声速的递增， 在水平不变区域中仅有两阶可传播模态.由于介质参数与y 无关， 本地本征函数即为实际本地简正波.声源为从 x =0 m 处向右入射的分布源 pi=εφ1(y;0)=sin(1.5πy/200) ， 即第二阶本地简正波.图5(a)为双向CMM ((10)式)计算得到的声压幅值(Pa)分布， 图5(b)为声压的模态系数 | Pn(x)| 的分布， 图5(c)为归一化能流E(x)/E(0).从图5(c)可以看出， 能流为常数， 代表数值解符合能量守恒.由于衰逝模态不传播能量，图5(b)中仅画出可传播模态对应的模态系数.入射波为第二阶模态， 即 P1，+(0)=1 ， 图中显示仅有第二阶模态在波导中传播， 表明声速的水平变化没有导致模态间的耦合.| P1(0)|=1.0005 ， 说明后向散射场的能量很小.从反射矩阵也可知，||RP||2～O(10-4)≪1.另外， 密度均匀时， 能流E(x)∝Re声速随传播距离增大， | P1(x)| 也随之增大， 符合能量守恒定律，与数值计算结果相符.

接着在图5 密度均匀、声速水平变化波导的基础上加入边界变化.假设 ρ1=0 ， ρ2=0 ， c1=0.1 ，c2=0 ， h1=-0.03.H =200 m， xr=1600 m.频率为 10 H z.在 x =0 m 处波导中有三阶可传播模态， 在水平不变区域中有两阶可传播模态.声源依旧为从 x =0 处向右入射的分布源pi=εφ1(y;0)=sin(1.5πy/200)， 即第二阶本地简正波.图6(c)显示数值解满足能量守恒.图6(b)说明虽然入射波为第二阶简正波， 但环境水平变化导致相邻位置的本地简正波之间发生耦合， 在 x =0 处前三阶模态在波导中同时传播， 随着声波的传播， 第三阶本地模态被截止， 最终在水平不变区域中仅有前两阶模态传播.此外， | |RP||2=0.0277 ， 相比边界不变波导(图5)， 边界水平变化较易产生不可忽视的后向散射作用.

3.4 浅海孤立子内波

孤立子内波是浅海中常见的声学现象.内波可能带来大尺度的水体垂直位移， 导致低频声信号在传播中产生波动[27].由于每个孤立子内波以各自的相速度传播， 多个内波间的相互作用复杂且难以描述内波对模态耦合的影响.因此， 以单个孤立子内波环境为例， 分析模态间的耦合作用.

考虑图7 所示的包含单个孤立子内波的浅海波导， 海面为绝对软界面， 海底为绝对硬界面.水体深度为60 m， 依据声速分布被分为三层， 上层为等声速层， 声速为 c1=1530 m/s， 中间为内波层，声速为

下层为等声速层， 声速为 c1=1490 m/s.各层中密度均匀 ρ =1000 kg/m3.内波的几何形状为

图5 密度均匀、声速水平变化波导中的声传播(声源为第二阶模态， 频率为10 Hz) (a)声压幅值分布; (b)声压的模态系数的水平分布; (c)归一化能流Fig.5.Sound propagation in a waveguide with constant mass density and range-dependent sound speed (the source is the second local mode， and its frequency is 10 Hz): (a) Sound field; (b) modal amplitudes distribution; (c) normalized energy flux distribution.

图6 密度均匀、声速及边界水平变化波导中的声传播(声源为第二阶简正波， 频率为10 Hz) (a)声压幅值分布; (b)声压的模态系数的水平分布; (c)归一化能流Fig.6.Sound propagation in a waveguide with constant mass density and range-dependent sound speed and boundary (the source is the second local mode and its frequency is 10 Hz): (a) Sound field; (b) modal amplitudes distribution; (c) normalized energy flux distribution.

图7 孤立子内波波导模型Fig.7.Configuration of a waveguide with an internal solitary wave.

这种类驼峰的内波形状是Korteweg-de Vries 方程的形式解， 描述弱非线性内波的传播[19].

研究波导中的模态耦合作用时， 常见且实用的手段是考虑单阶简正波激发的声场中每个水平位置处各阶简正波的分量.虽然本文提出的耦合模态法采用的展开基函数((3)式)不是真正的本地简正波， 展开系数也不是真正的各阶简正波的分量， 然而， 由(12)式的辐射条件可知， 矩阵的特征向量对应各阶本地简正波在本地本征函数 φn上的模态系数， 因此， 本地简正波及声场在本地简正波上的分量可通过正交基变换获得.据此， 我们数值计算得到图7 模型中 x =0 处的真正本地简正波， 示于图8.

图9(a)为入射第一阶本地简正波(即图8(a))，频率为100 Hz 时， 采用双向CMM((10)式)计算得到的声压幅值分布(Pa).图9(b)为声压在各个水平位置处各阶本地简正波的分量.显然地， 在声速水平不变区域中， 简正波之间没有耦合.在声速水平变化区域， 简正波之间产生耦合， 能量从第一阶简正波转移到高阶简正波中.另外， 连续水平变化的声速没有导致强烈的反射， 故而单向传播模型同样适用于求解该环境中声场.

4 讨论部分

首先理论验证双向CMM 满足能量守恒.假设 p (x，y) 和 p*(x，y) 为Helmholtz 方 程 的 声 场 解，也就是

图8 图7 所示模型中 x =0 m 处的简正波Fig.8.Local mode shape of the waveguide at x =0 m shown in Fig.7.

图9 浅海孤立子内波波导中的声传播(声源为第一阶简正波， 频率为100 Hz) (a)声压幅值分布; (b)前五阶简正波分量Fig.9.Sound propagation in a waveguide with an internal solitary wave (the incident wave is the first mode at 100 Hz): (a) Sound field; (b) modal amplitudes distribution.

(35)式 × p*- (36)式 × p ， 得到

(37)式对y 积分， 得到

代入声压展开式(3)， 根据(7)式， 有

因此，

(39)式对应声场互易性公式， (40)式对应能流守恒公式.可见， 双向耦合模态方法满足能量守恒定律.

本文提出的双向耦合模态方法为半解析半数值解法， 在推导过程中没有引入任何近似假设.然而， 水平变化因素—声速、密度及边界几何函数均默认为连续且可导.对于水平变化区域中边界导数不存在的情况， 可以通过保角变换[28]将复杂形状边界下的边界条件转换为水平不变边界下的边界条件， 将均匀Helmholtz 方程变为含有变化折射率的Helmholtz 方程， 再采用耦合模态法处理.该双向耦合模态方法可以用于研究水平分层波导中的声传播问题， 即介质参数不连续变化的波导环境， 只需恰当地加入交界面处的连续性条件.对于下层半无穷波导， 如Pekeris 波导， 第3 节中指出声速水平连续缓慢变化时， 后向散射场几乎可忽略， 因此可以通过在下层空间的声速上加入随深度逐渐增大的吸收项[29]， 近似分析深度半无穷波导中的声传播问题.

5 结 论

本文提出一种基于多模态导纳法的耦合模态方法以研究水平变化波导中的声传播问题.根据传统耦合模态理论和多模态法， 将每个位置处的声压用一组正交本地本征函数展开， 对Helmholtz 方程在本地本征函数上作投影， 推导出模态系数满足的二阶耦合模态方程组.为解决二阶耦合方程组数值计算发散的问题， 引入导纳矩阵， 将二阶耦合模态方程组简化为两个一阶演化方程， 并利用Magnus数值积分方法数值求解.利用该耦合模态方法数值求解水平变化波导中的声场， 与COMSOL 参考解比较， 结果吻合， 表明该耦合模态方法能够精确求解水平变化波导中的分布源及点源声传播问题.单向近似和绝热近似耦合模态方法可以近似求解水平缓变波导的声场.双向耦合模态方法在推导过程中没有任何近似， 计算误差来源于数值实现， 并且Magnus 数值积分方法具有大步长的计算特性， 满足能量守恒且数值解无条件稳定.总之， 双向耦合模态法方法充分考虑了波导环境对模态耦合的作用， 能够精确有效求解水平变化波导中的声传播问题， 可以用于求解实际声传播问题.