海洋环境噪声是海洋中的固有声场，其研究始于二战期间，利用环境噪声实测数据总结出的Knudsen谱^[1]和Wenz曲线^[2]现在仍然是水声工程中经常使用的背景场参考曲线。随着水声学的发展以及声呐技术对背景场需求的日益增加，人们在持续开展海洋环境噪声实验研究^[3-8]的基础上，着手海洋环境噪声的理论建模以及利用环境噪声数据反演海面、海底参数等技术^[9-11]进行研究。按照适用于海洋环境来分，目前海洋环境噪声模型主要分为2类：1)适于分层介质海洋环境的二维噪声模型^[12-14]，将海水介质沿深度分若干层，水平方向视为均匀，仅考虑海洋信道多种环境因素作用下不同分层相关参数的变化：如声速分布、海面和海底的反射以及介质的吸收等；2)考虑环境参数随水平距离或方位变化情况的三维海洋环境噪声模型^[15-18]。实际海底地形往往非常复杂，包括浅海大陆架、向深海延伸的大陆坡，以及海山、海沟、岛礁等，这些地形的起伏对海洋波导中的声传播具有重要影响^[19]。海洋环境噪声理论建模是基于所关注海域的环境噪声源模型和声传播模型，因此实际环境中模型的仿真需要深入了解海山、岛礁等复杂海底地形对海洋环境噪声的影响，关注复杂海底地形引起的海洋环境噪声建模的难度和仿真分析计算量的增加。

本文依据现有的海洋地形数据库，针对真实的复杂海底地形—海山/岛礁海域，采用射线三维声传播的N×2D和3D算法，建立了适于复杂地形海域的风关海洋环境噪声理论模型，研究复杂地形海域的海洋环境噪声特性。为了突出海底地形对风关海洋环境噪声特性的影响，假设风关噪声源随机均匀分布且声速剖面不随距离和方位变化。

1 复杂地形海域风关海洋环境噪声模型

由于射线法对高频声传播计算具有独特优势，计算简洁，适于求解深、浅海，特别是深海与距离有关环境的声场，结果直观、物理意义清晰明确，因此模型选取了射线三维声传播模型。利用射线声传播模型的N×2D和3D算法：采用柱坐标系，假设垂直接收阵位于z轴，将三维海洋环境按水平方位角划为N个扇区(见图 1)，计算每个位于方位扇区中心角度上声源至接收点的声场，叠加后获得接收点的总声场。N×2D和3D算法的区别在于后者考虑了方位分区之间的声耦合(即2个相邻垂直扇区之间能量交换)，因此3D算法精度高，适于计算涉及地形剧烈起伏或海洋水团变化等引起的声线水平折射问题，但3D算法的计算时间明显高于N×2D。若认为海洋环境在水平方位变化不太剧烈，声传播在水平方位上的耦合甚小于深度方向耦合，则采用N×2D算法更便捷、经济。

假设相互独立的风关噪声源随机均匀地分布在海面以下z_s深度无限大平面上(即将z_s深度的风关噪声源视为无指向性单极子声源)，且认为水平方向声速不变。在半径为r的圆形平面内，海面离散噪声源在经历不同的海洋环境传播后到达接收点，不同声源对环境噪声场的贡献各不相同，因此先将计算海域按照一定的角度步长分为L个方位分区，再按一定的水平距离步长将计算半径划分成J个区，则计算区域内的海面被划分为与方位角θ_l=lΔθ(l=1, 2，…, L)和水平距离r_j=r₀+j·Δr(j=1, 2，…, J)有关的多个扇环区，每个扇环的面积为：S_j=r_jΔθΔr。假设均匀分布的海面噪声源单位面积的强度级为NSL_{j, l}(dB)，则单位面积的噪声源强度I_{j, l}=ns_{j, l}²=10^{NSL_{j, l}/10}。用P(z₁, r_j, z_s, θ_l)表示位于水平距离r_j、方位θ_l、深度z_s处的声源在接收位置z₁处产生的声压，根据声学互易原理，可将接收阵元位置视为计算声源点，噪声源点则视为接收点，则接收点z₁处的风关环境噪声场可以由所有海面噪声源对接收点的贡献叠加获得^[18]：

$ \varPhi ({z_1}) = \sum\limits_{l = 1}^L {{{\rm{e}}^{{\rm{i}}{\psi _l}}}} \sum\limits_{i = 1}^J {{{\rm{e}}^{{\rm{i}}{\psi _j}}}} n{s_{j,l}}\sqrt {{S_j}} P({z_1},{r_j},{z_s},{\theta _l}) $

(1)

式中ψ_j和ψ_l是0~2π均匀分布的随机数，代表距离和方位2个方向的随机相位信息。

噪声场空间分布特征通常用两点声场复共轭积的系综平均表示，称为噪声互谱密度，它代表了噪声场的空间特性，定义为：

$ C({z_1},{z_2}) = \langle \varPhi ({z_1}){\varPhi ^*}({z_2})\rangle $

(2)

式中：*表示复共轭；〈〉号表示系综平均。将式(2)展开，本文考虑水平均匀分布的噪声源(l=m和j=n部分)互不相关，即假设不同网格单元(l≠m和j≠n)的噪声是不相干的，则噪声互谱密度可近似为：

$ C({z_1},{z_2}) \cong \sum\limits_{l = 1}^L {\sum\limits_{j = 1}^J n } s_{j,l}^2{S_j}P({z_1},{r_j},{z_s},{\theta _l}){P^*}({z_2},{r_j},{z_s},{\theta _l}) $

(3)

接收点z₁的声压值可通过射线传播程序计算获得：

$ P({r_j},{z_1}) = \sum\limits_\xi {{A_\xi }} ({r_j},{z_1},{\alpha _r})\exp [{\rm{i}}{k_0}{\varphi _\xi }({r_j},{z_1},{\alpha _r})] $

(4)

式中：φ为相位；ξ为声线索引，表示掠射角为α_r的不同声线；A为声压振幅。

令z₁=z₂，对方程(3)两边取对数获得接收深度z₁处的环境噪声强度级：

$ SL({z_1}) = 10\lg \langle |\varPhi ({z_1}){|^2}\rangle $

(5)

考虑接收阵含有N_p个垂直阵元，各阵元的声场可用列向量的形式表示：

$ \mathit{\boldsymbol{\varPhi }} = {\rm{col}} [\varPhi ({z_n})],n = 1,2, \cdots ,{N_p}, $

(6)

式中：col[ ]表示列矢量，Φ(z_n)为每个接收阵元的复声压。入射平面波阵响应向量可表示为：

$ \mathit{\boldsymbol{q}}(\alpha ,\beta ) = {\rm{col}} [{\mathit{\boldsymbol{q}}_n}(\alpha ,\beta )] $

(7)

式中：q_n(α, β)=exp(-ikx_n)；k为平面波矢量；c₀是接收阵位置的声速；β和α分别为方位角和俯仰角。

噪声场垂直指向性用矩阵向量来表示：

$ |{\mathit{\boldsymbol{q}}^*}(\alpha ,\beta ) \cdot \mathit{\boldsymbol{\varPhi }}{|^2} = {\mathit{\boldsymbol{q}}^{\rm{H}}}(\mathit{\boldsymbol{\varPhi }}{\mathit{\boldsymbol{\varPhi }}^{\rm{H}}})\mathit{\boldsymbol{q}} $

(8)

噪声场垂直指向性的均值为：

$ \langle |{\mathit{\boldsymbol{q}}^*}(\alpha ,\beta ) \cdot \mathit{\boldsymbol{\varPhi }}{|^2}\rangle = {\mathit{\boldsymbol{q}}^{\rm{H}}}{\mathit{\boldsymbol{C}}_{{\rm{noise}}}}\mathit{\boldsymbol{q}} $

(9)

式中C_noise=〈ΦΦ^H〉为平均的噪声互谱密度矩阵。

2 复杂地形海域风关海洋环境噪声仿真

本文选取近海山/岛礁海域，对复杂地形海域风关海洋环境噪声进行数值仿真，主要对接收点全方位接收的海洋环境噪声级随深度的变化、全方位海洋环境噪声垂直指向性等特性进行讨论。在仿真海洋环境噪声之前，首先分析几种海底地形方位的传播损失。

2.1 仿真计算海域环境

依据海底地形数据，选取某海山/岛礁海域用于数值仿真计算海洋环境噪声特性，其海底地形见图 2，垂直接收阵所在位置用★表示，该处海深为1 368 m；考虑50 km范围(在图中使用圆形标注)内的海面噪声源进行模型仿真计算，接收点正北方向的方位角为0°，按1°的步长将360°方位角沿顺时针方向划分为不同的扇区；计算频率为1 kHz，声速剖面采用深海声道形式(见图 3)，在本文计算范围所对应的海区，声速剖面基本呈现负梯度形式。

2.2 不同方位的声传播损失

各水平方位分区的海底地形不同，海底引起的反射也不同，导致各方位的声传播损失不同，因而在相同噪声源情形下，地形不同的方位区，对接收点海洋环境噪声的贡献不同。本文首先计算3种显著不同海底地形所在水平方位的声传播损失。

图 4是不同方位区间的海底地形变化图，方位范围分别为144°~160°、200°~220°、300°~320°，对应方位区间内最高海山的高度分别约为1 300、600和150 m。

图 5分别是N×2D算法计算的不同方位角区域内声源深度为200、800 m时，3条不同掠射角度的声线轨迹。图中的黑粗实线是海底深度线，实线、虚线、点划线分别代表掠射角0°、10°、20°的声线轨迹。由于计算海深范围内声速基本呈负梯度分布，声线均会与海底接触。在海底是海山的情况下，声线与海山斜坡每接触一次，其掠射角就会发生改变，改变量与斜坡的坡度成正比。在3种海底地形中，图 5(a)、(b)的海山高、坡度大，声线掠射角改变最大，声线轨迹变化也最大，同时相对于图 5(a)中深度200 m、图 5(b)中深度800 m的0°掠射角声线能够越过海山顶部继续向前传播，表明负梯度声速剖面和高海山情况下，更远距离的海面噪声能够以0°掠射角声线到达大深度的接收点；图 5(c)、(d)中的海山坡度相对较缓，声线掠射角改变量相对较少，其轨迹变化没有图 5(a)、(b)中的剧烈，因此在深度200 m的接收点，更多的海面噪声能够以上述3种掠射角声线的形式到达，而在800 m深度，海面噪声只有掠射角20°的声线能够到达，表明在负梯度声速剖面和小坡度海山情况下，200 m深接收点的海面噪声能量要高于深度800 m的接收点；图 5(e)、(f)的海底相对平坦，声线轨迹只发生微小的改变，同样在负梯度声速剖面情况下，200 m深接收点的海面噪声能量略高于深度800 m的接收点。上述声传播损失计算结果表明用N×2D算法在负梯度声速剖面条件下，对于坡度大的海底山，远处海面噪声源更容易到达深度深的接收点；而对于坡度小的海山和近似平坦海底，远处海面噪声源则更易到达接收深度浅的位置。

对于存在海山/岛礁海底地形的海域，位于不同位置的接收点，受到海底山/岛礁声反射的作用不同，导致不同位置的接收点环境噪声特性明显不同，下面针对图 2所示复杂海底地形，对接收点距离海底山/岛礁较远和较近2种情形，仿真计算讨论风关海洋环境噪声级和垂直指向性。

2.3 接收点距离海山较远的仿真结果

取图 2中的海底地形实际数据，接收点位于海山/岛礁群中间，距离海山较远，与其周边最近的海山/岛礁主峰的距离达20~40 km。图 6给出接收点海洋环境噪声级均值(1 kHz)随接收深度变化的数值仿真结果。2种算法的接收点噪声级在近海面处均随深度小幅增加，此后降低，不同的是随着接收深度的进一步增加，N×2D的噪声级结果基本趋于稳定，而3D的结果则一直持续小幅度降低，600 m海深内3D计算值稍高于N×2D，600 m以下海深3D计算值稍低于N×2D，二者在相同深度的差值均在4 dB以内。

从海底地形图中可以观察到接收点的东北、东南方向均存在海底山，图 7是N×2D和3D算法计算得到的各方位扇区(1°)内海洋环境噪声级随深度和方位角的变化情况，可以看出，二者的噪声级随方位变化均有起伏，N×2D的噪声级较大值位于海山高度较高的方位区域(如150°方位附近)，并且在此方位区内，噪声级较大的接收点位置较深；而海山较低的方位区域(如200°方位附近)，噪声级较大的接收点位置较浅。3D算法计算得到的水平方位分区噪声级，量级区间要大于N×2D，噪声级较大的接收点和深度，主要在海面以下几百米的深度范围内，且随水平方位角呈基本均匀分布(与N×2D结果明显不同)。从图 7可以看出，3D的扇区噪声级随方位变化的趋势(即高、低噪声强度对应的方位区间)与N×2D基本相似，只是二者的量值不同。应该说明的是，本例接收点位于海底山/岛礁群围中，四周边海底山/岛礁的反射效应，致使该点处形成近海面数百米深度处噪声级较高，3D计算模型考虑了各扇区声场的耦合，即计及水平分区之间的能量交换，基本反映了上述各种声场效应作用的综合结果。

图 8是N×2D和3D算法计算得到的各方位海洋环境噪声垂直指向性变化情况(1 kHz)，垂直阵元间隔0.7 m，首阵元深度100 m，末阵元深度121.7 m。从图中可以看出2种算法的指向性结果变化趋势基本相似：在海底较平坦的方位角区域(如310°附近)，噪声垂直指向性的结构较为稳定，存在常见的海洋环境噪声凹槽；而在存在海山的方位区，由于海山高度的不同，噪声垂直指向性形状的变化较大，环境噪声凹槽的深度得到不同程度的填补，其宽度也随之发生变化；但N×2D结果在各方位的最大值一般位于水平凹槽的左右两端或者俯仰角接近零度方向，3D结果在全方位都会存在靠近海面方向的较大值。

2.4 接收点距离海山较近的仿真结果

改变垂直接收阵的位置，如图 9所示，使其与东面的岛礁群的距离更近(约数km)，该接收点海深约240 m，其东南方位计算范围内也存在海山，计算时采用的声速剖面仍如图 3所示。图 9画出了距接收点50 km范围(图中圆圈)内的海山/岛礁情况。

图 10是上述接收点四周不同方位区间的海底地形变化图，其中图 10(a)、(b)、(c)对应的方位区间分别为：44°~60°，62°~100°，250°~280°；10(a)、(b)方位内分布着不同高度的海山，10(c)方位内深度随接收距离增大而逐渐变深。

图 11是用N×2D和3D算法计算该接收点位置的海洋环境噪声级(1 kHz)随深度变化，N×2D与3D算法的结果变化趋势相似，基本都是随深度缓慢增加，二者在同样深度的差值最大接近3 dB。

图 12是N×2D和3D算法得到的各方位扇区风关噪声级随接收深度和方位角的变化，两者的仿真结果变化趋势较为相似，存在海山/岛礁方位区域内的噪声级均较高，并且在此类方位内，3D计算噪声级略高于N×2D结果。

图 13是N×2D和3D算法计算的各方位扇区海洋环境噪声垂直指向性(1 kHz)随方位角变化情况，相邻垂直阵元间隔0.7 m，首阵元深度100 m，末阵元深度121.7 m。接收点的东面方位近距离存在海山/岛礁，方位环境噪声垂直指向性的能量角度结构发生明显的改变，近水平方向噪声凹槽的深度被填补；其他方位的噪声凹槽则比较明显。

归纳起来，上述仿真分析表明：岛礁海域风关海洋环境噪声场和接收点深度、接收点与岛礁的相对位置(包括接收点是否处于岛礁群中间)、岛礁海域海山高度、大小等密切相关。在互不相关风关噪声源均匀分布的假设和负梯度声速剖面条件下，用N×2D和3D模型仿真计算了1 kHz的环境噪声特性，虽然2种算法在不同方位的扇区噪声级存在差异，但趋势相同，且不同深度处，二者全方位的总噪声级差值不超过4 dB。3D算法由于考虑了方位间的耦合，理论上其结果比N×2D结果准确，但在计算时间上，3D算法耗时接近N×2D算法的6倍。

2.5 不同水平方位分区间隔对噪声级的影响

图 14是不同的水平方位分区步长(扇区大小)对N×2D和3D这2种算法的噪声级影响，14(a)、(b)分别是上述2个接收点距离海山/岛礁较远和较近算例在水平方位分区步长(扇区)分别是10°、5°和1°的结果。从图中可以看出，在本文采用的海洋环境下，无论是接收点距离海山/岛礁远近，不同的方位分区步长对N×2D算法的计算结果影响非常小，3D算法的噪声级计算结果会随着步长的减小而增大，也就是说对海底地形变化海域，理论上海底地形分辨率越高，水平分区越精细，3D算法的计算结果应该越准确，而水平方位分区越精细，则计算时间就会相应的增加，所以建模需要权衡计算精度和计算时间问题。

3 结论

1) 接收点位于岛礁附近，考虑方位扇区的风关噪声级和垂直指向性，存在海山/岛礁方位扇区的计算值明显不同于平坦海底方位扇区结果，如在海山/岛礁方位扇区的垂直指向性靠近水平方向(小俯仰角)的凹槽深度被填补，而其他方位扇区的噪声凹槽则比较明显。

2) 接收点与周边岛礁/海山相距较远时，考虑方位扇区的海洋环境噪声级，海面以下几百米海深内3D算法的结果明显高于N×2D算法，导致几百米海深内全方位总噪声级3D计算值也高于N×2D，在更深的接收深度3D计算值稍低于N×2D；接收点距离海山/岛礁较近时，负梯度声速剖面和海山/岛礁地形致使存在海山/岛礁的方位扇区噪声级均较高，3D计算的各扇区噪声级略高于N×2D；本文算例中2种算法的差值均在4 dB以内。

总体来说，基于风关海洋环境噪声模型的仿真计算结果表明，岛礁海域的风关噪声场和接收点的深度、接收点与岛礁的相对位置(包括接收点是否处于岛礁群中间)、岛礁海域海山的高度和大小等因素密切相关。在随机均匀分布、互不相关风关噪声源假设和负梯度声速剖面条件下，仿真计算1 kHz的风关噪声特性，虽然N×2D近似和3D这2种算法在不同方位的扇区噪声级计算结果存在差异但趋势相近，且不同深度处，二者全方位的总噪声级差值不超过4 dB，理论上本例3D算法结果比N×2D更准确，但相同水平方位分区条件下3D算法耗时约为N×2D算法的6倍，因此在水平方向声速没有剧烈变化情况下，当对计算时间有限制时，可选用N×2D算法来计算风关噪声场特性。