浅海水平变化波导下低频声能量传输特性

祝捍皓，肖 瑞，朱 军，唐 骏

(1.浙江海洋大学 海洋科学与技术学院，浙江 舟山 316022；2.浙江大学 浙江省海洋观测-成像试验区重点实验室，浙江 舟山 316021；3.中国科学院声学研究所 声场声信息国家重点实验室，北京 100190；4.浙江海洋大学 船舶与机电工程学院，浙江 舟山 316022；5.挪威科技大学 声学组，挪威 特隆赫姆 7491)

浅海中的声传播问题是水声领域研究的热点之一，是理解、预测及应用各种浅海声学现象的基础[1-3].随着潜艇隐身技术的发展和对各类物理海洋现象的监测需求，各类浅海水下声学探测手段逐渐转向低频(小于1 kHz)，这显然与低频信号在浅海中的传播问题紧密联系，因而对浅海低频声传播问题的研究获得了越来越广泛的关注[4-6].但目前在研究浅海环境下的声传播问题时，多数只关注海水层中声场的计算建模和传播规律，且通常假设海底为水平分层的液态介质，而有关考虑海底为弹性介质下地形变化复杂的浅海低频声场建模的研究则相对较少[7].

浅海海底结构一般可视为由近孔隙介质的沉积层和近弹性介质的基岩层构成，沉积层厚度一般在10 m左右[8-10].由于低频水声信号波长长且有较强的穿透力，其在浅海环境中传播时能穿透沉积层到达基岩层，所以，在研究浅海低频声传播问题时将海底视为弹性介质更为合理.此外，实际浅海海底受内陆河流冲刷及大陆板块挤压等因素影响，其海底地形、海底底质均存在较大差异.以海底地形为例，有的海域较平整，有的则成斜坡状，有的甚至为陡峭的山状强切割地形.若仍采用传统水平分层的浅海波导模型将无法准确描述低频声信号在浅海中的传播特性，更无法为浅海水声工程应用提供合理的理论参考.因此，在对浅海低频声传播问题的研究中，亟待开展在浅海水平变换波导下的讨论.

任何声传播问题的研究基础是合理的声场计算方法，经过近五、六十年的发展，目前在水声传播研究领域已建立起多种声场计算方法，如简正波方法、射线方法、抛物方程方法、快速场方法以及一系列衍生计算方法[11].但上述方法都是在对波动方程和海洋环境作出一些假设和近似的基础上推导得到，其通用性均受到了一定的限制，尤其是对浅海水平变化波导下低频声传播问题的计算[12].随着对浅海水平变换波导下低频声传播问题的日益关注，积极建立通用性更强的浅海声场计算方法日益重要.有限元法(FEM)[13]把环境划分为离散的单元进行求解，因而更适用于描述声场中的各类变化，但由于其计算量大，以往除用于提供基准解外，一般在海洋声传播问题中很少使用.近年来，计算机硬件水平的快速提高为FEM在海洋声传播问题中的应用提供了有效的支持.

基于上述研究现状，针对浅海水平波导下低频声信号的传播特性的研究需求，本文采用COMSOL Multiphysics(以下简称COMSOL)[14]有限元软件平台建立了浅海低频声传播的有限元计算模型.针对当前浅海声场建模研究中多数只关注海水层中声场的计算建模和传播规律且通常假设海底为液态介质的研究现状，以声能流为研究对象[15]，同步讨论低频声波在浅海水体/弹性海底中的传播变化规律，以期为低频声信号在浅海中的探测应用提供理论基础.

1 理论概述

1.1 声场的有限元形式

针对本文所关注的浅海环境下的低频声传播问题，在三维柱坐标系下建立了符合浅海环境特点的波导模型.如图1(a)所示，模型中深度为zs的简谐点声源位于柱坐标对称轴上，海水层与海底层分别近似为均匀各向同性的流体介质和弹性体介质；ρw、cw分别为流体层中密度和声速；ρb、cp及cs分别为海底层中的密度、纵波声速和横波声速；αp、αs分别为海底层中的纵波声速衰减和横波声速衰减；海水层深度H(r)随水平距离r变化；Sw为海水层，在模型中表示为均匀流体；Sb为海底层，在模型中表示为均匀弹性介质.由于柱坐标系的轴对称性，可以将三维问题转化为二维(r,z)平面上求解[14].

图1 浅海声场模型Fig.1 Acoustic field model in shallow water

有限元法求解过程中的定义域如图1(b)所示.图中：L1边界为柱坐标中心对称轴边界；L2为海水层与浅海海表交互作用边界；L3为海水层与弹性海底层的交互作用边界；模型四周为完美匹配层(Perfectly Matched Layer，PML)[13].

对Helmholtz方程进行权重积分并结合高斯理论，Sw中的声学有限元方程以及Sb中的有限元振动方程可表示为[13-14]

(Kw+iωCw-ω2Mw)pi=[Fwi]

(1)

(Kb+iωCb-ω2Mb)si=[Fbi]

(2)

式中：Kw、Cw及Mw分别为刚度矩阵、阻尼矩阵及质量矩阵；ω为模型中简谐点声源发出声波的角频率；pi为海水层节点i处的声压；Fwi为声学激励；Kb、Cb及Mb分别为弹性海底层中上没有受到边界约束部分的刚度矩阵、阻尼矩阵及质量矩阵；si为海底层节点i处的位移；Fbi为结构体上的激励载荷.

在计算时，设定L1为Dirichlet边界条件，即边界L1上声压p满足恒等于0的物理假设，其表达式为

p=0,在L1上

(3)

在海水层与弹性海底层的相互耦合边界L4上建立声-弹耦合方程，即满足法向位移、法向应力连续，切向应力为0，其表达式为

(4)

(5)

(6)

为模拟声信号在无限大的海洋环境中的传播，在模型四周增加完美匹配层(Perfectly Matched Layer，PML)[13]，通过在有限元方程中增加吸收系数转化为PML方程.

通过耦合上述有限元法建立的物理场方程、流-弹耦合方程和PML方程，便可实现对浅海环境下海水层中声压p和海底层中位移场s的计算，通过p、s与各层其他物理场间的对应关系，可进一步得到模型中如振速场v、应力场T等其它物理场数值[15].在本文研究中，上述计算过程将在有限元软件COMSOL平台上实现.

1.2 波导中声能流的计算

在传统标量声场的建模中，对声场能量的研究多是基于传播过程中声压的变化特性展开.相比声压，声能流更适合揭示声波能量传播中的一般性“流动”特征且更适合用于对水体/海底中声能量特性的统一讨论.因此本文研究中，也将借助声能流密度这一矢量对低频声信号在浅海环境下的传播特性展开讨论.在本文所定义的(r,z)平面内，频域上海水层与弹性海底层中的声能流Iw、Ib定义为

(7)

(8)

(9)

(10)

下文以In为研究对象，对浅海水平变化波导下的低频声信号传播特性进行研究.

2 有限元法声场计算结果验证

为验证FEM对海洋声传播问题计算结果的正确性，本节首先使用FEM仿真计算了水平[15]及楔形上坡[11]两类海底地形下的声压传播特性，并与已有声场计算程序所得结果进行对比.所设定两类类海洋环境模型与参数示意图如图2所示.

两类模型中，均设定水体中cw=1 500 m/s，ρw=1 000 kg/m3.图2(a)为一类水平分层海底模型，水深在r=0至r0=4 km (r0为模型最远距离)处恒为100 m.图2(b)上坡海底海洋模型设定为美国声学学会(Acousticcal Society of America，ASA)提出的楔形海域检验模型，即在水深在r=0至r0=4 km处由200 m线性减小为0，上坡倾斜角度约为α1=2.68°.两类模型中均考虑海底为弹性介质，图2(a)中ρb=1 500 kg/m3，cp=2 000 m/s，cs=1 000 m/s，αpλ=αsλ=0.1 dB (λ为波长)；图2(b)中ρb=1 500 kg/m3，cp=1 700 m/s，cs=800 m/s，αpλ=λαs=0.5 dB.图2(a)中zs=20 m，频率100 Hz.图2(b)中zs=100 m，频率25 Hz.

图3所示为分别采用FEM以及现有声场计算方法模拟得到的图2中两类海洋环境下的单接收点(z=30 m)声压传播损失(Transmission Loss，TL)对比图[13]，TL计算公式为

(11)

式中：pref│r=1 m为距离声源1 m处的参考声压[1].

图3中虚线为FEM计算结果，实线对应现有声场计算方法计算结果.计算中，根据环境模型与现有声场计算方法适用性，分别选择快速场程序(Scotter)和抛物方程程序(RAMs)进行求解[16].

可以看出，FEM与两种声场计算程序仿真计算得到的传播损失曲线完全一致，验证了本文所研究FEM方法在计算水声场结果的准确性.从仿真结果也可以看到，在声源频率及声源深度一致的前提下，海底地形、海底声学参数的变化将对波导中声能量的传播产生重要影响，而海底地形在水平方向上起伏、变化又广泛存在于实际浅海环境中.因此，深入研究海底地形变化对低频声信号的影响规律，才能准确将低频声传播特性用于水声工程实践应用中.

针对上述目标，在本节研究基础上，下一节将通过数值模拟进一步讨论海底地形变化对低频声传播的影响规律，以期为低频声信号的应用提供理论基础.考虑FEM在使用中对波动方程不做任何假设和近似，可以处理各种几何形状下的物理场问题，因此更适用于对浅海实际水平复杂变化波导环境下的计算，下文研究中在COMSOL平台上利用FEM方法开展仿真研究.

图2 两类海底地形下海洋环境模型示意图Fig.2 Schematic diagram for marine environment models of two types of seabeds

图3 接收深度z=30 m时两类海洋环境下TL曲线对比图Fig.3 TL versus r at z=30 m in two marine environments

3 海底地形变化对声信号传播特性的影响

3.1 楔形海底

图4～6分别为当海底声学参数不变，海底地形分别为水平、上坡及下坡时水下20 m处100 Hz点声源激发起声能流在波导中的传播特征，声能流强度用声强级SIL表示，其计算式为

Iref=6.76×10-19(W/m2)

(12)

式中：Iref为参考声强[1].

以下图中各点声能流的大小和传播方向分别用箭头的长度和方向来表示.除地形变化外，图4～6其余仿真条件均与图2(a)中一致，ΔH1为水平方向4 km处海面与海底的垂直距离.

图4 水平海底仿真环境下声能流传播图(ΔH1=100 m)Fig.4 Distribution of acoustic energy flux in simulation of horizontal seabed (ΔH1=100 m)

图5 楔形上坡海底仿真环境下声能流传播图Fig.5 Distribution of acoustic energy flux in simulation of up-sloping wedge-shaped seabed

图6 楔形下坡海底仿真环境下声能流传播图Fig.6 Distribution of acoustic energy flux in simulation of down-sloping wedge-shaped seabed

在图2(a)仿真条件下，100 Hz声波在其水体将激发起共9阶波导简正波，其对应在水平海底界面上的掠射角θ∈[4°,40°][1].对比分析图4和5的仿真结果可以看到，相比水平海底，由于受α1的影响，水体层中i′次反射后声波与海底接触时的掠射角将增加(2i′-1)α1，原本小掠射角对应的低阶简正波传播特性不断向大掠射角、易衰减的高阶简正波耦合，声能量向海底的“泄漏”效应显著增强且起伏跨度缩短.α1越大，低阶简正波向高阶简正波的耦合越多，高阶声能量向海底“泄漏”效应越强，对应水体层中总声能量相应的衰减也越快，传输至远场后所保留的简正波阶数也越少，干涉特性越简单.在图5仿真条件下，随着ΔH1从100 m变化至0，水层中100 Hz声能量传输至r=4 km处衰减为20 dB左右.

对比分析图4和6可以看到，由于有α2的存在，声源位于浅水处，随着水平传播距离的增大水深不断增加，水体层中i′次反射后声波与海底接触时的掠射角将减小(2i′-1)α2，与楔形上坡海底的影响规律相反，原本大掠射角对应的高阶简正波传播特性不断向小掠射角的低阶简正波耦合，声能量向海底的“泄漏”效应明显减弱且起伏跨度增大.α2越大，高阶简正波向低阶简正波的耦合越多，声能量向远处传播时的衰减越慢，声场的干涉特性也越简单，即楔形下坡海底的存在对声场形成了增强效应.

3.2 小型海底山，小型海沟

图7 海底山、海沟地形下海洋环境模型示意图Fig.7 Schematic diagram of marine environment model of seamount and trench

图8和9分别为当海底声学参数不变，海底存在小型海底山与小型海沟时100 Hz点声源激发起声能流在波导中的传播特征，仿真中设定海底山与海沟中心位置位于水平方向2 km处.除地形变化外，图8和9中其余仿真条件均与图2(a)中一致.

图8中设定在海底水平距离1.5～2.5 km处存在一小型海底山.

图8 海底山仿真环境下声能流传播图Fig.8 Distribution of acoustic energy flux in simulation of seamount

对比分析图4和8的仿真结果可以看到，相比图4水平海底环境下的声能量传播特性，由于小型海底山的“阻挡”，整个仿真传播过程中声能量的衰减增加，且海底山越高，声能量的衰减越大.从图8的对比中可以看出，当海底存在海底山影响时，海底山前、后声能量的传播规律不相同.由于仿真条件下的海底山地形具有底边短、变化快的起伏特点，低频声波受其影响主要集中于海山前后，在点声源激发起的声能量到达海底山前受海底山地形的影响并不明显.随着声波不断临近海底山，受海底山前沿的反射作用，使部分原沿水平方向传播的声能量向海底山上方水层反射，声能量在海底山上方升高，即对应z<ΔH2处的声能量受海底山的“阻挡”作用并不显著.而同时在此过程中，部分声能量由海底山上坡处透射进入海底，该部分声能量主要对应掠射角较大，易衰减的中、高阶简正波，因此当声能量在“翻越”海底山的过程中其包含的中、高阶简正波能量更易快速衰减，在“翻越”海底山后只余少量低阶简正波存在，因此海底山后声能量起伏稳定且衰减较慢；而当对接收点深度z>ΔH2时，海底山对声能量传播产生的明显的“阻挡”作用，声能量衰减明显.

与图8仿真结果相对应，图9所示为当海底存在海沟时的声波传播特性的仿真图.仿真中设定在海底r=1.5～2.5 km处存在一海沟.

图9 海沟仿真环境下声能流传播图Fig.9 Distribution of acoustic energy flux in simulation of trenches on seabed

4 结论

本文基于有限元法，结合具体仿真算例讨论了浅海水平变化波导下声能量的传播变化特性及机理，以期为浅海水声工程设备的研发，测试及使用等提供一些理论参考.对比浅海水平分层波导下的声能量传播规律，结果表明：

(1)利用有限元法可以准确实现对各类复杂海洋环境下声场分布的求解.

(2)无论是楔形上坡或是下坡海底，二者均会对浅海中声能量的传播产生重要影响.相比水平海底，上坡海底倾角α1越大，高阶声能量向海底“泄漏”效应越强，声能量向远处传播时的衰减越快；而楔形下坡海底地形的影响则正好相反，倾斜角度α2越大，声能量向远处传播时的衰减越慢.正波的转移越多，声能量向远处传播时的衰减越慢，即楔形下坡海底的存在对声场形成了增强效应.

(3)当海底存在小型海底山影响时，由于海底山上坡面的反射作用，点声源激发起的声能量在“翻越”海底山的过程中，将在海底山上方产生提升，海底山后声能量起伏稳定且衰减较慢.

(4)与海底山的影响类似，在点声源激发起的声能量到达海沟前其受地形的影响并不明显，当海沟深度较浅，声能量在“跨越”海沟过程中将产生部分声能量的额外衰减.随着海沟深度的增加，当存在夹角大于波导中部分低阶简正波的掠射角时，声能量在海沟上方水层产生提升，随着海沟深度的继续增加，“跨越”后声能量的起伏特性与水平海底基本一致.

致谢在此感谢哈尔滨工程大学水声工程学院朴胜春教授课题组为本论文的撰写提供有限元软件支持.