杨 磊 金诸斌 刘祖斌 卢奂采

(浙江工业大学特种装备制造与先进加工技术教育部重点实验室 杭州310014)

0 引言

随着人类对海洋探索的日益重视,水下声学定位系统应用日益广泛,目前已研制成功的有长基线定位系统、短基线定位系统、超短基线定位系统,以及近年来提出的变基线定位系统[1]等。但是,声速受到温度、盐度、深度变化的影响导致它是非均匀分布的,因此声线在海水中的传播路径通常也是弯曲的并且由于海面和海底的反射以及折射使得水下声线传播产生多径效应[2]。如果在水声定位系统中使用平均声速,则会导致系统定位精度较低,特别是当首达声线为非直达声情形时,会严重降低系统的定位精度。针对上述问题,近年来科研人员提出了众多的声速修正方法。

基于统计学的声速修正方法得益于人工智能算法的进步,近年得到了较大的发展。文献[3]针对长基线定位系统,利用粒子群算法求解最优参数计算有效声速。针对短基线定位系统,文献[4]将区域划分自适应粒子群优化算法用于寻找有效声速,结果表明该算法可有效提升系统性能。文献[5]提出利用期望最大化方法来估计水下两点间的有效声速以提高单信标系统定位精度。文献[6]首先将声速剖面进行分段线性处理,并默认声线以直线形式传播,随后通过数学解析式来描述声线的传播过程,进而实现声速修正。文献[7]在文献[6]研究的基础上认为可以结合声速反演的方法来进一步提高声线修正的精度。基于射线模型的声线修正方法由于更加符合声线在水下传播的实际过程,应用更加广泛。考虑到水中存在的噪声影响以及时变的水声环境,文献[8]对基于射线模型下直达声区域内水平距离和时延估算的克拉美罗界展开了研究,这对于基于射线模型的声速修正方法精度对比有重要意义。文献[9]将水下声线的传播形式进行分类后,提出通过射线追踪的方式,建立一个包含不同声传播路径的时延和水平距离信息表格供水下定位使用。Bellhop 是水声中重要计算工具,文献[10]直接通过使用该工具计算首达声线后进行声速修正。文献[11]利用二分迭代法计算直达声区域的首达声线后,再得到对应的有效声速,并通过遗传算法对区域有效声速表进行稀疏后提高实际使用效率。文献[12]在非直达声区域,利用二分迭代法计算两点间所有本征声线后,以最短时延得到首达声线后再计算出对应的有效声速。利用迭代法解算声线初始掠射角是一种有效的方法,但是由于声线在水下传播路径的复杂性,常常导致迭代法失效[13],这也是本文讨论的重点。

受声速剖面结构和声源与接收点间距离影响,两点间首达声线会存在直达、反射和折射3 种形式,而不再拘泥于直达声,此时需要通过计算两点间所有本征声线后以最短时延为标准确定首达声线。但是,在面对声线传播经历多跨度和包含反转点时,利用迭代法解算声线初始掠射角失效的现象在文献中却鲜有涉及。为了更加精准解算本征声线初始掠射角,本文提出依据反转点和多跨度对初始掠射角区间进行分类,并在包含反转点的单个跨度角度区间内采取先预处理后迭代的方式。最后,利用该算法针对南海某地负梯度声速剖面展开仿真研究,计算结果表明该算法在非直达声下依然可以精准计算任意一点首达声线和有效声速值。

1 声线修正方法

1.1 有效声速法

有效声速法[14]是基于水下声学射线模型提出的一类声速修正方法,它利用等效直线代替声线传播的曲线轨迹,不仅降低了声速修正的计算复杂度且有着较高的精度,因此被广泛应用于水下定位系统中。

在水下声学射线模型中,声线的传播路径如图1所示,其过程满足斯涅耳定律:

图1 声线传播路径示意图

其中,c(z) 是深度为z处的声速值,θ(z) 是深度为z处的掠射角,zS和zP分别是声源和接收点处深度,k0是满足斯涅耳定律的声线常数。声线在水下传播的水平距离h和时延t的一般表达式分别为

在不考虑远距离条件下,水下声速变化可以采取竖直分层介质模型,而忽略水平变化的影响[15]。为了精确逼近声线传播的轨迹,将声速剖面沿深度方向进行分层,层内默认为恒定声速梯度,则每层的声线轨迹为一段圆弧。对声源和接收点间声速剖面分为N层,层内声速梯度为

其声线传播的水平距离和时延公式可写为

有效声速可由首达声线传播的斜距与时延的比值计算得到:

其中,R表示两点间的真实斜距,t表示利用射线模型计算的首达声线传播时延。为了精准计算首达声线的传播时延,必须要通过搜索两点间所有本征声线后,以最短时延为标准获取首达声线的时延信息[16]。

1.2 多跨度

随着水平距离的不断增大,声线会经历多次反射或反转形成多跨度后到达接收点。图2 是南海某区域深海声道,首达声线经过海底反射后到达接收点,此时跨度数是0。图3 中声速剖面是文献[17]提及的一类双线性声速剖面,此时首达声线到达形式是经过2 次海面和海底反射后到达接收点,经历的跨度数为2。

图2 首达声线是经过海底反射的形式(跨度为0)

图3 首达声线是经过海面和海底反射的形式(跨度为2)

图4 是声源深度在zs=950 m 和接收点深度在zd=380 m 时,声线经历不同跨度时传播的水平距离和初始掠射角的变化曲线。图中跨度为0、1、2 的角度区间分别对应[14.31 °,16.73 °]、[16.74 °,28.56°]、[28.57°,38.71°],且每个区间内的变化曲线是单调的,但是相邻2 个变化曲线是间断的,因此采用二分迭代法计算声线初始掠射角时,若没有依据不同的跨度对应的角度区间选取迭代初值,可能会导致迭代区间受限而无法得到正确的结果。在图4 中,当声源和接收点间的真实水平距离为13 300 m时,跨度为0 对应的区间内并不存在本征声线;水平距离为10 000 m 时,3 个跨度都存在对应的本征声线,但需要分类计算。

图4 跨度为0、1、2 时声线初始掠射角和水平距离关系

1.3 反转点

水下声线传播形式中的折射会导致声线反转并产生反转点[18],它的特点是声线在此处的掠射角为0 °,且反转点深度随初始掠射角的改变而变化。由斯涅耳定律可得:

式中z0、θ0为声源处的深度和掠射角。如图5 是采用图3 中的声速剖面,计算得到的首达声线是经历2 次反转和海面反射后到达接收点;图6 是采用图4中的声速剖面,计算得到的首达声线是经历反转和海底反射后到达接受点。

图5 首达声线是经过两次反转和海面反射的形式

图6 首达声线是经过反转和海底反射的形式

由斯涅耳定律可知,声线从声源处出发经过不同路径传播到达接收点,若存在c(z)＞c(z0),且初始掠射角在合适范围内就会出现反转点,此处把包含反转点的反转区间临界角称为临界反转初始掠射角。图7 和8 是以南海和北极某地的声速剖面为例,通过改变声源深度计算得到的临界反转初始掠射角。由图可知,当声源位于声速剖面最大点对应的深度处(即c(z0)= max[c(z)])时,临界反转初始掠射角为0°,这就表明此时并不存在反转区间,无论初始掠射角多大,都不会出现声线反转的情形。因此,临界反转初始掠射角与声源位置处声速呈负相关变化关系。

图7 南海某地声速剖面下临界反转初始掠射角随声源深度的变化图

图8 北极某地声速剖面下临界反转初始掠射角随声源深度的变化图

反转点对计算本征声线的影响主要有如下2点。

(1)当声线中存在上(下)反转点时,声线在深度方向所能达到的最高(低)点就是反转点。当初始掠射角在较小区间时声线始终无法到达接收点层,但是随着角度的增大则存在一个最小的能到达接收点层的初始掠射角,如图9 中的θL,将其定义为最小到达初始掠射角。因此实际计算本征声线时,可以忽略小于最小到达初始掠射角的角度区间。

图9 最小到达初始掠射角示意图

(2)对于包含反转点的声线传播形式而言,结合式(5)和图9 可知,由于声线的反转点深度随初始掠射角的变化而变化,使得声线传播水平距离不仅是初始掠射角θ0的函数,还受声速梯度的影响。图10 是采用图4 中的声速剖面,当声源zs=900 m、接收点zd= 380 m、跨度为0 时,包含反转点(角度区间[6.71 °,14.30 °])和不包含反转点(角度区间[14.31 °,16.73 °])的声线传播水平距离随初始掠射角变化曲线。由图可知,包含反转点的变化曲线不仅呈现出非单调变化趋势,同时曲线也是非光滑的(这是由于对声速剖面分层后默认层内声速梯度是线性所导致的,文献[19]也对分层效应带来的误差进行了分析)。当两点间水平距离如图中所示为10 000 m 时,若直接使用二分迭代法解算本征声线初始掠射角,会出现算法不收敛的情形。由于在水下声线传播过程中反转点只出现在部分初始掠射角范围内,为了快速、精准地计算本征声线,只对包含反转点的声线计算过程进行预处理即可。因此,有必要将声线掠射角范围按照声线传播形式中是否包含反转点进行分类。

图10 反转点对声线初始掠射角与水平传播距离变化关系的影响

2 改进的非直达声下声速修正方法

2.1 本征声线时延和水平距离的计算

如图11 所示,声线传播形式可由子跨度进行描述[20],其中S1表示从声源出发经过海面反射或者上反转点折射后,回到声源层的子跨度;S2表示从声源层到达接收点层的子跨度;S3表示从接收点层出射后经过海底反射或下反转点折射后,回到接收点的子跨度。根据具体到达形式,本征声线可分为如下4 种到达形式:第1 类是直达声;第2 类是经过海面反射或者反转点折射后到达接受点;第3 类是经过海底反射或者反转点折射后到达接收点;第4类是经过海面反射或者反转点折射后,再经过海底反射或者反转点折射后到达接收点。

图11 声线传播的子跨度

结合式(5)可计算4 种到达形式的声线传播水平距离:

一个完整的跨度声线所传播的水平距离为

对于存在反转点的声线到达形式,应该首先确定反转点所在的深度位置。当声源和接收点位置确定后,利用式(9)可以确定不同初始掠射角的声线对应的反转点处的声速c(z反转)。然后在声速剖面中寻找c(z反转) 所处的层间位置,结合式(9)可得反转点深度:

综上,任意类型的本征声线的水平传播距离可以表示为

传播时延的计算可以类比声线传播水平距离的计算方法,并结合式(6)得到式(14)。

其中,k表示声线所经历的完整跨度数,T表示完整跨度所传播的时延,ti表示其中一种到达形式下去除完整跨度后部分路径所传播的时延。

2.2 非直达声下有效声速计算

非直达声下有效声速的计算前提是获取首达声线,因此整个算法从搜索本征声线开始。声线初始掠射角向上和向下的角度区间分别是[-90 °,0]和[0,90 °],斯涅耳定律中使用的是余弦函数,因此在算法中统一使用[0,90 °]计算即可。在给定声速剖面和声源位置后,应首先判断声线传播路径是否存在反转点。依据1.3 节可知,当声源位于声速剖面最大值对应深度时(c(z0)= max[c(z)]),临界反转初始掠射角为0 °,那么声线不存在反转点;当声源不在声速剖面最大值对应的深度处时,则声线可能包含反转点。此时,应假设以接收点层为反转点深度,计算得到最小到达初始掠射角θ1:

那么在[0,θ1] 初始掠射角区间内发出的声线都无法到达接收点,计算本征声线时可以忽略。再以声速剖面最大值max[c(z)] 层为反转点所在层,得到临界反转初始掠射角θ2:

那么在[θ2,90°] 初始掠射角区间内必然不会再次出现反转点。综上所述,可以将声线初始掠射角范围初步分为3 个区间:[0,θ1]、[θ1,θ2]、[θ2,90°]。[0,θ1] 是不存在本征声线的初始掠射角区间,[θ1,θ2] 是声线存在反转点的区间,[θ2,90°] 区间内的声线不会存在反转点。

图12 总结了该有效声速算法的逻辑流程。其中有反转点算法(算法1)和无反转点算法分别是针对声线包含反转点和不包含反转点情形下提出的解算本征声线初始掠射角方法。有反转点算法对初始掠射角完成两次分类后,还需要按照步进角Δθ遍历计算单个跨度内每个初始掠射角的声线水平传播距离(步骤1～6),再将两点间真实水平距离减去该结果,寻找到差值的零点位置后,将该零点的角度区间作为迭代角度区间解算本征声线初始掠射角(步骤7～22)。此外,二分迭代法中并没有确定迭代的方向(算法中默认为递增),因此当差值实际变化为递减时,需要调换迭代法中的上下限(步骤21)。无反转点算法可以直接省去有反转点算法中的步骤5～22,在步骤4 确定声线初始掠射角范围后,直接利用二分迭代法(算法2)计算不同情况即可。算法中步进角Δθ的选择可以依据声速剖面分层的情况来确定,通常当声速剖面的相对梯度较大时,建议选择较小的步长。

图12 有效声速计算流程图

3 算例验证

本文算例采用了Argo 系统中南海某处1 月平均声速剖面,将声源放置在不同位置来计算区域有效声速,并对有效声速空间分布展开分析,验证本算法的准确性。算例的水平距离范围是[0,2000] m、深度距离范围为[0,1000] m,计算步距为Δr=Δz=10 m。

图13 至图16 展示了声源在zs=1025 m 下的区域有效声速分布和区域某点首达声线到达形式。结合图13 和图14 可知,有效声速的空间分布随着水平距离增大可以分为3 个区域,首先经历平稳区间,然后经历增大区间,最后经历跃变区间。需要特别说明的是,随着水平距离的增大,上述的变化历程会不断重复。在平稳区间内,首达声线均为直达形式,此时采用传统的迭代法即可计算出首达声线;在增大区间内,首达声线在大声速区域内经历了反转,此时若没有考虑反转点对迭代法解算声线掠射角的影响,则可能导致无法得到正确结果;在跃变区间,首达声线为经过反转和海底反射后到达,且经历的跨度数为1,这表明在跨度为0 的角度区间并不存在首达声线,此时不仅需要考虑反转点的影响,还需要将初始掠射角依据跨度进行分类,分别计算后得到首达声线。

图13 有效声速空间分布图

图14 不同深度下有效声速随水平距离变化示意图

图15 接收点在zr=200 m、R=800 m 时,两点间的本征声线图(跨度为0)

图16 接收点在zr=200 m、R=18 000 m 时,两点间的本征声线图(跨度为1)

图17 和18 是声源位于水面(zs=0 m)时计算的区域有效声速分布图。此时,由于声源位于声速剖面最大值处,那么根据第2 节分析可知声线必然不会出现反转现象。但是,随着水平距离的增大,首达声线从直达声变化为经海底反射后到达(图19),再变化为经海底反射和海面反射后到达的形式(图20),声线经历的跨度数从0 增加到1,需要根据跨度数对初始掠射角进行分类后再解算本征声线的初始掠射角。

图17 有效声速空间分布图

图18 不同深度下有效声速随水平距离变化示意图

图19 接收点在zr=900 m、R=10 000 m 时两点间本征声线图(跨度为0)

图20 接收点在zr=900 m、R=12 000 m 时两点间本征声线图(跨度为1)

上述2 个案例计算中,首达声线都是从直达声开始,随着水平距离的增大,如果满足存在反转点条件,则会经过反转后到达;如果不满足,则会经历反射后到达;最后都会经历多次反射或者反转后到达的情况。计算结果表明,该算法可以准确地解算声线存在反转点和多跨度时的本征声线初始掠射角,并且在非直达声情形下得到任意一点的有效声速值。

4 结论

基于射线模型的有效声速法是通过迭代法解算两点间所有本征声线后,再利用首达声线信息计算出有效声速值的方式实现声速修正的。但是,当声线存在多跨度和反转点时,由于声线水平传播距离随初始掠射角变化曲线存在分段函数和非单调的现象,使得直接使用迭代法无法得到需要的结果。为了解决上述问题,本文提出了一种非直达声下的声速修正方法。

首先,文中指出了反转点对解算本征声线的影响,提出了最小到达掠射角和临界反转初始掠射角2 个概念。前者指的是可以到达接收点层的最小初始掠射角,小于该角度的区间在实际计算时可以直接忽略;后者指的是在小于该角度的区间内,声线均包含反转点。并且依据这2 个角度对初始掠射角区间([0,90 °])进行了初步的分类。与此同时,考虑到声线包含反转点时即使在单个跨度对应的角度区间内,声线传播水平距离随初始掠射角变化曲线为非单调的,算法首先将初始掠射角区间按照Δθ步进角计算所有角度对应的水平传播距离,再与两点间真实水平距离做差后寻找曲线的零点,并且根据实际做差的符号变化来调整迭代法中的上下限。通过这种先预处理再迭代的方式,有效保证了迭代法在声线存在反转点时依然可以准确地计算出首达声线。此外,依据在不同跨度下声线传播水平距离随初始掠射角变化曲线的特点,为了避免迭代区间受到限制从而无法获取准确的结果,有必要将初始掠射角区间再依据不同的跨度数进行二次分类,随后在不同跨度对应的角度区间内分别使用迭代法计算。

最后2 个案例的计算结果表明,本文提出的算法可以在声线存在反转点和多跨度条件下准确解算出首达声线,并得到对应的有效声速值。该方法具有解算速度快且精度较高的优点,对于提高水声定位系统的技术指标有一定的实际价值,因此后续会进行真实的水下定位实验,并对不同海洋水下环境下的应用进行对比研究。