利用海洋环境噪声空间特性估计浅海海底分层结构及地声参数*

李赫 郭新毅 马力

1) (中国科学院声学研究所,中国科学院水声环境特性重点实验室,北京 100190)

2) (中国科学院大学,北京 100049)

海洋环境噪声场中包含了海洋中的诸多信息,海底地声参数是影响海洋环境噪声场空间分布的主要因素之一.对于不同的海底分层结构,海底反射损失会根据沉积层厚度和各层声速呈现出不同的临界角和干涉条纹结构.本文利用Harrison能流理论,从理想反射系数出发,分别考虑了声速、密度、衰减系数、沉积层厚度等几种参数对无沉积层和单层沉积层中反射系数的影响,并对单层沉积层海底的反射系数进行了化简,结合互易原理解释了反射损失条纹结构的形成机理.中国黄海某海区试验结果表明,利用海洋环境噪声空间方向谱获得的海底反射损失,可以提取海底反射临界角和干涉条纹信息,由此可估计出海底分层结构、声速和沉积层厚度等海底参数信息.

1 引 言

海洋环境噪声是海洋中永恒存在的声场,主要由风浪、降雨、舰船、海洋生物、工业等因素形成,其中风成噪声在各个海域普遍存在,其频段覆盖几百至几十千赫兹[1].一方面海洋环境噪声是影响声呐工作性能的主要因素之一,只有对其进行充分的了解研究才能充分提高水下设备的工作性能,所以国内外学者纷纷对海洋环境噪声进行建模研究[2−5].另一方面,海洋环境噪声场中也包含了诸多水体、海底、海面等环境信息[6],因此利用海洋环境噪声进行参数反演,获取海水、海底等环境信息也成为时下的热点问题之一[7,8].

海底地声参数反演主要分为主动反演和被动反演两部分.Zeng等[9]利用爆炸声源对中国黄海海域海底衰减系数进行反演;李梦竹等[10]提出了一种适用于低声速沉积层的海底参数声学反演方法.相比于主动反演,利用海洋环境噪声进行被动反演无需主动发射声源信号,在实验中仅需数据接收装置,大幅度节省了实验消耗以及工作任务,并且对海洋生物没有任何影响.周建波等[11]利用矢量水听器通过海洋环境噪声提取了声场格林函数;江鹏飞等[12]根据海洋环境噪声的空间指向性在不同掠射角范围内对海底参数的敏感度不同,提出了分步反演的方法;骆文于[13]利用中国东海环境噪声数据反演了液态半空间海底声速、密度、衰减系数;早在1982年,布列霍夫斯基[14]就曾利用海洋环境噪声场的各向异性,来获取海底反射信息;Harrison在2002年进一步系统地研究了海洋环境噪声中的能流理论,利用上行波和下行波的比值提取了海底反射损失,并对海底进行了反演[15],并利用希尔伯特变换补全反射系数的虚部,通过计算垂直方向的冲击响应利用漂流浮标对海底进行成像[16];Siderius等[17]利用海洋环境噪声提取了格林函数用来表征两点间的冲击响应,对海底分层进行成像;Muzi等[18,19]在Harrison的基础上利用互谱密度矩阵的特性对阵元进行了理论扩展处理,提高了利用海洋环境噪声提取的海底损失曲线的分辨率.

本文在Harrison能流理论基础上,从理想环境反射系数出发,解释了不同海底参数对海底反射损失(bottom−loss,BL)曲线的影响,并利用实测海洋环境噪声数据对海底分层结构、声速及沉积层厚度进行估计.本文首先引入Harrison能流理论估计海底反射损失的方法,分析了不同海底分层情况下的曲线特点;然后利用简正波理论,讨论了浅海波导中近场连续谱和远场离散谱在噪声估计海底反射损失中的贡献;并从理想反射系数出发,分别讨论了不同海底参数对无沉积层、单一海底沉积层两种情况的影响,结合互易原理解释了反射损失条纹结构的产生机理,并根据噪声获得的海底反射损失曲线给出海底声速、沉积层厚度的估计方法;最后介绍了 2016年冬中国黄海某海区海上试验,对实测海洋环境噪声数据进行处理,估计其海底分层结构、海底声速和沉积层厚度.

2 利用海洋环境噪声垂直方向谱估计海底反射损失

从浅海理想波导出发,假设噪声源均匀分布在海面,利用垂直阵接收环境噪声信号,其空间结构如图1所示.波导中的声波可分为两部分： 接收阵上方的来波称为下行波,从噪声源直达接收阵;下方的称为上行波,经过海底一次反射后到达接收阵.在此近场情况下不考虑多次反射,后文将给出解释.

图1 浅海中利用垂直阵接收噪声空间结构Fig.1.Spatial structure of receiving ocean ambient noise using vertical line array in shallow water.

声速剖面为等声速的理想情况下,表面声源的出射角θs、海底掠射角θb与接收角θr均相等,在此均用掠射角θb表示.对于某一角频率ω下,海底反射损失BL有如下定义：

其中B为波束形成得到的不同角度和角频率下的接收能量,即垂直方向能量谱,

其中E代表期望;H为共轭转置;w是波束形成的加权向量;p为对应角频率ω的接收数据,在M阵元接收时,为M维向量;Cω为噪声互谱密度矩阵.

为说明此方法计算得到BL的有效性,先考虑极端情况,假设阵元间隔为0.2 m的垂直阵遍布在声速为1500 m/s、水深为40 m的全海深中,忽略海底横波的影响,海底参数如表1所列.

表1 海水及海底声学参数Table 1.Parameters of ocean and sub−bottom.

互谱密度矩阵直接利用OASES软件中的OASN模块计算获得[20],w设置为常规波束形成加权向量.图2(a)为根据理想情况下真实反射系数公式求得的BL[21],后文将对该公式进行讨论;图2(b)是根据噪声垂直空间指向性利用OASN模块仿真得到的;图2(c)提取出了1800 Hz下两种方法的掠射角BL曲线.图2(a)和图2(b)中的条纹与图2(c)中的曲线峰值成对应关系,阵元遍布全海深的理想状态下,图2(c)中两条曲线的峰值一一对应,仅在幅度上存在一定偏差,因此利用Harrison能流理论估计海底反射损失是可行的.

图2 (a) 真实反射损失BL;(b) 根据噪声垂直空间指向性利用OASN模块仿真得到的 ;(c) 1800 Hz下两种方法的比较,实线为BL,虚线为Fig.2.(a) True BL;(b) computed by vertical direc−tionality of ocean ambient noise using OASN;(c) two meth−ods compare under 1800 Hz,the full line is BL,the imagin−ary line is.

而在实际海上试验中,阵元个数是有限的,很难达到上述遍布全海深的理想状态.有限的阵元个数会带来不同掠射角方向分辨率的下降,这样会导致图2(b)中条纹信息一定程度上的缺失,同时也使图2(c)中峰值的模糊,但仍然有很多特征可以利用.图3给出了两种不同海底分层结构的,保留表1中的参数,图中实线为间隔0.2 m的水听器遍布全海深的结果,虚线为间隔0.2 m的42元水听器阵,第一个阵元设置在水深30 m处,频率均为1800 Hz.海底模型如图4所示,其中图4(a)为无限大液态声学半空间海底,图4(b)为带有一层沉积层的海底.

图3 1800 Hz下不同海底分层结构下的 (实线为阵元遍布全海深,点划线为42阵元,间隔均为0.2 m) (a) 海底为无限大液体声学半空间;(b) 海底为单层沉积层Fig.3.The of different structure of sub−bottom under 1800 Hz： (a) Infinite acoustic half space;(b) 1 layer of sedi−ment.The full line corresponds to the condition that the elements set across the sea.The imaginary line corresponds to the condition that 42 elements set at the depth of 30 m.

图4 海底分层模型 (a)液体无限大声学半空间海底;(b) 存在一层沉积层海底Fig.4.Model of sub−bottom stratification： (a) Infinite acoustic half space;(b) 1 layer of sediment.

3 不同方向噪声对 曲线影响分析

在浅海波导中,点声源在大掠射角范围内发出的声波很难进行远距离传播,因为这一部分经过海底反射存在能量损失,反射波携带着一定海底信息,这一特性对于均匀分布在表面的噪声源同样适用.在分层介质中,点源声场的波函数的积分表达式为

为了描述准确,这一部分在K/I模型下进行讨论[2].Kuperman给出海面噪声源均匀情况下噪声互谱密度的波数积分表达式

图5 Pekeris分支割线、极点和积分的复波数平面Fig.5.Pekeris branch cut map.

和简正波表达式

式中g为点源格林函数,Ψm为模式函数.由于接收器为垂直阵,式中r=0.对于近场连续谱部分采用(5)式波数积分表达式进行计算,而远场离散谱部分则采用(6)式简正波表达式,将其分别代入(3)式中求得不同谱域内的.图6给出了无限大液态声学半空间海底情况下,不同谱域内噪声场的垂直指向性和对应的,图6(a)和图6 (b)是利用(5)式计算近场连续谱部分,即 0＜kr＜k2,这一部分的贡献主要来自大掠射角,其对应的是大于临界角的部分,而这部分经过海底后存在能量泄漏,返回到接收器阵的声波可以表达出海底反射的信息;图6(c)和图6(d)是(6)式计算得到的远场离散谱部分,此时k2＜kr＜k1,其贡献主要来自于小掠射角,对应小于临界角的部分,这一角度内发生全内反射,海底没有声能吸收,对的大掠射角条纹部分没有贡献;图6(e)和图6 (f)为全谱域结果,与仅考虑近场条件下的基本一致.在近场条件下大掠射角部分反射损失较高,所以本文不考虑多次海底反射的情况.

图6 不同谱域内噪声场的指向性和 (a),(b)连续谱部分噪声场;(c),(d) 离散谱部分噪声场;(e),(f)全谱域噪声场Fig.6.Vertical directivity and of ocean ambient noise in different spectral domain： (a),(b) Near field continuous spectrum;(c),(d) far field discrete spectrum;(e),(f) full spectrum domain.

4 BL曲线条纹结构产生机理及海底参数估计

4.1 不同海底参数对无沉积层BL的影响

影响BL的海底地声参数有3个,分别是海底衰减系数、声速c和密度.图7别为不同海底参数下的BL.由图7(a)可见,仅在掠射角小于临界角时,对BL有影响,在建立浅海远程传播模型中起到了至关重要的作用[21].但第3节中,通过对近场连续谱和远场离散谱的分解已经证明： 噪声估计的主要来自于近场连续谱的贡献,在近场条件下,衰减系数又是一个不敏感参数,因此可以忽略衰减系数带来的影响.在分界面上掠射角θ1,θ2利用Snell定律建立关系

图7 不同海底参数下的BL (a)衰减系数 ;(b)声速c;(c)密度Fig.7.The BL under different parameters of sub−bottom：(a) Attenuation coefficient ;(b) sound speed c;(c) dens−ity.

透射波水平射出,即θ2=0°时存在临界角效应,此时,cosθ1=c1/c2,在海水声速c1已知的情况下,反射系数的临界角θ1完全由c2决定.除此之外,海底声速对于反射系数的幅度也有一定影响.图7(b)为不同海底声速下的BL.海底密度ρ仅对BL的幅值起到影响,与临界角不存在任何关系,如图7(c).

4.2 单层沉积层BL条纹结构产生机理

下面考虑存在一层沉积层时的反射情况,反射系数的表达式为

其中,R12,R23分别为两个分界面上的反射系数,φ=k2hsinθ2表示声波穿过沉积层的垂直相移.(8)式的分子和分母中均含有 exp(2jφ) 一项,难以讨论其各个参数对R的影响,特别是图2反射系数中影响条纹的因素,为了让各个参数在两式中的作用明显表达,需要对其进行一定程度上的化简,下面从单一沉积层的反射模型开始考虑,如图8.

图8 单层沉积层海底反射模型Fig.8.Reflection model of sub−bottom with 1 layer of sedi−ment.

声波在沉积层中传播,会形成多次返回海水中的反射波,反射系数的原始表达为各个反射波的叠加

式中,T12表示海水进入沉积层的透射系数,T21代表沉积层进入海水中的透射系数.近场条件中,海面风成噪声源源级较低,可忽略海底多次反射后返回海水中反射波带来的影响,仅考虑图8中前两条返回声线,这种情况下反射系数化简为

与(8)式相比,(10)式非常简洁,便于讨论不同参数与R的对应关系.图9分别为(8)和(10)两式计算的海底反射损失BL和1800 Hz下的BL曲线,依然沿用表1中参数.从图9可以看出,两式计算的结果基本一致,仅在幅值上存在细微差异.从级数的角度解释,大掠射角近场条件下,R为一纯小数,(9)式中第三项之后均为第二项的高阶小量,因此仅在幅值上存在一定影响.综上,可以利用(10)式化简后的反射系数R代替(8)式,对单层沉积层海底反射进行讨论.

用θc1和θc2表示两个临界角,当掠射角θ＜θc1时,没有透射波,(10)式仅存在R12一项,此时R为模近似1的复数;当θc1＜θ＜θc2时,(10)式中R12,T12,T21均为实数,R23为模近似1的复数,R12为随θ单调递增的负数,其模单调递减,T12=1+R12为正且单调递增,T21同理,因此对于此部分振荡完全由 exp(2jφ) 一项引起;当θ＞θc2时,(10)式中除 exp(2jφ) 一项以外均为随θ单调改变的实数,其本身不存在极值点,因此这一部分与θc1＜θ＜θc2时一样,BL曲线中的振荡均由exp(2jφ)一项引起.在计算BL过程中,取R的模,所以exp(2jφ)项中仅有实部对BL曲线有贡献,其实部通过Euler公式可知

因此BL曲线中的振荡由声速和沉积层厚度h共同决定.在海底声速确定的情况下,对于某一掠射角θ0,BL曲线对于频率f的振荡周期可由(11)式确定,把频率f看作变量,(11)式表示为

其振荡周期为c2/2hsinθ2,在某一频率fc范围内振荡次数为 2hsinθ2/c2·fc,振荡次数即BL的极值点数,这与图9(a)和图9(b)中BL的条纹个数是对应的.

根据互易原理,可以看作不同距离的表面噪声源到达同一接收点的路径,图8中前两条反射声线看作入射声线,原来的入射声线看作汇聚后的出射声线.两条入射声线的相位差等于声线在沉积层中走过路线的相位变化,即垂直相位变化 2φ.在相位相差 2π 的整数倍时发生相长干涉,

图9 (a) (8)式计算的真实BL;(b) (10)式计算化简后的BL;(c) 1800 Hz下两条BL曲线Fig.9.(a) True BL computed by equation (8);(b) BL after simplified computed by equation (10);(c) the curves of BL of two methods under 1800 Hz.

(13)式f0即为(12)式中f的峰值点,与其振荡周期相对应,同样也解释了图9中BL条纹的形成机理.

而在某一频率确定的情况下,(12)式中sinθ2可用(7)式确定,把 sinθ2作为变量,(12)式的振荡周期为c2/(2fh) ,从图10可以看出BL在大掠射角范围内关于 sinθ2均匀振荡.

图10 (10)式计算化简后的BL (横轴利用(7)式转为 sinθ2)Fig.10.Curve of BL after simplified computed by equation(10) with the horizontal axis converted to sinθ2 by equa−tion (7).

4.3 利用海洋环境噪声获得的估计海底地声参数

BL曲线的临界角效应和周期振荡这两种特征可以分别用来确定单沉积层海底的分层声速和沉积层深度.在利用噪声获得时,由于阵元的个数限制,会引起条纹的模糊,图11(a)仿真了42阵元,第一阵元设置在30 m处的,其他参数依旧沿用表1,与图2阵元遍布全海深的相比,条纹显得模糊很多.图11(b)单独提取出了1800 Hz下的曲线,该曲线中依然可以看到较为明显的临界角信息,即幅值发生跃变的点,通过两处临界角θc1,θc2可以利用上文所述的Snell定律,在折射角为0时确定沉积层声速

和基底声速

而在θc2后,存在一个明显“凸点”,这个“凸点”不一定是这条曲线本身的极值点,但实际上对应的是理想情况下反射损失BL的一个极值点,由于有限阵元带来的模糊,使得这个极值点变为了“凸点”.假设“凸点”的位置对应为θt,那么利用θt与θc2便可以得到(12)式在某一固定频率下关于sinθ2的半周期

图11 (a) 42阵元设置在水下30 m处 ;(b) 1800 Hz下 曲线Fig.11.(a) with 42 elements at the depth of 30 m;(b) the curve of under 1800 Hz.

从图11(b)可以看出,θt约为46°,θc2约为41°,利用(17)式估计出沉积层厚度约为2.84 m,与实际仿真设置的3 m相比偏小,这是因为有限阵元带来的模糊,即空间采样点不足导致曲线中的极值点之间的分辨率下降.也可以利用理论阵元合成的方法对噪声的互谱密度矩阵进行扩展,以提高曲线的分辨率[17,18].

5 利用噪声能流理论获取海底参数海上试验验证

实验数据于2016年12月在中国黄海某海区采集,当地实际水深35 m,冬季海水近似看作等声速,声速为1487 m/s,实测声速剖面如图12所示.实验采用阵元间隔为0.2 m的25元船挂阵采集噪声数据,将实测数据分为20组,每组时长5 min,对每组数据分别进行处理.图13和图14分别给出了第9组和第11组数据处理结果,1600 Hz下两图的曲线中明显存在两次突变,因此存在两个临界角,符合单层沉积层海底模型,可以利用上文所述的方法对海底分层结构及地声参数进行估计.图中已标注出θc1,θc2,θt的位置,θc1,θc2的选取依据主要是根据曲线中幅值突变的程度,选取其突变的初始角度,θt选取的是θc2后的第一个“凸点”,然后利用(14),(15),(17)三式对该海区海底声速和沉积层厚度进行估计.图15为20组数据计算结果的统计分布直方图,其中图15(a)为沉积层声速,图15(b)为基底声速,图15(c)为沉积层厚度.沉积层声速的估计平均值为1544 m/s,这与同海区的反演结果1531 m/s相近[22];基底声速的估计平均值为1935 m/s,相比沉积层声速其波动范围略大,主要集中在1900－1960 m/s之间,目前暂无其他方法对这一结果进行对比;沉积层厚度的估计平均值为2 m,考虑到有限阵元带来的模糊影响,实际沉积层厚度可能略高于估计值.

图12 实验实测海水声速剖面Fig.12.Sound speed profile of the experiment.

图13 第9组数据处理结果 (a) ;(b) 1600 Hz下的 曲线Fig.13.Results of the 9th set of data： (a) ;(b) curve of under 1600 Hz.

图14 第11组数据处理结果 (a) ;(b) 1600 Hz下 曲线Fig.14.Results of the 11th set of data： (a) ;(b) curve of under 1600 Hz.

图15 数据计算结果统计分布直方图 (a)沉积层声速;(b)基底声速;(c)沉积层厚度Fig.15.Statistical distribution histogram of 20 sets of data：(a) Sound speed of the sediment;(b) sound speed of the basement;(c) thickness of the sediment.

为验证所估计沉积层厚度的可靠性,试验船搭载了浅地层剖面仪(sub−bottom profiler),该地层剖面仪发射调频信号,利用声波在水中和水下沉积层内传播和反射的特性来探测海底地层结构.图16给出了从试验点向东北方向扫描的海底地形,35 m深处为海底表层,在其下方存在沉积层,图16(b)中用黑色标出,厚度约为2.2 m,噪声估计出的沉积层厚度与此结果基本符合.

图16 浅地层剖面仪探测试验海区海底结构(图(b)中用黑点标明海底分界面)Fig.16.Sub−bottom structure of the experimental area de−tected by sub−bottom profiler.The interface is marked by black spot in panel (b).

6 结 论

本文根据Harrison能流理论,假设噪声源均匀分布在海面,利用垂直阵接收的海洋环境噪声数据对海底反射损失进行估计,讨论了无限大液体声学半空间海底和存在一层沉积层两种情况的曲线特点,并将近场连续谱和远场离散谱两部分的贡献分开考虑,从理想反射系数出发,推导了各个海底参数对BL的影响.