基于多层水平变化浅海海底模型的低频反演方法*

康娟 彭朝晖 何利 李晟昊 于小涛

1) (中国科学院声学研究所,声场声信息国家重点实验室,北京 100190)

2) (中国科学院大学物理学院,北京 101407)

准确获取海底声学参数对声场分析、声纳应用等具有重要意义.反演是获取海底声学参数的重要手段之一.为了避免反演中的多值问题,地声反演往往尽量减少待反演参数,多采用单层或双层水平不变地声模型,反演后的声学参数仅能在有限频段应用,难以同时应用于较宽频段.本文以实验中浅地层剖面仪测量的海底沉积层分层结构为基础建立地声模型,采用多物理量分步联合手段反演各沉积层声学参数,并通过匹配场定位、后验概率分析、采样等手段对反演结果进行验证.将反演结果与单层等效地声模型反演结果进行对比分析,结果说明,本文建立的多层水平变化海底声学参数更接近真实情况,可以应用在较宽频段范围,解释不同声学现象,且精度更高.

1 引言

在浅海,海底作为海洋信道的下边界,其地声参数会对声波传播产生显著影响.反演是一种有效获取浅海地声参数的方法,可以快速和低成本地估计出一定范围内的海底沉积物声速、密度和衰减系数等声学参数,受到水声学界的广泛关注.国内外已经发展了多种反演方法,包括匹配场反演法[1,2]、海底反射系数反演法[3]、简正波到达结构反演法[4]、传播损失反演法[5,6]、脉冲时间到达结构反演法[7]、声场频率-距离干涉结构反演法[8]等等.采用单个物理量构成代价函数进行多个物理量反演时,往往要面对反演结果的多值性问题,因此有学者利用不同代价函数对海底参数敏感度不同,提出了多个代价函数分步反演地声参数的方法.2004 年,Li 和Zhang[9]提出的多种反演方法联合反演方案,经多个海区进行实验验证,一定程度上解决了反演中的多值性问题.2019 年,李梦竹等[10]先利用海底反射损失得到层厚以减少待反演参数维度,再利用匹配场反演其他参数,有效地降低了反演结果的不确定性.

在对浅海海底参数反演时,一般将海底等效为如图1 所示三种水平不变海底模型[11].对于较窄频率范围的声学应用,采用这三种模型通常可以达到比较好的效果.但是实际应用时,往往需要更接近实际的多层地声模型.崔智强等[12]利用传播损失作为代价函数,仿真表明在假设三层沉积层的海底环境下,能够提高地声参数的反演精度.鹿力成等[13]通过研究浅海海底模型对于声传播的影响,发现较高频率声场主要受海底表层影响,而影响低频声场的海底深度能够达到几十米.此时,要想同时满足较宽频段的声传播,需要采用多层海底模型.

图1 三种较为常见海底模型Fig.1.Three typical bottom models.

实际地声反演中,除了要考虑海底分层结构外,往往还需要考虑海底水平变化情况.祝捍皓等[14]针对浅海水平变化波导下的低频声能量传输问题,基于有限元方法,讨论了三种不同海底变化下带来的声能量差异.邱海宾和杨坤德[15]在利用拖线阵对海底参数进行反演时,建立了水平变化海底模型,在声源频率增多、拖线阵变长的情况下,依旧得到了收敛性好的匹配场结果.Kim等[16]考虑到海底地形在空间上的差异,将反演区域分为几段进行反演,反演结果很好地解释了实验场地随水平变化的声学现象.Bonnel等[17]在利用啁啾声反射手段得到多层结构的基础上,提出一种在模态号数未知情况下利用频散曲线反演随距离变化多层海底参数的方法,成功得到沉积层的地声参数.

2022 年8 月,中国科学院声学研究所声场声信息实验室在我国南海北部某大陆架海区进行了一次海底声学反演实验.实验中利用浅地层剖面仪对海底进行了测量,如图2 所示,可见海底具有明显的分层结构及水平变化.本文利用本次实验测量的浅地层剖面仪、海洋环境及声学数据,开展水平变化分层海底地声参数反演方法研究,建立起更接近真实海底并适用于较宽频带(50—500 Hz)声场计算的多层水平变化海底地声模型.

图2 浅地层剖面仪获取的海底分层结构Fig.2.Seabed layering obtained from sub-bottom profiler.

2 海上实验及海底地声模型

海上实验选用没有多次气泡脉动干扰的气枪作为声源,其标称深度为7.5 m.实验期间获取的气枪信号声源谱级曲线如图3 所示.使用潜标垂直接收阵接收记录声信号,垂直阵共有15 个阵元,深度分布在5—80 m 之间.本文反演过程中将选用10 km 范围内的气枪信号用于计算海底声速、密度和厚度,20 km 范围内的气枪信号用于计算衰减系数.

图3 气枪信号声源谱级Fig.3.The frequency response curve of source level of air gun.

垂直阵不同深度处共绑扎了16 个自容式温度压力计(TD)组成温度链,实验期间获取的声速剖面如图4 所示.可见是典型的浅海夏季温跃层声速剖面,存在线性内波,但声速剖面起伏不大.因此,本文不考虑水文环境的水平变化,声速剖面取信号发射期间的平均声速剖面,如图4 中黑线所示.

图4 声速剖面Fig.4.Sound-speed profile.

实验期间通过浅地层剖面仪获取了海底沉积层结构,如图2 所示.对该海底分层结构进行简化,得到地声模型如图5 所示.图5 中,接收点位处为0 km,传播路径上的海底主要由3 层沉积层和1 层基底层组成.假设各沉积层内的海底声速、密度及衰减系数不变,则待确定的参数主要包括各层的厚度、海底声速、密度及衰减系数,共15 个参数.

反演前,可根据先验信息减少待反演参数数量.首先,利用浅地层剖面仪得到的各沉积层声波传播时间,将各沉积层的厚度变换为海底声速的函数:

式中,tm+1-tm为m层两个边界的往返时间差;cbm为m层的海底声速.

其次,考虑到海底密度敏感性远小于海底声速,海底密度可以通过Hamilton 经验公式由海底声速表示[18]:

第三,利用较近距离声信号反演海底声学参数,一般代价函数对海底衰减系数不敏感.反演时,可以先将海底衰减系数设定为定值,先对图5 中各层声速cb1—cb4进行反演,然后将反演后的各层声速作为已知条件,利用不同距离和频率的声传播损失来反演海底衰减系数.

反演中,需要选择合适的声场计算模型.在对浅海低频声传播进行研究时通常采用简正波理论,由于所考虑的海底水平变化较弱,且将各层海底参数视为随距离不变,因此在计算声场时不考虑简正波模态之间的耦合,采用绝热简正波进行表示[19]:

其中krm为简正波特征值;φm为特征函数.

3 反演方法及过程

3.1 反演流程

基于上述分析,对分层海底地声模型的声学参数采取分步反演的方法,反演流程如图6 所示.首先,利用浅地层剖面仪获得海底分层结构,选取简正波频散结构对各层声速进行反演.再利用匹配场反演方法对表层声速进一步反演.最后将利用实测的传播损失反演海底衰减系数.在反演过程中采用多种方式对反演结果进行验证.

图6 反演流程图Fig.6.Flowchart of inversion.

为避免待反演参数陷入局部最优,采用并行多种群遗传算法[20]进行优化搜索.具体操作为:编码采取20 位二进制数;种群数目选择为3;每个种群有24 个个体;多个种群交叉概率在0.7—0.9 间随机产生;多个种群变异概率在0.001—0.05 间随机产生;精英种群代价函数连续20 代不再变化时停止搜索.

3.2 简正波频散结构反演各层声速

在浅海,接收的声压信号p可以用一系列简正波的叠加来表示[21]:

式中,An(t),Φn(t) 分别代表在t时刻、第n号简正波的幅度和相位.

利用简正波频散结构的反演海底声速,其基本原理是:由于简正波存在频散现象,不同阶数不同频率的简正波群速度不同,对一定距离外接收的脉冲声信号进行时频分析,提取出来频散曲线与简正波群速度密切相关,可用于海底声速反演,采用如下代价函数[22]:

其中M为待反演的参数向量; ΔTn(f) 为第n号简正波不同频率信号到达时间差; ΔTnm(f) 为频率f下第n号简正波和第m号简正波信号到达时间差.上标 e 代表实验提取数据,c 代表拷贝场数据,这里通过简正波模型KEAKENC[23]计算拷贝场.

在距离为r的情况下,不同频率同号简正波到达时间差 ΔTn(f)为

同一频率不同号简正波信号到达时间差ΔTnm(f)为

其中水平变化海底群速度近似为

利用简正波频散结构进行反演,关键在于如何精确地提取出各号简正波的到达时间差(f)及(f) .本文采用warping 变换结合Wigner-Ville 分布的方法提取简正波频散结构并计算各号简正波到达时间差.

Warping 变换是一种可逆的变换,可以通过逆变换算子对warping 变换后的信号采取时频域滤波再恢复的方式分离出各号简正波.分离得到的各号简正波可以利用Wigner-Ville分布做时频分析,提取出其频散结构,Wigner-Ville分布的表达式为[25]

式中,*代表复数共轭,它是一种对信号瞬时自相关函数的傅里叶变换.

在实际反演过程中,首先选取距接收点右侧约6 km 处的气枪信号如图7(a),提取接收深度处于75.93 m 的阵元接收信号的频散曲线.如图7(b)和图7(c)为warping 变换后的频域分布图和时频分布图,其中第一阶简正波受噪声干扰程度较大且与第二阶简正波形成干涉现象,故选取三到八阶简正波频散结构用于反演,图7(d)为提取出的三到八阶简正波频散曲线.

图7 采用warping 变换提取频散曲线 (a) 接收信号;(b) warping 变化后的频域信号(图中红色虚线为各号简正波截止频率);(c) warping 变换后信号时频图;(d) 接收信号时频图(图中白色虚线为提取出的频散曲线)Fig.7.Estimating dispersion curves using warping:(a) Received signal;(b) spectrum of the warped signal,and the cutoff frequency;(c) spectrogram of the warped signal;(d) spectrogram of the signal after warping and Wigner-Ville,and estimated dispersion curves.

采用后验概率分析方法可以评估反演方法的有效性[26].在遗传算法反演进化过程中,保存所有的搜索值Mk以及对应的目标函数 cost(Mk),对反演结果在搜索空间Ω内进行概率统计,将向量M按照其代价函数的大小进行排列,按照Boltzmann 分布进行加权,得到反演结果的第k组向量后验概率分布:

其中Nobs是保存下来的模式向量的个数;T是温度控制参数,一般取反演过程中的50 个最佳代价函数的均值较为合适;向量M中第n个参量值为m的边缘概率密度为

针对频散曲线图7(d)进行反演时,保留样本的边缘概率密度如图8 所示,其收敛程度依次为cb2,cb3,r,cb1,cb4.其中cb4的收敛 性最 差,且收敛于最高值2000 m/s 处,说明它对于代价函数敏感性较弱,后续将cb4固定为2000 m/s,不再对其进行反演.而相比于收敛性较好、收敛结果近似为单值的cb2和cb3,cb1存在一定的离散性.

图8 反演参数的边缘概率密度Fig.8.The marginal probability distribution of inversion parameters.

为了进一步检验反演方法的有效性,将反演得到的各海底参数作为输入,计算距离约为6.2 km、深度处于75.93 m 的阵元接收到的信号对应的简正波频散结构,并与实验数据提取的频散曲线进行比较,如图9 所示,红线为模型计算的频散曲线,蓝线为实验提取的频散曲线,可以看出两者符合较好.再将反演的海底参数用于匹配场定位,结果如图10 所示,图10 中用红色圆圈标出匹配场定位的距离为5.9 km、模糊深度中值为7.75 m,与实际气枪声源的距离5.98 km、标定深度7.5 m 相比,误差较小.

图9 频散曲线拟合图Fig.9.Dispersion curves as estimated from the data with warping (blue) and predicted after inversion (red).

图10 定位结果Fig.10.Localization result of Bartlett processor.

进一步对不同距离多个气枪信号进行反演,得到的结果如表1 所列.统计反演结果发现,cb1产生了40 m/s 范围的浮动,cb2产生了12 m/s 范围的浮动,cb3产生了33 m/s 范围的浮动,结合前面的收敛性分析,也表明了利用简正波频散结构反演二三层声速较第一层更为准确,因此需要利用其他方式完成对于第一层声速的进一步探究反演.

表1 不同距离下海底参数反演结果Table 1.Inversion results at different distances.

3.3 匹配场反演表层声速

以宽带相关Bartlett 匹配场处理器作为代价函数:

M为待反演参数向量;N为接收通道个数;pa(fi)为第a个垂直阵元接收到的复声压,e 代表实际数据; c 代表拷贝场数据; * 代表复数共轭.该代价函数φ位于0—1 之间,越靠近零越接近于真实场.

匹配场反演之前,考虑到不同频率声波在传播过程中作用于不同海底深度,首先通过仿真分析匹配场在不同频率时对海底参数的敏感性,以获得更适用于反演表层声速cb1的方案.

假设真实的海底结构如图5 所示,共分四层,各层声速cb1—cb4分别为1580,1630,1700,1800 m/s.与实验情况相同,声源深度为7.5 m,接收器位于5.9 km 处,深度位于5—80 m,共15 个水听器,仿真计算垂直阵接收的信号用于模拟实验获取的声信号,计算频率为50—500 Hz.

首先,采用半无限大的单层海底,在不同频段范围内(50 Hz 带宽)进行等效声速cbe的反演,声速搜索范围为1500—2000 m/s,步长为5 m/s,结果如表2 所列.可以发现,50 Hz 频率下反演得到的cbe为1725 m/s,位于真实海底结构中下层声速cb3和cb4之间,随着频率增加,cbe逐渐减小,当频率大于250 Hz 时,cbe接近于真实海底的表层声速1580 m/s,说明较高频率下声波主要受表层作用,此时匹配场反演对底层声速较不敏感,因此在利用匹配场进行多层海底地声反演时,应选取较低频率进行.而在较低频率上下层海底共同产生影响时,需要对各层声速进行敏感性分析,分析时选取100—150 Hz 频段.

表2 仿真中不同频率下反演获取的单层等效声速Table 2.Single-layer equivalent sound speeds at different frequencies.

下面,利用100—150 Hz 较低频段对各层声速进行敏感性分析,通过改变其中一个声速,保持其他声速不变的情况下,来测试φ(M) 的敏感性,结果如图11 所示,圆圈代表其真实值.由图11 可以看出,待反演的参数中第一层声速cb1敏感性远远大于cb2和cb3.因此,可以利用高低频匹配场结合的方式对表层声速进行反演.

图11 各层声速敏感性分析Fig.11.One-dimension cost function value of sound speeds in each layer.

接下来,再利用实验数据进行分析.首先将海底视为一个半无限大的单层海底,利用匹配场方法反演不同频率下的单层等效声速cbe,结果如图12所示.图12 中,红点代表各频段中代价函数最小值对应的cbe.图12 中可以看出,与上面的仿真结果类似,频率低于250 Hz 时反演得到的cbe应该受较深沉积层声速影响较大,而频率大于250 Hz 时反演的cbe应该只能反映表层声速cb1的变化情况.尽管图12 中不同距离下得到的cbe随频率的变化关系存在一定的差异,但仍然可以总结出如下规律:250—300 Hz的cbe在1600—1630 m/s之间,与3.2 节利用简正波频散结构反演得到的cb1较为接近;在350 Hz处cbe存在一个明显的声速跃变,此频率下的cbe在1560—1580 m/s 之间;大于350 Hz后随着频率增加,cbe逐渐减小.因此可以根据随着频率增加,声线可穿透的海底厚度减小,猜想表层可能存在一个明显的分层,该分层以上声速呈正梯度分布,平均声速在1560—1580 m/s 之间;而该分层以下声速较大且较为稳定.

图12 不同频率下的单层等效声速 (a) r=3.5 km ;(b) r=3.7 km ;(c) r=3.9 km ;(d)r=4.1 kmFig.12.Single-layer equivalent sound speeds at different frequencies:(a) r=3.5 km ; (b) r=3.7 km ; (c) r=3.9 km ;(d) r=4.1 km .

为此,对浅地层剖面仪数据的表层分层结果进一步分析,建立更细化的沉积层分层结构,如图13所示.设上表层存在cb11—cb12的正声速梯度,平均厚度为h11,下表层声速为cb13,平均厚度为h12.在反演过程中,将简正波频散结构反演得到的二三层平均声速(cb2=1572 m/s,cb3=1620 m/s)作为已知,利用100—150 Hz 频带下的信号通过匹配场方式反演表层声速变化情 况(即cb11,cb12,cb13),得到的结果如表3,其中cba为cb11和cb12的平均声速.和图12 进行对比,可以看出表层情况同单层等效声速变化较为一致,上表层声速呈正梯度分布,平均声速在1562.5—1585 m/s 之间,与350 Hz下的cbe较为一致,下表层声速较为稳定在1670 m/s附近,且h11和h12的和近似于表1中h1的结果.

表3 在100—150 Hz 频段表层参数的反演结果Table 3.Inversion results of the surface layer at 100—150 Hz.

图13 细化的海底分层结构Fig.13.Finer-grained seabed layering along the consider track.

实验期间进行了海底重力取样,如表4 为5 km 位点处采样样品分析结果(100 kHz 换能器测量).由表4 可以看出,0—4 m 深度内整体采样趋势为随深度增加,声速逐渐增大,其中的平均声速为1560 m/s,另外,0 km 处海底取样样品分析0—2.35 m 深度内的平均声速为1576 m/s,10 km处海底取样样品分析0—4 m 深度内的平均声速为1574 m/s.均与本文反演得到的cba较为接近.

表4 5 km 站位采样样品分析结果Table 4.Sample results at 5 km site.

3.4 传播损失反演衰减系数

一定距离内的传播损失曲线是反演低频海底衰减系数的有效途径.将实验获得的声传播损失用于海底衰减系数反演,构造代价函数如下:

式中,N为声传播距离上的采样点数;M为阵元个 数; e 代表实 际数据; c 代表拷 贝场数 据; TL 表示以f0为中心频率的1/3 倍频程带宽内的平均声传播损失.

将简正波频散结构反演得到的二三层声速及其密度和匹配场反演得到的表层声速及其密度的均值作为已知,将各层衰减系数视作相等,利用100—500 Hz 各阵元在20 km 内共计39 个距离处的传播损失反演衰减系数,结果如表5,平均误差为2 dB.图14 所示为深度8 和54 m 处阵元在225和375 Hz 处传播损失的对比结果,其中圆圈代表实验数据,直线代表利用反演结果计算的声传播损失曲线,可以看出两者符合较好.

表5 不同频率下反演的衰减系数及反演结果预报声传播损失的平均误差Table 5.The inversion result of attenuation coefficient and their errors.

图14 不同频率的传播损失 (a) 225 Hz;(b) 375 HzFig.14.Comparison of theoretical and experimental transmission-loss at different frequencies:(a) 225 Hz;(b) 375 Hz.

将100—500 Hz 间衰减系数随频率的变化关系进行拟合,如图15 所示,为α=0.36f1.92dB/m,其中,f的单位是kHz,α与文献[27]总结中的砂质海底条件下衰减系数与频率的非线性关系α=(0.38±0.08)f(1.87±0.13)dB/m 基本一致.

图15 不同频率的衰减系数Fig.15.Inverted attenuation coefficient at different frequencies.

为了进一步对比本文反演结果的效果,建立传统单层海底模型,根据图12 中不同频段下等效声速的反演结果,对于声速较为稳定的100—300 Hz频段范围,利用传播损失反演衰减系数,反演结果及对应的声传播损失预报误差如表6 所示.比较表5 和表6,可以看出相同频率下多层海底模型反演结果预报传播损失的误差小于单层等效海底反演结果预报传播损失的误差,且多层海底声速可以同时适用于较宽频段.

表6 单层海底模型下反演结果Table 6.Inversion results for single-layer bottom model.

4 结论

2022 年8 月在南海北部大陆架斜坡海区开展的海上实验中,浅地层剖面仪获取的海底呈现出明显的水平变化分层结构.本文据此建立的多层水平变化海底地声模型,提出一种低频联合反演方法,即先利用简正波频散结构反演各层海底声速、再利用匹配场方法反演对简正波频散结构较不敏感的表层海底声速,最后将得到的结果用到传播损失中反演衰减系数.形成结论如下.

1) 在复杂海底环境中,利用浅地层剖面仪可以预先获取海底分层结构,在此基础上进行反演,可以有效地减少待反演参数,获得更接近实际海底情况的反演结果.

2) 根据敏感性不同,利用简正波频散结构和匹配场联合反演的方式,可以提高各层声速反演结果的可信性,反演结果可以验证多种声学现象.

3) 相比于等效单层海底模型,多层水平变化海底地声模型的反演结果在较低频传播损失上的精度得到明显提升,且多层海底声速可以同时适用于较宽频段.

感谢自然资源部第一海洋研究所提供的浅地层剖面仪数据,感谢中山大学海洋与技术学院提供的采样数据,感谢南方海洋科学与工程广东省实验室(广州)组织实施的“关键海区底质声学特性调查航次”的全体海上工作人员,是他们的辛勤劳动为本论文提供了可靠的数据支撑.